# Improve speed in Sum of probabilities ~ Normal Distribution

I am having some performance problems, and I need some advice. I understand that my code is not well written

l = 86.;
m = 76.;
SD = 14.;

function1a[x_] :=
CDF[NormalDistribution[l - 0.5, SD], x] -
CDF[NormalDistribution[l + 0.5, SD], x];

factor1 = Sum[function1a[i], {i, 0., 200.}];

function1[x_] := If[x < 0., 0., function1a[x]*1./factor1];

function2a[x_] :=
CDF[NormalDistribution[m - 0.5, SD], x] -
CDF[NormalDistribution[m + 0.5, SD], x];

factor2 = Sum[function2a[i], {i, 0., 200.}];

function2[x_] := If[x < 0., 0., function2a[x]*1./factor2];

index[u_] := If[u < 0, -1., 1.];
chi[x_] := If[x < 0., Ceiling[Abs[x]], Floor[x]];
function[x_, y_] := function1[x]*function2[y];
Sum[function[i, j], {i, 0., 200.}, {j, 0.,
i + index[-10.5]*chi[-10.5]}] // Timing


What I do is, having two normal distributed independent variables, I am doing a continuity correction, then I make atransformation in negative values and the I calculate the joint PDF. Finally I need to calculate the given sum, but it takes more than 2 seconds to do so. ANy idea how to optimize it?

• For one, you can use the fact that CDF works on lists. So instead of doing Sum[function1a[i], {i, 0., 200.}], you can simply do Total[function1a[Range[0., 200.]]], which is significantly faster. – Sjoerd Smit Feb 8 '18 at 14:57

You compute sums over j and function evaluations of function1 and function2 over and over again. We may exploit that CDF works on lists (thanks to @SjoerdSmit) and we may get rid of some expensive Ifs by using UnitStep.

Try this:

RepeatedTiming[
iran = Range[0., 200.];
cut = Floor[index[-10.5]*chi[-10.5]];
vals1 = UnitStep[iran] function1a[iran];
vals2 = UnitStep[iran] function2a[iran];
Divide[
vals1[[-cut + 1 ;;]].Accumulate[vals2[[1 ;; cut - 1]]],
Total[vals1] Total[vals2]
]
]


{0.00043, 0.489927}

Here I use

$$\sum_{i=1}^{n} \sum_{j=1}^{m+i} (a_i \, b_j) = \sum_{i=1}^n a_i S_i = \langle a, S \rangle$$

with $S_i = \sum_{j=1}^{m+i} b_j$. The vector S is efficiently computed by Accumulate, probably by the recursive formula $S_{j+1} = S_j + b_j$.

Sum[function[i, j], {i, 0., 200.}, {j, 0., i + index[-10.5]*chi[-10.5]}] // RepeatedTiming


{1.81, 0.489927}

• Very elegant. I have a lot of new things to digest here. Appreciate it. I was just asking about Unitstep, but you already edited your answer. – Fierce82 Feb 8 '18 at 15:48
• Tom, you're welcome. If you are interested in UnitStep you might also have a look at Clip for later use. – Henrik Schumacher Feb 8 '18 at 15:50
• can u explain the logic behind last line? – Fierce82 Feb 8 '18 at 15:58
• Better? Or do you mean the occurence of Divide? – Henrik Schumacher Feb 8 '18 at 16:06
• vals1[[-cut + 1 ;;]] is just the part of vals1 that runs from indices -cut+1 to the end. See Span (that is the InputForm of ;;) in the docs. Similarly, vals2[[1 ;; cut - 1]] reads a part from vals. Since cut-1 is negative (-12) it counts the indices from the end; in out case 1 ;; cut - 1 is equivalent to Span[1,201+cut-1]. The . (itsInputForm is Dot) computes the inner product between two vectors (and also matrix-vector and matrix-matrix products). – Henrik Schumacher Feb 8 '18 at 16:26