Thiết kế website giá rẻ

Question

I have a Markov chain transition matrix and I’m attempting to simulate client distribution starting from year 1. Clients transition between classes according to probabilities outlined in the transition matrix. Here’s the transition matrix:

transition_matrix <- matrix(c(1-p[1],p[1],0,0,0,0,0,0,0,0,
                              1-p[1],0,p[1],0,0,0,0,0,0,0,
                              1-p[1]-p[2],p[2],0,p[1],0,0,0,0,0,0,
                              1-p[1]-p[2]-p[3],p[3],p[2],0,p[1],0,0,0,0,0,
                              1-p[1]-p[2]-p[3]-p[4],p[4],p[3],p[2],0,p[1],0,0,0,0,
                              1-p[1]-p[2]-p[3]-p[4]-p[5],p[5],p[4],p[3],p[2],0,p[1],0,0,0,
                              1-p[1]-p[2]-p[3]-p[4]-p[5]-p[6],p[6],p[5],p[4],p[3],p[2],0,p[1],0,0,
                              1-p[1]-p[2]-p[3]-p[4]-p[5]-p[6]-p[7],p[7],p[6],p[5],p[4],p[3],p[2],0,p[1],0,
                              1-p[1]-p[2]-p[3]-p[4]-p[5]-p[6]-p[7]-p[8],p[8],p[7],p[6],p[5],p[4],p[3],p[2],0,p[1],
                              1-p[1]-p[2]-p[3]-p[4]-p[5]-p[6]-p[7]-p[8]-p[9],p[9],p[8],p[7],p[6],p[5],p[4],p[3],p[2],p[1]), nrow = 10,ncol = 10,byrow = TRUE)

Maybe you can write a function that simulates the distribution of customers after many years if I use different transition matrices.
I tried to simulate distribution manually with 82990 customers:

n <- 100
x9 <- numeric(n)
x8 <- numeric(n)
x7 <- numeric(n)
x6 <- numeric(n)
x5 <- numeric(n)
x4 <- numeric(n)
x3 <- numeric(n)
x2 <- numeric(n)
x1 <- numeric(n)
x0 <- numeric(n)
x9[1] = 82990 * (1 - p[1])
x8[1] = 82990 * p[1]
print(sum(c(x9[1], x8[1])))
x9[2] = x9[1] * (1 - p[1]) + x8[1] * (1 - p[1])
x8[2] = x9[1] * p[1]
x7[2] = x8[1] * p[1]
print(sum(c(x9[2], x8[2], x7[2])))
x9[3] = x9[2] * (1 - p[1]) + x8[2] * (1 - p[1]) + x7[2] * (1 - p[1] - p[2])
x8[3] = x9[2] * p[1] + x7[2] * p[2]
x7[3] = x8[2] * p[1]
x6[3] = x7[2] * p[1]
print(sum(c(x9[3], x8[3], x7[3], x6[3])))
x9[4] = x9[3] * (1 - p[1]) + x8[3] * (1 - p[1]) + x7[3] * (1 - p[1] - p[2]) + x6[3] * (1 - p[1] - p[2] - p[3])
x8[4] = x9[3] * p[1] + x7[3] * p[2] + x6[3] * p[3]
x7[4] = x8[3] * p[1] + x6[3] * p[2]
x6[4] = x7[3] * p[1] 
x5[4] = x6[3] * p[1]
print(sum(c(x9[4], x8[4], x7[4], x6[4], x5[4])))
x9[5] = x9[4] * (1 - p[1]) + x8[4] * (1 - p[1]) + x7[4] * (1 - p[1] - p[2]) + x6[4] * (1 - p[1] - p[2] - p[3]) + x5[4] * (1 - p[1] - p[2] - p[3] - p[4])
x8[5] = x9[4] * p[1] + x7[4] * p[2] + x6[4] * p[3] + x5[4] * p[4]
x7[5] = x8[4] * p[1] + x6[4] * p[2] + x5[4] * p[3]
x6[5] = x7[4] * p[1] + x5[4] * p[2]
x5[5] = x6[4] * p[1]
x4[5] = x5[4] * p[1]
print(sum(c(x9[5], x8[5], x7[5], x6[5], x5[5], x4[5])))
x9[6] = x9[5] * (1 - p[1]) + x8[5] * (1 - p[1]) + x7[5] * (1 - p[1] - p[2]) + x6[5] * (1 - p[1] - p[2] - p[3]) + x5[5] * (1 - p[1] - p[2] - p[3] - p[4]) + x4[5] * (1 - p[1] - p[2] - p[3] - p[4] - p[5])
x8[6] = x9[5] * p[1] + x7[5] * p[2] + x6[5] * p[3] + x5[5] * p[4] + x4[5] * p[5]
x7[6] = x8[5] * p[1] + x6[5] * p[2] + x5[5] * p[3] + x4[5] * p[4]
x6[6] = x7[5] * p[1] + x5[5] * p[2] + x4[5] * p[3]
x5[6] = x6[5] * p[1] + x4[5] * p[2]
x4[6] = x5[5] * p[1]
x3[6] = x4[5] * p[1] 
print(sum(c(x9[6], x8[6], x7[6], x6[6], x5[6], x4[6], x3[6]))) 
x9[7] = x9[6] * (1 - p[1]) + x8[6] * (1 - p[1]) + x7[6] * (1 - p[1] - p[2]) + x6[6] * (1 - p[1] - p[2] - p[3]) + x5[6] * (1 - p[1] - p[2] - p[3] - p[4]) + x4[6] * (1 - p[1] - p[2] - p[3] - p[4] - p[5]) + x3[6] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6])
x8[7] = x9[6] * p[1] + x7[6] * p[2] + x6[6] * p[3] + x5[6] * p[4] + x4[6] * p[5] + x3[6] * p[6]
x7[7] = x8[6] * p[1] + x6[6] * p[2] + x5[6] * p[3] + x4[6] * p[4] + x3[6] * p[5]  
x6[7] = x7[6] * p[1] + x5[6] * p[2] + x4[6] * p[3] + x3[6] * p[4]
x5[7] = x6[6] * p[1] + x4[6] * p[2] + x3[6] * p[3]
x4[7] = x5[6] * p[1] + x3[6] * p[2]
x3[7] = x4[6] * p[1]
x2[7] = x3[6] * p[1]
print(sum(c(x9[7], x8[7], x7[7], x6[7], x5[7], x4[7], x3[7], x2[7])))
x9[8] = x9[7] * (1 - p[1]) + x8[7] * (1 - p[1]) + x7[7] * (1 - p[1] - p[2]) + x6[7] * (1 - p[1] - p[2] - p[3]) + x5[7] * (1 - p[1] - p[2] - p[3] - p[4]) + x4[7] * (1 - p[1] - p[2] - p[3] - p[4] - p[5]) + x3[7] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6]) + x2[7] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6] - p[7])
x8[8] = x9[7] * p[1] + x7[7] * p[2] + x6[7] * p[3] + x5[7] * p[4] + x4[7] * p[5] + x3[7] * p[6] + x2[7] * p[7] 
x7[8] = x8[7] * p[1] + x6[7] * p[2] + x5[7] * p[3] + x4[7] * p[4] + x3[7] * p[5] + x2[7] * p[6] 
x6[8] = x7[7] * p[1] + x5[7] * p[2] + x4[7] * p[3] + x3[7] * p[4] + x2[7] * p[5]
x5[8] = x6[7] * p[1] + x4[7] * p[2] + x3[7] * p[3] + x2[7] * p[4] 
x4[8] = x5[7] * p[1] + x3[7] * p[2] + x2[7] * p[3] 
x3[8] = x4[7] * p[1] + x2[7] * p[2] 
x2[8] = x3[7] * p[1]
x1[8] = x2[7] * p[1]
print(sum(c(x9[8], x8[8], x7[8], x6[8], x5[8], x4[8], x3[8], x2[8], x1[8]))) 
x9[9] = x9[8] * (1 - p[1]) + x8[8] * (1 - p[1]) + x7[8] * (1 - p[1] - p[2]) + x6[8] * (1 - p[1] - p[2] - p[3]) + x5[8] * (1 - p[1] - p[2] - p[3] - p[4]) + x4[8] * (1 - p[1] - p[2] - p[3] - p[4] - p[5]) + x3[8] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6]) + x2[8] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6] - p[7]) + x1[8] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6] - p[7] - p[8])
x8[9] = x9[8] * p[1] + x7[8] * p[2] + x6[8] * p[3] + x5[8] * p[4] + x4[8] * p[5] + x3[8] * p[6] + x2[8] * p[7] + x1[8] * p[8]
x7[9] = x8[8] * p[1] + x6[8] * p[2] + x5[8] * p[3] + x4[8] * p[4] + x3[8] * p[5] + x2[8] * p[6] + x1[8] * p[7]
x6[9] = x7[8] * p[1] + x5[8] * p[2] + x4[8] * p[3] + x3[8] * p[4] + x2[8] * p[5] + x1[8] * p[6]
x5[9] = x6[8] * p[1] + x4[8] * p[2] + x3[8] * p[3] + x2[8] * p[4] + x1[8] * p[5]
x4[9] = x5[8] * p[1] + x3[8] * p[2] + x2[8] * p[3] + x1[8] * p[4]
x3[9] = x4[8] * p[1] + x2[8] * p[2] + x1[8] * p[3]
x2[9] = x3[8] * p[1] + x1[8] * p[2]
x1[9] = x2[8] * p[1]
x0[9] = x1[8] * p[1]
print(sum(c(x9[9], x8[9], x7[9], x6[9], x5[9], x4[9], x3[9], x2[9], x1[9], x0[9])))
x9[10] = x9[9] * (1 - p[1]) + x8[9] * (1 - p[1]) + x7[9] * (1 - p[1] - p[2]) + x6[9] * (1 - p[1] - p[2] - p[3]) + x5[9] * (1 - p[1] - p[2] - p[3] - p[4]) + x4[9] * (1 - p[1] - p[2] - p[3] - p[4] - p[5]) + x3[9] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6]) + x2[9] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6] - p[7]) + x1[9] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6] - p[7] - p[8]) + x0[9] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6] - p[7] - p[8] - p[9])
x8[10] = x9[9] * p[1] + x7[9] * p[2] + x6[9] * p[3] + x5[9] * p[4] + x4[9] * p[5] + x3[9] * p[6] + x2[9] * p[7] + x1[9] * p[8] + x0[9] * p[9] 
x7[10] = x8[9] * p[1] + x6[9] * p[2] + x5[9] * p[3] + x4[9] * p[4] + x3[9] * p[5] + x2[9] * p[6] + x1[9] * p[7] + x0[9] * p[8]
x6[10] = x7[9] * p[1] + x5[9] * p[2] + x4[9] * p[3] + x3[9] * p[4] + x2[9] * p[5] + x1[9] * p[6] + x0[9] * p[7]
x5[10] = x6[9] * p[1] + x4[9] * p[2] + x3[9] * p[3] + x2[9] * p[4] + x1[9] * p[5] + x0[9] * p[6]
x4[10] = x5[9] * p[1] + x3[9] * p[2] + x2[9] * p[3] + x1[9] * p[4] + x0[9] * p[5]
x3[10] = x4[9 ]* p[1] + x2[9] * p[2] + x1[9] * p[3] + x0[9] * p[4]
x2[10] = x3[9] * p[1] + x1[9] * p[2] + x0[9] * p[3]
x1[10] = x2[9] * p[1] + x0[9] * p[2]
x0[10] = x1[9] * p[1] + x0[9] * p[1]
print(sum(c(x9[10], x8[10], x7[10], x6[10], x5[10], x4[10], x3[10], x2[10], x1[10], x0[10])))
x9[11] = x9[10] * (1 - p[1]) + x8[10] * (1 - p[1]) + x7[10] * (1 - p[1] - p[2]) + x6[10] * (1 - p[1] - p[2] - p[3]) + x5[10] * (1 - p[1] - p[2] - p[3] - p[4]) + x4[10] * (1 - p[1] - p[2] - p[3] - p[4] - p[5]) + x3[10] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6]) + x2[10] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6] - p[7]) + x1[10] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6] - p[7] - p[8]) + x0[10] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6] - p[7] - p[8] - p[9])
x8[11] = x9[10] * p[1] + x7[10] * p[2] + x6[10] * p[3] + x5[10] * p[4] + x4[10] * p[5] + x3[10] * p[6] + x2[10] * p[7] + x1[10] * p[8] + x0[10] * p[9] 
x7[11] = x8[10] * p[1] + x6[10] * p[2] + x5[10] * p[3] + x4[10] * p[4] + x3[10] * p[5] + x2[10] * p[6] + x1[10] * p[7] + x0[10] * p[8]
x6[11] = x7[10] * p[1] + x5[10] * p[2] + x4[10] * p[3] + x3[10] * p[4] + x2[10] * p[5] + x1[10] * p[6] + x0[10] * p[7]
x5[11] = x6[10] * p[1] + x4[10] * p[2] + x3[10] * p[3] + x2[10] * p[4] + x1[10] * p[5] + x0[10] * p[6]
x4[11] = x5[10] * p[1] + x3[10] * p[2] + x2[10] * p[3] + x1[10] * p[4] + x0[10] * p[5]
x3[11] = x4[10] * p[1] + x2[10] * p[2] + x1[10] * p[3] + x0[10] * p[4]
x2[11] = x3[10] * p[1] + x1[10] * p[2] + x0[10] * p[3]
x1[11] = x2[10] * p[1] + x0[10] * p[2]
x0[11] = x1[10] * p[1] + x0[10] * p[1]
print(sum(c(x9[11], x8[11], x7[11], x6[11], x5[11], x4[11], x3[11], x2[11], x1[11], x0[11])))

for (i in 12:n) {
  x9[i] = 
    x9[i-1] * (1 - p[1]) + 
    x8[i-1] * (1 - p[1]) + 
    x7[i-1] * (1 - p[1] - p[2]) + 
    x6[i-1] * (1 - p[1] - p[2] - p[3]) + 
    x5[i-1] * (1 - p[1] - p[2] - p[3] - p[4]) + 
    x4[i-1] * (1 - p[1] - p[2] - p[3] - p[4] - p[5]) + 
    x3[i-1] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6]) + 
    x2[i-1] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6] - p[7]) + 
    x1[i-1] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6] - p[7] - p[8]) + 
    x0[i-1] * (1 - p[1] - p[2] - p[3] - p[4] - p[5] - p[6] - p[7] - p[8] - p[9])
  x8[i] = x9[i-1] * p[1] + x7[i-1] * p[2] + x6[i-1] * p[3] + x5[i-1] * p[4] + x4[i-1] * p[5] + x3[i-1] * p[6] + x2[i-1] * p[7] + x1[i-1] * p[8] + x0[i-1] * p[9] 
  x7[i] = x8[i-1] * p[1] + x6[i-1] * p[2] + x5[i-1] * p[3] + x4[i-1] * p[4] + x3[i-1] * p[5] + x2[i-1] * p[6] + x1[i-1] * p[7] + x0[i-1] * p[8]
  x6[i] = x7[i-1] * p[1] + x5[i-1] * p[2] + x4[i-1] * p[3] + x3[i-1] * p[4] + x2[i-1] * p[5] + x1[i-1] * p[6] + x0[i-1] * p[7]
  x5[i] = x6[i-1] * p[1] + x4[i-1] * p[2] + x3[i-1] * p[3] + x2[i-1] * p[4] + x1[i-1] * p[5] + x0[i-1] * p[6]
  x4[i] = x5[i-1] * p[1] + x3[i-1] * p[2] + x2[i-1] * p[3] + x1[i-1] * p[4] + x0[i-1] * p[5]
  x3[i] = x4[i-1] * p[1] + x2[i-1] * p[2] + x1[i-1] * p[3] + x0[i-1] * p[4]
  x2[i] = x3[i-1] * p[1] + x1[i-1] * p[2] + x0[i-1] * p[3]
  x1[i] = x2[i-1] * p[1] + x0[i-1] * p[2]
  x0[i] = x1[i-1] * p[1] + x0[i-1] * p[1]
  }
for (i in 1:n) {
  print(round(c(i, x9[i], x8[i], x7[i], x6[i], x5[i], x4[i], x3[i], x2[i], x1[i], x0[i]), digits = c(0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)))
}

Thiết kế website giá rẻ

Danh mục

simulation of customer distribution in r