3
$\begingroup$

Are there any implementations of Inverse Stieltjes Transform in Mathematica?

For instance, in the example below, I'm hoping to recover Beta(1/2,1) density $\frac{1}{2 \sqrt{x}}$ from its Stieltjes Transform $\frac{\coth ^{-1}\left(\sqrt{s}\right)}{\sqrt{s}}$.

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*)

This would give an alternative way to we recover density from its moment generating function. Instead of inverse Laplace Transform of the MGF (Mathematica gets stuck), we would do forward Laplace, followed by the inverse Stieltjes.

This answer of JimB shows that it's possible, I'm looking for a robust standalone implementation.

Collection of other inverse Laplace Transform tricks: laplace-transform-tricks.nb

Improve this question
$\endgroup$
0

1 Answer 1

Reset to default
8
+500
$\begingroup$

On way is:

f[s_] := ArcCoth[Sqrt[s]]/Sqrt[s]

InverseLaplaceTransform[FullSimplify[
InverseMellinTransform[f[s] // TrigToExp, s, w, 
GenerateConditions -> False] /. w -> Exp[t], 
Assumptions -> t > 0], t, x]

(* HeavisideTheta[1 - x]/(2 Sqrt[x]) *)

Simplify[%, Assumptions -> 0 < x < 1]

(* 1/(2 Sqrt[x]) *)

Another way:

Simplify[InverseMellinTransform[Simplify[
InverseMellinTransform[f[s] // TrigToExp, s, Exp[t], 
GenerateConditions -> False], Assumptions -> t > 0], t, Exp[-x], 
GenerateConditions -> False], Assumptions -> x > 0]

(* HeavisideTheta[1 - x]/(2 Sqrt[x]) *)

EDITED 12.10.2024

I managed to find the Inverse Stieltjes transform is similar like Hilbert transform:

 InverseStieltjesTransform = 
 1/\[Pi]^2*Integrate[1/(s - x)*f[s], {x, -Infinity, Infinity}, 
 PrincipalValue -> True] == 
 Simplify[-1/Pi*
 InverseFourierTransform[-I Sign[z] FourierTransform[f[s], s, z, 
 FourierParameters -> {1, -1}], z, x, 
 FourierParameters -> {1, -1}], Assumptions -> x > 0] == 
 Simplify[
 1/Pi*InverseFourierTransform[
 I*(2*HeavisideTheta[z] - 1)*
 FourierTransform[f[s], s, z, FourierParameters -> {1, -1}], z, 
 x, FourierParameters -> {1, -1}], Assumptions -> x > 0]

Stieltjes transform of:

$\begin{cases} \frac{1}{2 \sqrt{x}} & 0<x<1 \\ 0 & \text{True} \end{cases}$ is: $\begin{cases} \frac{\tanh ^{-1}\left(\frac{1}{\sqrt{s}}\right)}{\sqrt{s}} & \Re(s)\geq 1\lor \Re(s)\leq 0\lor \Im(s)>0\lor \Im(s)<0 \\ \frac{\tanh ^{-1}\left(\sqrt{s}\right)}{\sqrt{s}} & \text{True} \end{cases}$ NOT $\frac{\coth ^{-1}\left(\sqrt{s}\right)}{\sqrt{s}}$ !

let:

$Version
(*"14.1.0 for Microsoft Windows (64-bit) (July 16, 2024)"*)

func = {Sin[x], Cos[x], Exp[I x], 1/(2 Sqrt[x]), Piecewise[{{1/(2*Sqrt[x]), 0 < x < 1}}, 0]};(*Function examples.*)

A = Integrate[1/(s - x)*(#), {x, -Infinity, Infinity}, 
PrincipalValue -> True] & /@ func(*Stieltjes Transform*)

A1 = Integrate[1/(s - x)*(#), {x, -Infinity, Infinity}, 
PrincipalValue -> True, GenerateConditions -> False] & /@ 
func(*Stieltjes Transform*)

B = 1/\[Pi]^2*Integrate[1/(s - x)*(#), {s, -Infinity, Infinity}, 
PrincipalValue -> True] & /@ A (*Inverse Stieltjes Transform*)
Simplify[B, Assumptions -> x > 0](*OK, can't calculate *)

B1 = 1/\[Pi]^2*Integrate[1/(s - x)*(#) // Expand, {s, -Infinity, Infinity}, 
 PrincipalValue -> True, Assumptions -> x > 0] & /@ A1[[1 ;; 4]](*OK! Inverse Stieltjes Transform*)

 (*{Sin[x], Cos[x], E^(I x), 1/(2 Sqrt[x])}*)

Hard example to invert :

 1/\[Pi]^2*Integrate[
 1/(s - x)*(-(\[Pi]/(2 Sqrt[-s]))), {s, -Infinity, Infinity}, 
 PrincipalValue -> True, 
 Assumptions -> x > 0](*Inverse Stieltjes Transform. Ok Works fine*)  


 1/\[Pi]^2*Integrate[1/(s - x)*(1/
 4 (-2 I \[Pi] Sqrt[1/s] - (-2 Log[-(1/s)] + 
    2 (Log[1/s] + Log[-(1/Sqrt[s])]) + Log[s])/Sqrt[
   s])), {s, -Infinity, Infinity}, PrincipalValue -> True, 
 Assumptions -> x > 0](*Inverse Stieltjes Transform. Ok Works fine*)

 1/\[Pi]^2*Integrate[1/(s - x)*(Piecewise[{{ArcTanh[1/Sqrt[s]]/Sqrt[s], 
   Re[s] >= 1 || Re[s] <= 0 || Im[s] > 0 || Im[s] < 0}}, 
  ArcTanh[Sqrt[s]]/Sqrt[s]]), {s, -Infinity, Infinity}, 
 PrincipalValue -> True, 
 Assumptions -> 
 x > 0](*Inverse Stieltjes Transform.I Aborted to long time to compute!*)

  1/\[Pi]^2*Integrate[1/(s - x)*((2 ArcCoth[Sqrt[s]] + Log[-(1/s)] - Log[1/s])/(
  2 Sqrt[s])), {s, -Infinity, Infinity}, PrincipalValue -> True, 
  Assumptions -> 
  x > 0](*Inverse Stieltjes Transform.I Aborted to long time to compute!*)

  Simplify[-1/Pi*InverseFourierTransform[-I Sign[z] *FourierTransform[
  Piecewise[{{ArcTanh[1/Sqrt[s]]/Sqrt[s], 
    Re[s] >= 1 || Re[s] <= 0 || Im[s] > 0 || Im[s] < 0}}, 
  ArcTanh[Sqrt[s]]/Sqrt[s]], s, z, FourierParameters -> {1,-1}], 
  z, x, FourierParameters -> {1, -1}], 
  Assumptions -> 
  x > 0](*Inverse Stieltjes Transform.I Aborted to long time to compute!*)

Workaround for one hard example,I use a Trick:

$$\frac{\int_{-\infty }^{\infty } \frac{ \begin{array}{cc} \{ & \begin{array}{cc} \frac{\tanh ^{-1}\left(\sqrt{s}\right)}{\sqrt{s}} & \Re(s)\geq 1\lor \Re(s)\leq 0\lor \Im(s)>0\lor \Im(s)<0 \\ \frac{\tanh ^{-1}\left(\sqrt{s}\right)}{\sqrt{s}} & \text{True} \\ \end{array} \\ \end{array} }{s-x} \, ds}{\pi ^2}=\int_0^1 \frac{\int_{-\infty }^{\infty } \frac{\frac{\partial }{\partial Q}\left( \begin{array}{cc} \{ & \begin{array}{cc} \frac{\tanh ^{-1}\left(\frac{Q}{\sqrt{s}}\right)}{\sqrt{s}} & \Re(s)\geq 1\lor \Re(s)\leq 0\lor \Im(s)>0\lor \Im(s)<0 \\ \frac{\tanh ^{-1}\left(Q \sqrt{s}\right)}{\sqrt{s}} & \text{True} \\ \end{array} \\ \end{array} \right)}{s-x} \, ds}{\pi ^2} \, dQ$$

L = Assuming[x \[Element] Reals, 1/\[Pi]^2*
Integrate[
1/(s - x)*(D[
   Piecewise[{{ArcTanh[Q/Sqrt[s]]/Sqrt[s], 
      Re[s] >= 1 || Re[s] <= 0 || Im[s] > 0 || Im[s] < 0}}, 
    ArcTanh[Q*Sqrt[s]]/Sqrt[s]], Q]), {s, -Infinity, Infinity}, 
 PrincipalValue -> True]]

L1 = FullSimplify[L, Assumptions -> {0 < Q < 1, x \[Element] Reals}]
Integrate[L1, {Q, 0, 1}, Assumptions -> x \[Element] Reals]
(*To long time to compute !!!*)

we compute numerically:

R[x_] := NIntegrate[
Piecewise[{{Log[-(((-1 + Q^2)*x)/(-1 + x))]/(-1 + Q^2*x) - 
   Log[-(((-1 + Q^2)*x)/(Q^2*(-1 + x)))]/(Q^2 - x), 
  x > 1}, {((-I)*Pi + Log[((-1 + Q^2)*x)/(-1 + x)])/(-1 + Q^2*x) -
    Log[-(((-1 + Q^2)*x)/(Q^2*(-1 + x)))]/(Q^2 - x), x <= 0}}, 
Log[((-1 + Q^2)*x)/(-1 + x)]/(-1 + Q^2*x) - 
 Log[((-1 + Q^2)*x)/(Q^2*(-1 + x))]/(Q^2 - x)]/Pi^2, {Q, 0, 1}];

Plot[{R[x], Piecewise[{{1/(2*Sqrt[x]), 0 < x < 1}}, 0]}, {x, 0, 2}, 
PlotStyle -> {Red, {Dashed, Black}},PlotLabels -> {"numeric", "analytically"}]

enter image description here

Computing Inverse Stieltjes transform numerically:

N[Piecewise[{{1/(2*Sqrt[x]), 0 < x < 1}}, 0] /. x -> 1/2, 20]
(* 0.70710678118654752440 *)

g[s_]:=(Piecewise[{{ArcTanh[1/Sqrt[s]]/Sqrt[s], 
Re[s] >= 1 || Re[s] <= 0 || Im[s] > 0 || Im[s] < 0}}, 
ArcTanh[Sqrt[s]]/Sqrt[s]])

1/\[Pi]^2*NIntegrate[1/(s - x)*g[s] /. x -> 1/2, {s, -Infinity, 1/2, 
Infinity}, Method -> "PrincipalValue", WorkingPrecision -> 20] 

(*0.70710678118654752494*)

Another example:

X = LaplaceTransform[MomentGeneratingFunction[NormalDistribution[\[Mu], \[Sigma]], t], t, s]
X1 = X /. \[Mu] -> 1/2 /. \[Sigma] -> 1/3
1/\[Pi]^2*NIntegrate[1/(s - x)*X1 /. x -> 1, {s, -Infinity, 1, Infinity}, 
Method -> "PrincipalValue", WorkingPrecision -> 20]

(* 0.38855278699767518268 + 0.71823717421850285373 I *)

(*and here's the puzzle, I don't know why I have to take the real part *)

Re[%]
(* 0.38855278699767518268 *)

N[PDF[NormalDistribution[\[Mu], \[Sigma]], x] /. \[Mu] -> 
 1/2 /. \[Sigma] -> 1/3 /. x -> 1, 20]

(*0.38855278699767518284*)

EDITED 13.10.2024

If we don't use the principal-value integral can be written explicitly:

 B = {-\[Pi] Cos[s], \[Pi] Sin[s], -I E^(I s) \[Pi], -(\[Pi]/(2 Sqrt[-s]))};
 g[s_] = B[[n]] // Quiet

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

 {Sin[x], Cos[x], E^(I x), 1/(2 Sqrt[x])}
Improve this answer
$\endgroup$
7
  • $\begingroup$ Added ilaplaceMellinTrig function that uses this trick here, right now ilaplaceMellinTrig, ilaplaceMellin and ilaplaceFeynmann utils in that file seem to cover most of my use-cases (thanks mostly to your help) $\endgroup$
    – Yaroslav Bulatov
    Commented Oct 11 at 5:47
  • $\begingroup$ I expected formula from wikipedia to work, but it produces 0 S[x_] := ArcTanh[1/Sqrt[x]]/Sqrt[x]; Limit[(S[x - I eps] - S[x + I eps])/(2 I Pi), eps -> 0, Direction -> "FromAbove"] $\endgroup$
    – Yaroslav Bulatov
    Commented Oct 11 at 6:00
  • $\begingroup$ @YaroslavBulatov I updated my answer. $\endgroup$
    – Mariusz Iwaniuk
    Commented Oct 12 at 22:09
  • $\begingroup$ Thanks, this is useful reference! (will link in an overview of techniques on wolfram community soon). PS, I found this simple case of Stieltjes transform works out of the box (a little confused by - sign in the result, whether it's a bug) S[z_] = (2 ArcCoth[Sqrt[z]])/Sqrt[z]; expr = 1/Pi Limit[ Im@ComplexExpand@S[x + I \[CurlyEpsilon]], \[CurlyEpsilon] -> 0]; Assuming[{0 < x < 1}, expr // FullSimplify] $\endgroup$
    – Yaroslav Bulatov
    Commented Oct 12 at 22:21
  • 1
    $\begingroup$ @YaroslavBulatov I updated my answer for more one definition of Inverse Stieltjes Transform. $\endgroup$
    – Mariusz Iwaniuk
    Commented Oct 13 at 9:49

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.