# SOLUTION TO STA 247, ASSIGNMENT #4. # The matrix of transition probabilities. Row i, column j in this matrix # is the probability of transitioning to state j when in state i. Trans <- matrix( c(0.85, 0.05, 0.05, 0.00, 0.05, 0.20, 0.60, 0.10, 0.00, 0.10, 0.15, 0.05, 0.50, 0.25, 0.05, 0.30, 0.05, 0.05, 0.00, 0.60, 0, 0, 0, 0, 1), 5, 5, byrow=TRUE) # Initial state probabilties. The chain always starts in state 1. init <- c (1, 0, 0, 0, 0) # FUNCTIONS TO FIND THE MARGINAL DISTRIBUTION OF THE STATE AT A GIVEN TIME. # First version. Due to R's funny idea of what : does, it won't work # when i is zero. marg1 <- function (i) { r <- init for (t in 1:i) r <- r %*% Trans r } # Second version. Note the fudge needed to get this to work when i is one. marg2 <- function (i) { Ti <- Trans if (i>1) for (t in 2:i) Ti <- Ti %*% Trans init %*% Ti } # PRINT THE MARGINAL DISTRIBUTIONS REQUESTED. cat("First version:\n\n") for (i in c(1,2,10,20)) cat ("Marginal at time",i,"is",marg1(i),"\n") # Try the second version. cat("\n\nSecond version:\n\n") for (i in c(1,2,10,20)) cat ("Marginal at time",i,"is",marg2(i),"\n") # BONUS QUESTION - FIND THE EXPECTED NUMBER OF PURCHASES. # # Let H be the number of purchases made, which is the number of times that # the customer is in state 4. # # We find this by finding all of H_1, H_2, H_3, and H_4, where H_i is the # expected number of purchases if they start in state i, for i = 1, 2, 3, 4. # Since customers actually start in state 1, we know that H = H_1. Of course, # the customers don't make any purchases after entering state 5. # # We can use the formula E(H|X_0=i) = E(E(H|X_1)|X_0=i) + d, where d is 1 if # i=4 and 0 otherwise. (This might seem undefined, since P(X_0=i) is zero # except for i=1, but it's meaningful if we imagine that the other values of # X_0 have some tiny non-zero probability, which won't affect the numerical # results.) This just says that the expected number of purchases starting # in state i is the average of the expected number starting in other states, # weighting the following state by how probable it is to move to that state # from state i, plus one if i=4, since in that case the customer is making a # purchase right then. # # Since E(H|X_1=j) = H_j, this gives a set of simultaneous linear equations # for H_1, H_2, H_3, and H_4. # # H_1 = 0.85 H_1 + 0.05 H_2 + 0.05 H_3 + 0.00 H_4 + 0 # H_2 = 0.20 H_1 + 0.60 H_2 + 0.10 H_3 + 0.00 H_4 + 0 # H_3 = 0.15 H_1 + 0.05 H_2 + 0.50 H_3 + 0.25 H_4 + 0 # H_4 = 0.30 H_1 + 0.05 H_2 + 0.05 H_3 + 0.00 H_4 + 1 # # The following matrix operations solve this system of equations. (Note # that diag(4) produces a 4x4 identity matrix.) T4 <- Trans[1:4,1:4] E4 <- solve(diag(4)-T4) %*% c(0,0,0,1) cat("\n\nExpected number of purchases:",E4[1],"\n") # Here's a brute force approximate check on the above answer... Epur <- 0 for (i in 1:1000) Epur <- Epur + marg1(i)[4] cat("Expected number of purchases (brute force check):",Epur,"\n")