# Visualizing density given by its Stieltjes transformation?

How can I use Mathematica to get $$f(y)$$, the PDF of a random variable given its Stieltjes transform $$g(s)$$?

$$g(s)=E\left[1/(s-X)\right]$$

I need to visualize probability density of a random variable $$X\in(0,1]$$ with the following Stieltjes transform:

stieltjes=(Sqrt[2] (\[Pi] - 2 ArcTan[Sqrt[s]/Sqrt[2]]))/(2 Sqrt[s] + Sqrt[2] s ArcTan[Sqrt[2]/Sqrt[s]])


Expecting this PDF to look similar to $$\operatorname{Beta}\left(\frac{1}{2},1\right)$$

Trying on some simpler distributions first:

Using Beta$$(1/2,1)$$ we have the following pdf $$f(y)$$ and $$g(s)$$

pdf = 1/2 y^(-1/2); (* BetaDistribution[1/2,1] *)
stieltjes = Integrate[pdf 1/(s - y), {y, 0, 1}]


$$f(y)=\frac{\frac{1}{\sqrt{y}}}{2} \leftrightarrow g(s)=\frac{\coth ^{-1}\left(\sqrt{s}\right)}{\sqrt{s}}$$

Also, using $$1/x^2$$ where $$x$$ is Zipf distributed as in the earlier question we have the following

distY = TransformedDistribution[1/z^2,
z \[Distributed] ZipfDistribution[1]];
Expectation[1/(s - i), i \[Distributed] distY]


$$f(y)=\left(1, \frac{6}{\pi ^2}\right),\left(\frac{1}{4}, \frac{3}{2 \pi ^2}\right), \left(\frac{1}{9} , \frac{2}{3 \pi ^2}\right),\ldots \leftrightarrow g(s)=-\frac{3 \left(\pi \cot \left(\frac{\pi }{\sqrt{s}}\right)-\sqrt{s}\right)}{\pi ^2 \sqrt{s}}$$

You can get Stieltjes Transforms by chaining two Laplace transforms together as seen below. Wondering if there's a way to use this to get symbolic or numeric inversion working for examples above.

dist = BetaDistribution[1/2, 1];
mgf = MomentGeneratingFunction[BetaDistribution[1/2, 1], t];
doubleLaplace = LaplaceTransform[mgf, t, s] // FullSimplify;
stieltjes = Expectation[1/(s - y), y \[Distributed] dist];
doubleLaplace == stieltjes (*True if s>1*)


(note you often need Feynman trick to do inverse Laplace transforms in Mathematica, documented here)

From the formula at Wiki Stieltjes Transformation, we should be able to get the density as

$${\displaystyle \rho (x)=\lim _{\varepsilon \to 0^{+}}{\frac {S_{\rho }(x-i\varepsilon )-S_{\rho }(x+i\varepsilon )}{2i\pi }}.}$$

where $$S_\rho$$ is equivalent to your $$g$$. I have not been successful at taking that limit with your stieltjes (as I keep getting 0 as the limit) but using very small values of $$\epsilon$$ I get the following approximation:

g[s_] := (Sqrt[2] (π - 2 ArcTan[Sqrt[s]/Sqrt[2]]))/(2 Sqrt[s] + Sqrt[2] s ArcTan[Sqrt[2]/Sqrt[s]])

f[x_, ϵ_] := (g[x - I ϵ] - g[x + I ϵ])/(2 π I) // N // Chop

Plot[f[x, 1/1000000], {x, -3, 1}, PlotRangeClipping -> False]


But that doesn't look like your expected beta distribution with parameters 1/2 and 1.

Update:

By setting $$x$$ to rational values, using Limit worked and a pattern became apparent. The pdf is given by

pdf[x_] := Piecewise[{{(4 Sqrt[2])/(Sqrt[-x] (8 - π^2 x +
4 Sqrt[2] Sqrt[x] ArcTan[(2 Sqrt[2] Sqrt[x])/(-2 + x)] +
x ArcTan[(2 Sqrt[2] Sqrt[x])/(-2 + x)]^2)), -2 < x < 0}}, 0]


$$\frac{4 \sqrt{2}}{\sqrt{-x} \left(-\pi ^2 x+x \tan ^{-1}\left(\frac{2 \sqrt{2} \sqrt{x}}{x-2}\right)^2+4 \sqrt{2} \sqrt{x} \tan ^{-1}\left(\frac{2 \sqrt{2} \sqrt{x}}{x-2}\right)+8\right)}$$

• Interesting!.... actually that looks close to what I expected, my expression was obtained using two Laplace transforms instead of MFG followed by Laplace, so that may explain the negative sign Commented Jun 25, 2023 at 2:29
• Does that formula work for the regular beta distribution? (Away from mathematica right now) Commented Jun 25, 2023 at 2:44
• Yes but I've only got it to work for rational values of $x$: stieltjes[s_] := ArcCoth[Sqrt[s]]/Sqrt[s]; Limit[(stieltjes[x - I \[Epsilon]] - stieltjes[x + I \[Epsilon]])/(2 \[Pi] I) /. x -> 7/8, \[Epsilon] -> 0, Direction -> "FromAbove"] gives the same result as PDF[BetaDistribution[1/2, 1], 7/8].
– JimB
Commented Jun 25, 2023 at 4:44
• Thanks for the update...this expression is also the one I hit a dead end trying to compute the inverse Laplace transform of, being able to get the explicit expression for inverse Stieltjes is promising Commented Jun 25, 2023 at 7:01
• I confirmed that this transformation is correct, limit + rational values seems like a useful trick to know! (background here) Commented Jun 25, 2023 at 18:19