Saya sudah mendapatkan kode berikut

P <- matrix(c(0.8,0.2,0.5,0.5),2,2,byrow=TRUE)
ev <- eigen(P)
rvec <- ev$vectors
lvec <- ginv(rvec)
pi_eig <-lvec[1,]/sum(lvec[1,])
powers <- 0:20
upowers <- lapply(powers, function(k) P %^% k)
A <- for (i in 1:20) {upowers[i]-pi_eig}

Saya ingin menghitung A yang berbeda dari kekuatan P yang berbeda dikurangi distribusi stasioner P, seperti yang ditunjukkan pada gambar.

enter image description here

Kode terakhir salah, dan saya tidak tahu bagaimana menyelesaikan masalah ini.

0
Berry 6 Maret 2019, 18:40

1 menjawab

Jawaban Terbaik

Berikut adalah solusi untuk menghitung A^k:

Misalkan Anda sudah memiliki matriks P dan P_infty yang akan saya rujuk sebagai Q. Anda kemudian dapat melakukan hal berikut:

# Data
set.seed(123)
P <- matrix(runif(9), 3L, 3L)
Q <- matrix(runif(9), 3L, 3L)
# Computes P^k (matrix multiplication, not R's P^k!)
matMultK <- function (k) {
  txt <- paste("P", paste(rep("%*% P", k), collapse = " "))
  eval(expr = parse(text = txt))
}
# Compute the results
lapply(1:3, function (k) matMultK(k) - Q)

Fungsi matMultK akan menghitung pangkat matriks P^k dari P bukan versi R (dalam R P^k dilakukan berdasarkan elemen). Perhatikan bahwa dengan beberapa modifikasi kecil, Anda juga dapat memiliki P menjadi dan argumen.

Terakhir, menggunakan lapply adalah cara yang mudah untuk menghitung barisan A^k. Di sini sekali lagi, ada banyak alternatif.

0
niko 6 Maret 2019, 16:55