Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm trying to reproduce the simulation with demographic stochasticity in Figure 1 from the paper entitled "Dynamical Resonance Can Account For Seasonality of Influenza Epidemics" (http://www.pnas.org/content/101/48/16915.full).

My code works just fine when $t$ is small $(0<t<0.02)$, but when I increase the value of $t$ (as below), the simulation doesn't finish compiling. Thus, I'm wondering if there's an issue with my code. Ideally, I would like my code to function when $(10<t<20)$ so I can compare my results to the paper. Any input would be greatly appreciated.

Clear[S, Infe, R, No, b0, b1, b, t, Ro, d, L, event, fun, main, list1, list2, data] 
No = 500000; 
b0 = 500; 
b1 = 0.04; 
d = 0.02; 
L = 4; 
list1 = {}; 
list2 = {}; 
b[t_] := b0 *(1 + b1* Cos[2 Pi t]); 
main[] := Module[{S = 499999, Infe = 1, R = 0, t = 0}, 
While[t < .03 && Infe > 0, {S, Infe, R, t} = fun[S, Infe, R, t];
AppendTo[list1, Infe]; AppendTo[list2, t] ]];

fun[S0_, Infe0_, R0_, t0_] := 
Module[ {S = S0, Infe = Infe0, R = R0, t = t0}, 
TransRate = (b[t]*S*Infe)/No; 
RecoveryRate = Infe/d; 
ImmunityLoss = (No - S - Infe)/L; 
TotalRate = TransRate + RecoveryRate + ImmunityLoss; 

t = t + RandomReal[ExponentialDistribution[TotalRate]]; 
event = RandomReal[UniformDistribution[]]; 

If[event < TransRate/TotalRate, 
S = S - 1 ; Infe = Infe + 1, 
If[event < (TransRate + RecoveryRate)/TotalRate, 
Infe = Infe - 1 ; R = R + 1, R = R - 1 ; S = S + 1]]; 
{S, Infe, R, t} 
   ]

main[] 
data = Transpose@{list2, list1};
ListPlot[data]
share|improve this question
1  
The value of t approaches your stopping limit of .03 from below slowly. The repeated Append is killing you. Change that operation to a Sow and Reap operation- that will speed things greatly (as in order of magnitude +). Then take a look at the top FAQ questions here to refactor this into a more functional/Mathematica style - the imperative style here is seldom efficient. Lastly, consider just doing it as a Markov chain perhaps - Mathematica has built-in capabilities for these that are quite efficient. –  rasher Jul 21 at 23:19
    
Also, consider pre-generating the random variates and indexing into them rather than the one-at-a-time way you're doing it now: in cases where there will be many iterations, the set-up/tear-down cost of making individual calls to Random... will get expensive. –  rasher Jul 21 at 23:26
    
A last thought (I don't have time today to re-write your example incorporating my thoughts above): for your desired range of 10<t<20 you'll on average iterate a quarter to a half billion times. Better to generate the variates in chunks, checking the running sum until your t condition is met, and vectorizing the remaining calculations where possible. –  rasher Jul 21 at 23:45
    
Thank you Rasher! Your advice helped tremendously. –  Omar Jul 23 at 2:12

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.