using DocumenterQuarto
using Empirikos
using Distributions
using JuMP
using LaTeXStrings
using Plots
using StatsPlots
using Hypatia
using RandomEmpirikos.jl: Reproducing REBayes Examples
Reproducing the REBayes R vignette
The goal of this notebook is to showcase the functionality of the Empirikos.jl package by (partially) reproducing results shown in the REBayes vignette.
gr()Plots.GRBackend()Gaussian Mixture models
Needles and haystacks
Here we demonstrate the classic “needles in a haystack” scenario from the REBayes vignette. This example simulates data where most observations come from a standard normal distribution (\(\mathrm{N}(0,1)\)), while a small fraction (10%) come from a shifted normal (\(\mathrm{N}(2,1)\)).
zs = rand(MersenneTwister(100), Normal(), 1000) .+ [fill(0,900); fill(2,100)]
zs = sort(zs);In Empirikos.jl, we represent each observation as a StandardNormalSample object. This differs from REBayes where you would directly use a vector of observations as input to the GLmix function.
The StandardNormalSample type encapsulates the observation along with knowledge that it follows a standard normal likelihood (i.e., \(Z_i \mid \mu_i \sim \mathrm{N}(\mu_i, 1)\)). This type-based approach is a key design feature of Empirikos.jl:
- It allows the package to know the likelihood structure without additional arguments;
- It enables multiple dispatch to handle different likelihood models;
- It maintains a consistent interface across all empirical Bayes problems.
normal_Zs = StandardNormalSample.(zs)
first(normal_Zs, 3)3-element Vector{StandardNormalSample{Float64}}:
N(-3.1746277553967115; μ, σ=1.0)
N(-2.990089450782722; μ, σ=1.0)
N(-2.8559587291548225; μ, σ=1.0)We now fit the nonparametric maximum likelihood estimator (NPMLE) using the Mosek optimizer. In REBayes, this would be done with a single call to GLmix(y), while in Empirikos we use the more general fit function with NPMLE specification:
normal_npmle = fit(NPMLE(DiscretePriorClass(), Hypatia.Optimizer), normal_Zs)
iter p_obj d_obj | abs_gap x_feas z_feas | tau kap mu | dir_res prox step alpha
0 3.2373e-02 -2.8080e+02 | 2.30e+03 2.05e+02 2.33e+00 | 1.00e+00 1.00e+00 1.00e+00 |
1 -3.9922e-01 -2.8538e+02 | 2.19e+03 2.08e+02 2.37e+00 | 9.36e-01 1.02e+00 9.50e-01 | 2.9e-12 8.7e-02 co-a 5.00e-02
2 -9.4365e-01 -2.8751e+02 | 2.08e+03 2.09e+02 2.38e+00 | 8.84e-01 1.02e+00 9.03e-01 | 1.3e-12 1.6e-01 co-a 5.00e-02
3 -2.4184e+00 -2.8975e+02 | 1.87e+03 2.10e+02 2.39e+00 | 7.93e-01 1.02e+00 8.12e-01 | 4.0e-13 4.5e-01 co-a 1.00e-01
4 -4.6345e+00 -2.9294e+02 | 1.68e+03 2.10e+02 2.40e+00 | 7.11e-01 1.03e+00 7.31e-01 | 3.9e-12 8.2e-02 co-a 1.00e-01
5 -7.7639e+00 -2.9620e+02 | 1.51e+03 2.11e+02 2.40e+00 | 6.39e-01 1.03e+00 6.58e-01 | 5.4e-12 3.3e-02 co-a 1.00e-01
6 -1.2275e+01 -3.0066e+02 | 1.36e+03 2.11e+02 2.40e+00 | 5.75e-01 1.03e+00 5.92e-01 | 1.6e-12 3.8e-02 co-a 1.00e-01
7 -1.8987e+01 -3.0723e+02 | 1.23e+03 2.11e+02 2.40e+00 | 5.17e-01 1.03e+00 5.33e-01 | 1.2e-12 5.1e-02 co-a 1.00e-01
8 -2.9477e+01 -3.1752e+02 | 1.10e+03 2.11e+02 2.40e+00 | 4.66e-01 1.03e+00 4.80e-01 | 2.8e-12 5.2e-02 co-a 1.00e-01
9 -4.5655e+01 -3.3343e+02 | 9.94e+02 2.11e+02 2.40e+00 | 4.19e-01 1.03e+00 4.32e-01 | 1.4e-11 4.8e-02 co-a 1.00e-01
10 -6.8694e+01 -3.5614e+02 | 8.94e+02 2.11e+02 2.40e+00 | 3.77e-01 1.03e+00 3.89e-01 | 2.8e-12 2.9e-02 co-a 1.00e-01
11 -9.9142e+01 -3.8621e+02 | 8.05e+02 2.11e+02 2.40e+00 | 3.40e-01 1.03e+00 3.50e-01 | 1.3e-11 2.6e-02 co-a 1.00e-01
12 -1.3729e+02 -4.2390e+02 | 7.24e+02 2.11e+02 2.40e+00 | 3.06e-01 1.03e+00 3.15e-01 | 1.2e-12 2.1e-02 co-a 1.00e-01
13 -3.2935e+02 -6.1386e+02 | 5.07e+02 2.10e+02 2.39e+00 | 2.15e-01 1.03e+00 2.20e-01 | 2.5e-12 9.5e-01 co-a 3.00e-01
14 -6.3843e+02 -9.1902e+02 | 3.55e+02 2.09e+02 2.38e+00 | 1.51e-01 1.02e+00 1.54e-01 | 3.2e-12 2.4e-01 co-a 3.00e-01
15 -7.5222e+02 -1.0296e+03 | 3.19e+02 2.07e+02 2.36e+00 | 1.37e-01 1.01e+00 1.39e-01 | 2.1e-11 4.5e-02 co-a 1.00e-01
16 -8.7220e+02 -1.1447e+03 | 2.88e+02 2.04e+02 2.32e+00 | 1.26e-01 9.96e-01 1.25e-01 | 2.7e-11 7.5e-02 co-a 1.00e-01
17 -9.6568e+02 -1.2242e+03 | 2.59e+02 1.93e+02 2.20e+00 | 1.19e-01 9.47e-01 1.12e-01 | 2.1e-11 2.9e-01 co-a 1.00e-01
18 -7.8195e+02 -9.7089e+02 | 2.33e+02 1.40e+02 1.60e+00 | 1.48e-01 6.35e-01 1.01e-01 | 4.8e-11 1.2e-01 co-a 1.00e-01
19 -6.0313e+02 -7.3789e+02 | 2.10e+02 9.97e+01 1.14e+00 | 1.87e-01 4.87e-01 9.11e-02 | 4.2e-11 7.9e-01 co-a 1.00e-01
20 -4.9603e+02 -5.9886e+02 | 1.89e+02 7.59e+01 8.64e-01 | 2.21e-01 3.73e-01 8.20e-02 | 3.1e-11 4.2e-01 co-a 1.00e-01
21 -4.2283e+02 -5.0496e+02 | 1.70e+02 6.05e+01 6.89e-01 | 2.49e-01 2.99e-01 7.38e-02 | 4.1e-11 1.3e-01 co-a 1.00e-01
22 -2.7948e+02 -3.2455e+02 | 1.19e+02 3.31e+01 3.77e-01 | 3.19e-01 1.81e-01 5.17e-02 | 2.4e-11 1.2e-01 co-a 3.00e-01
23 -2.0630e+02 -2.3488e+02 | 8.33e+01 2.10e+01 2.39e-01 | 3.53e-01 1.05e-01 3.62e-02 | 2.2e-11 6.6e-01 co-a 3.00e-01
24 -2.0088e+02 -2.2838e+02 | 7.91e+01 2.02e+01 2.30e-01 | 3.48e-01 9.87e-02 3.44e-02 | 6.2e-11 5.7e-01 co-a 5.00e-02
25 -1.8483e+02 -2.0908e+02 | 7.12e+01 1.78e+01 2.03e-01 | 3.56e-01 8.69e-02 3.09e-02 | 1.4e-11 3.8e-01 co-a 1.00e-01
26 -1.7155e+02 -1.9323e+02 | 6.40e+01 1.59e+01 1.81e-01 | 3.58e-01 7.77e-02 2.78e-02 | 1.5e-11 3.5e-01 co-a 1.00e-01
27 -1.5902e+02 -1.7835e+02 | 5.76e+01 1.42e+01 1.61e-01 | 3.61e-01 6.93e-02 2.50e-02 | 1.0e-11 2.5e-01 co-a 1.00e-01
28 -1.4812e+02 -1.6550e+02 | 5.18e+01 1.27e+01 1.45e-01 | 3.62e-01 6.23e-02 2.25e-02 | 8.2e-12 2.2e-01 co-a 1.00e-01
29 -1.3792e+02 -1.5353e+02 | 4.66e+01 1.14e+01 1.30e-01 | 3.62e-01 5.59e-02 2.03e-02 | 1.3e-11 1.7e-01 co-a 1.00e-01
30 -1.2878e+02 -1.4288e+02 | 4.20e+01 1.03e+01 1.18e-01 | 3.61e-01 5.05e-02 1.82e-02 | 1.0e-11 1.5e-01 co-a 1.00e-01
31 -1.2024e+02 -1.3297e+02 | 3.78e+01 9.33e+00 1.06e-01 | 3.60e-01 4.56e-02 1.64e-02 | 9.2e-12 1.3e-01 co-a 1.00e-01
32 -1.1245e+02 -1.2397e+02 | 3.40e+01 8.45e+00 9.62e-02 | 3.58e-01 4.13e-02 1.48e-02 | 5.9e-12 1.2e-01 co-a 1.00e-01
33 -1.0516e+02 -1.1560e+02 | 3.06e+01 7.65e+00 8.72e-02 | 3.56e-01 3.74e-02 1.33e-02 | 6.0e-12 1.1e-01 co-a 1.00e-01
34 -9.8435e+01 -1.0790e+02 | 2.75e+01 6.94e+00 7.91e-02 | 3.53e-01 3.39e-02 1.20e-02 | 7.2e-12 1.0e-01 co-a 1.00e-01
35 -9.2138e+01 -1.0073e+02 | 2.48e+01 6.30e+00 7.18e-02 | 3.50e-01 3.08e-02 1.08e-02 | 9.4e-12 9.9e-02 co-a 1.00e-01
36 -8.6285e+01 -9.4088e+01 | 2.23e+01 5.72e+00 6.52e-02 | 3.47e-01 2.79e-02 9.68e-03 | 5.4e-12 9.8e-02 co-a 1.00e-01
37 -8.0800e+01 -8.7890e+01 | 2.01e+01 5.20e+00 5.92e-02 | 3.43e-01 2.54e-02 8.72e-03 | 6.4e-12 9.7e-02 co-a 1.00e-01
38 -7.5680e+01 -8.2125e+01 | 1.80e+01 4.73e+00 5.39e-02 | 3.40e-01 2.31e-02 7.84e-03 | 5.1e-12 9.7e-02 co-a 1.00e-01
39 -7.0878e+01 -7.6738e+01 | 1.62e+01 4.30e+00 4.90e-02 | 3.36e-01 2.10e-02 7.06e-03 | 5.8e-12 9.7e-02 co-a 1.00e-01
40 -6.6385e+01 -7.1714e+01 | 1.46e+01 3.91e+00 4.45e-02 | 3.33e-01 1.91e-02 6.35e-03 | 5.7e-12 9.7e-02 co-a 1.00e-01
41 -6.2169e+01 -6.7015e+01 | 1.32e+01 3.56e+00 4.05e-02 | 3.29e-01 1.74e-02 5.72e-03 | 5.4e-12 9.7e-02 co-a 1.00e-01
42 -5.8219e+01 -6.2626e+01 | 1.18e+01 3.23e+00 3.68e-02 | 3.26e-01 1.58e-02 5.14e-03 | 2.9e-12 9.8e-02 co-a 1.00e-01
43 -5.4513e+01 -5.8519e+01 | 1.07e+01 2.94e+00 3.35e-02 | 3.23e-01 1.43e-02 4.63e-03 | 3.1e-12 9.8e-02 co-a 1.00e-01
44 -5.1037e+01 -5.4679e+01 | 9.59e+00 2.67e+00 3.05e-02 | 3.19e-01 1.30e-02 4.17e-03 | 2.9e-12 9.9e-02 co-a 1.00e-01
45 -4.7777e+01 -5.1087e+01 | 8.63e+00 2.43e+00 2.77e-02 | 3.16e-01 1.19e-02 3.75e-03 | 2.6e-12 9.9e-02 co-a 1.00e-01
46 -4.4720e+01 -4.7727e+01 | 7.76e+00 2.21e+00 2.52e-02 | 3.13e-01 1.08e-02 3.37e-03 | 2.2e-12 1.0e-01 co-a 1.00e-01
47 -4.1854e+01 -4.4584e+01 | 6.99e+00 2.01e+00 2.28e-02 | 3.10e-01 9.78e-03 3.04e-03 | 3.3e-12 1.0e-01 co-a 1.00e-01
48 -3.9167e+01 -4.1646e+01 | 6.29e+00 1.82e+00 2.07e-02 | 3.08e-01 8.87e-03 2.73e-03 | 2.5e-12 1.0e-01 co-a 1.00e-01
49 -3.6650e+01 -3.8898e+01 | 5.66e+00 1.65e+00 1.88e-02 | 3.05e-01 8.05e-03 2.46e-03 | 1.6e-12 1.0e-01 co-a 1.00e-01
50 -3.4293e+01 -3.6331e+01 | 5.09e+00 1.50e+00 1.71e-02 | 3.03e-01 7.30e-03 2.21e-03 | 2.6e-12 1.0e-01 co-a 1.00e-01
51 -3.2087e+01 -3.3933e+01 | 4.58e+00 1.36e+00 1.54e-02 | 3.01e-01 6.61e-03 1.99e-03 | 1.4e-12 1.0e-01 co-a 1.00e-01
52 -3.0023e+01 -3.1694e+01 | 4.12e+00 1.23e+00 1.40e-02 | 2.99e-01 5.98e-03 1.79e-03 | 1.9e-12 1.0e-01 co-a 1.00e-01
53 -2.8094e+01 -2.9606e+01 | 3.71e+00 1.11e+00 1.26e-02 | 2.98e-01 5.41e-03 1.61e-03 | 1.3e-12 1.0e-01 co-a 1.00e-01
54 -2.6293e+01 -2.7658e+01 | 3.34e+00 1.00e+00 1.14e-02 | 2.97e-01 4.89e-03 1.45e-03 | 1.1e-12 1.0e-01 co-a 1.00e-01
55 -2.4612e+01 -2.5844e+01 | 3.00e+00 9.06e-01 1.03e-02 | 2.96e-01 4.41e-03 1.31e-03 | 1.2e-12 1.0e-01 co-a 1.00e-01
56 -2.3044e+01 -2.4155e+01 | 2.70e+00 8.17e-01 9.30e-03 | 2.95e-01 3.98e-03 1.18e-03 | 1.2e-12 1.0e-01 co-a 1.00e-01
57 -2.1583e+01 -2.2584e+01 | 2.43e+00 7.36e-01 8.38e-03 | 2.95e-01 3.58e-03 1.06e-03 | 9.9e-13 1.0e-01 co-a 1.00e-01
58 -2.0223e+01 -2.1123e+01 | 2.19e+00 6.62e-01 7.54e-03 | 2.95e-01 3.22e-03 9.52e-04 | 9.8e-13 1.1e-01 co-a 1.00e-01
59 -1.8957e+01 -1.9765e+01 | 1.97e+00 5.94e-01 6.77e-03 | 2.96e-01 2.89e-03 8.56e-04 | 1.1e-12 1.1e-01 co-a 1.00e-01
60 -1.7779e+01 -1.8503e+01 | 1.77e+00 5.33e-01 6.07e-03 | 2.97e-01 2.59e-03 7.71e-04 | 7.5e-13 1.1e-01 co-a 1.00e-01
61 -1.6684e+01 -1.7333e+01 | 1.60e+00 4.77e-01 5.43e-03 | 2.99e-01 2.32e-03 6.93e-04 | 7.7e-13 1.1e-01 co-a 1.00e-01
62 -1.5671e+01 -1.6251e+01 | 1.44e+00 4.26e-01 4.85e-03 | 3.01e-01 2.07e-03 6.24e-04 | 1.0e-12 1.0e-01 co-a 1.00e-01
63 -1.4738e+01 -1.5255e+01 | 1.29e+00 3.80e-01 4.33e-03 | 3.04e-01 1.85e-03 5.62e-04 | 2.0e-12 9.5e-02 co-a 1.00e-01
64 -1.3883e+01 -1.4344e+01 | 1.16e+00 3.38e-01 3.85e-03 | 3.07e-01 1.65e-03 5.05e-04 | 1.1e-12 8.7e-02 co-a 1.00e-01
65 -1.3103e+01 -1.3513e+01 | 1.05e+00 3.01e-01 3.43e-03 | 3.10e-01 1.47e-03 4.55e-04 | 3.1e-12 7.8e-02 co-a 1.00e-01
66 -1.2395e+01 -1.2759e+01 | 9.42e-01 2.68e-01 3.05e-03 | 3.14e-01 1.30e-03 4.09e-04 | 1.6e-12 6.9e-02 co-a 1.00e-01
67 -1.1753e+01 -1.2077e+01 | 8.48e-01 2.38e-01 2.71e-03 | 3.18e-01 1.16e-03 3.68e-04 | 9.5e-13 6.1e-02 co-a 1.00e-01
68 -1.0024e+01 -1.0242e+01 | 5.93e-01 1.60e-01 1.82e-03 | 3.31e-01 7.76e-04 2.58e-04 | 1.5e-12 9.4e-01 co-a 3.00e-01
69 -9.9984e+00 -1.0215e+01 | 5.87e-01 1.59e-01 1.81e-03 | 3.30e-01 7.73e-04 2.55e-04 | 2.9e-10 8.7e-01 co-a 1.00e-02
70 -9.9475e+00 -1.0161e+01 | 5.81e-01 1.57e-01 1.79e-03 | 3.31e-01 7.62e-04 2.52e-04 | 3.0e-11 5.6e-01 co-a 1.00e-02
71 -9.5661e+00 -9.7569e+00 | 5.23e-01 1.40e-01 1.60e-03 | 3.34e-01 6.81e-04 2.27e-04 | 6.4e-12 8.7e-01 co-a 1.00e-01
72 -8.5075e+00 -8.6359e+00 | 3.65e-01 9.42e-02 1.07e-03 | 3.47e-01 4.56e-04 1.59e-04 | 2.1e-12 9.1e-01 co-a 3.00e-01
73 -8.4833e+00 -8.6103e+00 | 3.61e-01 9.32e-02 1.06e-03 | 3.48e-01 4.52e-04 1.57e-04 | 1.5e-10 7.3e-01 co-a 1.00e-02
74 -8.3679e+00 -8.4882e+00 | 3.43e-01 8.83e-02 1.01e-03 | 3.49e-01 4.28e-04 1.49e-04 | 1.2e-11 7.4e-01 co-a 5.00e-02
75 -7.7036e+00 -7.7856e+00 | 2.40e-01 6.01e-02 6.85e-04 | 3.58e-01 2.91e-04 1.04e-04 | 3.6e-12 9.6e-01 co-a 3.00e-01
76 -7.6274e+00 -7.7050e+00 | 2.27e-01 5.69e-02 6.48e-04 | 3.60e-01 2.75e-04 9.88e-05 | 2.0e-11 7.4e-01 co-a 5.00e-02
77 -7.4831e+00 -7.5524e+00 | 2.05e-01 5.08e-02 5.79e-04 | 3.63e-01 2.45e-04 8.89e-05 | 1.8e-12 4.2e-01 co-a 1.00e-01
78 -7.1040e+00 -7.1517e+00 | 1.43e-01 3.49e-02 3.98e-04 | 3.69e-01 1.68e-04 6.22e-05 | 4.5e-12 6.4e-01 co-a 3.00e-01
79 -7.0185e+00 -7.0613e+00 | 1.29e-01 3.14e-02 3.57e-04 | 3.70e-01 1.51e-04 5.59e-05 | 7.2e-12 9.4e-01 co-a 1.00e-01
80 -6.7875e+00 -6.8171e+00 | 8.99e-02 2.17e-02 2.47e-04 | 3.74e-01 1.05e-04 3.91e-05 | 4.0e-12 9.2e-01 co-a 3.00e-01
81 -6.7823e+00 -6.8116e+00 | 8.90e-02 2.15e-02 2.45e-04 | 3.75e-01 1.03e-04 3.87e-05 | 1.8e-10 7.5e-01 co-a 1.00e-02
82 -6.7563e+00 -6.7841e+00 | 8.45e-02 2.04e-02 2.32e-04 | 3.75e-01 9.81e-05 3.67e-05 | 2.0e-11 6.5e-01 co-a 5.00e-02
83 -6.6072e+00 -6.6265e+00 | 5.91e-02 1.42e-02 1.61e-04 | 3.77e-01 6.81e-05 2.57e-05 | 5.9e-12 9.2e-01 co-a 3.00e-01
84 -6.5902e+00 -6.6085e+00 | 5.61e-02 1.35e-02 1.53e-04 | 3.77e-01 6.47e-05 2.44e-05 | 5.5e-11 9.5e-01 co-a 5.00e-02
85 -6.5575e+00 -6.5740e+00 | 5.05e-02 1.21e-02 1.38e-04 | 3.78e-01 5.80e-05 2.19e-05 | 1.9e-11 7.2e-01 co-a 1.00e-01
86 -6.5284e+00 -6.5432e+00 | 4.54e-02 1.09e-02 1.24e-04 | 3.79e-01 5.21e-05 1.97e-05 | 3.9e-12 3.6e-01 co-a 1.00e-01
87 -6.4494e+00 -6.4598e+00 | 3.18e-02 7.57e-03 8.62e-05 | 3.80e-01 3.63e-05 1.38e-05 | 3.8e-12 7.8e-01 co-a 3.00e-01
88 -6.4403e+00 -6.4501e+00 | 3.02e-02 7.19e-03 8.19e-05 | 3.80e-01 3.45e-05 1.31e-05 | 2.8e-11 6.9e-01 co-a 5.00e-02
89 -6.4229e+00 -6.4317e+00 | 2.72e-02 6.47e-03 7.37e-05 | 3.81e-01 3.10e-05 1.18e-05 | 4.8e-12 4.6e-01 co-a 1.00e-01
90 -6.3760e+00 -6.3822e+00 | 1.90e-02 4.52e-03 5.14e-05 | 3.82e-01 2.16e-05 8.26e-06 | 2.5e-12 8.8e-01 co-a 3.00e-01
91 -6.3705e+00 -6.3764e+00 | 1.80e-02 4.29e-03 4.89e-05 | 3.82e-01 2.06e-05 7.84e-06 | 1.5e-10 9.0e-01 co-a 5.00e-02
92 -6.3601e+00 -6.3654e+00 | 1.62e-02 3.86e-03 4.40e-05 | 3.82e-01 1.85e-05 7.05e-06 | 1.8e-11 7.8e-01 co-a 1.00e-01
93 -6.3509e+00 -6.3556e+00 | 1.46e-02 3.47e-03 3.95e-05 | 3.82e-01 1.66e-05 6.35e-06 | 5.3e-12 4.2e-01 co-a 1.00e-01
94 -6.3256e+00 -6.3290e+00 | 1.02e-02 2.43e-03 2.76e-05 | 3.83e-01 1.16e-05 4.44e-06 | 3.8e-12 8.2e-01 co-a 3.00e-01
95 -6.3197e+00 -6.3227e+00 | 9.19e-03 2.18e-03 2.49e-05 | 3.83e-01 1.04e-05 3.99e-06 | 2.5e-11 9.5e-01 co-a 1.00e-01
96 -6.3145e+00 -6.3172e+00 | 8.27e-03 1.96e-03 2.24e-05 | 3.83e-01 9.39e-06 3.59e-06 | 6.2e-12 5.6e-01 co-a 1.00e-01
97 -6.3002e+00 -6.3021e+00 | 5.78e-03 1.37e-03 1.57e-05 | 3.83e-01 6.57e-06 2.51e-06 | 5.2e-12 9.4e-01 co-a 3.00e-01
98 -6.2986e+00 -6.3003e+00 | 5.49e-03 1.31e-03 1.49e-05 | 3.83e-01 6.24e-06 2.39e-06 | 5.8e-11 8.5e-01 co-a 5.00e-02
99 -6.2954e+00 -6.2970e+00 | 4.94e-03 1.18e-03 1.34e-05 | 3.83e-01 5.61e-06 2.15e-06 | 9.9e-12 7.0e-01 co-a 1.00e-01
100 -6.2926e+00 -6.2940e+00 | 4.45e-03 1.06e-03 1.20e-05 | 3.83e-01 5.05e-06 1.93e-06 | 8.1e-12 4.0e-01 co-a 1.00e-01
101 -6.2849e+00 -6.2859e+00 | 3.11e-03 7.40e-04 8.43e-06 | 3.83e-01 3.53e-06 1.35e-06 | 5.2e-12 7.4e-01 co-a 3.00e-01
102 -6.2831e+00 -6.2840e+00 | 2.80e-03 6.66e-04 7.59e-06 | 3.83e-01 3.18e-06 1.22e-06 | 1.7e-11 6.9e-01 co-a 1.00e-01
103 -6.2783e+00 -6.2789e+00 | 1.96e-03 4.66e-04 5.31e-06 | 3.83e-01 2.22e-06 8.51e-07 | 4.8e-12 9.8e-01 co-a 3.00e-01
104 -6.2777e+00 -6.2783e+00 | 1.86e-03 4.43e-04 5.05e-06 | 3.83e-01 2.11e-06 8.08e-07 | 1.3e-10 9.2e-01 co-a 5.00e-02
105 -6.2766e+00 -6.2772e+00 | 1.67e-03 3.99e-04 4.54e-06 | 3.83e-01 1.90e-06 7.27e-07 | 9.7e-11 8.8e-01 co-a 1.00e-01
106 -6.2757e+00 -6.2762e+00 | 1.50e-03 3.59e-04 4.09e-06 | 3.83e-01 1.71e-06 6.54e-07 | 1.7e-11 5.2e-01 co-a 1.00e-01
107 -6.2731e+00 -6.2734e+00 | 1.05e-03 2.51e-04 2.86e-06 | 3.83e-01 1.20e-06 4.58e-07 | 7.9e-12 7.9e-01 co-a 3.00e-01
108 -6.2725e+00 -6.2728e+00 | 9.47e-04 2.26e-04 2.58e-06 | 3.83e-01 1.08e-06 4.12e-07 | 6.2e-11 6.9e-01 co-a 1.00e-01
109 -6.2708e+00 -6.2711e+00 | 6.63e-04 1.58e-04 1.80e-06 | 3.83e-01 7.53e-07 2.88e-07 | 2.5e-11 8.7e-01 co-a 3.00e-01
110 -6.2705e+00 -6.2707e+00 | 5.96e-04 1.43e-04 1.62e-06 | 3.83e-01 6.77e-07 2.59e-07 | 1.7e-10 8.6e-01 co-a 1.00e-01
111 -6.2701e+00 -6.2703e+00 | 5.36e-04 1.28e-04 1.46e-06 | 3.83e-01 6.10e-07 2.33e-07 | 6.1e-11 4.3e-01 co-a 1.00e-01
112 -6.2692e+00 -6.2693e+00 | 3.75e-04 8.99e-05 1.02e-06 | 3.82e-01 4.27e-07 1.63e-07 | 2.4e-11 6.0e-01 co-a 3.00e-01
113 -6.2690e+00 -6.2691e+00 | 3.38e-04 8.09e-05 9.22e-07 | 3.82e-01 3.84e-07 1.47e-07 | 4.1e-11 4.4e-01 co-a 1.00e-01
114 -6.2684e+00 -6.2685e+00 | 2.36e-04 5.66e-05 6.45e-07 | 3.82e-01 2.69e-07 1.03e-07 | 1.1e-11 5.5e-01 co-a 3.00e-01
115 -6.2683e+00 -6.2684e+00 | 2.13e-04 5.10e-05 5.81e-07 | 3.82e-01 2.42e-07 9.24e-08 | 1.2e-10 4.2e-01 co-a 1.00e-01
116 -6.2679e+00 -6.2680e+00 | 1.49e-04 3.57e-05 4.07e-07 | 3.82e-01 1.69e-07 6.47e-08 | 2.0e-11 5.1e-01 co-a 3.00e-01
117 -6.2677e+00 -6.2677e+00 | 1.04e-04 2.50e-05 2.85e-07 | 3.82e-01 1.19e-07 4.52e-08 | 2.0e-10 9.7e-01 co-a 3.00e-01
118 -6.2676e+00 -6.2677e+00 | 9.37e-05 2.25e-05 2.56e-07 | 3.82e-01 1.07e-07 4.07e-08 | 3.5e-10 5.0e-01 co-a 1.00e-01
119 -6.2675e+00 -6.2675e+00 | 6.55e-05 1.58e-05 1.80e-07 | 3.82e-01 7.47e-08 2.85e-08 | 7.2e-11 5.6e-01 co-a 3.00e-01
120 -6.2674e+00 -6.2674e+00 | 4.58e-05 1.10e-05 1.26e-07 | 3.81e-01 5.23e-08 1.99e-08 | 3.5e-10 8.6e-01 co-a 3.00e-01
121 -6.2673e+00 -6.2673e+00 | 3.21e-05 7.74e-06 8.81e-08 | 3.81e-01 3.66e-08 1.39e-08 | 6.7e-10 9.4e-01 co-a 3.00e-01
122 -6.2673e+00 -6.2673e+00 | 2.88e-05 6.97e-06 7.94e-08 | 3.81e-01 3.29e-08 1.25e-08 | 1.8e-09 5.4e-01 co-a 1.00e-01
123 -6.2672e+00 -6.2672e+00 | 2.02e-05 4.88e-06 5.56e-08 | 3.81e-01 2.30e-08 8.77e-09 | 3.8e-10 4.7e-01 co-a 3.00e-01
124 -6.2672e+00 -6.2672e+00 | 1.41e-05 3.42e-06 3.89e-08 | 3.81e-01 1.61e-08 6.14e-09 | 3.1e-10 5.0e-01 co-a 3.00e-01
125 -6.2672e+00 -6.2672e+00 | 9.88e-06 2.39e-06 2.72e-08 | 3.80e-01 1.13e-08 4.29e-09 | 2.8e-10 4.0e-01 co-a 3.00e-01
126 -6.2672e+00 -6.2672e+00 | 6.91e-06 1.67e-06 1.91e-08 | 3.80e-01 7.90e-09 3.00e-09 | 2.5e-10 3.3e-01 co-a 3.00e-01
127 -6.2672e+00 -6.2672e+00 | 4.84e-06 1.17e-06 1.34e-08 | 3.80e-01 5.53e-09 2.10e-09 | 2.5e-10 2.3e-01 co-a 3.00e-01
128 -6.2671e+00 -6.2671e+00 | 2.42e-06 5.86e-07 6.68e-09 | 3.80e-01 2.77e-09 1.05e-09 | 1.4e-10 8.9e-01 co-a 5.00e-01
129 -6.2671e+00 -6.2671e+00 | 2.18e-06 5.28e-07 6.01e-09 | 3.80e-01 2.49e-09 9.46e-10 | 2.1e-08 4.8e-01 co-a 1.00e-01
130 -6.2671e+00 -6.2671e+00 | 1.52e-06 3.70e-07 4.21e-09 | 3.80e-01 1.74e-09 6.62e-10 | 2.5e-10 3.3e-01 co-a 3.00e-01
131 -6.2671e+00 -6.2671e+00 | 1.37e-06 3.33e-07 3.79e-09 | 3.80e-01 1.57e-09 5.96e-10 | 3.3e-10 3.8e-02 co-a 1.00e-01
132 -6.2671e+00 -6.2671e+00 | 1.23e-06 2.99e-07 3.41e-09 | 3.80e-01 1.41e-09 5.36e-10 | 3.7e-10 6.3e-03 co-a 1.00e-01
133 -6.2671e+00 -6.2671e+00 | 6.17e-07 1.50e-07 1.71e-09 | 3.80e-01 7.07e-10 2.68e-10 | 1.6e-08 5.8e-01 co-a 5.00e-01
134 -6.2671e+00 -6.2671e+00 | 3.09e-07 7.50e-08 8.54e-10 | 3.80e-01 3.54e-10 1.34e-10 | 3.5e-10 5.9e-01 co-a 5.00e-01
135 -6.2671e+00 -6.2671e+00 | 2.78e-07 6.75e-08 7.68e-10 | 3.80e-01 3.18e-10 1.21e-10 | 3.4e-10 7.9e-02 co-a 1.00e-01
136 -6.2671e+00 -6.2671e+00 | 1.94e-07 4.72e-08 5.38e-10 | 3.79e-01 2.23e-10 8.45e-11 | 2.3e-09 1.9e-01 co-a 3.00e-01
137 -6.2671e+00 -6.2671e+00 | 1.75e-07 4.25e-08 4.84e-10 | 3.80e-01 2.00e-10 7.61e-11 | 9.5e-10 3.4e-02 co-a 1.00e-01
138 -6.2671e+00 -6.2671e+00 | 1.22e-07 2.98e-08 3.39e-10 | 3.79e-01 1.40e-10 5.32e-11 | 1.3e-09 1.6e-01 co-a 3.00e-01
139 -6.2671e+00 -6.2671e+00 | 1.21e-07 2.95e-08 3.36e-10 | 3.79e-01 1.39e-10 5.27e-11 | 7.8e-10 3.3e-03 co-a 1.00e-02
140 -6.2671e+00 -6.2671e+00 | 1.09e-07 2.65e-08 3.02e-10 | 3.79e-01 1.25e-10 4.74e-11 | 7.9e-11 5.0e-03 co-a 1.00e-01
141 -6.2671e+00 -6.2671e+00 | 9.82e-08 2.39e-08 2.72e-10 | 3.79e-01 1.12e-10 4.27e-11 | 2.3e-10 2.8e-03 co-a 1.00e-01
142 -6.2671e+00 -6.2671e+00 | 3.93e-08 9.55e-09 1.09e-10 | 3.79e-01 4.50e-11 1.71e-11 | 7.9e-12 5.2e-01 co-a 6.00e-01
optimal solution found; terminating
status is Optimal after 142 iterations and 20.802 seconds
Fitted NPMLE with Hypatia.Optimizer and 𝒢:
DiscretePriorClass | support = -3.1747277553967117:0.030043037229013896:5.808140376078443The NPMLE estimates the mixing distribution nonparametrically, which in this case should identify two mass points - one near 0 and another near 2.
Next, we define our posterior targets - in this case, the posterior means:
normal_postmean_targets = PosteriorMean.(normal_Zs)1000-element Vector{PosteriorMean{StandardNormalSample{Float64}}}:
𝔼[μ | N(-3.1746277553967115; μ, σ=1.0)]
𝔼[μ | N(-2.990089450782722; μ, σ=1.0)]
𝔼[μ | N(-2.8559587291548225; μ, σ=1.0)]
𝔼[μ | N(-2.7380023441265133; μ, σ=1.0)]
𝔼[μ | N(-2.724951899176289; μ, σ=1.0)]
𝔼[μ | N(-2.6950346914636967; μ, σ=1.0)]
𝔼[μ | N(-2.5584323721075397; μ, σ=1.0)]
𝔼[μ | N(-2.5305668218766053; μ, σ=1.0)]
𝔼[μ | N(-2.3632683788658877; μ, σ=1.0)]
𝔼[μ | N(-2.361941909453861; μ, σ=1.0)]
𝔼[μ | N(-2.355281185624943; μ, σ=1.0)]
𝔼[μ | N(-2.342143690964997; μ, σ=1.0)]
𝔼[μ | N(-2.284751655993758; μ, σ=1.0)]
⋮
𝔼[μ | N(3.1330895294747294; μ, σ=1.0)]
𝔼[μ | N(3.371564222667284; μ, σ=1.0)]
𝔼[μ | N(3.3919915435058083; μ, σ=1.0)]
𝔼[μ | N(3.4247099212446166; μ, σ=1.0)]
𝔼[μ | N(3.42836573942505; μ, σ=1.0)]
𝔼[μ | N(3.835965091861546; μ, σ=1.0)]
𝔼[μ | N(4.011422275963797; μ, σ=1.0)]
𝔼[μ | N(4.129379367933984; μ, σ=1.0)]
𝔼[μ | N(4.1511807079199565; μ, σ=1.0)]
𝔼[μ | N(4.46857344617616; μ, σ=1.0)]
𝔼[μ | N(4.538449778609994; μ, σ=1.0)]
𝔼[μ | N(5.808040376078443; μ, σ=1.0)]The following code creates a three-panel plot showing: 1. The true mixture density 2. The estimated prior distribution (mixing distribution) 3. The posterior mean function (the shrinkage estimator)
normal_pl_marginal = plot(marginalize(NormalSample(1.0),
DiscreteNonParametric([0.0;2.0], [0.9;0.1])),
components=false, label=nothing, xguide="z",
yguide = "f(z)",
title = "True mixture")
normal_pl_npmle_prior = plot(support(normal_npmle.prior),
probs(normal_npmle.prior), seriestype=:sticks,
xguide = L"\mu", yguide=L"g(\mu)",
title = "NPMLE prior estimate", label=nothing)
normal_pl_posterior_mean = plot(response.(location.(normal_postmean_targets)),
normal_postmean_targets.(normal_npmle.prior),
label=nothing, title="NPMLE posterior mean",
xguide="z", yguide=L"E[\mu \mid Z=z]")
plot(normal_pl_marginal, normal_pl_npmle_prior, normal_pl_posterior_mean,
size=(1000,350), layout=(1,3))This plot is similar to what you would see in REBayes but is generated using Empirikos.jl’s object-oriented approach. The posterior mean function (right panel) shows the shrinkage effect: observations are pulled toward the mass points of the estimated prior distribution. The shrinkage is stronger for values close to the mass points and weaker in regions between them.
Binomial NPMLE
Now we’ll analyze the Tacks dataset using a binomial mixture model. This classic dataset from Beckett and Diaconis involves outcomes from tossing thumbtacks.
tacks_tbl = Empirikos.Tacks.load_table()320-element CSV.File:
CSV.Row: (x = 7, k = 9)
CSV.Row: (x = 4, k = 9)
CSV.Row: (x = 6, k = 9)
CSV.Row: (x = 6, k = 9)
CSV.Row: (x = 6, k = 9)
CSV.Row: (x = 6, k = 9)
CSV.Row: (x = 8, k = 9)
CSV.Row: (x = 6, k = 9)
CSV.Row: (x = 5, k = 9)
CSV.Row: (x = 8, k = 9)
CSV.Row: (x = 6, k = 9)
CSV.Row: (x = 3, k = 9)
CSV.Row: (x = 3, k = 9)
⋮
CSV.Row: (x = 8, k = 9)
CSV.Row: (x = 7, k = 9)
CSV.Row: (x = 6, k = 9)
CSV.Row: (x = 6, k = 9)
CSV.Row: (x = 8, k = 9)
CSV.Row: (x = 6, k = 9)
CSV.Row: (x = 9, k = 9)
CSV.Row: (x = 6, k = 9)
CSV.Row: (x = 7, k = 9)
CSV.Row: (x = 6, k = 9)
CSV.Row: (x = 6, k = 9)
CSV.Row: (x = 6, k = 9)We convert the table to a vector of BinomialSample objects, where each sample represents the number of successes out of a known number of trials:
tack_Zs = BinomialSample.(tacks_tbl.x, Int64.(tacks_tbl.k))320-element Vector{BinomialSample{Int64, Int64}}:
ℬ𝒾𝓃(7; p, n=9)
ℬ𝒾𝓃(4; p, n=9)
ℬ𝒾𝓃(6; p, n=9)
ℬ𝒾𝓃(6; p, n=9)
ℬ𝒾𝓃(6; p, n=9)
ℬ𝒾𝓃(6; p, n=9)
ℬ𝒾𝓃(8; p, n=9)
ℬ𝒾𝓃(6; p, n=9)
ℬ𝒾𝓃(5; p, n=9)
ℬ𝒾𝓃(8; p, n=9)
ℬ𝒾𝓃(6; p, n=9)
ℬ𝒾𝓃(3; p, n=9)
ℬ𝒾𝓃(3; p, n=9)
⋮
ℬ𝒾𝓃(8; p, n=9)
ℬ𝒾𝓃(7; p, n=9)
ℬ𝒾𝓃(6; p, n=9)
ℬ𝒾𝓃(6; p, n=9)
ℬ𝒾𝓃(8; p, n=9)
ℬ𝒾𝓃(6; p, n=9)
ℬ𝒾𝓃(9; p, n=9)
ℬ𝒾𝓃(6; p, n=9)
ℬ𝒾𝓃(7; p, n=9)
ℬ𝒾𝓃(6; p, n=9)
ℬ𝒾𝓃(6; p, n=9)
ℬ𝒾𝓃(6; p, n=9)In REBayes, a similar analysis would use the Bmix function, but here we use our common fit interface with NPMLE. The information about the binomial likelihood is provided by the BinomialSample type instead, which makes the code more extensible and consistent across different likelihood models.
tacks_npmle = fit(NPMLE(DiscretePriorClass(), Hypatia.Optimizer), tack_Zs);
iter p_obj d_obj | abs_gap x_feas z_feas | tau kap mu | dir_res prox step alpha
0 3.8484e-01 -5.0613e+00 | 3.19e+02 1.46e+00 1.01e+00 | 1.00e+00 1.00e+00 1.00e+00 |
1 -4.4676e+00 -6.8414e+00 | 9.56e+01 8.72e-01 6.01e-01 | 5.03e-01 7.39e-01 3.00e-01 | 1.8e-15 6.2e-01 co-a 7.00e-01
2 -2.0140e+00 -2.8978e+00 | 4.81e+01 2.01e-01 1.39e-01 | 1.09e+00 4.21e-03 1.50e-01 | 3.7e-15 9.7e-01 co-a 5.00e-01
3 -1.6879e+00 -2.0839e+00 | 3.35e+01 9.87e-02 6.81e-02 | 1.56e+00 6.11e-02 1.05e-01 | 7.1e-15 3.0e-01 co-a 3.00e-01
4 -1.4505e+00 -1.6458e+00 | 2.34e+01 4.72e-02 3.26e-02 | 2.28e+00 2.92e-02 7.34e-02 | 8.2e-15 3.3e-01 co-a 3.00e-01
5 -1.2610e+00 -1.3175e+00 | 1.17e+01 1.33e-02 9.16e-03 | 4.05e+00 8.18e-03 3.66e-02 | 1.1e-13 9.2e-01 co-a 5.00e-01
6 -1.1868e+00 -1.2061e+00 | 5.84e+00 4.50e-03 3.11e-03 | 5.97e+00 3.38e-03 1.83e-02 | 9.5e-13 9.2e-01 co-a 5.00e-01
7 -1.1831e+00 -1.2004e+00 | 5.25e+00 4.03e-03 2.78e-03 | 6.00e+00 2.75e-03 1.64e-02 | 1.1e-11 5.5e-01 co-a 1.00e-01
8 -1.1697e+00 -1.1810e+00 | 3.67e+00 2.63e-03 1.81e-03 | 6.45e+00 1.78e-03 1.15e-02 | 2.5e-12 5.1e-01 co-a 3.00e-01
9 -1.1609e+00 -1.1685e+00 | 2.56e+00 1.76e-03 1.21e-03 | 6.73e+00 1.19e-03 8.03e-03 | 4.6e-12 4.1e-01 co-a 3.00e-01
10 -1.1547e+00 -1.1599e+00 | 1.79e+00 1.20e-03 8.28e-04 | 6.91e+00 8.11e-04 5.61e-03 | 3.6e-12 4.0e-01 co-a 3.00e-01
11 -1.1503e+00 -1.1539e+00 | 1.25e+00 8.22e-04 5.67e-04 | 7.07e+00 5.54e-04 3.92e-03 | 1.2e-12 3.1e-01 co-a 3.00e-01
12 -1.1451e+00 -1.1468e+00 | 6.24e-01 3.98e-04 2.75e-04 | 7.29e+00 2.66e-04 1.95e-03 | 1.7e-12 7.0e-01 co-a 5.00e-01
13 -1.1436e+00 -1.1448e+00 | 4.35e-01 2.74e-04 1.89e-04 | 7.42e+00 1.83e-04 1.36e-03 | 1.5e-11 8.2e-01 co-a 3.00e-01
14 -1.1426e+00 -1.1434e+00 | 3.03e-01 1.88e-04 1.30e-04 | 7.57e+00 1.26e-04 9.50e-04 | 8.8e-12 4.9e-01 co-a 3.00e-01
15 -1.1419e+00 -1.1425e+00 | 2.12e-01 1.29e-04 8.91e-05 | 7.71e+00 8.62e-05 6.63e-04 | 1.9e-11 4.8e-01 co-a 3.00e-01
16 -1.1415e+00 -1.1418e+00 | 1.48e-01 8.91e-05 6.15e-05 | 7.82e+00 5.93e-05 4.63e-04 | 2.1e-11 4.9e-01 co-a 3.00e-01
17 -1.1412e+00 -1.1414e+00 | 1.03e-01 6.15e-05 4.24e-05 | 7.93e+00 4.09e-05 3.23e-04 | 9.8e-11 5.8e-01 co-a 3.00e-01
18 -1.1410e+00 -1.1411e+00 | 7.20e-02 4.26e-05 2.94e-05 | 8.01e+00 2.83e-05 2.26e-04 | 1.3e-10 6.3e-01 co-a 3.00e-01
19 -1.1408e+00 -1.1410e+00 | 5.02e-02 2.96e-05 2.04e-05 | 8.07e+00 1.96e-05 1.58e-04 | 1.1e-10 7.1e-01 co-a 3.00e-01
20 -1.1407e+00 -1.1408e+00 | 3.50e-02 2.07e-05 1.43e-05 | 8.10e+00 1.36e-05 1.10e-04 | 1.1e-10 7.0e-01 co-a 3.00e-01
21 -1.1407e+00 -1.1408e+00 | 2.45e-02 1.44e-05 9.96e-06 | 8.12e+00 9.48e-06 7.67e-05 | 9.5e-11 6.5e-01 co-a 3.00e-01
22 -1.1406e+00 -1.1407e+00 | 1.22e-02 7.22e-06 4.98e-06 | 8.12e+00 4.73e-06 3.81e-05 | 1.8e-10 8.3e-01 co-a 5.00e-01
23 -1.1406e+00 -1.1407e+00 | 8.46e-03 5.08e-06 3.51e-06 | 8.07e+00 3.31e-06 2.65e-05 | 8.7e-10 7.9e-01 co-a 3.00e-01
24 -1.1406e+00 -1.1406e+00 | 3.35e-03 2.05e-06 1.41e-06 | 8.00e+00 1.33e-06 1.05e-05 | 5.1e-10 6.7e-01 co-a 6.00e-01
25 -1.1406e+00 -1.1406e+00 | 1.31e-03 8.41e-07 5.80e-07 | 7.80e+00 5.45e-07 4.10e-06 | 3.7e-09 8.5e-01 co-a 6.00e-01
26 -1.1406e+00 -1.1406e+00 | 6.42e-04 4.29e-07 2.96e-07 | 7.65e+00 2.69e-07 2.01e-06 | 2.5e-09 8.4e-01 co-a 5.00e-01
27 -1.1406e+00 -1.1406e+00 | 3.18e-04 2.16e-07 1.49e-07 | 7.58e+00 1.33e-07 9.97e-07 | 1.1e-09 7.2e-01 co-a 5.00e-01
28 -1.1406e+00 -1.1406e+00 | 1.27e-04 8.69e-08 6.00e-08 | 7.55e+00 5.29e-08 3.97e-07 | 7.9e-10 8.6e-01 co-a 6.00e-01
29 -1.1406e+00 -1.1406e+00 | 2.53e-05 1.74e-08 1.20e-08 | 7.54e+00 1.05e-08 7.93e-08 | 3.1e-10 8.2e-01 co-a 8.00e-01
30 -1.1406e+00 -1.1406e+00 | 1.01e-05 7.00e-09 4.82e-09 | 7.51e+00 4.23e-09 3.16e-08 | 8.2e-10 2.6e-01 co-a 6.00e-01
31 -1.1406e+00 -1.1406e+00 | 9.08e-06 6.29e-09 4.34e-09 | 7.51e+00 3.79e-09 2.85e-08 | 3.5e-09 9.3e-02 co-a 1.00e-01
32 -1.1406e+00 -1.1406e+00 | 8.62e-06 5.99e-09 4.12e-09 | 7.51e+00 3.60e-09 2.70e-08 | 1.9e-08 2.4e-02 co-a 5.00e-02
33 -1.1406e+00 -1.1406e+00 | 8.53e-06 5.93e-09 4.08e-09 | 7.51e+00 3.56e-09 2.68e-08 | 5.5e-10 1.0e-05 co-a 1.00e-02
34 -1.1406e+00 -1.1406e+00 | 5.97e-06 4.15e-09 2.86e-09 | 7.51e+00 2.49e-09 1.87e-08 | 3.3e-10 6.1e-04 co-a 3.00e-01
35 -1.1406e+00 -1.1406e+00 | 5.91e-06 4.11e-09 2.83e-09 | 7.51e+00 2.47e-09 1.85e-08 | 3.0e-10 2.2e-05 co-a 1.00e-02
36 -1.1406e+00 -1.1406e+00 | 4.14e-06 2.86e-09 1.98e-09 | 7.51e+00 1.73e-09 1.30e-08 | 3.0e-09 9.2e-03 co-a 3.00e-01
optimal solution found; terminating
status is Optimal after 36 iterations and 0.178 seconds
Let’s visualize the estimated prior distribution of success probabilities:
plot(support(tacks_npmle.prior),
probs(tacks_npmle.prior),
seriestype=:sticks,
xguide = "p",
yguide = "g(p)",
size=(500,400),
label=nothing,
title="NPMLE prior for Tacks data")The plot shows the estimated distribution of success probabilities for the tacks. The NPMLE identifies several mass points, indicating groups of tacks with different tendencies to land point-up.
Poisson NPMLE
Finally, we’ll analyze the Norberg insurance claims dataset, which consists of claim counts and exposure times for different occupational groups.
norberg_tbl = Empirikos.Norberg.load_table()72-element CSV.File:
CSV.Row: (OccGroup = 1, Exposure = 3349.02, Death = 16)
CSV.Row: (OccGroup = 2, Exposure = 1394.0, Death = 5)
CSV.Row: (OccGroup = 3, Exposure = 479.81, Death = 0)
CSV.Row: (OccGroup = 4, Exposure = 23.98, Death = 0)
CSV.Row: (OccGroup = 5, Exposure = 11.3, Death = 0)
CSV.Row: (OccGroup = 6, Exposure = 2273.55, Death = 11)
CSV.Row: (OccGroup = 7, Exposure = 179.85, Death = 0)
CSV.Row: (OccGroup = 8, Exposure = 3947.44, Death = 24)
CSV.Row: (OccGroup = 9, Exposure = 4332.56, Death = 9)
CSV.Row: (OccGroup = 10, Exposure = 109.84, Death = 1)
CSV.Row: (OccGroup = 11, Exposure = 647.77, Death = 2)
CSV.Row: (OccGroup = 12, Exposure = 1308.36, Death = 1)
CSV.Row: (OccGroup = 13, Exposure = 154.81, Death = 4)
⋮
CSV.Row: (OccGroup = 61, Exposure = 1551.3, Death = 3)
CSV.Row: (OccGroup = 62, Exposure = 113.49, Death = 0)
CSV.Row: (OccGroup = 63, Exposure = 147.44, Death = 0)
CSV.Row: (OccGroup = 64, Exposure = 826.4, Death = 2)
CSV.Row: (OccGroup = 65, Exposure = 775.19, Death = 1)
CSV.Row: (OccGroup = 66, Exposure = 634.16, Death = 2)
CSV.Row: (OccGroup = 67, Exposure = 2293.45, Death = 6)
CSV.Row: (OccGroup = 68, Exposure = 942.03, Death = 4)
CSV.Row: (OccGroup = 69, Exposure = 232.31, Death = 2)
CSV.Row: (OccGroup = 70, Exposure = 37.08, Death = 0)
CSV.Row: (OccGroup = 71, Exposure = 233.25, Death = 0)
CSV.Row: (OccGroup = 72, Exposure = 39.52, Death = 0)We convert the data to PoissonSample objects. Each sample captures both the count (death) and the exposure time:
norberg_Zs = Empirikos.PoissonSample.(norberg_tbl.Death, norberg_tbl.Exposure ./ 344);Let’s examine the range of claim rates:
extrema(response.(norberg_Zs) ./ nuisance_parameter.(norberg_Zs))(0.0, 8.888314708352173)And check how many groups had zero claims:
sum(response.(norberg_Zs) .== 0) # how many zeros in the samples?23Now we fit the NPMLE with a finer grid of 1000 points:
norberg_npmle = fit(NPMLE(DiscretePriorClass(),
Hypatia.Optimizer;
prior_grid_length=1000),
norberg_Zs);
iter p_obj d_obj | abs_gap x_feas z_feas | tau kap mu | dir_res prox step alpha
0 1.2878e-01 -2.3832e+01 | 4.45e+02 1.29e+01 6.83e-01 | 1.00e+00 1.00e+00 1.00e+00 |
1 -4.8042e+00 -2.8956e+01 | 3.11e+02 1.32e+01 7.03e-01 | 6.80e-01 1.05e+00 7.00e-01 | 1.4e-14 2.9e-01 co-a 3.00e-01
2 -3.7142e+01 -5.6908e+01 | 1.25e+02 1.19e+01 6.30e-01 | 3.04e-01 9.84e-01 2.80e-01 | 1.4e-14 9.7e-01 co-a 6.00e-01
3 -3.9425e+01 -5.3262e+01 | 8.74e+01 8.14e+00 4.33e-01 | 3.10e-01 6.10e-01 1.96e-01 | 1.5e-14 6.6e-01 co-a 3.00e-01
4 -2.3077e+01 -2.8299e+01 | 4.36e+01 2.95e+00 1.57e-01 | 4.27e-01 2.15e-01 9.80e-02 | 5.0e-14 6.9e-01 co-a 5.00e-01
5 -1.7909e+01 -2.1253e+01 | 3.05e+01 1.88e+00 1.00e-01 | 4.69e-01 1.44e-01 6.84e-02 | 7.9e-14 6.1e-01 co-a 3.00e-01
6 -1.4033e+01 -1.6182e+01 | 2.13e+01 1.20e+00 6.36e-02 | 5.16e-01 9.08e-02 4.78e-02 | 3.3e-14 4.8e-01 co-a 3.00e-01
7 -1.1359e+01 -1.2814e+01 | 1.48e+01 8.09e-01 4.30e-02 | 5.34e-01 6.18e-02 3.33e-02 | 5.7e-14 6.0e-01 co-a 3.00e-01
8 -9.1893e+00 -1.0156e+01 | 1.03e+01 5.35e-01 2.84e-02 | 5.66e-01 4.04e-02 2.33e-02 | 4.3e-14 6.6e-01 co-a 3.00e-01
9 -7.5456e+00 -8.1914e+00 | 7.21e+00 3.56e-01 1.89e-02 | 5.95e-01 2.67e-02 1.62e-02 | 8.7e-14 9.6e-01 co-a 3.00e-01
10 -6.9671e+00 -7.5117e+00 | 6.48e+00 2.99e-01 1.59e-02 | 6.38e-01 2.26e-02 1.46e-02 | 1.9e-13 5.1e-01 co-a 1.00e-01
11 -5.7868e+00 -6.1320e+00 | 4.52e+00 1.88e-01 9.99e-03 | 7.10e-01 1.39e-02 1.02e-02 | 6.8e-14 6.1e-01 co-a 3.00e-01
12 -4.9867e+00 -5.2060e+00 | 3.15e+00 1.19e-01 6.31e-03 | 7.87e-01 8.82e-03 7.07e-03 | 1.1e-13 7.9e-01 co-a 3.00e-01
13 -4.4583e+00 -4.5972e+00 | 2.19e+00 7.48e-02 3.98e-03 | 8.74e-01 5.55e-03 4.93e-03 | 2.2e-13 6.3e-01 co-a 3.00e-01
14 -4.1229e+00 -4.2119e+00 | 1.53e+00 4.77e-02 2.54e-03 | 9.59e-01 3.55e-03 3.43e-03 | 2.4e-13 5.1e-01 co-a 3.00e-01
15 -3.9121e+00 -3.9696e+00 | 1.07e+00 3.08e-02 1.64e-03 | 1.04e+00 2.29e-03 2.40e-03 | 4.8e-13 3.6e-01 co-a 3.00e-01
16 -3.7835e+00 -3.8213e+00 | 7.45e-01 2.02e-02 1.07e-03 | 1.11e+00 1.50e-03 1.67e-03 | 4.5e-13 3.5e-01 co-a 3.00e-01
17 -3.7039e+00 -3.7291e+00 | 5.20e-01 1.34e-02 7.12e-04 | 1.17e+00 9.99e-04 1.17e-03 | 5.0e-13 4.1e-01 co-a 3.00e-01
18 -3.6540e+00 -3.6708e+00 | 3.63e-01 8.98e-03 4.77e-04 | 1.22e+00 6.68e-04 8.17e-04 | 1.4e-12 5.0e-01 co-a 3.00e-01
19 -3.6222e+00 -3.6336e+00 | 2.54e-01 6.06e-03 3.22e-04 | 1.27e+00 4.50e-04 5.70e-04 | 2.4e-12 6.2e-01 co-a 3.00e-01
20 -3.6019e+00 -3.6096e+00 | 1.77e-01 4.11e-03 2.18e-04 | 1.31e+00 3.05e-04 3.98e-04 | 4.5e-12 7.9e-01 co-a 3.00e-01
21 -3.5980e+00 -3.6049e+00 | 1.59e-01 3.67e-03 1.95e-04 | 1.32e+00 2.72e-04 3.58e-04 | 8.6e-12 4.5e-01 co-a 1.00e-01
22 -3.5867e+00 -3.5915e+00 | 1.11e-01 2.51e-03 1.33e-04 | 1.35e+00 1.85e-04 2.50e-04 | 5.9e-12 5.7e-01 co-a 3.00e-01
23 -3.5793e+00 -3.5825e+00 | 7.76e-02 1.71e-03 9.07e-05 | 1.39e+00 1.26e-04 1.74e-04 | 6.9e-12 8.5e-01 co-a 3.00e-01
24 -3.5780e+00 -3.5809e+00 | 6.98e-02 1.53e-03 8.14e-05 | 1.40e+00 1.13e-04 1.57e-04 | 1.1e-11 3.7e-01 co-a 1.00e-01
25 -3.5740e+00 -3.5760e+00 | 4.88e-02 1.05e-03 5.58e-05 | 1.43e+00 7.70e-05 1.10e-04 | 4.4e-12 5.6e-01 co-a 3.00e-01
26 -3.5715e+00 -3.5728e+00 | 3.41e-02 7.20e-04 3.83e-05 | 1.45e+00 5.28e-05 7.65e-05 | 1.5e-11 8.4e-01 co-a 3.00e-01
27 -3.5710e+00 -3.5723e+00 | 3.06e-02 6.46e-04 3.43e-05 | 1.46e+00 4.72e-05 6.88e-05 | 2.4e-11 3.6e-01 co-a 1.00e-01
28 -3.5697e+00 -3.5705e+00 | 2.14e-02 4.45e-04 2.37e-05 | 1.48e+00 3.25e-05 4.81e-05 | 1.3e-11 6.0e-01 co-a 3.00e-01
29 -3.5687e+00 -3.5693e+00 | 1.50e-02 3.07e-04 1.63e-05 | 1.51e+00 2.24e-05 3.36e-05 | 6.2e-11 9.5e-01 co-a 3.00e-01
30 -3.5686e+00 -3.5691e+00 | 1.34e-02 2.76e-04 1.46e-05 | 1.51e+00 2.01e-05 3.02e-05 | 3.9e-11 4.2e-01 co-a 1.00e-01
31 -3.5681e+00 -3.5684e+00 | 9.39e-03 1.91e-04 1.01e-05 | 1.53e+00 1.39e-05 2.11e-05 | 3.7e-11 7.3e-01 co-a 3.00e-01
32 -3.5680e+00 -3.5683e+00 | 8.44e-03 1.71e-04 9.09e-06 | 1.53e+00 1.24e-05 1.90e-05 | 1.0e-10 3.8e-01 co-a 1.00e-01
33 -3.5677e+00 -3.5679e+00 | 5.90e-03 1.19e-04 6.31e-06 | 1.54e+00 8.60e-06 1.33e-05 | 2.5e-11 7.0e-01 co-a 3.00e-01
34 -3.5676e+00 -3.5678e+00 | 5.31e-03 1.07e-04 5.66e-06 | 1.55e+00 7.70e-06 1.19e-05 | 9.6e-11 3.9e-01 co-a 1.00e-01
35 -3.5674e+00 -3.5675e+00 | 3.71e-03 7.40e-05 3.93e-06 | 1.56e+00 5.35e-06 8.33e-06 | 3.1e-11 7.3e-01 co-a 3.00e-01
36 -3.5674e+00 -3.5675e+00 | 3.34e-03 6.65e-05 3.53e-06 | 1.56e+00 4.80e-06 7.50e-06 | 1.3e-10 4.1e-01 co-a 1.00e-01
37 -3.5673e+00 -3.5673e+00 | 2.33e-03 4.63e-05 2.46e-06 | 1.57e+00 3.34e-06 5.24e-06 | 7.0e-11 7.5e-01 co-a 3.00e-01
38 -3.5672e+00 -3.5673e+00 | 2.10e-03 4.16e-05 2.21e-06 | 1.57e+00 2.99e-06 4.71e-06 | 1.9e-10 4.2e-01 co-a 1.00e-01
39 -3.5672e+00 -3.5672e+00 | 1.46e-03 2.90e-05 1.54e-06 | 1.58e+00 2.09e-06 3.29e-06 | 6.9e-11 7.4e-01 co-a 3.00e-01
40 -3.5671e+00 -3.5672e+00 | 1.32e-03 2.61e-05 1.39e-06 | 1.58e+00 1.87e-06 2.96e-06 | 3.9e-10 4.0e-01 co-a 1.00e-01
41 -3.5671e+00 -3.5671e+00 | 9.20e-04 1.82e-05 9.67e-07 | 1.59e+00 1.31e-06 2.07e-06 | 1.2e-10 6.6e-01 co-a 3.00e-01
42 -3.5671e+00 -3.5671e+00 | 6.43e-04 1.27e-05 6.76e-07 | 1.59e+00 9.11e-07 1.45e-06 | 3.9e-10 8.5e-01 co-a 3.00e-01
43 -3.5671e+00 -3.5671e+00 | 5.78e-04 1.15e-05 6.08e-07 | 1.59e+00 8.19e-07 1.30e-06 | 2.2e-10 3.6e-01 co-a 1.00e-01
44 -3.5670e+00 -3.5671e+00 | 4.04e-04 8.00e-06 4.25e-07 | 1.59e+00 5.72e-07 9.09e-07 | 1.5e-10 4.4e-01 co-a 3.00e-01
45 -3.5670e+00 -3.5670e+00 | 2.83e-04 5.60e-06 2.97e-07 | 1.59e+00 3.99e-07 6.36e-07 | 1.1e-10 4.7e-01 co-a 3.00e-01
46 -3.5670e+00 -3.5670e+00 | 1.41e-04 2.80e-06 1.49e-07 | 1.59e+00 2.00e-07 3.17e-07 | 1.4e-10 9.9e-01 co-a 5.00e-01
47 -3.5670e+00 -3.5670e+00 | 9.86e-05 1.96e-06 1.04e-07 | 1.59e+00 1.40e-07 2.22e-07 | 7.9e-10 7.4e-01 co-a 3.00e-01
48 -3.5670e+00 -3.5670e+00 | 4.92e-05 9.82e-07 5.22e-08 | 1.59e+00 6.98e-08 1.10e-07 | 3.7e-10 8.8e-01 co-a 5.00e-01
49 -3.5670e+00 -3.5670e+00 | 1.47e-05 2.96e-07 1.57e-08 | 1.58e+00 2.10e-08 3.30e-08 | 3.0e-10 9.0e-01 co-a 7.00e-01
50 -3.5670e+00 -3.5670e+00 | 7.28e-06 1.49e-07 7.93e-09 | 1.57e+00 1.05e-08 1.64e-08 | 2.2e-09 9.4e-01 co-a 5.00e-01
51 -3.5670e+00 -3.5670e+00 | 3.63e-06 7.49e-08 3.98e-09 | 1.56e+00 5.24e-09 8.15e-09 | 1.3e-09 7.6e-01 co-a 5.00e-01
52 -3.5670e+00 -3.5670e+00 | 1.14e-06 2.52e-08 1.23e-09 | 1.52e+00 1.69e-09 2.57e-09 | 3.7e-08 5.9e-01 co-a 7.00e-01
53 -3.5670e+00 -3.5670e+00 | 3.45e-07 7.60e-09 3.68e-10 | 1.52e+00 5.09e-10 7.75e-10 | 3.7e-09 1.1e-01 co-a 7.00e-01
optimal solution found; terminating
status is Optimal after 53 iterations and 0.273 seconds
In REBayes, this would be accomplished with the Pmix function with exposure arguments, but our approach provides a unified interface with the type system handling the exposure information.
Let’s visualize the results:
pl = plot(
plot(support(norberg_npmle.prior),
probs(norberg_npmle.prior)./sum(probs(norberg_npmle.prior)),
seriestype=:sticks,
xguide = L"\theta",
yguide = L"g(\theta)",
legend=:topright,
label="NPMLE",
title="Risk distribution"),
plot(support(norberg_npmle.prior),
(probs(norberg_npmle.prior) ./ sum(probs(norberg_npmle.prior))).^(1/3),
seriestype=:sticks,
xguide = L"\theta",
yguide = L"g(\theta)^{1/3}",
legend=:topright,
label="NPMLE",
title="Cube-root transformed"),
size=(800,400)
)The left panel shows the estimated risk distribution, while the right panel shows the cube-root transformation to better visualize smaller masses. Similar to the REBayes analysis, we can see evidence of groups with particularly high risk (around θ=8), which would be difficult to capture with a parametric approach.
Finally, we compute the posterior means for each group:
norgberg_postmean_npmle = (PosteriorMean.(norberg_Zs)).(norberg_npmle.prior)72-element Vector{Float64}:
1.516497463019897
1.1477304056884876
0.8900678400349428
1.1229084038155277
1.1343002295969666
1.4664836253077358
1.0164762017419169
1.6494798867585212
0.7913509956260509
1.2713016956237706
1.0779960615909798
0.7848319618996998
2.4249473789626075
⋮
0.8706989730386656
1.0562150646243584
1.035113175842552
0.9955905035611794
0.8996424331413105
1.085110916237371
0.9348278530015118
1.2292844609156781
1.3854937394831124
1.1117890803224781
0.9884559917026545
1.1097830775100936These values represent the empirical Bayes credibility estimates for each occupation group’s risk factor. In insurance terminology, these would be used to set premiums that reflect both the overall portfolio patterns and each group’s specific experience.
# Display the mean posterior estimate
mean(norgberg_postmean_npmle)1.14506213208152This mean value represents the overall risk level across all occupational groups after empirical Bayes adjustment.