Sampling from custom plausibility function

Modeling
1

I have fitted a Gaussian process model to a simulation of a thermal decomposition process. From this meta-model I can get the model prediction and variance functions. Given some new values the probability for a new observation is exp(-[(obs-model)^2/model variance]).

For this post I used a simplified model with constant variance. I also limited the unknowns to a single variable with a uniform distribution but the aim is to ultimately include prior distributions rather than point estimates for all the terms in the meta-model.

The custom function correctly returns a probability of 0.83 when all the true values are used.
I am a little stuck on how to get STAN to return the 95 % plausible range for the parameters.

“”"
Created on Thu Aug 23 10:09:41 2018

“”"
import stan_utility

Import data into STAN

stan_dict = {}
stan_dict[‘N’] = 1
stan_dict[‘lambda’] = 0.044820773600000002
stan_dict[‘Cp’] = 1303.5546505
stan_dict[‘rho’] = 301.41640818000002
stan_dict[‘EaPerR_true’] = 12544.486510999999
stan_dict[‘dH’] = -48024275.07
stan_dict[‘K0’] = 54799472.333000004
stan_dict[‘r’] = 0.21342190019999999
stan_dict[‘Toven’] = 340.43299999999999
stan_dict[‘alpha’] = 9.6725974666999992
stan_dict[‘y’] = 1

Define Model

model_code = “”"

functions {
real sampled_log(real lambda, real Cp,real rho, real dH, real K0,real r,real Toven,real alpha,real EaPerR){
return exp(-1* ((Toven - (460.592092717159 + 7.25706884765198 * ((lambda
- 0.22) / 0.18) + -15.038050709874 * ((rho - 850) / 750)
+7.3771764342968 * ((dH - (-30000000)) / 20000000) +
-23.5986485095225 * ((K0 - 50005000) / 49995000) +
-13.5323572491023 * ((r - 0.1375) / 0.1125) + 3.0360815919351 * ((
alpha - 10) / 5) + ((lambda - 0.22) / 0.18) * (((r - 0.1375) /
0.1125) * 9.50770545201133) + ((rho - 850) / 750) * (((dH - (
-30000000)) / 20000000) * -5.28112528903226) + ((rho - 850) / 750)
* (((r - 0.1375) / 0.1125) * 15.1745700813374) + 277.960976564701
* ((EaPerR - 15500) / 9500) + ((rho - 850) / 750) * (((EaPerR
-15500) / 9500) * -21.1397159704067) + ((EaPerR - 15500) / 9500)
* (((K0 - 50005000) / 49995000) * 24.314407428981) + ((EaPerR
-15500) / 9500) * (((r - 0.1375) / 0.1125) * -24.5774862093203))) ^ 2 )/17);
}
}
data {
int<lower=0> N;
real lambda;
real Cp;
real rho;
real EaPerR_true;
real dH;
real K0;
real r;
real Toven;
real alpha;
real y;
}
parameters {
real<lower=10000, upper=20000> EaPerR;

}
model {
y ~ sampled_log(lambda,Cp,rho,dH,K0,r,Toven,alpha);
}
generated quantities {
real expected;
expected = sampled_log(lambda,Cp,rho,dH,K0,r,Toven,alpha,EaPerR_true);
}
“”"

#Model settings
modelName = ‘PredictEaPerR’
nIter=2500
nChains=4

Load model or create new

sm = stan_utility.StanModel_cache(model_code=model_code,model_name=modelName)

fit = sm.sampling(data=stan_dict,chains=nChains,iter=nIter)

print(fit)

2

Instead of trying to calculate the probability directly with the function what I needed to do was:
obs ~ normal(model,model sigma).

so the code in my model block becomes:

Toven ~ normal(model(lambda,Cp,rho,dH,K0,r,alpha,EaPerR),4.123);

and the function is simply the model equation.