Cannot evaluate a vector for some ProductDistribution

Specific Domains Probabilistic programming
1

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[1])
        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[1], 0.) + logpdf(d1.marginals[2], 0.5) + logpdf(d1.marginals[3], 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.

2

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[1])
        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