STAT 131 Final Lab
STAT | 统计代写 | assignment – 这个项目是STAT代写，assignment代写的代写题目
Central Limit Theory and Markov Chain
In this section, we will simulate the behavior of central limit theory and Markov chain.
Visualizations of central limit theory
Suppose our distribution of interest isUnif(0,1), we are interested in the distribution of X 10 , i.e., n = 10is the selected sample size. We create a matrix with 10 columns, each row represents one sample.
n = 10 # number of samples nsim = 104 # repeat the sampling for 10000 times x = matrix(runif(nsim*n), nrow=nsim, ncol=n) # each row is a replicate of 10 samples
To calculate X 10 , we simply take the average of each row of the matrixx, we can do this with therowMeans orapplyfunction:
xbar = rowMeans(x)
Finally, we create a histogram of X 10 :
hist(xbar)
Histogram of xbar
xbar
Frequency
0.2 0.3 0.4 0.5 0.6 0.7 0.
0
500
1000
1500
2000
Figure 1: Histogram of the sample mean.
Markov chain: Chutes and Ladders
In the game called Chutes and Ladders , also referred as Snakes and Ladders, players try to be first to reach a ceritain distribution on a board. The board is a grid of squares, numbered from 1 to the number of squares.
The board has some chutes and some ladders, each of which connects a pair of squares. Here we will
consider the one player version of the game.
Figure 2: Chutes and Ladders
We consider a simplified version of this game on the 3 3 board shown in figure 2. The player starts out at square 1, and wants to get to square 9. On each move, a fair coin is flipped, and the player gets to advance 1 square if the coint lands Heads and 2 squares if the coin lands Tails. There are 2 ladders (shown as upward-point arrow) and 2 chutes (shown as downward-pointing arrows) on the board. If the resulting square is the base of a ladder, e.g., square 5, the player gets to climb the ladder, instantly arriving at square 7. If the resulting square is the top of a chute, e.g., square 8, the player instantly slides down to the bottom of the chute, which is 1. In fact, Chutes and Ladders can be represented as a Markov chain with transition matrix Q as follows:
Q =
0 0. 5 0. 5 0 0
0 0 0. 5 0. 5 0
0. 5 0 0 0. 5 0
0. 5 0 0 0 0. 5
0 0 0 0 1
.
The first state is the initial state (starting point), and the final state is called an absorbing state because
from this state the chain can never leave.
Q = matrix(c(0, 0.5, 0.5, 0, 0,
0, 0, 0.5, 0.5, 0,
0.5, 0, 0, 0.5, 0,
0.5, 0, 0, 0, 0.5,
0, 0, 0, 0, 1),
nrow=5,ncol=5,byrow=TRUE)
Note that we typed in the entries by row so we need to setbyrow=TRUE, and we need to tell R the number of
rows (nrow) and the number of columns (ncol) in the matrix.
To obtain higher order transition probabilities, we can multiply Q by itself repeatedly. The matrix multiplica-
tion command in R is%*%( not just*). So
Q2 = Q %*% Q
Q3 = Q2 %*% Q
Q4 = Q3 %*% Q
Q5 = Q4 %*% Q
produces Q^2 to Q^5. If we want to know the probability of going from state 4 to state in exactly 5 steps, we
can extract the(3 , 4)entry of Q^5 :
Q5[4,5]
## [1] 0.
To compute a power Qn without directly doing repeated matrix multiplications, we can use the command
Q %% nafter installing and loading theexpmpackage. For example,Q %% 42yields Q^42. By exploring
the behavior of Qn as n grows, we can see Qn converges to a matrix in which each row is the stationary
distribution s (theorem 11.3.6), and get a sense of how long it takes for the chain to get very close to its
stationary distribution.
library(expm)
Q %% 100
## [,1] [,2] [,3] [,4] [,5]
## [1,] 4.423778e-07 2.534345e-07 3.986250e-07 3.735593e-07 0.
## [2,] 3.735593e-07 2.140090e-07 3.366129e-07 3.154466e-07 0.
## [3,] 3.986250e-07 2.283689e-07 3.591995e-07 3.366129e-07 0.
## [4,] 2.534345e-07 1.451905e-07 2.283689e-07 2.140090e-07 0.
## [5,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.
Another way to obtain the stationary distribution numerically is to useeigen(t(Q))is compute the eigenvalues
and eigenvectors of the transpose of Q ; then the eigenvector corresponding to the eigenvalue 1 can be selected
and normalized so that the components sum to 1.
decomp = eigen(t(Q))
eigval = decomp$values
eigvec = decomp$vectors
s = eigvec[,eigval==1] # stationary distribution
Re(s) # show as real values
## [1] 0 0 0 0 1
We call state 5 the absorbing state, since from this state the chain can never leave.
Simulating from a finite-state Markov chain
To simulate a 4-state Markov chain, first we need to define a 4 4 transition matrix Q = ( qij ), where
qij = P ( Xn = j | Xn 1 = i ):
M = 4 # the number of states
Q = matrix(c(1/3, 1/3, 1/3, 0,
0, 0, 1/2, 1/2,
0, 1, 0, 0,
1/2, 0, 0, 1/2),
nrow=M, ncol=M,
byrow=TRUE)
We store the realized values of this Markov chain to a vectorx. x[1]is the initial value that randomly chooses a value from 1 to M. We observe this Markov chain forntimetotal time points (including the initial point). ntime = 10 x = rep(0,ntime) x[1] = sample(1:M,1)
At timet,x[t]randomly choose a value j { 1 ,... , M }based on a conditional probability{ P ( xt = j | xt 1 =
i ) , i = 1 ,... , M }, stored atQ[x[t-1],].
for (t in c(2:ntime)) {
x[t] = sample(1:M,1,prob=Q[x[t-1],])
}
We can plot the tail of the simulated chain to see if it reaches stationarity.
plot(tail(x,200),type="l",xlab="",ylab="x")
0 50 100 150 200
1.
1.
2.
2.
3.
3.
4.
x
Figure 3: Simulated Markov chain after reaching stationarity. Here we plot the last 200 simulated samples.
Since we setntimeto a large number, it may be reasonable to believe that the chain is close to stationarity
during the latter portion of the simulation. To check this, we eliminate the first half of the simulations to
give the chain time to reach stationarity:
x = x[-( 1:(ntime/2) )]
We then use thetablecommand to calculate the number of times the chain visited each state; dividing bylength(x)converts the counts into proportions. The result is an approximation to the stationarity distribution. table(x)/length(x)
## x
## 1 2 3 4
## 0.2208 0.2840 0.2034 0.
The Lab
The assignment in this case is to run the attached code, and answer the questions included at the end. In
order to have all the points from the assignment you need to:
(1) Include the output of the code as proof that the code ran correctly.
(2) Write a complete and correct solution to the questions included afterward.
(3) Upload this as a pdf on Canvas.
Questions:
- Visualizations of central limit theory. Compare the shape of the histograms when sample size n=5,10,30,100, what do you observe? Also tryBin(10 , 0_._ 6)(rbinom) andBeta(0_._ 5 , 0_._ 1)(rbeta), does the central limit theory still holds? Comment on why it does or doesnt. Note that in total there should be 12 plots, and three comments, one per distribution (Uniform, Binomial and Beta). (4pts)
- Markov chain: Chutes and Ladders.
- Explain why, despite the fact that there are 9 squares, we can represent the game using the 5 5 transition matrix Q as shown above. (1pt)
- Find the expected the number of times the player reaches state 4, starting from state 1, denoted by n 14. (1pt)
- Find the expected duration of the game, where duration is the number of coin flips to reach state 5, denoted by N. Hint: before reaching state 5, we visit all other states, 2 , 3 , 4 , as many times as necessary,{ n 12 , n 13 , n 14 }. (1pt)
- Simulating from a finite-state Markov chain. Solve question 24 in Chapter 11, page 532. You can solve it by hand or with R. (3pts) – Extra credits: plot the simulated Markov chain after it reaches stationarity as shown in the above figure 3. (2pts)