Skip to content

Commit 7f02e13

Browse files
committed
feat: add density estimation examples to docs
- Fix mathjax warning in `make.jl` - Add `CairoMakie` and `Distributions` as dependencies.
1 parent 03ca500 commit 7f02e13

File tree

5 files changed

+96
-14
lines changed

5 files changed

+96
-14
lines changed

docs/Project.toml

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,6 @@
11
[deps]
2+
CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
3+
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
24
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
35
DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
46
DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"

docs/make.jl

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -6,10 +6,11 @@ using Revise
66
using Documenter
77
using Documenter.Remotes
88
using DocumenterCitations
9+
using Distributions
910
using DualPerspective
1011

1112
# Define macros here.
12-
const MATHJAX3_MACROS = Dict(
13+
mathjax3_macros = Dict(
1314
:tex => Dict(
1415
:inlineMath => [["\$","\$"], ["\\(","\\)"]],
1516
:tagSide => "left",
@@ -44,14 +45,15 @@ makedocs(
4445
canonical = "https://MPF-Optimization-Laboratory.github.io/DualPerspective.jl/stable",
4546
repolink = "https://github.com/MPF-Optimization-Laboratory/DualPerspective.jl",
4647
# mathengine = Documenter.KaTeX(KATEX_MACROS)
47-
mathengine = Documenter.MathJax3(MATHJAX3_MACROS)
48+
mathengine = Documenter.MathJax3(mathjax3_macros)
4849
),
4950
modules = [DualPerspective],
5051
authors = "Michael P. Friedlander and contributors",
5152
pages = [
5253
"Home" => "index.md",
53-
"User Guide" => "guide.md",
54+
# "User Guide" => "guide.md",
5455
"Theory" => "theory.md",
56+
"Density Estimation" => "density.md",
5557
"Examples" => "examples.md",
5658
"API Reference" => "api.md",
5759
"Development" => "dev.md",

docs/src/density.md

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
```@meta
2+
CurrentModule = DualPerspective.DensityEstimation
3+
CollapsedDocStrings = true
4+
```
5+
6+
# Density Estimation from Moments
7+
8+
The **density-from-moments** problem seeks to reconstruct an unknown probability density from a empirical moment measurements. When the density is supported on the unit interval, the problem is known as the _Haussdorf moment problem_ [gzyl_HausdorffMomentProblem_2010, gzyl_SuperresolutionMaximumEntropy_2017, gzyl_SearchBestApproximant_2007](@cite).
9+
10+
The approach implemented in `DualPerspective.jl` and its submodule `DensityEstimation` is based on the **maximum entropy principle**: among all densities matching the observed moments, select the one with maximal entropy. This yields a robust, principled estimator that converges to the true density as the number of moments increases [borwein_ConvergenceBestEntropy_1991](@cite).
11+
12+
## Formulation
13+
14+
Here we describe the discrete case, where the density is supported on a finite set of locations $x = (x_1, ..., x_n)$. Given an unknown density $p = (p_1, ..., p_n)$ (with $p_j \geq 0$, $\sum_j p_j = 1$), the first $m$ moments are given by
15+
16+
```math
17+
\mu_i = \sum_{j=1}^n x_j^i p_j, \quad i = 1,\ldots, m.
18+
```
19+
20+
In practice, we don't observe the true moments $\mu_i$, but instead we might estimate these from samples $\{X^{(k)}\}_{k=1}^N$:
21+
22+
```math
23+
\hat{\mu}_i = \frac{1}{N} \sum_{k=1}^N \left(X^{(k)}\right)^i, \quad i = 1,\ldots, m.
24+
```
25+
26+
We then solve the **maximum entropy** problem
27+
28+
```math
29+
\max_{p\in\Delta} \left\{
30+
\textstyle\sum_j p_j \log p_j/q_j
31+
\mid
32+
\ A p \approx b
33+
\right\}
34+
```
35+
36+
where $\Delta$ is the set of discrete densities of length $n$, $A$ is the **moment operator** with entries $A_{ij} = x_j^i$, the $m$-vector $b$ collects the empirical moments $\hat{\mu}_i$, and $q\in\Delta$ is a reference density, i.e., a _prior_ on the unknown density $p$.
37+
38+
The function `reconstruct` solves this problem and returns the estimated density.
39+
40+
```@docs
41+
reconstruct
42+
```
43+
44+
## Example
45+
46+
```@example density
47+
using Distributions
48+
using CairoMakie
49+
using DualPerspective.DensityEstimation
50+
51+
# Generate a mixture of two Gaussians
52+
means, vars, weights = [1.0, 5.0], [0.7, 2.0], [0.4, 0.6]
53+
fMix = MixtureModel([Normal(means[i], vars[i]) for i in 1:2], weights)
54+
55+
# Sample the mixture and compute the empirical moments
56+
samples = rand(fMix, 10000)
57+
b = [mean(samples .^ i) for i in 1:6]
58+
59+
# Setup a grid
60+
bnds = [minimum(samples), maximum(samples)]
61+
Δx = 0.01
62+
x = range(bnds[1], bnds[2], step=Δx)
63+
64+
# Reconstruct the density and rescale to convert from a probability mass to density
65+
p = reconstruct(b, x; λ=1e-8, rtol=1e-10)
66+
p = p / Δx
67+
68+
# Plot the true and estimated densities
69+
fig = Figure()
70+
ax = Axis(fig[1, 1], title="Mixture of two Gaussians")
71+
lines!(ax, bnds[1]..bnds[2], x -> pdf(fMix, x), label="true density")
72+
scatter!(ax, x, p, label="estimated density", markersize=7, color=:red)
73+
axislegend(ax)
74+
fig
75+
```
76+
77+
## Additional Examples
78+
79+
See `experiments/julia/density-moments.jl` for end-to-end examples, including:
80+
81+
- Biased die (discrete)
82+
- Gaussian and Gaussian mixtures (continuous)
83+
- Empirical moment estimation

docs/src/refs.bib

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,6 @@ @book{nesterovLecturesConvexOptimization2018
77
publisher = {Springer International Publishing},
88
address = {Cham},
99
doi = {10.1007/978-3-319-91578-4},
10-
urldate = {2025-04-29},
1110
copyright = {http://www.springer.com/tdm},
1211
isbn = {978-3-319-91577-7 978-3-319-91578-4},
1312
langid = {english},
@@ -33,7 +32,6 @@ @book{beck_FirstOrderMethodsOptimization_2017
3332
address = {Philadelphia, PA},
3433
doi = {10.1137/1.9781611974997},
3534
url = {https://doi.org/10.1137/1.9781611974997},
36-
urldate = {2024-05-30},
3735
isbn = {978-1-61197-498-0 978-1-61197-499-7},
3836
langid = {english},
3937
}
@@ -57,7 +55,6 @@ @article{borwein_ConvergenceBestEntropy_1991
5755
pages = {191--205},
5856
issn = {1052-6234, 1095-7189},
5957
doi = {10.1137/0801014},
60-
urldate = {2024-11-05},
6158
abstract = {Given a finite number of moments of an unknown density on a finite measure space, the best entropy estimate--that nonnegative density x with the given moments which minimizes the BoltzmannShannon entropy I(x):= x log x--is considered. A direct proof is given that I has the Kadec property in L--if Yn converges weakly to 37 and I(yn) converges to I(37), then yn converges to 37 in norm. As a corollary, it is obtained that, as the number of given moments increases, the best entropy estimates converge in LI norm to the best entropy estimate of the limiting problem, which is simply in the determined case. Furthermore, for classical moment problems on intervals with strictly positive and sufficiently smooth, error bounds and uniform convergence are actually obtained.},
6259
langid = {english},
6360
}
@@ -73,7 +70,6 @@ @article{gzyl_HausdorffMomentProblem_2010
7370
pages = {3319--3328},
7471
issn = {0096-3003},
7572
doi = {10.1016/j.amc.2010.04.059},
76-
urldate = {2024-11-16},
7773
abstract = {Hausdorff moment problem is considered and a solution, consisting of the use of fractional moments, is proposed. More precisely, in this work a stable algorithm to obtain centered moments from integer moments is found. The algorithm transforms a direct method into an iterative Jacobi method which converges in a finite number of steps, as the iteration Jacobi matrix has null spectral radius. The centered moments are needed to calculate fractional moments from integer moments. As an application few fractional moments are used to solve finite Hausdorff moment problem via maximum entropy technique. Fractional moments represent a remedy to ill-conditioning coming from an high number of integer moments involved in recovering procedure.},
7874
keywords = {Fractional moments,Hankel matrices,Hausdorff moment problem,Maximum entropy,Moments space},
7975
}
@@ -89,7 +85,6 @@ @article{gzyl_SearchBestApproximant_2007
8985
pages = {652--661},
9086
issn = {00963003},
9187
doi = {10.1016/j.amc.2006.11.125},
92-
urldate = {2024-11-16},
9388
abstract = {In this work we want to show that, in general, a universal-best approximation does not exist, and that the choice of the approximant should be related to the context in which the approximation is needed. As we will clearly show, a good choice of approximant to compute expected values, does not necessarily remain the best choice when the context is changed, for example to approximately compute hazard rates.},
9489
copyright = {https://www.elsevier.com/tdm/userlicense/1.0/},
9590
langid = {english},
@@ -107,6 +102,5 @@ @article{gzyl_SuperresolutionMaximumEntropy_2017
107102
publisher = {Taylor \& Francis},
108103
issn = {1741-5977},
109104
doi = {10.1080/17415977.2016.1273918},
110-
urldate = {2024-11-05},
111105
abstract = {The method of maximum entropy has proven to be a rather powerful way to solve the inverse problem consisting of determining a probability density fS(s) on [0,{$\infty$}) from the knowledge of the expected value of a few generalized moments, that is, of functions gi(S) of the variable S. A version of this problem, of utmost relevance for banking, insurance, engineering and the physical sciences, corresponds to the case in which S{$\geq$}0 and gi(s)=exp(-{$\alpha$}is), the expected values E[exp(-{$\alpha$}iS)] are the values of the Laplace transform of S the points {$\alpha$}i on the real line. Since inverting the Laplace transform is an ill-posed problem, to devise numerical techniques that are efficient is of importance for many applications, especially in cases where all we know is the value of the transform at a few points along the real axis. A simple change of variables transforms the Laplace inversion problem into a fractional moment problem on [0,~1]. In this note, we examine why the maximum entropy method provides a good approximation to the density from a few fractional moments, or correspondingly, why a density can be recovered from a few values of its Laplace transform along the real axis.},
112106
}

src/DensityEstimation.jl

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,4 @@
11
module DensityEstimation
2-
32
using ..DualPerspective
43
using QuadGK
54
using LinearAlgebra
@@ -26,17 +25,19 @@ function moment_operator(x, m)
2625
return A
2726
end
2827

29-
"""
30-
reconstruct(μvec::Vector, x_grid::Vector)
28+
@doc raw"""
29+
reconstruct(μvec::Vector, x_grid::Vector; λ=1e-6, kwargs...)
3130
32-
Reconstruct a density from the moments μ_1, μ_2, ..., μ_m contained in the m-vector `μvec`.
31+
Reconstruct a density from the moments $μ_1, μ_2,\ldots, μ_m$ contained in the m-vector `μvec`.
3332
3433
# Arguments
3534
- `μvec`: m-vector of moments.
3635
- `x_grid`: grid of points at which to evaluate moment operator.
36+
- `λ`: regularization parameter (default: 1e-6).
37+
- `kwargs`: passed to the solver. See [`DualPerspective.solve!`](@ref) for details.
3738
3839
# Returns
39-
- `density`: function that esti
40+
- `density`: function that estimates the density at the points in `xgrid`.
4041
"""
4142
function reconstruct(μvec::Vector, xgrid; λ=1e-6, kwargs...)
4243
m = length(μvec)

0 commit comments

Comments
 (0)