2014-01-17 4 views
1

Мне сложно определить стационарное распределение марковской модели. Я начинаю понимать теорию и связи: Учитывая стохастическую матрицу, для dermine стационарного распределения мы должны найти собственный вектор для наибольшего собственного значения (то есть 1)Вычисление собственных значений/собственных векторов стохастической матрицы

Я начал с созданием стохастической матрицы

set.seed(6534) 
stoma <- matrix(abs(rnorm(25)), nrow=5, ncol=5) 
stoma <- (stoma)/rowSums(stoma) # that should make it a stochastic matrix rowSums(stoma) == 1 

Затем я использую R eigen функцию

ew <- eigen(stoma) 

Но я не понимаю, результаты

> ew 
$values 
[1] 1.000000e+00+0.000000e+00i -6.038961e-02+0.000000e+00i -3.991160e-17+0.000000e+00i 
[4] -1.900754e-17+1.345763e-17i -1.900754e-17-1.345763e-17i 

$vectors 
       [,1]   [,2]   [,3]     [,4]     [,5] 
[1,] -0.4472136+0i 0.81018968+0i 0.3647755+0i -0.0112889+0.1658253i -0.0112889-0.1658253i 
[2,] -0.4472136+0i 0.45927081+0i -0.7687393+0i 0.5314923-0.1790588i 0.5314923+0.1790588i 
[3,] -0.4472136+0i 0.16233945+0i 0.2128250+0i -0.7093859+0.0000000i -0.7093859+0.0000000i 
[4,] -0.4472136+0i -0.09217315+0i 0.4214660+0i -0.1305497-0.1261247i -0.1305497+0.1261247i 
[5,] -0.4472136+0i -0.31275073+0i -0.2303272+0i 0.3197321+0.1393583i 0.3197321-0.1393583i 

Вектор для наибольшего значения (1) имеет все те же значения компонента «-0.4472136». Даже если я изменю семя, чтобы нарисовать различное число, я снова получаю те же значения. Что мне не хватает? Почему компоненты собственного вектора все эквалы? Почему они не суммируются до 1 - так как это должно быть стационарное распределение?

Благодарим за помощь!

+0

если размер корреляционной матрицы меньше, то количество периодов, то матрица вырождена – Qbik

ответ

4

Это путаница между левыми и правыми собственными векторами. Попробуйте eigen(t(stoma)).

Я нашел, что это полезно, чтобы попробовать пример из wikipedia article on stochastic matrices:

p <- matrix(c(0,0,1/4,0,0,0,0,1/4,0,0, 
     1/2,1,0,1/2,0,0,0,1/4,0,0,1/2,0,1/4,1/2,1), 
    ncol=5) 
p 
    [,1] [,2] [,3] [,4] [,5] 
## [1,] 0.00 0.00 0.5 0.00 0.50 
## [2,] 0.00 0.00 1.0 0.00 0.00 
## [3,] 0.25 0.25 0.0 0.25 0.25 
## [4,] 0.00 0.00 0.5 0.00 0.50 
## [5,] 0.00 0.00 0.0 0.00 1.00 

Проверьте это стохастическая матрица:

rowSums(p) 
## [1] 1 1 1 1 1 

Собственные значения:

zapsmall(eigen(p)$values) 
## [1] 1.0000000 0.7071068 -0.7071068 0.0000000 0.0000000 

собственных векторах:

print(zapsmall(eigen(p)$vectors),digits=3) 

##  [,1] [,2] [,3] [,4] [,5] 
## [1,] 0.447 0.354 0.354 -0.802 -0.609 
## [2,] 0.447 0.707 0.707 0.535 -0.167 
## [3,] 0.447 0.500 -0.500 0.000 0.000 
## [4,] 0.447 0.354 0.354 0.267 0.776 
## [5,] 0.447 0.000 0.000 0.000 0.000 

(Результаты как ваши, но знак произвольного переворачивается. R масштабирует собственные векторы, так что сумма квадратов каждого столбца равна 1: sqrt(1/5) составляет ок. 0,447 ...)

Вы ищете другой (левой) собственных векторов:

print(zapsmall(eigen(t(p))$vectors),digits=3) 
##  [,1] [,2] [,3] [,4] [,5] 
## [1,] 0 -0.149 0.3011 -0.5 0.707 
## [2,] 0 -0.149 0.3011 0.5 0.000 
## [3,] 0 -0.422 -0.8517 0.0 0.000 
## [4,] 0 -0.149 0.3011 -0.5 -0.707 
## [5,] 1 0.869 -0.0517 0.5 0.000 
+0

Спасибо! Я не знал разницы между левыми и правыми собственными векторами. – Drey

3

стационарное распределение (а.к.а. устойчивое состояние) наиболее легко рассчитывается с markovchain пакета.

library(markovchain) 
mc <- new("markovchain", transitionMatrix = stoma) 
steadyStates(mc) 

Это дает тот же ответ, как

ev <- eigen(t(stoma)) 
ev$vectors[, 1]/sum(ev$vectors[, 1]) 
+0

Спасибо! Насколько я понимаю из поста Бена, мне нужно было найти собственные векторы для транспонированной матрицы. Просто добавим небольшую поправку: 'Re (ev $ vectors [, 1]/sum (ev $ vectors [, 1]))' нам нужно будет нормализовать сумму. – Drey

+0

@Drey: хорошо пятнистый; теперь исправлено. –

Смежные вопросы