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.

In CrossValidated, the answer to a question on prediction of number of viewers of Gangnam Style is given in R. I have mapped the answer directly to Mathematica:

S0 = 10000000; a = 0.00000005; b = 0.01; k = 1.2;
views = {{0, 0}}; S = S0; inf = 1;
For[t = 1, t < 365, t++, dS = -a*inf*S;
  dinf = a*inf*S - b*inf;
  S = S + dS;
  inf = inf + dinf;
  noviews = Last[Last[views]] + k inf; 
  views = Append[views, {t, noviews}]];
ListPlot[views]

Result from function

However, this doesn't really feel or look like good Mathematica code. I tried to follow the advice in Alternatives to procedural loops and iterating over lists in Mathematica. But I get stuck in trying to create the function that the FoldList should use. It seems like the solution would be to have FoldList create a number of lists, one for each "state" variable.

So my question is: how do I rewrite the For loop that has a number of "state" variables and where the output depends on both these variables and earlier values of the output?

share|improve this question
add comment

3 Answers 3

up vote 12 down vote accepted

I wouldn't try to force FoldList, if there are better suited built-ins like RecurrenceTable

We could try to simplify your equations first

neweqs = {
     views[n + 1] == views[n] + k inf[n + 1], 
     inf[n + 1] == inf[n] + dinf, 
     dinf == a inf[n] S[n] - b inf[n], 
     S[n + 1] == S[n] + dS, 
     dS == -a inf[n] S[n]
     }~Eliminate~{dS, dinf} /. And -> List // Simplify

to get

(* {b inf[n] + inf[1 + n] == inf[n] (1 + a S[n]), 
   a inf[n] S[n] + S[1 + n] == S[n], 
   k inf[1 + n] + views[n] == views[1 + n]} *)

So now we can write

Block[{S0 = 10000000, a = 0.00000005, b = 0.01, k = 1.2,
  view0 = {{0, 0}}, inf0 = 1, max = 365, views, inf, S, n},

 view2 = RecurrenceTable[neweqs~Join~{
       views[0] == Last@Last@view0, inf[0] == inf0, S[0] == S0},
     {views[n]}, {n, 1, max},
     DependentVariables -> {inf, S, views}] /. 
    vals_ :> Transpose@{Range[max], Flatten@vals};
 ]
ListPlot @ view2

and get the same output as your example

For a maximum of 30000 I get in my machine an absolute timing of

your version: 5.6 s

this version: 0.017 s

So, at least for that size, thats a 330× speedup.

share|improve this answer
add comment

UPDATE: With FoldList and NestList one can eliminate the need for Join at every iteration and get much simpler and faster forms:

 foldListVws = Transpose[{Range[0, 364],
  (FoldList[{#[[1]] + k #[[3]] (1 + a #[[2]] - b),
     #[[2]] (1 - a #[[3]]), #[[3]] (1 + a #[[2]] - b)} &,
   {0, S0, 1}, Range[0, 364]][[All, 1]] //  Most)}]; 
 foldListVws == views
 (* True *)

 nestListVws =  NestList[{#[[1]] + 1, #[[2]] + k #[[4]] (1 + a #[[3]] - b),
     #[[3]] (1 - a #[[4]]), #[[4]] (1 + a #[[3]] - b)} &,
    {0, 0, S0, 1}, 364][[All, ;; 2]]; 
 nestListVws == views
 (* True *)

Using max= 30 000 as in @Rojo's post, the timings I get for the five methods (on a Windows Vista 64bit laptop with Intel Core2 Duo 2.80GHz CPU and 8G memory) are:

enter image description here

More readable update functions for FoldList and NestList:

ClearAll[nstLstUF, fldLstUF];
nstLstUF =Function[{t,x,y,z}, {t+1, x+k z (1 + a y - b), y (1 - a z), z (1 + a y - b)}];
fldLstUF = Function[{x,y,z}, {x + k z (1 + a y - b), y (1 - a z), z (1 + a y - b)}]; 

... and almost-one liners:

fldLstVws2 = Transpose[{Range[0, 364],
  (FoldList[fldLstUF[Sequence@@#]&, {0, S0, 1}, Range[0, 364]][[All, 1]] // Most)}];
nstLstVws2 = NestList[nstLstUF[Sequence @@ #]&,{0, 0, S0, 1}, 364][[All, ;; 2]];
fldLstVws2 == nstLstVws2 == views
(* True *)

Original post:

 foldVws =First@Fold[{Join[#1[[1]],
    {{#2, Last@Last[#1[[1]]] + k #1[[3]] (1 + a #1[[2]] - b)}}],
   #1[[2]] (1 - a #1[[3]]), 
   #1[[3]] (1 + a #1[[2]] - b)} &,
 {{{0, 0}}, S0, 1},
 Range[364]];
 foldVws==views
 (* True *)

Alternatively, you can use Nest:

 nestVws = Nest[{Join[#1[[1]],
  {{#1[[4]], Last@Last[#1[[1]]] + k #1[[3]] (1 + a #1[[2]] - b)}}],
  #1[[2]] (1 - a #1[[3]]), 
  #1[[3]] (1 + a #1[[2]] - b), #1[[4]] + 1} &,
  {{{0, 0}}, S0, 1, 1}, 364] // First;
 nestVws==views
 (* True *)
share|improve this answer
add comment

I suggest using Fold instead of FoldList. The conversion from For loop to Fold loop is algorithmic, could be done even programmatically. The definiton of Fold is:

Fold[f, x, {a, b, ...}] returns the last of the list {x, f[x,a], f[f[x, a], b], ...} (that is, it returns the last element of FoldList[f, x, {a, b, ...}]).

First, identify the function f that you want to use in Fold.

(* parameters *)
S0 = 10000000; a = 0.00000005; b = 0.01; k = 1.2;
(* initial values *)
{views, S, inf} = {{{0, 0}}, S0, 1};
update[{S_, inf_, views_}, t_] := Module[{dS, dinf, noviews},
   dS = -a*inf*S;
   dinf = a*inf*S - b*inf;
   noviews = Last[Last[views]] + k (inf + dinf);
   (* update will return the updated values in a list, as it accepts the same first argument *)
   {S + dS, inf + dinf, Append[views, {t, noviews}]}
   ];
For[t = 1, t < 365, t++,
  {S, inf, views} = update[{S, inf, views}, t]
  ];

Now it is easy to convert it to Fold. Define t as a range of values that the function will iterate through.

{views, S, inf} = {{{0, 0}}, S0, 1}; (* these must be reset *)
ListPlot@Last@Fold[update, {S, inf, views}, Range[1, 364]]

Mathematica graphics

share|improve this answer
add comment

Your Answer

 
discard

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.