# Best approximation of a compactly supported density by a single Gaussian

Note: This is a follow-up question inspired by a previous (more difficult) question I asked on MathOverflow.

Let $$f:\mathbb{R}\to\mathbb{R}$$ be a (sufficiently regular, e.g. smooth) probability density supported on the (compact) interval $$[a,b]$$. Let $$\Vert\cdot\Vert$$ be some norm (e.g. $$L^2$$, total variation). What is the best approximation to $$f$$ by a single, univariate Gaussian? In other words, solve $$\text{argmin}_{m\in\mathbb{R},v^2\ge 0}\Vert f-g_{m,v^2}\Vert$$

where $$g_{m,v^2}$$ is the Gaussian PDF with mean $$m$$ and variance $$v^2$$. As I mention above, this question is related to a previous question; I realized while trying to solve that problem that even the simple case of a single Gaussian seems nontrivial.

Hint: I have verified numerically, for both $$L^2$$ and TV norm, that the answer is not $$m=\mathbb{E}_{X\sim f} X$$ or $$v^2=\text{var}_{X\sim f} (X)$$.

Note: The details of which norm (or even metric) or regularity assumptions are made are not important here. I am mostly curious to see if this computation can be carried out for any "reasonable" distance and regularity. In fact, even the case of minimizing over $$m$$ alone (holding $$v^2$$ fixed) appears difficult. Note that for the KL-divergence (which is neither a norm nor a metric), this can be solved in closed form quite easily.

• The natural thing to try, for say $L^2$, is to differentiate (under the integral sign) with respect to $m$ and $v$, set them equal to zero, and try to solve. Does anything useful happen when you do that? – Nate Eldredge Mar 28 '19 at 16:39
• @NateEldredge Yes, I have tried that and the resulting integral equation is daunting. – JohnA Mar 28 '19 at 16:41
• You certainly shouldn't expect to get a closed-form solution to the equations (say for $L^2$), but numerical methods should work reasonably well. – Robert Israel Mar 28 '19 at 17:57
• Among the other possibilities for norms to use: en.wikipedia.org/wiki/Statistical_distance#Examples – Matt F. Mar 28 '19 at 18:20
• @MattF Absolutely. The reason I focused the question on norms (and $L^2$ in particular) is that I suspect this should be easier to solve in a Hilbert space. Any probability metric (e.g. Wasserstein, Hellinger) would be of interest as well. – JohnA Mar 28 '19 at 20:03

(i) You can try to minimize the squared Hellinger distance $$H(m,v)^2:=\frac12\,\int(\sqrt f-\sqrt{g_{m,v^2}})^2 =1-\int \sqrt f \sqrt{g_{m,v^2}},$$ which amounts to the maximization of the comparatively simple expression $$J(m,v):=\int \sqrt{f(x)}\frac1{\sqrt v}\,\exp\Big\{-\frac{(x-m)^2}{4v^2}\Big\}dx$$ in $$m,v$$.

In particular, if $$f$$ is log concave, then the integrand in the latter integral is log concave in $$(m,x)$$. So, by the Prékopa--Leindler theorem (see e.g. Corollary 3.5), $$J(m,v)$$ is log concave in $$m$$ and hence will usually have a unique maximum in $$m$$. If, moreover, $$f$$ is symmetric about the midpoint $$(a+b)/2$$ of the interval $$[a,b]$$, then the maximum of $$J(m,v)$$ in $$m$$ is at $$m=(a+b)/2$$, and then it remains to maximize $$J((a+b)/2,v)$$ in $$v>0$$, which should usually be easy to do.

E.g., if $$f$$ is the density of the beta distribution Beta$$(p,p)$$ with parameters $$\alpha=p$$ and $$\beta=p$$ with $$p\ge1$$, then $$J(m,v)$$ is log concave in $$m$$ and its maximum in $$m$$ is at $$m=1/2$$. Here is a graph of the ratio of the maximizer of $$J(1/2,v)$$ in $$v>0$$ to the true standard deviation of the Beta$$(p,p)$$ distribution as the function of $$p\ge1$$:

We see that this ratio is pretty close to $$1$$ already for $$p=1$$.

(ii) Alternatively, you can try to minimize $$\int(\ln f-\ln{g_{m,v^2}})^2d\mu =\int\Big(\ln f(x)+\ln\sqrt{2\pi}+\ln v+\frac{(x-m)^2}{2v^2}\Big)^2\,\mu(dx)$$ for some measure $$\mu$$. Expanding the latter integrand and then integrating term-wise, we see that, for a given measure $$\mu$$, this reduces to the minimization in $$v>0$$ and $$m$$ of a linear combination of the form $$\sum_{j=1}^{15}c_j w_j(m,v),$$ where, for each $$j$$, $$c_j$$ is a known real number and $$w_j(m,v)$$ is a rather simple elementary function of $$m,v$$.

• Why the downvote? This approach seems to lead to a tractable case. Also, is the normalization in your maximization expression correct? – lcv Mar 28 '19 at 20:14
• @lcv : Thank you for your comment. I have now corrected that normalization factor. – Iosif Pinelis Mar 28 '19 at 20:32
• For i), would you suggest integration over $[a,b]$ or over $\mathbb{R}$? For ii), the only sensible bounds are $[a,b]$. – Matt F. Mar 30 '19 at 3:43
• @MattF. : All the integrals are over $\mathbb R$. Note that in (ii) the integral is with respect to a measure $\mu$, to be appropriately chosen. – Iosif Pinelis Mar 31 '19 at 0:55
• I have added some details. – Iosif Pinelis Apr 1 '19 at 14:58