Markov chain simulation of a sample path - to get empirical probabilities

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Markov chain simulation of a sample path - to get empirical probabilities

chantel7777777
This post was updated on .
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.
Reply | Threaded
Open this post in threaded view
|

Re: 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.

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Markoc chain simulation of a sample path - to get empirical probabilities

Jim Lemon
In reply to this post by chantel7777777
On Thu, 15 May 2014 10:03:03 AM chantel7777777 wrote:

Dear Chantel,
The problem:
A perplexed student decides to post her/his (note non-sexist usage)
homework to the R help list on a Friday evening. On Monday morning,
she finds that she/he has received five different answers. Realizing that
some might be incorrect, she/he persuades her/his smart friend to
check them. Sadly, the friend, while very smart, is also very lazy and
only does one now and then. The probabilities of the number of the
answers judged incorrect each day are:

Incorrect Probability
0 0.25
1 0.3
2 0.25
3 0.2

Whenever she/he gets this information, she/he strikes off those
answers from the list. If she/he has one or no remaining answers at
the end of each day, she/he posts the question to the R help list that
evening, and receives five more answers the next morning. What is
the most probable value for the number of remaining answers that
she/he will have on Friday morning?
 First you need the transition matrix. Assuming that the steps are from
morning to morning:

pmat<-
 matrix(c(0.25,0,0,0.45,0.3,
 0.3,0.25,0,0.2,0.25,
 0.25,0.3,0.25,0,0.2,
 0.2,0.25,0.3,0.25,0,
 0,0.2,0.25,0.3,0.25),
 nrow=5,ncol=5,byrow=TRUE)
 rownames(pmat)<-colnames(pmat)<-2:6
pmat

Since she/he began with five answers, we begin with state 4. Notice
here that since each step is from morning to morning and the R help
list is amazingly reliable, there are no 0 or 1 states.

initstate<-c(0,0,0,1,0)

One method of calculating the distribution after a number of steps is
to multiply the initial state by the transition matrix raised to the power
of the number of steps, so:

initstate%*%pmat%*%pmat%*%pmat%*%pmat

While this produces the correct answer, it is not the correct method.
You want to collect a number of probabilistic outcomes and use this to
estimate the most likely value of the number of potential answers
remaining on Friday morning.

outcomes<-rep(NA,100)
for(round in 1:100) {
 start=4
 for(i in 1:4) start<-sample(1:5,1,prob=pmat[start,])
 outcomes[round]<-start
}
table(outcomes)

So our hypothetical student might be able to work out from this how to
correct her/his code.

Jim

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.