2
$\begingroup$

I have this simple code, but it takes ages to plot the last function: FS12... how to speed it up? Thank you

Clear["Global`*"];

(*Initial Input data*)

u11 = 1;
(*s11=u11*0.05;*)

s11 = u11*1;

d1 = 50; p1 = 0.1;

u22 = 1;
(*s22=u22*0.25;*)

s22 = u22*0.4;

d2 = 50; p2 = 0.1;

t0 = 50;

F = 0.99;

fS11[r_] := PDF[NormalDistribution[u11, s11], r];

Quiet[Plot[fS11[r], {r, 0, 2}, AspectRatio -> 0.2, PlotRange -> Full]]

FS11[r_] := CDF[NormalDistribution[u11, s11], r];
Quiet[Plot[FS11[r], {r, 0, 4}, AspectRatio -> 0.2, PlotRange -> Full]]

FS1[x_] := (1 - p1)*UnitStep[x] + p1*FS11[x];
Quiet[plt1 = 
  Plot[FS1[x], {x, -2, 3}, AspectRatio -> 0.2, PlotRange -> Full, 
   PlotStyle -> Blue]]

Quiet[dataS1 = 
   Table[x /. FindRoot[FS1[x] == y, {x, -1, 5}], {y, 0, 0.99, 0.001}]];

FS1n[x_] := CDF[SmoothKernelDistribution[dataS1], x]

Quiet[plt2 = 
  Plot[FS1n[x], {x, -2, 4}, AspectRatio -> 0.2, PlotRange -> Full, 
   PlotStyle -> Red]]

Show[plt1, plt2]

fS1n[x_] := PDF[SmoothKernelDistribution[dataS1], x];

Quiet[Plot[fS1n[x], {x, -1, 2}, AspectRatio -> 0.2, PlotRange -> Full]]

fS22[x_] := PDF[NormalDistribution[u22, s22], x] 

Quiet[Plot[fS22[x], {x, -3, 3}, AspectRatio -> 0.2, PlotRange -> Full]]

FS22[x_] := CDF[NormalDistribution[u22, s22], x]

Quiet[Plot[FS22[x], {x, 0, 3}, AspectRatio -> 0.2, PlotRange -> Full]]

FS2[x_] := ((1 - p2)*UnitStep[x] + p2*FS22[x])^(d1/d2)

Quiet[Plot[FS2[x], {x, -2, 4}, AspectRatio -> 0.2, PlotRange -> Full]]

FS12[x_] := Evaluate[NIntegrate[fS1n[y]*FS2[x - y], {y, -4, +4}]];

Quiet[Plot[FS12[x], {x, 0, 2}, AspectRatio -> 0.2, PlotRange -> Full]]
$\endgroup$
12
  • $\begingroup$ Have you considered possibly rephrasing your problem using tools such as (N)Probability, (N)Expectation, Distributed and Conditioned? $\endgroup$
    – kirma
    Mar 4 '20 at 17:37
  • 1
    $\begingroup$ The definition of FS12 throws the following message in v12 on Mac: NIntegrate::inumr: The integrand ([Piecewise] <<1>>) (0.05 Erfc[1.76777 (1-x+y)]+0.9 UnitStep[x-y]) has evaluated to non-numerical values for all sampling points in the region with boundaries {{-4,-3.99987}}. I'm not saying this is the problem but perhaps someone who understands this better might. $\endgroup$
    – Mark R
    Mar 4 '20 at 17:57
  • $\begingroup$ @kirma I don't see how those functions help me in getting FS12 (which is a convolution, although Convolute didn't work for me because fS1n is not a distribution per se).... $\endgroup$
    – jpcgandre
    Mar 4 '20 at 17:58
  • 2
    $\begingroup$ The more I think about this is that the distrbiution from which dataS1 was sampled is just not smooth. So using SmoothKernelDistribution does not make any sense. The same holds for the convolution against its PDF: It just does not exist. $\endgroup$ Mar 5 '20 at 7:31
  • 1
    $\begingroup$ Hmm. You actually want to calculate the distribution of random variable X+Y, right? You can do that by convolving their PDFs, if existent. Here, one one of the two random variables has a density, so you may apply this to determine the CDF of X+Y instead. $\endgroup$ Mar 5 '20 at 8:39
3
$\begingroup$

This is one of the cases where Set should be preferred over SetDelayed in a function definition:

FS1n[x_] := CDF[SmoothKernelDistribution[dataS1], x];

Quiet[plt2 = 
     Plot[FS1n[x], {x, -2, 4}, AspectRatio -> 0.2, PlotRange -> Full, 
      PlotStyle -> Red]]; // AbsoluteTiming // First

Block[{x},
  FS1n[x_] = CDF[SmoothKernelDistribution[dataS1], x];
  ];

Quiet[plt2 = 
     Plot[FS1n[x], {x, -2, 4}, AspectRatio -> 0.2, PlotRange -> Full, 
      PlotStyle -> Red]]; // AbsoluteTiming // First

7.21294

0.047187

So a more than 100-fold speed up. Why? SmoothKernelDistribution[dataS1] has to build an interpolation function first by doing some statistical tricks. If you use SetDelayed, this computation has to be performed _every time you call FS1n. With Set, this has to be done only once, precisely at the time of definition of FS1n.

This does not completely resolve the slowness of the last plot, though. That is due to some problem in NIntegrate. Have to think about it...

$\endgroup$
9
  • $\begingroup$ Actually if I useBlock[{x},FS12[x_]=NIntegrate[fS1n[y]*FS2[x-y],{y,-4,+4}];]; it takes 3sec longer that if I use FS12[x_]:=NIntegrate[fS1n[y]*FS2[x-y],{y,-4,+4}];... $\endgroup$
    – jpcgandre
    Mar 5 '20 at 10:39
  • $\begingroup$ Apparently... when I use fS1[x_]:=(1-p1)*DiracDelta[x]+p1*fS11 it works super fast... $\endgroup$
    – jpcgandre
    Mar 5 '20 at 11:08
  • 1
    $\begingroup$ NIntegrate is a case where you do not want to use Set (or SetDelayed+Evaluate)... $\endgroup$ Mar 5 '20 at 11:10
  • 1
    $\begingroup$ And towards the second comment: Yes, now you are convolutin the correct CDF against the correct PDF... =) Convolution against a DiracDelta is pretty easy in the end ;) $\endgroup$ Mar 5 '20 at 11:11
  • 1
    $\begingroup$ Sorry, I am quite confused and actually also busy, you know... But if it works for you, it is fine for me. =) $\endgroup$ Mar 5 '20 at 11:28
2
$\begingroup$

Perhaps not very satisfying, but I was able to get this to complete in a reasonable amount of time (1 minute, 16 seconds):

plotValues = ParallelTable[{x, FS12[x]}, {x, 0.01, 2, .01}];
ListPlot[plotValues, AspectRatio -> 0.2, PlotRange -> Full]

My machine is a Mac running v12 and it has 6 kernels.

Here is the plot: ListPlot of FS12

UPDATE:
I went back and started at 0 and notice that there is a discontinuity between 0 and .01: ListPlot of FS12 starting at 0 This doesn't explain the problem since even starting Plot at .01 takes a long time.

$\endgroup$
5
  • $\begingroup$ That's good. I'll try in my machine. Thank you! $\endgroup$
    – jpcgandre
    Mar 4 '20 at 18:16
  • 1
    $\begingroup$ I did notice that the interpolating function used by FS12 has some discontinuities with very small numbers. Here are the first few values with a step size of .001: {0., 130.44}, {0.001, 110.191}, {0.002, 89.9411}, {0.003, 69.6915}, \ {0.004, 49.442}, {0.005, 29.1925}, {0.006, 8.94294}, {0.007, \ 5.28131}, {0.008, 4.32973}, {0.009, 3.37815}, {0.01, 2.42657}, \ {0.011, 1.47499}, {0.012, 0.523408}, {0.013, 1.88788*10^-17}, {0.014, 1.5711*10^-17}, {0.015, 1.25433*10^-17}, {0.016, 9.3755*10^-18}, {0.017, 6.20772*10^-18}, {0.018, 3.03995*10^-18}, {0.019, 0.}, {0.02, 0.}, {0.021, 0.} $\endgroup$
    – Mark R
    Mar 4 '20 at 19:46
  • $\begingroup$ I think you are referring to dataS1, well yes, FS1 is a mixture distribution.... $\endgroup$
    – jpcgandre
    Mar 4 '20 at 19:56
  • 1
    $\begingroup$ Yes, I mean that fS1n evaluates to an interpolating function and that is what has the discontinuities. $\endgroup$
    – Mark R
    Mar 4 '20 at 20:24
  • $\begingroup$ I had to do it because the way FS1 is built.... I couldn't get mathematica to work when determining FS12 if I determined fS1n by the derivative of FS1 or if I defined fS1[x_]:=(1-p1)*DiracDelta[x]+p1*fS11. But it works with SmoothKernelDistribution or EmpiricalDistribution... Just awful slowly $\endgroup$
    – jpcgandre
    Mar 5 '20 at 0:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.