Thanks for the information. But still it seems like it skips the likelihood function, just samples from the prior distribution.

The whole problem is summarized like: `theta_1,theta_2`

are the control parameters, given by lognormal prior distribution, they somehow determine the observations `f_1, f_2`

through a function. Given a pair of theta, f_1 and f_2 could be calculated.

Now, I want to infer theta given by an observation of fm_1,fm_2. Likelihood function is set up by a custom function.

Based on your advise, I tried both of them, I don’t know where is wrong:

(1) The effect of the observation has already been taken into account. The likelihood function just simply loop through all the parameters. pm.Potential was used.

```
with pm.Model() as model:
theta_1 = pm.Lognormal('theta_1', mu=mu_1, sd=sigma_1)
theta_2 = pm.Lognormal('theta_2', mu=mu_2, sd=sigma_2)
def likelihood(theta_1,theta_2):
a = m_1*m_2
b = m_1*k_2*theta_2+m_2*k_2*theta_2+m_2*k_1*theta_1
c = theta_1*theta_2*k_1*k_2
# Controlled parameter
i0 = 9
w_1_square = (b-np.sqrt(b**2-4*a*c))/(2*a)
w_2_square = (b+np.sqrt(b**2-4*a*c))/(2*a)
f_1_square = w_1_square/(4*np.pi**2)
f_2_square = w_2_square/(4*np.pi**2)
return np.exp(-2**(i0-2)*((f_1_square/fm_1**2-1)**2+
(f_2_square/fm_2**2-1)**2))
like = pm.Potential('like',likelihood(theta_1,theta_2))
start = pm.find_MAP(model = model)
step = pm.Metropolis()
trace = pm.sample(50000, step = step,start=start)
pm.traceplot(trace)
print(pm.summary(trace))
```

(2) The second method is to consider the observation separately: using pm.DensityDist

```
with pm.Model() as model:
theta_1 = pm.Lognormal('theta_1', mu=mu_1, sd=sigma_1)
theta_2 = pm.Lognormal('theta_2', mu=mu_2, sd=sigma_2)
def likelihood(f_1, f_2):
a = m_1*m_2
b = m_1*k_2*theta_2+m_2*k_2*theta_2+m_2*k_1*theta_1
c = theta_1*theta_2*k_1*k_2
# Controlled parameter
i0 = 9
f_1_square = (b-np.sqrt(b**2-4*a*c))/(2*a)/(4*np.pi**2)
f_2_square = (b+np.sqrt(b**2-4*a*c))/(2*a)/(4*np.pi**2)
return np.exp(-2**(i0-2)*((f_1_square/f_1**2-1)**2+
(f_2_square/f_2**2-1)**2))
like = pm.DensityDist('like',likelihood, observed={'f_1':fm_1,'f_2':fm_2})
start = pm.find_MAP(model = model)
step = pm.Metropolis()
trace = pm.sample(10000, step = step,start=start)
```

Apparently, they are sampling from the prior distribution. I just start to use PyMc3 this week and fail to figure out what went wrong.