4
$\begingroup$

I am looking to recreate a derivation of a PDF using a transform of random variables. I did the original derivation in Mathematica with the Mathstatica package, but when I updated Mathematica recently my Mathstatica installation no-longer works and I am getting no response from Mathstatica support.

Thus, I want to try and do the same thing in plain Mathematica, but am not having any luck.

Here is my Mathematica + Mathstatica code

(* Define our Gaussian distribution *)
f = 1/(\[Sigma] Sqrt[2 \[Pi]])
    Exp[-((x - vergMean)^2/(2 \[Sigma]^2))];   

(* Domain and assumptions *)
domain[f] = {x, 0, 
    Pi/2}   && {vergMean \[Element] Reals, \[Sigma] \[Element] 
     Reals, \[Sigma] > 0, h > 0, Re[h] != 0};

(* Transformation of random variables *)
g = Transform[y == h/Tan[x], f]


This gives the result


(E^(-((vergMean - ArcCot[y/h])^2/(2 \[Sigma]^2))) h)/(Sqrt[
 2 \[Pi]] (h^2 + y^2) \[Sigma])

I am assuming that I should be able to use TransformDistribution to do the same, but I cannot get it to work at all.

For example,


(* Define our Gaussian distribution *)
f = 1/(\[Sigma] Sqrt[2 \[Pi]]) Exp[-((x - vergMean)^2/(2 \[Sigma]^2))];

(* Transformation of random variables *)
transF = 
 TransformedDistribution[y = h/Tan[x], x \[Distributed] f,
  Assumptions -> {vergMean \[Element] Reals, \[Sigma] \[Element] 
     Reals, \[Sigma] > 0, h > 0, Re[h] != 0, x > 0, x < Pi/2}]

PDF[transF, x]

This code doesn't work and just outputs a duplicate of the original code.

I think I am probably doing something dumb / missing something obvious, but I cannot see what it is.

Any help would be gratefully received.

$\endgroup$
3
  • $\begingroup$ The correct syntax would be: f = NormalDistribution[vergMean, σ]; transF = TransformedDistribution[h/Tan[x], x \[Distributed] f, Assumptions -> Element[vergMean, Reals] && σ > 0 && h > 0]; PDF[transF, x] , but unfortunately Mathematica isn't able to find the PDF. $\endgroup$
    – flinty
    Jun 30, 2021 at 14:46
  • $\begingroup$ mathStatica is a great. But it does assume that you determine that the function is one-to-one to be able to use its Transform function. The chapter on transformations is excellent and states "As per Theorem 1, Transform and TransformExtremum should only be used on transformations that are one-to-one." And it even tells you what to do when the transformation is not one-to-one. $\endgroup$
    – JimB
    Jun 30, 2021 at 20:51
  • $\begingroup$ @JimB, thanks. Yes, the it is clear the function has a one-to-one mapping. I am really hoping that I can get Mathstatica up and running again. Updating Mathematica caused some issues and after initial contact via support, I am having trouble getting hold of them. From what they said it is a simple fix as well. $\endgroup$
    – flyingmind
    Jul 1, 2021 at 7:33

2 Answers 2

6
$\begingroup$

EDIT

It seems to work now. The key is to use the TruncatedDistribution so $h \cot(x)$ is monotonic.

transF=TransformedDistribution[h Cot[x],
  x\[Distributed]TruncatedDistribution[{0,π/2},NormalDistribution[μ,σ]],
  Assumptions->μ∈Reals&&σ>0&&h>0];
FullSimplify[PDF[transF,y],{y>0,h>0,μ∈Reals,σ>0}]

gives

(E^(-((μ + 2 ArcTan[(y - Sqrt[h^2 + y^2])/h])^2/(2 σ^2))) h Sqrt[2/π])/((h^2 + y^2) σ *
  (Erf[μ/(Sqrt[2] σ)] - Erf[(-(π/2) + μ)/(Sqrt[2] σ)]))

which is similar to what you had before but properly normalized.

Original

I've been trying to coax Mathematica into doing this but it is still refusing. As a workaround, you can use the definition of a transformed distribution directly. From Wikipedia

$$ f_Y(y) =f_X\big(g^{-1}(y)\big) \left| \frac{d}{dy} \big(g^{-1}(y)\big) \right| $$

Which Mathematica actually handles fine

f=Exp[-((#-μ)^2/(2 σ^2))]/(σ Sqrt[2 π])&;
g=h Cot[#]&;
gi=Quiet@InverseFunction[g];
Simplify[f[gi[y]]Abs[D[gi[y],y]],{y>0,h>0}]

gives

(E^(-((μ-ArcCot[y/h])^2/(2 σ^2))) h)/(Sqrt[2 π] (h^2+y^2) σ)

as desired.

I know that's not ideal, but at least it's a way forward. Still a bit surprised TransformedDistribution with Assumptions or TruncatedDistribution can't seem to handle it.

$\endgroup$
11
  • 1
    $\begingroup$ Is $(0,\pi/2)$ associated with $x$ or $y$? That resulting pdf for $y$ doesn't integrate to 1 over $(0,\pi/2)$ or $(0,\infty)$. For that transformation formula to work, the function needs to be monotonic but only if $x$ is restricted to say $(0,\pi/2)$ or say $(0,\pi)$). I probably haven't had enough coffee this morning so I'm probably missing something. $\endgroup$
    – JimB
    Jun 30, 2021 at 17:18
  • $\begingroup$ Ah, that's fantastic. Yes, by definition the domain of x is (0, Pi/2). In reality for the physical problem it will never get near the the limits. But, geometrically that is the case. I was also surprise vanilla Mathematica could not handle it. Even when putting in all assumptions possible. How would the code need to be adjusted in light of the limits of x? Many thanks again. $\endgroup$
    – flyingmind
    Jun 30, 2021 at 17:31
  • 2
    $\begingroup$ I didn't do any sanity checks since OP gave the desired result in the question. That said, it should be $0<x<\pi/2$ which restricts $y>0$. I have updated the Simplify line to reflect this. Indeed the resulting PDF doesn't integrate to 1 because of the restriction of $0<x<\pi/2$ which makes the original distribution a TruncatedDistribution which itself does not integrate to 1 (they're only considering a slice of the full normal distribution). $\endgroup$
    – bRost03
    Jun 30, 2021 at 17:34
  • 1
    $\begingroup$ Thanks, that's great. I am being a bit slow today. You are correct, the distribution over x is truncated, so the integral will not be 1. I think the etiquette is to give others a chance to answer prior to accepting. But this seems to replicate the behaviour of Mathstatica. Thankyou. $\endgroup$
    – flyingmind
    Jun 30, 2021 at 17:43
  • 1
    $\begingroup$ @flyingmind can't DM on this site, but feel free to ask another question if you want more help! $\endgroup$
    – bRost03
    Jun 30, 2021 at 18:19
4
$\begingroup$

Amplifying on answer by @bRost03

distX = TruncatedDistribution[{0, Pi/2}, NormalDistribution[μ, σ]];

distY = TransformedDistribution[h/Tan[x], x \[Distributed] distX, 
   Assumptions -> h > 0];

g = Assuming[{y > 0, h > 0, σ > 0, μ ∈ Reals},
  PDF[distY, y] // FullSimplify]

(* (E^(-((μ + 2 ArcTan[(y - Sqrt[h^2 + y^2])/h])^2/(
  2 σ^2))) h Sqrt[2/π])/((h^2 + 
   y^2) σ (Erf[μ/(Sqrt[2] σ)] - 
   Erf[(-(π/2) + μ)/(Sqrt[2] σ)])) *)

Verifying that total probability is unity

Assuming[{h > 0, σ > 0, μ ∈ Reals},
 Integrate[g, {y, 0, Infinity}]]

(* 1 *)

Plot3D[g /. {μ -> 0, σ -> 1}, {y, 0, 5}, {h, 0, 5},
 AxesLabel -> (Style[#, 12, Bold] & /@ {y, h, PDF}),
 PlotRange -> {0, 1},
 ClippingStyle -> None]

enter image description here

G = Assuming[
  {y > 0, h > 0, σ > 0, μ ∈ Reals},
  CDF[distY, y] // FullSimplify]

(Erf[(π - 2 μ)/(2 Sqrt[2] σ)] + 
 Erf[(μ + 2 ArcTan[(y - Sqrt[h^2 + y^2])/h])/(Sqrt[2] σ)])/(
Erf[μ/(Sqrt[2] σ)] - Erf[(-(π/2) + μ)/(Sqrt[2] σ)])

Limit[G, y -> Infinity]

(* 1 *)

Plot3D[G /. {μ -> 0, σ -> 1}, {y, 0, 20}, {h, 0, 5},
 AxesLabel -> (Style[#, 12, Bold] & /@ {y, h, CDF}),
 PlotRange -> {0, 1},
 WorkingPrecision -> 15,
 MaxRecursion -> 3]

enter image description here

$\endgroup$
1
  • $\begingroup$ Thanks this is really great. I should have thought about using TruncatedDistribution. I've been using Mathematica for a few years now, but have barely scratched the surface of what it is capable of. $\endgroup$
    – flyingmind
    Jul 1, 2021 at 7:37

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.