2
$\begingroup$

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?

$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – Sjoerd Smit Feb 8 '18 at 14:57
3
$\begingroup$

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}

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Fierce82 Feb 8 '18 at 15:48
  • 1
    $\begingroup$ Tom, you're welcome. If you are interested in UnitStep you might also have a look at Clip for later use. $\endgroup$ – Henrik Schumacher Feb 8 '18 at 15:50
  • $\begingroup$ can u explain the logic behind last line? $\endgroup$ – Fierce82 Feb 8 '18 at 15:58
  • 1
    $\begingroup$ Better? Or do you mean the occurence of Divide? $\endgroup$ – Henrik Schumacher Feb 8 '18 at 16:06
  • 1
    $\begingroup$ 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). $\endgroup$ – Henrik Schumacher Feb 8 '18 at 16:26

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.