Hi,
here is the important part of the code:
def mcmc_inference_diffusion(Data, Surface, Thickness, A_fact):
#### Setting up the Data and Timeframe
time = np.array(Data['Time'])[:i+1]*1e-9 #s
y_combined = np.zeros((len(time),len(Surface)))
for a in range(len(Surface)):
y_combined[:,a] = np.array(Data[str(a)])[:i+1]
with pm.Model() as model:
#### Defining Model Parameters
## Diffusion
S_ratio_power = pm.Normal("S_ratio_power",0,3, initval=1)
S_ratio_model = pm.Deterministic('S_ratio_model', 10**(S_ratio_power))
Beta_val_power = pm.HalfNormal("Beta_val_power",1, initval=1)
Beta_val_model = pm.Deterministic("Beta_val_model", np.exp(-Beta_val_power))
S_substrate_power_offset = pm.HalfNormal('S_substrate_power_offset',1, initval=0.1)
S_substrate_model = pm.Deterministic('S_substrate_model', 10**(S_substrate_power_offset*3))
## Recombination
k1_power_offset = pm.HalfNormal('k1_power_offset', 1)
k1_model = pm.Deterministic('k1_model', 10**(3 + k1_power_offset*2))
p0_power_offset = pm.HalfNormal('p0_power_offset', 1)
p0_model = pm.Deterministic('p0_model', 10**(13 + p0_power_offset*2))
#### Simulation of Time-Resolved PL
N_calc = diffusion(S_ratio_model, Beta_val_model, S_substrate_model, k1_model, p0_model, shared(time), shared(Surface),Thickness, shared(A_fact))
## Likelihood Function (Normal Distribution)
Y_obs = pm.Normal("Y_obs_late", mu=at.sqrt(y_combined[1:i,:]/N_calc[1:i,:]), sigma=1, observed=np.ones(shape=np.shape(y_combined[1:i,:]))) #sigma=np.sqrt((i-max_locator)/5000)
#### Draw Samples from the Posterior Distribution
trace = pm.sample_smc(progressbar=False)
return trace
The function diffusion() is:
def diffusion(S_ratio_model, Beta_val_model, S1_model, k1_model, p0_model, time, Surface, thickness, A_fact):
### Bulk Recombination Parameters
k1 = k1_model
p0 = p0_model
beta_est = beta_estimator(thickness, S1_model, S_ratio_model, Beta_val_model) # Polynomial Function
Diffusion = at.switch(at.gt(beta_diffusion,0.5),Diff_calc_high(S1,S2,beta_diffusion,thickness),Diff_calc_low(S1,S2,beta_diffusion,thickness)) #Diff_calc_high and _low are polynomial functions
beta_model = beta_est*np.pi/(thickness*1e-7)
S_front = at.switch(at.eq(Surface,1),S1_model,S1_model*S_ratio_model)
### Equation (24)
S_front_4d = S_front.dimshuffle(0,'x','x','x')
S_front_4d.broadcastable
(False, True, True, True)
Diffusion_4d = Diffusion#.dimshuffle(0,'x','x','x')
beta_4d = beta_model.dimshuffle('x', 0 ,'x','x')
beta_4d.broadcastable
(True, False, True, True)
### Defining the spacial domains
x = np.arange(20)
stretch = 2.5 # Factor that stretches the tanh-grid towards the edges
thickness_tanh_spacing = thickness/2*(1-np.tanh(stretch*(1-2*x/len(x)))/np.tanh(stretch))
z_array = at.as_tensor_variable(thickness_tanh_spacing*1e-7)
z_array_4d = z_array.dimshuffle('x','x',0, 'x')
z_array_4d.broadcastable
(True, True, False, True)
U_z = at.cos(beta_4d*z_array_4d)+S_front_4d/(Diffusion_4d*beta_4d)*at.sin(beta_4d*z_array_4d)
U_z.broadcastable
(False, False, False, True)
### Equation (25)
A_fact_4d= (A_fact).dimshuffle(0,'x','x','x')
A_fact_4d.broadcastable
(False, True, True, True)
A_param = ((at.exp(-A_fact_4d*z_array_4d)*U_z).sum(axis=2)/(at.power(U_z,2)).sum(axis=2))
A_param_4d = A_param.dimshuffle(0,1,2,'x')
A_param_4d.broadcastable
(False, False, True, True)
time_4d = time.dimshuffle('x','x','x',0)
time_4d.broadcastable
(True, True, True, False)
n_tz0 = (A_param_4d*U_z * at.exp(-(Diffusion_4d*at.power(beta_4d,2))*time_4d)).sum(axis=1) # sum over all beta values
## Recombination
time_3d = time.dimshuffle('x','x',0)
time_3d.broadcastable
(True, True, False)
n_tz1 = (n_tz0*at.exp(-(k1*(n_tz0+2*p0)/(n_tz0+p0))*time_3d))
n_tz = n_tz1.sum(axis=1)# sum over thickness
### Equation (13) & (32)
I_t = (2*p0*n_tz+n_tz**2)
### Normalizing to compare with the data
PL_value_norm = at.mul(I_t.T,1/I_t.max(axis=1))
return PL_value_norm
Sorry for the long code example.
Any help is greatly appreciated!