# Plotting function takes forever

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

• Have you considered possibly rephrasing your problem using tools such as (N)Probability, (N)Expectation, Distributed and Conditioned? Mar 4, 2020 at 17:37
• 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. Mar 4, 2020 at 17:57
• @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).... Mar 4, 2020 at 17:58
• 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. Mar 5, 2020 at 7:31
• 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. Mar 5, 2020 at 8:39

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

• 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}];... Mar 5, 2020 at 10:39
• Apparently... when I use fS1[x_]:=(1-p1)*DiracDelta[x]+p1*fS11 it works super fast... Mar 5, 2020 at 11:08
• NIntegrate is a case where you do not want to use Set (or SetDelayed+Evaluate)... Mar 5, 2020 at 11:10
• 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 ;) Mar 5, 2020 at 11:11
• Sorry, I am quite confused and actually also busy, you know... But if it works for you, it is fine for me. =) Mar 5, 2020 at 11:28

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:

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

• That's good. I'll try in my machine. Thank you! Mar 4, 2020 at 18:16
• 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.} Mar 4, 2020 at 19:46
• I think you are referring to dataS1`, well yes, FS1 is a mixture distribution.... Mar 4, 2020 at 19:56
• Yes, I mean that fS1n evaluates to an interpolating function and that is what has the discontinuities. Mar 4, 2020 at 20:24
• 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 Mar 5, 2020 at 0:10