How to fit a normal approximation to data in Julia

Specific Domains Statistics
5

Or the third one (which is, frankly, probably overkill here):

f = [12, 199, 3125, 29616, 175426, 662401, 1608671, 2477850, 2441350,
     1580290, 709625, 235426, 60899, 12643, 2126, 296, 41, 4, 0, 0]
x = 0:(length(counts)-1)

using DynamicHMC, TransformVariables, Distributions, UnPack, LogDensityProblems
import ForwardDiff, Random

struct NormalModel{TF,TX}
    f::TF
    x::TX
end

function (model::NormalModel)(θ)
    @unpack μ, σ = θ
    @unpack f, x = model
    D = Normal(μ, σ)
    sum(f .* logpdf.(D, x))
end

p = NormalModel(f, x)
t = as((μ = asℝ, σ = asℝ₊))
ℓ = TransformedLogDensity(t, p)
∇ℓ = ADgradient(:ForwardDiff, ℓ)
results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇ℓ, 1000)
posterior = t.(results.chain)
quantile(map(p -> p.μ, posterior), [0.25, 0.75]) # μ
3 Likes