[R] Markoc chain simulation of a sample path - to get empirical probabilities

[Rolf Turner]

(1) This list is *NOT* for helping people with homework.

(2) Your code is kludgy, inefficient, and apparently wrong.

(3) Think about what is going to happen in your very first call to 
sample() --- it will be:

     sample(6,1,prob=P[path[5],])

(Note that "replace=TRUE" is unnecessary and nonsensical here.)

But path[5] has not been set yet --- so it is NA.  Crash!!!

Think through what you are doing more carefully.

(4) Talk to your instructor!  I am *fairly* confident that he or she 
will not bite.

cheers,

Rolf Turner

On 16/05/14 05:03, chantel7777777 wrote:
> The question:
> A computer store sells iPods. If at the end of the day they have 0 or 1 unit
> of stock, they order enough new units so that their total number of units on
> hand is 5. New merchandise arrives before the store opens the next day.
> Demand each day is random with the following distribution
>
> Demand     Probability
>       0                0.3
>       1                0.4
>       2                0.2
>       3                0.1
>
> a) What is the most likely number of units in stock on Friday given that the
> store opened with 5 units of stock on Monday?
>
> My assignment is to write an R code that takes as input a transition matrix
> and an initial distribution, and then simulates data from which I can
> empirically calculate probabilities...
>
>
> I am confident that my transition matrix is correct, it is
>          0.0  0.0  0.1  0.2  0.4  0.3
>          0.0  0.0  0.1  0.2  0.4  0.3
> P=    0.3  0.4  0.3  0.0  0.0  0.0
>          0.1  0.2  0.4  0.3  0.0  0.0
>          0.0  0.1  0.2  0.4  0.3  0.0
>          0.0  0.0  0.1  0.2  0.4  0.3
> where the row1=o units of stock
>                 row2= 1 unit of stock....
>                 row3= 2 units of stock etc, and the columns also follow the
> same order
>
> (Manually I know how to do this and I have done so, but I cannot do it in R)
> I have the following code:
>
>> set.seed(1413974749)
>> mChainSimulation=function(P,pathLength, initDist)
> + {
> +   path=numeric(pathLength)
> +   path[1]=6 # Because on Monday we open with 5 units of stock (which is
> 6th state in Markov chain)
> +   for(i in 6:pathLength)
> +   {
> +     path[i]=sample(6, 1 , replace=TRUE,prob=P[path[i-1],])
> +   }
> +   return(path)
> + }
>> numStates=6
>> initDist=c(1/6,1/6,1/6,1/6,1/6,1/6)
>> P=matrix(c(0, 0, 0.1, 0.2, 0.4, 0.3, 0, 0, 0.1, 0.2, 0.4, 0.3, 0.3, 0.4,
>> 0.3, 0, 0, 0, 0.1, 0.2, 0.4, 0.3, 0, 0, 0, 0.1, 0.2, 0.4, 0.3, 0, 0, 0,
>> 0.1, 0.2, 0.4, 0.3),6,6, byrow=T)
>> P
>       [,1] [,2] [,3] [,4] [,5] [,6]
> [1,]  0.0  0.0  0.1  0.2  0.4  0.3
> [2,]  0.0  0.0  0.1  0.2  0.4  0.3
> [3,]  0.3  0.4  0.3  0.0  0.0  0.0
> [4,]  0.1  0.2  0.4  0.3  0.0  0.0
> [5,]  0.0  0.1  0.2  0.4  0.3  0.0
> [6,]  0.0  0.0  0.1  0.2  0.4  0.3
>> pathLength=4
>> simNum=10000
>> stock0Num=0 # number of times in the simulation that it ends on 0 units of
>> stock, i.e. in simulations it ends on 1
>> stock1Num=0 # number of times in the simulation that it ends on 1 unit of
>> stock, i.e in simulations it ends on 2
>> stock2Num=0 # number of times in the simulation that it ends on 2 units of
>> stock, i.e in simulations it ends on 3
>> stock3Num=0 # number of times in the simulation that it ends on 3 units of
>> stock, i.e in simulations it ends on 4
>> stock4Num=0 # number of times in the simulation that it ends on 4 units of
>> stock, i.e in simulations it ends on 5
>> stock5Num=0 # num
>
> !!!!!!!!!!!!!!!!!!! # I get an error here in the next few
> lines!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>
>> for(i in 1:simNum)
> + {
> +   pathSimulation=mChainSimulation(P,pathLength, initDist)
> +   pathEnd=pathSimulation[pathLength]
> +
> +   if(pathEnd==1)
> +   {
> +     stock0Num=stock0Num+1
> +   }
> +   if(pathEnd==2)
> +   {
> +     stock1Num=stock1Num+1
> +   }
> +   if(pathEnd==3)
> +   {
> +     stock2Num=stock2Num+1
> +   }
> +   if(pathEnd==4)
> +   {
> +     stock3Num=stock3Num+1
> +   }
> +   if(pathEnd==5)
> +   {
> +     stock4Num=stock4Num+1
> +   }
> +   if(pathEnd==6)
> +   {
> +     stock5Num=stock5Num+1
> +   }
> + }
>   Hide Traceback
>
>   Rerun with Debug
>   Error in sample.int(length(x), size, replace, prob) :
>    incorrect number of probabilities
> 3 sample.int(length(x), size, replace, prob)
> 2 sample(1:numStates, 1, prob = P[path[i - 1], 6])
> 1 mChainSimulation(P, pathLength, initDist)
>
>
> I don't know how to fix my error and would greatly appreciate any help.