chantel7777777
2014-May-15 17:03 UTC
[R] Markoc chain simulation of a sample path - to get empirical probabilities
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.
--
View this message in context:
http://r.789695.n4.nabble.com/Markoc-chain-simulation-of-a-sample-path-to-get-empirical-probabilities-tp4690658.html
Sent from the R help mailing list archive at Nabble.com.
Rolf Turner
2014-May-16 00:01 UTC
[R] Markoc chain simulation of a sample path - to get empirical probabilities
(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.
Jim Lemon
2014-May-16 09:25 UTC
[R] Markoc chain simulation of a sample path - to get empirical probabilities
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