I have a simple example working with a univariate normal loglikelihood, but I am having trouble extending it to multivariate normal. I am implementing simulated moments, with the asymptotic Gaussian distribution used to form the likelihood function. Below are two code examples, the first of which works (univariate case), and the second of which is a simple modification of the first, adding a second statistic to the problem. The second example does not run, giving the error
ERROR: LoadError: MethodError: no method matching loglikelihood(::MvNormal{Float64,PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}}, ::Array{Array{ForwardDiff.Dual{...
which makes me think that the AD is not able to handle the MvNormal loglikelihood, but perhaps I’m wrong about that. Here are the examples.
Any pointers would be very welcome!
This code, which just uses the sample mean to form the univariate statistic, is working;
# simple example of method of simulated moments estimated by
# NUTS MCMC
# load packages
using TransformVariables, LogDensityProblems, DynamicHMC, MCMCDiagnostics,
Parameters, Statistics, Distributions, ForwardDiff
"""
DGP is``y ∼ μ + ϵ``, where ``ϵ ∼ N(0, 1)`` IID,
and the statistic is the sample mean.
Prior for μ is flat
"""
function dgp(μ)
y = randn(100) .+ μ
m = mean(y)
end
# generate data
μ_true = 3.0
m = dgp(μ_true)
# Define a structure for the problem
# Should hold the data and the parameters of prior distributions.
struct MSM_Problem{Tm <: Real}
"statistic"
m::Tm
end
# Make the type callable with the parameters *as a single argument*.
function (problem::MSM_Problem)(θ)
@unpack m = problem # extract the data
@unpack μ = θ # extract parameters (only one here)
S = 100
ms = 0.0
for s = 1:S
ms += dgp(μ)/S
end
loglikelihood(Normal(0, 1), [m .- ms])
end
# original problem, without transformation of parameters
p = MSM_Problem(m)
# define the transformation of parameters (in this case, an identity)
problem_transformation(p::MSM_Problem) = as((μ=as(Array,1),))
# Wrap the problem with the transformation
t = problem_transformation(p)
P = TransformedLogDensity(t, p)
# use AD (Flux) for the gradient
∇P = ADgradient(:ForwardDiff, P);
# Sample from the posterior. `chain` holds the chain (positions and
# diagnostic information), while the second returned value is the tuned sampler
# which would allow continuation of sampling.
chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000);
# We use the transformation to obtain the posterior from the chain.
posterior = transform.(Ref(t), get_position.(chain));
However, when I try to add a second statistic, I am not able to get the AD to work. The code:
# simple example of method of simulated moments estimated by
# NUTS MCMC
# load packages
using TransformVariables, LogDensityProblems, DynamicHMC, MCMCDiagnostics,
Parameters, Statistics, Distributions, ForwardDiff
"""
DGP is``y ∼ μ + ϵ``, where ``ϵ ∼ N(0, 1)`` IID,
and the statistic is the sample mean and sample std. dev. (which is useless in this case)
Prior for μ is flat
"""
function dgp(μ)
y = randn(100) .+ μ
m = [mean(y);std(y)]
end
# generate data
μ_true = 3.0
m = dgp(μ_true)
# Define a structure for the problem
# Should hold the data and the parameters of prior distributions.
struct MSM_Problem{Tm <: Vector}
"statistic"
m::Tm
end
# Make the type callable with the parameters *as a single argument*.
function (problem::MSM_Problem)(θ)
@unpack m = problem # extract the data
@unpack μ = θ # extract parameters (only one here)
S = 100
ms = zeros(2)
for s = 1:S
ms += dgp(μ)/S
end
loglikelihood(MvNormal(zeros(2), eye(2)), [m .- ms])
end
# original problem, without transformation of parameters
p = MSM_Problem(m)
# define the transformation of parameters (in this case, an identity)
problem_transformation(p::MSM_Problem) = as((μ=as(Array,1),))
# Wrap the problem with the transformation
t = problem_transformation(p)
P = TransformedLogDensity(t, p)
# use AD (Flux) for the gradient
∇P = ADgradient(:ForwardDiff, P);
# Sample from the posterior. `chain` holds the chain (positions and
# diagnostic information), while the second returned value is the tuned sampler
# which would allow continuation of sampling.
chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000);
Consider something like