# Custom posterior distribution in Turing
**URL:** https://discourse.julialang.org/t/custom-posterior-distribution-in-turing/55457
**Category:** General Usage
**Tags:** question, turing, distributions
**Created:** [February 17, 2021, 12:45pm UTC](https://discourse.julialang.org/t/custom-posterior-distribution-in-turing/55457 "2021-02-17T12:45:21Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@mcql](https://discourse.julialang.org/u/mcql)
#### Post date: [February 17, 2021, 12:45pm UTC](https://discourse.julialang.org/t/custom-posterior-distribution-in-turing/55457/1 "2021-02-17T12:45:21Z")
Hello,
I have the following simple Turing model:
```julia
@model function opt(y, λ, sa, sb)
α ~ Exponential(λ)
β ~ Exponential(λ)
end
```
I would like to sample from a posterior custom distribution similar to:
```julia
(beta(α, β) ^ - length(y)) * exp(- α * (λ - sa)) * exp(- β * (λ -sb))
```
(where beta is not the Beta distribution but the beta function in SpecialFunctions.jl)
I have tried to use the for loop on the posterior function:
```julia
for i in 1:length(y)
y[i] ~ (beta(α, β) ^ - length(y)) * exp(- α * (λ - sa)) * exp(- β * (λ -sb))
end
```
But doesn’t work because the function is not of type Distribution. I think I should define a custom posterior distribution.
I have alrady read the documentation at [Advanced Usage](https://turing.ml/dev/docs/using-turing/advanced) . But the explained example of custom distribution is very different from the distribution I would like to obtain. Could you suggest me any {solution/something to read/rtfm} to define a custom distribution similar to that function? In general I would like to undestrand how define a custom distribution having a statistical model definition.
Thank you.
---
### Author:  [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)
#### Post date: [February 17, 2021, 5:46pm UTC](https://discourse.julialang.org/t/custom-posterior-distribution-in-turing/55457/2 "2021-02-17T17:46:30Z")
I think you may have your terminology a bit confused–the _posterior_ is what you estimate by running MCMC on a specific instance of your model. In that sense, the Turing model _is_ the “custom posterior.” I think you’re actually looking for a way to define a custom distribution to use for your data _likelihood_. Am I interpreting it correctly that the formula you included above is a probability density function for `y`, conditional on `α`, `β`, `λ`, `sa`, and `sb`?
The Turing docs you linked to use a very simple distribution as an example, but show pretty much what you need to do. For your case, if I’m understanding it correctly, the implementation would look something like:
```julia
using Distributions
struct MyDistribution{T} <: ContinuousUnivariateDistribution
α::T
β::T
sa::T
sb::T
end
function Distributions.logpdf(d::MyDistribution, y)
α, β, λ, sa, sb = d.α, d.β, d.λ, d.sa, d.sb
# code to calculate log-PDF of `y`, given `α`, `β`, `λ`, `sa`, and `sb`
end
function Distributions.rand(rng::AbstractRNG, d::MyDistribution)
# code to draw random sample from `d`, given `α`, `β`, `λ`, `sa`, and `sb`
end
```
---
### Author:  [@mcql](https://discourse.julialang.org/u/mcql)
#### Post date: [February 18, 2021, 8:10am UTC](https://discourse.julialang.org/t/custom-posterior-distribution-in-turing/55457/3 "2021-02-18T08:10:26Z")
Hi,
> think you’re actually looking for a way to define a custom distribution to use for your data _likelihood_ .
Yes, it’s correct.
> that the formula you included above is a probability density function for `y` , conditional on `α` , `β` , `λ` , `sa` , and `sb`
In this case the formula is p(α,β | y).
λ, sa and sb are numbers, not distributions. Sa, and sb can be built using y.
I tried to define logpdf() using the formula above but i think it’s not correct.
```julia
function Distributions.logpdf(d::MyDist, y::Real)
α, β, λ = d.α, d.β, d.λ
# code to calculate log-PDF of `y`, given `α`, `β`, `λ`, `sa`, and `sb`
sa = sum([log(i) for i in y])
sb = sum([log(1 - i) for i in y])
return (beta(α, β) ^ - length(y)) * exp(- α * (λ - sa)) * exp(- β * (λ - sb))
end
```
```julia
@model function mutation(y, λ)
α ~ Exponential(λ)
β ~ Exponential(λ)
for i in 1 : length(y)
y[i] ~ MyDist(α, β, λ)
end
end
```
When i try to plot the sample I always obtain this plot:
I also tried to insert trivial istructions in logpdf and randr like:
```julia
return 10
```
but the plot is always the same.
I have also set function Distributions.logpdf(d::MyDist, y::AbstractVector{\<:Real})
The samples are generated with Gibbs(MH(:α), MH(:β))
Thank you.
---
### Author:  [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)
#### Post date: [February 18, 2021, 7:51pm UTC](https://discourse.julialang.org/t/custom-posterior-distribution-in-turing/55457/4 "2021-02-18T19:51:46Z")
Okay, a couple of follow up questions. It looks like your custom distribution actually a multivariate distribution, is that the case? It kind of looks like that, since the formula depends on all the elements of `y`…what’s the name of the distribution? What kind of process is it trying to capture?
Also, is the formula you gave for the PDF, or the log(PDF)? If it’s actually for the PDF, your custom `logpdf` function will probably be returning a tiny number on the wrong measurement scale, which would explain why your MCMC chains just look like samples from the Exponential priors.
---
### Author:  [@mcql](https://discourse.julialang.org/u/mcql)
#### Post date: [February 18, 2021, 9:45pm UTC](https://discourse.julialang.org/t/custom-posterior-distribution-in-turing/55457/5 "2021-02-18T21:45:22Z")
Hi,
> [@ElOceanografo](#):
>
> It looks like your custom distribution actually a multivariate distribution, is that the case?
The function I wrote is p(α,β | y) (what I call “posterior”) but I think I should use p(y|α,β) (the likelihood, as you wrote).
```julia
for i in 1:length(y)
y[i] ~ p(y|α,β)
end
```
Where my p(y|α,β) is a Beta distribution. Is it the correct way to sample?
---
### Author:  [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)
#### Post date: [February 18, 2021, 11:10pm UTC](https://discourse.julialang.org/t/custom-posterior-distribution-in-turing/55457/6 "2021-02-18T23:10:52Z")
The posterior distribution, p(α, β | y), is what you estimate using Bayesian inference (typically MCMC). If you’ve already got a closed-form expression for p(α, β | y), you don’t need to use Turing at all, you can just plug the values of y, α, and β in and calculate the posterior directly.
If I’m understanding correctly, though, you have some data `y` which are i.i.d. and restricted to the interval between 0 and 1, which you’re modeling using a beta distribution. Is that right? If so, the following does what you want “out of the box:”
```julia
@model function mutation(y, λ)
α ~ Exponential(λ)
β ~ Exponential(λ)
for i in 1 : length(y)
y[i] ~ Beta(α, β)
end
end
```