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 am trying to construct a recursion, g[t], where

  • g[t] = a[t] if g[t-1] = 1
  • g[t] = b[t] if g[t-1] = 2.

Additionally a and b are defined as:

a[t_] := If[0 <= RandomReal[] <= 0.94, 1, 2]
b[t_] := If[0 <= RandomReal[] <= 0.9, 2, 1]

So then I want (I've also tried it where I replace a and b with their original defintions):

g[t_]:= Piecewise[{{a[t], g[t - 1] = 1}, {b[t], g[t - 1] = 2}}]
g[1] = 1

This, however, is not working.

share|improve this question
3  
Welcome here. You are aware that g[t-1]=1 does not mean is equal to? Therefore, you first two bullet points don't make sense to me, because after the if I expect a true/false. –  halirutan Sep 30 '13 at 12:27
    
My first two bullet points are not actually written into my code; I was just trying to give a general idea as to what I'm trying to do. The parts in boxes are what I'm actually working with. But do you mean that my conditions, g[t-1] = 1 or 2, defined in my piecewise function do not make sense? I've noticed that if i choose a simple condition, like t > 1, then my piecewise function works. –  Gigi Galwah Sep 30 '13 at 15:36
    
Possible duplicate: 18697. –  Jacob Akkerboom Oct 1 '13 at 15:48
    
In my opinion, that one is indeed duplicate. However, maybe it is not too localized, as I guess this is a great place to learn for the first time the difference between == and =. I suggest we close this Q&A (but I'm still glad I got to write an answer!) –  Jacob Akkerboom Oct 1 '13 at 15:50
add comment

2 Answers

up vote 2 down vote accepted

As halirutan has pointed out, the main problem with your code is the equals sign that in the second elements (conditions) of the lists in the list in the Piecewise function. This should be ==.

The t in a[t] and b[t] suggests to me that you want to store values of these random sequences, rather than generating them over and over.

Also you can use RandomChoice, rather than your construction with RandomReal.

To store your values in the code below, I have used something that looks like a technique called memoization. For a and b, this changes the outcome, as things are generated only once, so it is not really memoization. But given a and b, calculating g from its stored previous value rather than going all the way back to t=1 is proper memoization. Note that in g[t_]:=g[t]=body the brackets should be read to be g[t_]:= (g[t]=body). Search this site for memoization if you are not sure how it works :).

The code then becomes

ClearAll[a, b, g]
a[t_] := a[t] = RandomChoice[{0.94, 0.06} -> {1, 2}]
b[t_] := b[t] = RandomChoice[{0.9, 0.1} -> {1, 2}]
g[t_] := g[t] = 
  Piecewise[{{a[t], g[t - 1] == 1}, {b[t], g[t - 1] == 2}}]
g[1] = 1;

Note that more improvements can be made. For example, it is much faster to first generate say 100000 random variables all at once, than to generate them on the fly. But maybe that is something for a future question.

share|improve this answer
    
I did not search for duplicates, which I am sure there are. Do not hesitate to close this Q&A or whatever, do not let the presence of my answer or any upvotes for it discourage you to close/delete :P. –  Jacob Akkerboom Sep 30 '13 at 16:44
add comment

I originally posted an incorrect post. This may have been better as a comment and does not address the code correct. However, I note that if the aim is a simulation, this two state discrete Markov process staring in state 1 could be coded using Mathematica builtin functions:

markov = DiscreteMarkovProcess[{1, 0}, {{0.94, 0.06}, {0.1, 0.9}}]

Underlying graph:

el = EdgeList[Graph[markov]];
prop = PropertyValue[{Graph[markov], #}, "Probability"] & /@ el;
Graph[markov, EdgeLabels -> Thread[el -> prop], EdgeLabelStyle -> 12]

enter image description here

Random function (e.g 10 samples):

ranf=RandomFunction[markov,{0,10}]

This creates a TemporalData object

Do probability, probability in state 2 after 100 time steps

 Probability[x[100] == 1, x \[Distributed] markov]

(=0.625)

Compute Stationary Distribution:

stat = StationaryDistribution[markov];
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.