# Obtaining the posterior distribution of parameters in Turing using the MCMC method
**URL:** https://discourse.julialang.org/t/obtaining-the-posterior-distribution-of-parameters-in-turing-using-the-mcmc-method/89111
**Category:** Probabilistic Programming
**Tags:** turing
**Created:** [October 22, 2022, 4:43pm UTC](https://discourse.julialang.org/t/obtaining-the-posterior-distribution-of-parameters-in-turing-using-the-mcmc-method/89111 "2022-10-22T16:43:33Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@Aru](https://discourse.julialang.org/u/Aru)
#### Post date: [October 22, 2022, 4:43pm UTC](https://discourse.julialang.org/t/obtaining-the-posterior-distribution-of-parameters-in-turing-using-the-mcmc-method/89111/1 "2022-10-22T16:43:33Z")
Hello everyone. I am a student majoring in physics with no solid background in statistics. I recently tried to use Turing, but I found that there are some hard-to-understand problems that make it difficult for me to use Turing in my work. To express my problem more clearly, I will show some code below.
My ultimate goal is to obtain the posterior distributions of the parameters with MCMC method in some physical models. Let’s say there is a machine which can measure the time (`t`), observational value (`obsVal`), and observational uncertainty (`uncertainty`) at each observation. The model is quite simple (just as an example):
```julia
begin
a = 1.5
b = 3.2
c = 0.08
f(x) = (a * x + sinpi(x / c)) * cospi(b)
end
```
where `a`, `b` and `c` are three parameters and `f(x)` is the model. The ground-truth values for the three parameters are given above. Then we measure the machine for 30 times, and we have 30 values for `time`, `trueVals`, `obsVals`, and `uncertainties`:
```julia
function observations(time::Array)
uncertainties = 2 * rand(Float64, length(time)) # Uncertainties ∈ [0, 2]
errors = rand.(Normal.(0, uncertainties)) # Assume error ~ Normal(0, uncertainty)
trueVals = f.(time)
obsVals = trueVals .+ errors
return trueVals, obsVals, uncertainties
end
```
The plot of the model is somehow like
```julia
begin
p = plot(xlabel="Time", ylabel="Value")
time = 20 * rand(Float64, 30) # Observate 30 times, time ∈ [0, 20]
trueVals, obsVals, uncertainties = observations(time)
plot!(p, time, obsVals, yerror=uncertainties, seriestype=:scatter, label="Obs")
plot!(p, time, trueVals, seriestype=:scatter, label="True")
end
```
Now I want to define a Turing model, and build MCMC chains to obtain the posterior distributions of `a`, `b`, and `c`:
```julia
@model function MyModel(obsVals, uncertainties)
# Prior
a ~ Uniform(1, 2)
b ~ Uniform(2.5, 3.5)
c ~ Normal(0, 1)
# Codes below may be wrong, I'm not sure...
obsVals .~ Normal.(f.(obsVals), uncertainties)
end
```
My questions are summarized as below:
1. How should I define a Turing model, which can take `obsVals`, `uncertainties`, and the custom defined model `f(x)` into consideration simultaneously, then build and run MCMC chains with suitable algorithm (like `HMC`) to get the posterior distributions of parameters?
2. Previously I used `emcee` to do similar tasks in Python. And it is the programmer’s obligation to do some burdensome tasks like calculating likelihood and checking marginal conditions (e.g., see functions `log_prior` and `log_probability` in [Fitting a model to data — emcee](https://emcee.readthedocs.io/en/stable/tutorials/line/#marginalization-uncertainty-estimation)). Nevertheless, similar functions are not found in Turing. Does this mean that Turing automatically take these tasks?
Thanks for any helpful suggestions! Really rudimentary questions… 🙂
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [October 22, 2022, 4:57pm UTC](https://discourse.julialang.org/t/obtaining-the-posterior-distribution-of-parameters-in-turing-using-the-mcmc-method/89111/2 "2022-10-22T16:57:21Z")
> [@Aru](#):
>
> `obsVals .~ Normal.(f.(obsVals), uncertainties)`
I’m confused here why f is predicting obsVals from the obsVals rather than from some covariate values. Like I’d expect your observations are predicted from the time and some parameters or something similar.
Beyond that, you have the right basic idea. Once you have your model defined you instantiate it with the particular data you’ve collected, then call sample on it. The Turing docs do a good job of giving those examples
---
### Author:  [@Aru](https://discourse.julialang.org/u/Aru)
#### Post date: [October 22, 2022, 5:45pm UTC](https://discourse.julialang.org/t/obtaining-the-posterior-distribution-of-parameters-in-turing-using-the-mcmc-method/89111/3 "2022-10-22T17:45:20Z")
Thanks for your reply. Yes, this is exactly what I feel confusing. The thing here is, at any time `t` I can obtain a measurement `obsVal ± uncetainty`. Then how can I express this function with the randomly sampled parameters `a`, `b`, and `c` in the Turing model, then measure the corresponding likelihood?
I rewrite the function as
```julia
f(time, a, b, c) = (a * x + sinpi(x / c)) * cospi(b)
```
and I define the Turing model as
```julia
@model function MyModel(obsVals, time, uncertainties)
# Prior
a ~ Uniform(1, 2)
b ~ Uniform(2.5, 3.5)
c ~ Normal(0, 1)
obsVals .~ Normal.(f.(time, a, b, c), uncertainties)
end
```
It seems like the model could (literally) work now. Hope this code is correct but I still cannot feel certain on its correctness… Maybe I should have a try on some more sophisticated models to check whether I get the right idea 😐 Anyway, thank you again!
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [October 22, 2022, 6:35pm UTC](https://discourse.julialang.org/t/obtaining-the-posterior-distribution-of-parameters-in-turing-using-the-mcmc-method/89111/4 "2022-10-22T18:35:37Z")
This looks more likely to be right. Give it a try and call sample see if you get something reasonable