stanhf normsys example

data.json

{
    "_metadata": {
        "fix_mu": "from POI.stan_data_card [config.py:L155]",
        "fixed_mu": "from POI.stan_data_card [config.py:L155]",
        "lu_channel_signal_normsys_mu": "from NormSys.stan_data_card [modifier.py:L293]",
        "lu_mu": "from POI.stan_data_card [config.py:L155]",
        "nominal_channel_background": "from Sample.stan_data_card [sample.py:L73]",
        "nominal_channel_signal": "from Sample.stan_data_card [sample.py:L73]",
        "observed_channel": "from Channel.stan_data_card [channel.py:L66]"
    },
    "fix_mu": 0,
    "fixed_mu": 1.0,
    "lu_channel_signal_normsys_mu": {
        "1": 0.9,
        "2": 1.1
    },
    "lu_mu": {
        "1": 0.0,
        "2": 10.0
    },
    "nominal_channel_background": [
        50.0,
        60.0
    ],
    "nominal_channel_signal": [
        5.0,
        10.0
    ],
    "observed_channel": [
        44,
        62
    ]
}

main.stan

// histfactory json examples/normsys.json
// histfactory spec version 1.0.0
// converted with stanhf 1.1.0

// no patch applied

functions {
  // [.local/share/pipx/venvs/stanhf/lib/python3.13/site-packages/stanhf/stanhf.stanfunctions]
  real poisson_real_lpdf(data vector k, vector lambda) {
    return sum(k .* log(lambda) - lambda - lgamma(k + 1.));
  }
  
  vector term_interp(real alpha, data vector x, data tuple(vector, vector) lu) {
    if (alpha > 1.) {
      return alpha * (lu.2 - x);
    }
    
    if (alpha < -1.) {
      return alpha * (x - lu.1);
    }
    
    vector[size(x)] s = 0.5 * (lu.2 - lu.1);
    vector[size(x)] a = 0.0625 * (lu.2 + lu.1 - 2. * x);
    real alpha_square = square(alpha);
    real r = alpha_square * (alpha_square * (alpha_square * 3. - 10.) + 15.);
    
    return r * a + alpha * s;
  }
  
  real factor_interp(real alpha, data tuple(real, real) lu) {
    if (alpha > 1.) {
      return lu.2 ^ alpha;
    }
    
    if (alpha < -1.) {
      return lu.1 ^ (-alpha);
    }
    
    real log_l = log(lu.1);
    real log_u = log(lu.2);
    
    vector[6] b = [lu.2 - 1., lu.1 - 1., log_u * lu.2, -log_l * lu.1,
                   square(log_u) * lu.2, square(log_l) * lu.1]';
    
    matrix[6, 6] m = [[0.9375, -0.9375, -0.4375, -0.4375, 0.0625, -0.0625],
                      [1.5, 1.5, -0.5625, 0.5625, 0.0625, 0.0625],
                      [-0.625, 0.625, 0.625, 0.625, -0.125, 0.125],
                      [-1.5, -1.5, 0.875, -0.875, -0.125, -0.125],
                      [0.1875, -0.1875, -0.1875, -0.1875, 0.0625, -0.0625],
                      [0.5, 0.5, -0.3125, 0.3125, 0.0625, 0.0625]];
    
    return 1. + alpha ^ linspaced_row_vector(6, 1, 6) * m * b;
  }
}
data {
  vector[2] nominal_channel_background; // from Sample.stan_data [sample.py:L66]
  vector[2] nominal_channel_signal; // from Sample.stan_data [sample.py:L66]
  int<lower=0, upper=1> fix_mu; // from POI.stan_data [config.py:L144]
  real fixed_mu; // from POI.stan_data [config.py:L144]
  tuple(real, real) lu_mu; // from POI.stan_data [config.py:L144]
  tuple(real, real) lu_channel_signal_normsys_mu; // from NormSys.stan_data [modifier.py:L286]
  array[2] int observed_channel; // from Channel.stan_data [channel.py:L59]
}
parameters {
  array[1 - fix_mu] real<lower=lu_mu.1, upper=lu_mu.2> free_mu; // from POI.stan_pars [config.py:L122]
}
transformed parameters {
  vector[2] expected_channel_background = nominal_channel_background; // from Sample.stan_trans_pars [sample.py:L59]
  vector[2] expected_channel_signal = nominal_channel_signal; // from Sample.stan_trans_pars [sample.py:L59]
  real mu = fix_mu ? fixed_mu : free_mu[1]; // from POI.stan_trans_pars [config.py:L130]
  expected_channel_signal *= factor_interp(mu, lu_channel_signal_normsys_mu); // from NormSys.stan_trans_pars [modifier.py:L300]
  vector[2] expected_channel = expected_channel_background
                               + expected_channel_signal; // from Channel.stan_trans_pars [channel.py:L43]
}
model {
  observed_channel ~ poisson(expected_channel); // from Channel.stan_model [channel.py:L51]
  mu ~ std_normal(); // from StandardNormal.stan_model [modifier.py:L319]
}
generated quantities {
  array[2] int rv_expected_channel = poisson_rng(expected_channel); // from Channel.stan_gen_quant [channel.py:L73]
}