Finding a continuous distribution that fits the empirical density of dataset $1,1/4,1/9,1/16\ldots$

Question

I have a list of coefs of the form $1,1/4,1/9,1/16,\ldots,1/d^2$ sampled with relative sampling frequencies $1,1/4,1/9,1/16,\ldots,1/d^2$.

How do I find a nice continuous density whose CDF closely matches the CDF of the distribution above?

The code below finds a close match to a simpler distribution, which uses sampling frequencies $1,1,1,1,\ldots$ instead of $1,1/4,1/9,1/16,\ldots$, I was looking for an elegant way to extend this to the distribution above

ClearAll["Global`*"];
d = 20; p = 2; $Assumptions = {0 < x < 1, 0 < y < 1};
coefs = Table[i^-p, {i, 1, d}];
{min, max} = MinMax[coefs];
dist = EmpiricalDistribution[coefs];
pdf = y^-(1 + 1/p) Boole[min <= y <= max];
pdf = pdf/Integrate[pdf, {y, 0, 1}];
Print["continuous pdf: ", pdf];
fittedCDF = Integrate[pdf, {y, 0, x}];
LogLogPlot[{CDF[dist, x], fittedCDF}, {x, min, 1/3}, 
 ExclusionsStyle -> Dashed, PlotLegends -> {"empirical", "fitted"}, 
 PlotLabel -> "CDF"]

For a more general setting, related question by Vitaly Kaurov is here.

---

Motivation: moment generating function of this density gives the loss curve observed when minizing $x_1^2+\frac{1}{4}x_2^2+\frac{1}{9}x_3^2+\ldots$ from a random starting point, see this community post

Summarizing briefly

1. loss after $t$ steps of GD minimizing $h_1 x_1^2 + h_2 x_2^2 \ldots + h_d x_d^2$ using step size $1$ and starting with $\mathbf{x0}$ is $$\sum_i h_i(1-h_i)^{2t} \mathbf{x0}_i^2$$

2. using a certain kind of random initial point $\mathbf{x0}$ we get the following in expectation

$$\sum_i h_i(1-h_i)^{2t}$$

3. approximating $(1-h_i)$ with $\exp(-h_i)$ we get $$\sum_i h_i \exp(-2 t h_i)$$

This approximation is valid when $\max_i h_i < 1$ and the problem is "high-dimensional", it has a high effective rank R defined in Bartlett paper

4. Now writing the previous sum in terms of expectation we get the following for some constant $c$ and expectation taken over empirical density of $h$

$$c E_h[h \exp(-2t h)]$$

5. Now modify the empirical density of $h$ by sampling each $h_i$ with relative frequency proportional to $h_i$, call this density $\hat{h}$. We can write our loss as expectation below $$c E_\hat{h}[\exp(-2t h)]$$

6. Observe the close relation of formula above to the moment generating function of $\hat{h}$

Answer 1

Note: Edited to focus on the CDF rather than the moment generating function.

If $X$ is a discrete random variable taking on values $1,1/2^2,1/3^2,1/4^2,\ldots$ with relative frequencies $1,1/2^2, 1/3^2, 1/4^2,\ldots$, then the probability that $X=z^2$ is $6/(z\pi)^2$. (Multiplying by $6/\pi^2$ causes the sum of the probabilities to be 1.)

The objective (I think) is to find a continuous function that either approximates the CDF or reproduces the CDF exactly.

We can find the distribution of $X$ in the following manner:

dist = ProbabilityDistribution[6/(z^2 π^2), {z, 1, ∞, 1}]
distX = TransformedDistribution[1/z^2, z \[Distributed] dist]

The CDF of $X$ is

CDF[distX, x]

So a continuous function that exactly produces the same CDF values at $1,1/4,1/9,1/16,\ldots$ is

Update:

Now that there is a finite population of numbers ($1,1/4,1/9,\ldots,1/d^2)$, the commands to obtain the CDF of X are:

dist = ProbabilityDistribution[6/(z^2 \[Pi]^2), {z, 1, d, 1}, Method -> "Normalize"];
distX = TransformedDistribution[1/z^2, z \[Distributed] dist];
CDF[distX, x] /. Ceiling[1/Sqrt[x]] -> 1/Sqrt[x]

Answer 2

If I understand your idea correctly, you want a function $f(x)$ (your CDF) that has the property that for $n \in \mathbb{N}$, $$ f\left(\frac{1}{n^2}\right) = A \sum_{m = n}^{d} \frac{1}{m^2} $$ (the right-hand side being the total probability of selecting an outcome $1/m^2$ with $m \geq n$ and therefore $1/m^2 < 1/n^2 = x$.) Above, $A$ is some normalization factor chosen such that $f(1) = 1$. Mathematica can solve this in terms of the PolyGamma function:

f[x_, d_] = a Sum[1/m^2, {m, 1/Sqrt[x], d}]
(* a (-PolyGamma[1, 1 + d] + PolyGamma[1, 1/Sqrt[x]]) *)

The normalization factor can then be found by requiring that

normrule = Solve[f[1, d] == 1, a]
(* {{a -> 6/(\[Pi]^2 - 6 PolyGamma[1, 1 + d])}} *)

and so

newcdf[x_, d_] = f[x, d] /. First[normrule]
(* (6 (-PolyGamma[1, 1 + d] + PolyGamma[1, 1/Sqrt[x]]))/(\[Pi]^2 - 6 PolyGamma[1, 1 + d]) *)

i.e., $$ f(x) = \frac{6 \left(\psi ^{(1)}\left(\frac{1}{\sqrt{x}}\right)-\psi ^{(1)}(d+1)\right)}{\pi ^2-6 \psi ^{(1)}(d+1)} $$ Plotting it:

With[{d = 20}, LogLogPlot[newcdf[x, d], {x, 0, 1}]]

... looks kinda dull, actually, after all that work. But it does have the desired property that by construction, its value is exactly equal to the value of the discrete CDF at all values of $x = 1/n^2$ for some integer $n \leq d$.

Finally, note that in the limit $d \to \infty$, a couple of the polygamma functions vanish, and we have

newcdf[x,Infinity]

$$\frac{6}{\pi ^2} \psi ^{(1)}\left(\frac{1}{\sqrt{x}}\right)$$ which looks a little nicer but unfortunately still involves polygamma functions.

Answer 3

You can pick a probability distribution $p(x)\propto x^{-1/2}$:

p[x_] = x^(-1/2);

In this way, the probability of picking a number in the range $\left[\left(n+\frac12\right)^{-2},\left(n-\frac12\right)^{-2}\right]$ (bracketing the number $1/n^2$) is

Assuming[n >= 1,
  Integrate[p[x], {x, 1/(n + 1/2)^2, 1/(n - 1/2)^2}]]
(*    8/(-1 + 4 n^2)    *)

which, for large $n$, is approximately proportional to $1/n^2$. We can interpret this to mean that numbers close to $1/n^2$ have a probability proportional to $1/n^2$, at least for $n\gg1$.

To get a better formula, we can set $p(x)=\sum_{i=0}^{\infty}c_ix^{i-1/2}$, series-expand the integral, and set all terms to zero in order to find the coefficients $c_i$:

p[x_] = 1/2 x^{-1/2} - 1/8 x^{1/2} + 7/96 x^(3/2) - 31/384 x^(5/2) +
        381/2560 x^(7/2) - 2555/6144 x^(9/2) + 1414477/860160 x^(11/2) -
        286685/32768 x^(13/2) + 118518239/1966080 x^(15/2);
s[n_] = Assuming[n >= 1, 
  Integrate[p[x], {x, 1/(n + 1/2)^2, 1/(n - 1/2)^2}] // FullSimplify];
Series[s[n], {n, ∞, 19}]
(*    1/n^2 + O[1/n]^20    *)

This formula works well, except for $n=1$:

Table[s[n], {n, 1, 10}] // N
(*    {893263., 0.255352, 0.111112, 0.0625, 0.04, 0.0277778,
       0.0204082, 0.015625, 0.0123457, 0.01}                    *)