Visualizing density given by its Stieltjes transformation?

Question

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)

Answer 1

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)}$$