2

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.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Trevor
  • 53
  • 3

4 Answers4

3

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}
J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
bill s
  • 68,936
  • 4
  • 101
  • 191
  • 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. – Trevor Sep 15 '17 at 11:53
  • OK -- I've changed the placement of the randomness, and used RandomChoice instead of RandomInteger. – bill s Sep 15 '17 at 12:34
1

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

ubpdqn
  • 60,617
  • 3
  • 59
  • 148
  • 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. – Trevor Sep 15 '17 at 11:55
  • 1
    @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) – Szabolcs Sep 15 '17 at 13:04
1

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).

Szabolcs
  • 234,956
  • 30
  • 623
  • 1,263
0

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]]]
J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574