5
$\begingroup$

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)

Improve this question
$\endgroup$
1
  • $\begingroup$ @MariuszIwaniuk 's answer below made me realize I should have noticed that the density function in my answer is positive in $(-2,0)$ and integrates to 1 but you mentioned that the stieltjes function is desired for $(0,1]$. I can only conclude that either the stieltjes function is not exactly what you want or that I have made a mistake or both. $\endgroup$
    – JimB
    Commented Oct 14 at 16:49

2 Answers 2

Reset to default
4
$\begingroup$

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]

Estimate of probability density function

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

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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 $\endgroup$
    – Yaroslav Bulatov
    Commented Jun 25, 2023 at 2:29
  • $\begingroup$ Does that formula work for the regular beta distribution? (Away from mathematica right now) $\endgroup$
    – Yaroslav Bulatov
    Commented Jun 25, 2023 at 2:44
  • $\begingroup$ 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]. $\endgroup$
    – JimB
    Commented Jun 25, 2023 at 4:44
  • $\begingroup$ 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 $\endgroup$
    – Yaroslav Bulatov
    Commented Jun 25, 2023 at 7:01
  • $\begingroup$ I confirmed that this transformation is correct, limit + rational values seems like a useful trick to know! (background here) $\endgroup$
    – Yaroslav Bulatov
    Commented Jun 25, 2023 at 18:19
3
$\begingroup$

Using Numeric Inverse Stieltjes Transform from here:

ClearAll["`*"]; Remove["`*"];

f[s_] := (Sqrt[2] (\[Pi] - 2 ArcTan[Sqrt[s]/Sqrt[2]]))/(2 Sqrt[s] + 
Sqrt[2] s ArcTan[Sqrt[2]/Sqrt[s]])

g1[x_] := 1/\[Pi]^2*NIntegrate[1/(s - x)*f[s], {s, -Infinity, x, Infinity}, 
Method -> "PrincipalValue", WorkingPrecision -> 20, 
MaxRecursion -> 20](*Method One*)

A = ListLinePlot[Table[{x, Re[g1[x]]}, {x, -3 + 10^-15, 3, 1/50}], 
PlotStyle -> Red] // Quiet

enter image description here

g2[x_] := -(2/Pi^2)*NIntegrate[(f[x - s] - f[x + s])/(2 s), {s, 10^-15, 
Infinity}](*Method Two*)
B = ListLinePlot[Table[{x, Re[g2[x]]}, {x, -3 + 10^-15, 3, 1/50}], 
PlotStyle -> {Black, Dashed}] // Quiet

enter image description here

Show[{A, B}](*Two Methods on one Plot*)

enter image description here

But that doesn't plot look like @JimB answer!.

EDITED 17.10.2024

From comments @Yaroslav Bulatov Harder example, to invert:

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

(*(3 (Sqrt[s] - \[Pi] Cot[\[Pi]/Sqrt[s]]))/(\[Pi]^2 Sqrt[s])*)

If I use Inverse Stieltjes Transform to:$\frac{3 \left(\sqrt{s}-\pi \cot \left(\frac{\pi }{\sqrt{s}}\right)\right)}{\pi ^2 \sqrt{s}}$ i should get:$\begin{cases} \frac{6 \left(\psi ^{(1)}\left(\left\lceil \frac{1}{\sqrt{x}}\right\rceil \right)-\psi ^{(1)}\left(\left\lfloor \frac{1}{\sqrt{x}}\right\rfloor +1\right)\right)}{\pi ^2} & 0<x\leq 1 \\ 0 & \text{True} \end{cases}$.

Let try:

 PDF[TransformedDistribution[1/z^2, 
 z \[Distributed] ZipfDistribution[1]], x]

 (*Piecewise[{{(6*(PolyGamma[1, Ceiling[1/Sqrt[x]]] - 
 PolyGamma[1, 1 + Floor[1/Sqrt[x]]]))/Pi^2, 
 Inequality[0, Less, x, LessEqual, 1]}}, 0]*)

 Plot[%, {x, -2, 2}, PlotRange -> Full]
 (*I get Zero for all range*)

Using: $$\sum _{n=1}^{\infty } \frac{6}{\pi ^2 \left(-1+n^2 s\right)}=\frac{3 \left(\sqrt{s}-\pi \cot \left(\frac{\pi }{\sqrt{s}}\right)\right)}{\pi ^2 \sqrt{s}}$$ Then:

f[s_] := 6/(\[Pi]^2 (-1 + n^2 s));

Sum[InverseLaplaceTransform[
InverseLaplaceTransform[f[s], s, t] /. t -> -t, t, x], {n, 1, 
Infinity}](*Inverse Stieltjes Transform -method 1*)

(*Sum[(6*DiracDelta[-(1/n^2) + x])/(n^2*Pi^2), {n, 1, Infinity}]*)

Sum[InverseLaplaceTransform[FullSimplify[
InverseMellinTransform[f[s], s, w, GenerateConditions -> False] 
/. w -> Exp[t], Assumptions -> t > 0], t, x], {n, 1, Infinity}] 
(*Inverse Stieltjes Transform -method 2*)
(* 0 *)

Sum[Simplify[InverseMellinTransform[
Simplify[
InverseMellinTransform[f[s], s, Exp[t], 
 GenerateConditions -> False], Assumptions -> t > 0], t, Exp[-x], 
GenerateConditions -> False], Assumptions -> x > 0], {n, 1, 
Infinity}](*Inverse Stieltjes Transform -method 3*)

(* 0 from InverseMellinTransform *)

Sum[FullSimplify[-1/Pi*InverseFourierTransform[-I Sign[z] FourierTransform[f[s], s, z, 
FourierParameters -> {1, -1}], z, x, 
FourierParameters -> {1, -1}], Assumptions -> x > 0], {n, 1, Infinity}]
(*Inverse Stieltjes Transform -method 4*)
(*Sum[(12*DiracDelta[-(1/n^2) + x])/(n^2*Pi^2), {n, 1, Infinity}] two times bigger than method 1*)

Sum[Simplify[1/Pi*InverseFourierTransform[
I*(2*HeavisideTheta[z] - 1)*
FourierTransform[f[s], s, z, FourierParameters -> {1, -1}, 
Assumptions -> n > 0], z, x, FourierParameters -> {1, -1}], 
Assumptions -> x > 0], {n, 1, Infinity}]
(*Inverse Stieltjes Transform -method 5*)
(*Sum[(12*DiracDelta[-(1/n^2) + x])/(n^2*Pi^2), {n, 1, Infinity}] two times bigger than method 1*)

Sum[1/\[Pi]^2*Integrate[1/(s - x)*f[s], {x, -Infinity, Infinity}, 
PrincipalValue -> True, Assumptions -> x > 0], {n, 1, 
Infinity}](*Inverse Stieltjes Transform -method 6*)
(*Integral does not converge*)

Sum[-(2/Pi^2)*Limit[Integrate[(f[x - s] - f[x + s])/(2 s), {s, eps, Infinity}, 
Assumptions -> {eps > 0, x > 0}, GenerateConditions -> False], 
eps -> 0, Assumptions -> x > 0, Direction -> -1], {n, 1, 
Infinity}](*Inverse Stieltjes Transform -method 7*)
(*Sum[-((6*(Log[1/n^2 - x] - Log[-(1/n^2) + x]))/(Pi^4*(-1 + n^2*x))), {n, 1, Infinity}]*)

h[x_] := Sum[-((6 (Log[1/n^2 - x] - Log[-(1/n^2) + x]))/(\[Pi]^4 (-1 + n^2 x))), {n, 1, 100}]
ReImPlot[h[x], {x, -2, 2}](*For Real part h[x] is 0 *)

Computing with Numerics:

R[s_] := (3 (Sqrt[s] - \[Pi] Cot[\[Pi]/Sqrt[s]]))/(\[Pi]^2 Sqrt[s]);

g1[x_] := 1/\[Pi]^2*
NIntegrate[1/(s - x)*R[s], {s, -Infinity, x, Infinity}, 
Method -> "PrincipalValue", WorkingPrecision -> 20, 
MaxRecursion -> 20](*Method One*)

dane = Table[{x, Re[g1[x]]}, {x, -2, 2, 1/20}] // Quiet;
ListLinePlot[dane, PlotStyle -> Red, PlotRange -> All] // Quiet
(*Random noise*)

g2[x_] := -(2/Pi^2)*NIntegrate[(R[x - s] - R[x + s])/(2 s), {s, 10^-15, Infinity}, 
WorkingPrecision -> 20](*Method Two*)
B = ListLinePlot[Table[{x, Re[g2[x]]}, {x, -2, 2, 1/20}]] // Quiet
(*Random noise*)

Could say that the solution is distribution nature that's why numeric's methods fail.

Improve this answer
$\endgroup$
7
  • $\begingroup$ +1 My answer is not only different but I failed to point out that the density was requested to be positive over $(0,1]$ and my result has the density in that region being zero. And your result has negative values for that region. So either the given Stieltjes transform is wrong or I've made a pretty big mistake or both. $\endgroup$
    – JimB
    Commented Oct 14 at 14:55
  • $\begingroup$ Hm....can this be specified to return a probability distribution? I kind of expected inverse Stieltjes transform of a probability distribution to return a probability distribution. The unusual part is that this formula comes from "discrete stieltjes transform", but we are applying continuous inverse, I thought it would be a "smoothed" approximation of distY in question $\endgroup$
    – Yaroslav Bulatov
    Commented Oct 16 at 5:10
  • $\begingroup$ @YaroslavBulatov You can be more specific about what you want to calculate, you mean : InverseLaplaceTransform[ InverseStieltjesTransform[x], x, t] ? or: PDF[InverseStieltjesTransform[x], x] ? $\endgroup$
    – Mariusz Iwaniuk
    Commented Oct 16 at 14:06
  • $\begingroup$ I'm trying to compute Stieltjes transform and inverse Stieltjes for various random variables z. Your previous recipes work for Beta[1/2,1]. But here's a harder example, can you invert it for random variable distY? distY = TransformedDistribution[1/z^2, z \[Distributed] ZipfDistribution[1]]; Expectation[1/(s - i), i \[Distributed] distY] $\endgroup$
    – Yaroslav Bulatov
    Commented Oct 16 at 17:50
  • $\begingroup$ @MariuszIwaniuk I've forked this into a separate question which lists all the Stieltjes transforms I'm interested in -- mathematica.stackexchange.com/questions/307837/… $\endgroup$
    – Yaroslav Bulatov
    Commented Oct 16 at 18:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.