Skip to content

Instantly share code, notes, and snippets.

@andrewfowlie
Last active March 28, 2025 06:36
Show Gist options
  • Save andrewfowlie/8734f8bf0b2841300f46ff2d4b0bba5e to your computer and use it in GitHub Desktop.
Save andrewfowlie/8734f8bf0b2841300f46ff2d4b0bba5e to your computer and use it in GitHub Desktop.
stanhf example
{
"_metadata": {
"fix_mu": "from POI.stan_data_card [config.py:L155]",
"fixed_mu": "from POI.stan_data_card [config.py:L155]",
"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_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
]
}
// histfactory json examples/model.json
// histfactory spec version 1.0.0
// converted with stanhf 1.0.0
// no patch applied
functions {
// [.local/lib/python3.12/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]
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 *= mu; // from Factor.stan_trans_pars [modifier.py:L85]
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]
}
generated quantities {
array[2] int rv_expected_channel = poisson_rng(expected_channel); // from Channel.stan_gen_quant [channel.py:L73]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment