0
$\begingroup$

I want to evaluate the following sum at lots of different values of $R_0$ and $N$ (roughly 100 of each, R0 ranging from 1 to 6, N ranging from 1000 to 10000):

$ \langle Y \rangle = \frac{\sum_{Y=1}^N \frac{(N-1)!}{(N-Y)!} \left( \frac{R_0}{N}\right)^{Y-1}}{\sum_{Y=1}^N \frac{(N-1)!}{Y(N-Y)!} \left(\frac{R_0}{N}\right)^{Y-1}}$.

I'd been using the following to do this:

solveSIS[R0_, NN_] :=
 Module[{pp1, meanpp},
  pp1 = 1/Sum[((NN - 1)!)*1/(Y*(NN - Y)!)*(R0/NN)^(Y - 1), {Y, NN}];
  meanpp = pp1*Sum[((NN - 1)!)*1/((NN - Y)!)*(R0/NN)^(Y - 1), {Y, NN}]]

Parallelize[Table[solveSIS[105/100 + i*5/100, 1000 + 100*j], {i, 99}, {j, 99}]]

Is there a faster way for me to go about doing this, as this took about two days without finishing before my laptop accidentally ran out of battery.

For example, the top and bottom sums are very similar, and a lot of the calculation may be wasted if it calculates them both separately. Is this something that mathematica would optimise by itself? If not, how could I go about optimising this?

$\endgroup$
2
$\begingroup$

I hope you don't need the exact fractions that come out of these calculations, because the numbers become so huge that Mathematica spends a lot of time working with them. I'll give a solution to machine precision.

The fastest way is to make use of vectorized operations (functions that act on whole lists at a time, rather than looping through individual elements) and combining these in suitable ways. Here is my proposal:

r0Vec = Table[105/100 + i*5/100, {i, 99}];
nVec = Table[1000 + 100*j, {j, 99}];

Table[
 Block[{yVec = N@Range[n], nFactors, terms},
  nFactors = (n - 1)!/(n - yVec)!*n^(1 - yVec);
  Table[
   terms = nFactors*r0^(yVec - 1);
   Total[terms]/Total[terms/yVec]
   , {r0, r0Vec}
   ]
  ]
 , {n, nVec}
 ]

I did not time the whole thing to complete, but if I let j run to 10 instead of 99, the whole thing completes in about 4 seconds on my 2013 MacBook Pro.

$\endgroup$

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.