2
$\begingroup$

I want to generate the random sequence defined by $a_n = a_{n-2} + \beta(n)a_{n-1}$ where $\beta$ takes the values $\pm1$. My attempt is:

 RecurrenceTable[{a[n] == a[n - 2] + (-1)^(RandomInteger[{0, 1}])*a[n - 1], 
 a[0] == 1, a[1] == 1}, a, {n, 0, 10}]

But of course, RandomInteger only computes on the first loop, and then saves that value for all subsequence computations. I found this solution here:

 r[n_] := RandomInteger[{0, 1}];
 rfib[0] = 1;
 rfib[1] = 1;
 rfib[n_] := rfib[n] = rfib[n - 2] + (-1)^r[n]*rfib[n - 1];
 Table[rfib[i], {i, 0, 10}]

However, I am wondering if my original intuition can be salvaged; i.e., is there a way to use RandomInteger inside RecurrenceTable and get a new integer each time?

Edit:

After posting this and looking for something else, I found this post. From this, I can write the code as 

 rr[n_?NumericQ] := RandomInteger[{0, 1}];
 RecurrenceTable[{a[n] == a[n - 2] + (-1)^(rr[n])*a[n - 1], 
 a[0] == 1, a[1] == 1}, a, {n, 0, 10}]

as I originally wanted to.

Improve this question
$\endgroup$
1
  • $\begingroup$ Where does this come from? Is it meant to model some particular physical problem? $\endgroup$
    – bill s
    Sep 15, 2017 at 12:34

4 Answers 4

Reset to default
3
$\begingroup$

One way to approach this is to directly implement the recursion:

Clear[a];
a[n_] := a[n] = RandomChoice[{-1, 1}] a[n - 1] + a[n - 2];
a[0] = a[1] = 1;

For example, the first 20 terms might be:

a /@ Range[20]
{1, 2, 1, -1, -2, -1, 1, 2, 1, 3, 2, -1, -3, -4, -1, 3, 4, 7, 3, -4}
Improve this answer
$\endgroup$
2
  • $\begingroup$ Doesn't RandomInteger[{-1,1}] allow for the choice of zero? Other than that question, I see how this works. As a side note, I want my subtractions to go left to right, so I would put the random $\pm1$ with the $n-1$ term. $\endgroup$
    – Trevor
    Sep 15, 2017 at 11:53
  • $\begingroup$ OK -- I've changed the placement of the randomness, and used RandomChoice instead of RandomInteger. $\endgroup$
    – bill s
    Sep 15, 2017 at 12:34
1
$\begingroup$

Just another way:

func[a0_, a1_, n_] := Module[{r = RandomChoice[{-1, 1}, n]},
   FoldList[{#1[[2]], #1[[1]] + #2 #1[[2]]} &, {a0, a1}, r]] [[All, 
   1]]

Examples:

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ I was looking at FoldList but my ultimate goal was to generate the arcsine distribution that occurs when one varies the heads/tails probability of random Fibonacci sequences (when the coin is fair, you get Viswanath's number as the growth rate, and when the coin is 100% heads, you get the Golden Ratio,) so I needed a clean/fast way to get very large terms. FoldList just didn't feel right to me for that. Your example is quite nice, though. $\endgroup$
    – Trevor
    Sep 15, 2017 at 11:55
  • 1
    $\begingroup$ @Trevor "Doesn't feel nice" isn't really an argument. :P Is there an actual practical problem you expect? Performance? Then I would make a compiled function, but still retain this basic design. I did not pay attention when I posted my answer (as you can see it's very similar), but now that I read ubpdqn's carefully, I think his should perform better as it generates the random number in bulk (instead of calling RandomChoice in each iteration) $\endgroup$
    – Szabolcs
    Sep 15, 2017 at 13:04
1
$\begingroup$

I would use NestList or Nest:

NestList[
  {Last[#], {1, RandomChoice[{-1, 1}]}.#} &,
  {1, 1}, (* two initial values in the sequence *)
  100
][[All, 2]]

This is directly compilable:

cf = Compile[{{a0, _Integer}, {a1, _Integer}, {n, _Integer}},
  NestList[
    {Last[#], {1, RandomChoice[{-1, 1}]}.#} &,
    {a0, a1},
    n
    ][[All, 2]]
  ]

There is no MainEvaluate in the compiled function, but it will switch back to standard evaluation as soon as the values exceed $2^{63}-1$ (on a 64-bit machine).

Improve this answer
$\endgroup$
0
$\begingroup$

One can exploit the equivalence of evaluating three-term recurrence relations with repeated multiplication of $2\times 2$ matrices for this task. Just like in ubpdqn's solution, I use RandomChoice[] to generate a bunch of $\pm1$ multipliers all at once:

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"]; (* for reproducibility *)
            With[{n = 20}, 
                 FoldList[Dot, {{0, 1}, {1, 1}}, 
                          Transpose[{{ConstantArray[0, n], ConstantArray[1, n]}, 
                                     {ConstantArray[1, n], RandomChoice[{-1, 1}, n]}},
                                    {2, 3, 1}]][[All, 2, 2]]]]
   {1, 2, -1, 3, -4, 7, -11, -4, -15, -19, -34,
    -53, -87, 34, -121, 155, -276, -121, -397, -518, 121}

Compare this with using RecurrenceTable[]:

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"];
            With[{n = 20}, pm = RandomChoice[{-1, 1}, n];
                 Quiet[RecurrenceTable[{y[k] == Indexed[pm, k] y[k - 1] + y[k - 2], 
                                        y[-1] == 1, y[0] == 1},
                                       y, {k, 0, n}], Indexed::partw]]]
Improve this answer
$\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.