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 `rand`

, `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])
```

Peter