0
$\begingroup$

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.

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

1 Answer 1

1
$\begingroup$

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

Plot of desired function using ListPlot

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

Plot of desired function using Plot

$\endgroup$

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.