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.
______________________________________________
R-help@r-project.org 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.
