# Cannot evaluate a vector for some ProductDistribution

Dear all,

I wish to specify a scaled inverse Wishart as formulated in this pdf, bottom of page 3. I understand how to implement it and have the following code so far (focused on the distributions). I’m focused on the case where my variance covariance matrix is only a diagonal matrix, whose elements are from a scaled inverse chi squared distribution. My MWE has the inverse chi squared:

``````using Distributions
using BayesianTools

function inverse_chi_squared(ν::Int)
""" v = the degrees of freedom parameter"""

rv        = 1 ./ rand(Chisq(ν), 1_000_000)
rv_cdf    = ecdf(rv)
obs       = unique(rv)
# Percentiles
cdf_obs = map(x -> rv_cdf(x), sort(obs))

# Generate probability mass function for next part
pmf = Float64[]
for i in 1:length(obs)
if i == 1
append!(pmf,cdf_obs)
else
append!(pmf, cdf_obs[i]-cdf_obs[i-1])
end
end

# Generate the distribution object
d = DiscreteNonParametric(sort(obs), pmf)
return d
end

d(n) = ProductDistribution(LogNormal(0, 1), InverseChiSquared(n + 1))
``````

Objective: To evaluate the log-pdf of this ProductDistribution given a vector in the support of this ProductDistribution.

Issue: I can only evaluate for each marginal separately, but not the distribution itself i.e.,

``````d1 = ProductDistribution(Normal(0,1), Uniform(0,1), InverseGamma(0.01, 0.01))
insupport(d1, [0., 1., 1.])

logpdf(d1, [0., 0.5, 0.5]) .== logpdf(d1.marginals, 0.) + logpdf(d1.marginals, 0.5) + logpdf(d1.marginals, 0.5)
``````

EDIT: I just realized I can reparameterize this into using an Inverse Gamma. I think I will try to just implement the entire distribution since it could be helpful going forward and it’s simply the product of two distributions.

So, I implemented the following workaround:

``````# So first we draw from the product distribution
function LogNormal_SqrtScaledInverseChiSquared(ν, tau²)
dgf(n)   = (n + 1) / 2
pd       = ProductDistribution(LogNormal(0, 1), InverseGamma(dgf(ν), dgf(ν) * tau²))
rv       = Random.rand!(pd, zeros(2, 1000000))
rv[2, :] = sqrt.(rv[2, :])
rv       = vec(prod(rv, dims=1))

rv_cdf = ecdf(rv)
obs    = unique(rv)
# Percentiles
cdf_obs = map(x -> rv_cdf(x), sort(obs))

# Generate probability mass function for next part
pmf = Float64[]
for i in 1:length(obs)
if i == 1
append!(pmf,cdf_obs)
else
append!(pmf, cdf_obs[i]-cdf_obs[i-1])
end
end
# Generate, using MLE, the distribution object
d      = DiscreteNonParametric(sort(obs), pmf)
sample = rand(d, 1000000)

# Make continuous
U = kde(sample)
return U
end
``````
1 Like