Custom posterior distribution in Turing

General Usage
, ,
1

Hello,
I have the following simple Turing model:

@model function opt(y, λ, sa, sb)
    α ~ Exponential(λ)
    β ~ Exponential(λ)
end

I would like to sample from a posterior custom distribution similar to:

  (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:

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 . 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.

1 Like
2

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:

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
1 Like
3

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.

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


@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:

nuovib
nuovib1000×500 56.9 KB

I also tried to insert trivial istructions in logpdf and randr like:

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.

4

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.

5

Hi,

The function I wrote is p(α,β | y) (what I call “posterior”) but I think I should use p(y|α,β) (the likelihood, as you wrote).

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?

6

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:”

@model function mutation(y, λ)
    α ~ Exponential(λ)
    β ~ Exponential(λ)
    for i in 1 : length(y)
        y[i] ~ Beta(α, β)
    end
end