Plotting CDFs of transformed normal random variables

Question

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.

Answer 1

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}]]