Using a posterior from a previous sample as a prior

Specific Domains Probabilistic programming
1

Sorry if this is a question that has already been asked, or if this is the wrong place for something like this - I’ve done a bit of searching and haven’t found anything that works for the Turing package. In short, what I want to do is, given some model with priors & a likelihood function, run an MCMC algorithm and get my chain/posterior output, then use that posterior as the new prior for the same likelihood function (if you’re familiar with data cloning imagine something like this). I’m not aware of a way to do this in Turing, and the closest thing I can come up with is to fit some sort of well-known (e.g., gamma) distribution to the posteriors and use that as a prior. Any thoughts?

2

Try wrapping the multivariate kernel density estimator from https://github.com/JuliaRobotics/KernelDensityEstimate.jl in a custom distribution with rand, logpdf, etc. (follow the Turing docs for custom distributions). If you only need univariate kernel density estimators, there is also https://github.com/JuliaStats/KernelDensity.jl. Once you do this, you can use this custom distribution in Turing as a prior to “update” the posterior sample. This would be a valuable contribution so if you do it, a PR and/or tutorial would be welcome.

2 Likes
3

Good idea! I’ve done my best to do this for a very simple example (copied directly from the Turing docs). I use the KernelDensityEstimate package (KernelDensity as far as I can tell only offers a bivariate kernel as opposed to a kernel that can work in 3 or more dimensions) to get the KDE from the chain, which works fine; the only part I am struggling with is the logpdf function. This might be because I’m a tad new to Julia, but is there a way to extract that value from the density? If you have any ideas, feel free to let me know!

# Import libraries.
using Turing, StatsPlots, Random

# Set the true probability of heads in a coin.
p_true = 0.5

# Iterate from having seen 0 observations to 100 observations.
Ns = 0:100

# Draw data from a Bernoulli distribution, i.e. draw heads or tails.
Random.seed!(12)
data = rand(Bernoulli(p_true), last(Ns))

# Declare our Turing model.
@model function coinflip(y)
    # Our prior belief about the probability of heads in a coin.
    p ~ Beta(1, 1)

    # The number of observations.
    N = length(y)
    for n in 1:N
        # Heads or tails of a coin are drawn from a Bernoulli distribution.
        y[n] ~ Bernoulli(p)
    end
end

# Settings of the Hamiltonian Monte Carlo (HMC) sampler.
iterations = 1000
ϵ = 0.05
τ = 10

# Start sampling.
chain = sample(coinflip(data), HMC(ϵ, τ), iterations)

# Plot a summary of the sampling process for the parameter p, i.e. the probability of heads in a coin.
histogram(chain[:p])

#### NOW WE SAMPLE WITH A SECOND PRIOR

using KernelDensityEstimate, LinearAlgebra

chain_matrix = convert(Array, get(chain, :p).p)
p_kde = kde!(transpose!(zeros(1, iterations), chain_matrix))

struct PosteriorDist <: ContinuousUnivariateDistribution
end

# Need to include rand for the prior
Distributions.rand(d::PosteriorDist, n::Int) = rand(d, n)

# Missing logpdf!

# Declare our Turing model.
@model function coinflip_newprior(y)
    # Our prior belief about the probability of heads in a coin.
    p_2 ~ PosteriorDist()

    # The number of observations.
    N = length(y)
    for n in 1:N
        # Heads or tails of a coin are drawn from a Bernoulli distribution.
        y[n] ~ Bernoulli(p)
    end
end

# Settings of the Hamiltonian Monte Carlo (HMC) sampler.
iterations = 1000
ϵ = 0.05
τ = 10

# Start sampling.
chain_2 = sample(coinflip_newprior(data), HMC(ϵ, τ), iterations)

# Plot a summary of the sampling process for the parameter p, i.e. the probability of heads in a coin.
histogram(chain_2[:p_2])

Peter

1 Like