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?
Try wrapping the multivariate kernel density estimator from GitHub - JuliaRobotics/KernelDensityEstimate.jl: Kernel Density Estimate with product approximation using multiscale Gibbs sampling in a custom distribution with
logpdf, etc. (follow the Turing docs for custom distributions). If you only need univariate kernel density estimators, there is also GitHub - JuliaStats/KernelDensity.jl: Kernel density estimators for Julia. 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.
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])