Hi there, I’ve attached some images to demonstrate the setup of my model. I have a problem learning the sig_s parameter shown in the traces - the vertical line indicates the expected value in the plot. Is there some obvious problem in the set up of the model in PyMC that might be causing me some issues? I am relatively new to PyMC.
Where do the “expected values” come from?
Hey, thanks for the reply. The expected values are the mean values used in normal distributions in generating the ‘observed’ data used in the model . Effectively I am trying to relearn these parameters using the hierarchical model.
If you have code to simulate data, could you post it? That might help figure out what’s going on.
Hi Cluhmann, please see below:
from ansys.mapdl.core import launch_mapdl
from ansys.mapdl.core.convert import convert_apdl_block
import numpy as np
import pyvista
import matplotlib.pyplot as plt
import os
import pymc as pm
import arviz as az
import graphviz
from scipy.stats import halfnorm
import pytensor
import pytensor.tensor as pt
from linear_population import linear_model
# Keep it reproducable
np.random.seed(0)
##### Generate data Y
# set up our data
K = 10 # number of structures
N = 20 # number of data points
data_description = "_"+ str(T) + "T_" + str(N) + "N"
# # True global parameters (these will need priors on in the inference). We use these to draw local parameters
µ_s0_true = 20e6
µ_s1_true = 2e6
σ_s0_true = 2e5
σ_s1_true = 1e4
# parameters for local distributions
# Draw K local means
s_k0 = np.random.normal(µ_s0_true, µ_s1_true, size=K)
# Draw K local variances
s_k1 = np.random.normal(σ_s0_true, σ_s1_true, size=K)
# ---- draw realisations for each K
s_k = []
turbines = []
# balanced bit
for i in range(K-1):
s_k.append(np.random.normal(s_k0[i], s_k1[i], size=N)) # size N
for j in range(N):
turbines.append(i)
s_k = np.array(s_k).flatten()
turbines = np.array(turbines).flatten()
# imbalanced data for final K
s_k = np.append(s_k, np.random.normal(s_k0[-1], s_k1[-1], size=2))
turbines = np.append(turbines, np.array([K-1, K-1]))
sigma0_true = 0.0001 # measurement noise
### True model
truemodel = .. this bit calls an external model (excluded here for privacy)
# Create true data
measurement_noise = np.random.normal(0, sigma0_true, size=100)
data = measurement_noise + truemodel # add some gaussian noise to the data, scaled by sigma
np.savetxt("data" + data_description + ".csv", data, delimiter=",") # save
np.savetxt("s_k0" + data_description + ".csv", s_k0, delimiter=",") # save
np.savetxt("s_k1" + data_description + ".csv", s_k1, delimiter=",")
np.savetxt("s_k" + data_description + ".csv", s_k, delimiter=",") # save
np.savetxt("turbines" + data_description + ".csv", turbines, delimiter=",") # save
np.savetxt("y_true" + data_description + ".csv", truemodel, delimiter=",") # save