I am trying to output a posterior probability distribution of parameters using 10000 points extracted from the posterior of a GP that was trained on a small set of simulation data. The prior of the parameters given the outputs is assumed to be non-informative uniform for each parameter. There are 13 parameters and 13 output properties in the model. I am using a multivariate normal likelihood function that is derived as the following:
where Yobs is the experimental value for the property calculated in the simulation, muGPM is the mean GP prediction and sigma is the diagonal covariance matrix with each element given by the sum of the variance of the GP prediction and a noise variance where the noise standard deviation prior is assumed to be Half Normal.
Could someone take a look at my code and see if I am in the right direction? I am new to data science and to PyMC3 so any advice would be extremely helpful. I established the 13 parameter priors separately, but I will only include a couple because they are similar and I am not sure if this is the best way to guide this. I establish testvals based on optimized parameter outputs from the initial simulation. The chain lengths are at 10000 just to get a decent initial look at the parameter distributions, but 80000 is the approximated length for convergence.
with pm.Model() as basic_model: var1 = pm.Uniform('var1', lower=16.0, upper=18.0, testval=16.94) var12 = pm.Uniform('var12', lower=0.4, upper=0.8, testval=0.679) var13 = pm.Uniform('var13', lower=6.0, upper= 10.2, testval=8.16) sigNoise = pm.HalfNormal('sigNoise') cov_noise = pm.Deterministic("cov_noise",(cov + sigNoise**2) * np.eye(13)) likelihood = pm.MvNormal('likelihood', mu=means, cov=cov_noise, observed=yobs) trace = pm.sample(10000, chains=2,cores=1)`