# Plotting CDFs of transformed normal random variables

I'm using the following transformed random variables

hX = TransformedDistribution[(z + 1/3)^2,
z \[Distributed] NormalDistribution[0, 1]]
hY = TransformedDistribution[(z + 1/3)^2,
z \[Distributed] NormalDistribution[1, 2]]


First I tried to plot the PDF of hY by running

Plot[PDF[hY, x], {x,0,1}]


But the plotting fails as some branches of the PDF are not well-defined (even though in the interval (0,1) everything is of course okey). I found a work around on this with Refine[..., x>0] but I would like to know what causes Mathematica to not be able to plot it without it? Any other ways to do it?

Next I wanted to plot the composition function

1 - Refine[CDF[hY, InverseCDF[hX, 1 - x]], InverseCDF[hX, 1 - x] > 0]


But when I run the code

Plot[1 - Refine[CDF[hY, InverseCDF[hX, 1 - x]],
InverseCDF[hX, 1 - x] > 0], {x, 0, 1}, ImageSize -> 600]


It is quite slow (something between 30-60 seconds). Any way to make this faster? I know that it is not that slow but if the distributions were even more complex I bet it would slow down even more. Any way to make it faster? Maybe there are some commands that I don't know of make it faster? Or maybe there is some much faster approximate method.

• tried Plot[Evaluate@PDF[hY, x], {x, 0, 1}]?
– kglr
Apr 10, 2020 at 8:39
• @kglr Great. That solves the first problem but not the second. Apr 10, 2020 at 13:27

## 1 Answer

The first thing to do is to find out where the most time is taken in

Plot[1 - Refine[CDF[hY, InverseCDF[hX, 1 - x]],
InverseCDF[hX, 1 - x] > 0], {x, 0, 1}, ImageSize -> 600]


We can break up the Plot into two parts:

(* The InverseCDF[hX, 1 - x] part evaluated at 102 points *)
xVals = Flatten[{0, 0.001, Range[99]/100., 0.999}];
AbsoluteTiming[y = Table[InverseCDF[hX, 1 - x], {x, xVals}];]
(* {3.50871, Null} *)

(* The CDF[hY, InverseCDF[hX, 1 - x] part *)
AbsoluteTiming[u = 1 - CDF[hY, y];]
(* {16.6575, Null} *)


So the CDF[hY, InverseCDF[hX, 1 - x] part is taking up the most time. But fortunately for your particular function there is an explicit expression for that cdf:

AbsoluteTiming[u = 1 - 1/2 (-Erf[(4 - 3 Sqrt[y])/(6 Sqrt[2])] + Erf[(4 + 3 Sqrt[y])/(6 Sqrt[2])]);]
(* {0.0004323, Null} *)


Putting this altogether on can produce the plot using a ListPlot:

AbsoluteTiming[xVals = Flatten[{0, 0.001, Range[99]/100., 0.999}];
y = Table[InverseCDF[hX, 1 - x], {x, xVals}];
u = 1 - 1/2 (-Erf[(4 - 3 Sqrt[y])/(6 Sqrt[2])] + Erf[(4 + 3 Sqrt[y])/(6 Sqrt[2])]);
ListPlot[Transpose[{xVals, u}], Joined -> True]]


If you want to use Plot rather than ListPlot, then there isn't as much time savings:

AbsoluteTiming[Plot[1 - 1/2 (-Erf[(4 - 3 Sqrt[InverseCDF[hX, 1 - x]])/(6 Sqrt[2])] +
Erf[(4 + 3 Sqrt[InverseCDF[hX, 1 - x]])/(6 Sqrt[2])]), {x, 0, 1}]]