Использование optim () для создания ковариационной матрицы в R - PullRequest
Новая
       45

Использование optim () для создания ковариационной матрицы в R

0 голосов
/ 18 ноября 2018

Я хотел бы создать ковариационную матрицу в R, используя функцию optim (). В частности, ковариационная матрица выглядит так: 

[ (sigma_a)^2               rho*sigma_a*sigma_b     0 ]
[ rho*sigma_a*sigma_b       (sigma_b)^2             0 ]
[      0                         0                  0 ]

, где (sigma_a)^2 и (sigma_b)^2 - дисперсии, а rho*sigma_a*sigma_b - ковариационная часть (rho - корреляция).

Я добавил свой код ниже, в частности, я использую пакет KFAS в R, поэтому большая часть кода не относится к моему вопросу - это часть updatefn, которая имеет ковариационную матрицу Q который я хотел бы, чтобы кто-нибудь проверил, если это возможно, так как я не уверен в этой части: 

# This section sets up everything
library('KFAS')
set.seed(1)
xx <- rnorm(1000)
Zt <- matrix(c(1,1,0),1,3)
Ht <- matrix(0)
Tt <- matrix(c(1,0,0,0,NA,NA,0,1,0),3,3,byrow=TRUE)
Rt <- matrix(c(1,0,0,0,1,0,0,0,0),3,3,byrow=TRUE)
Qt <- matrix(c(NA,NA,0,NA,NA,0,0,0,0),3,3,byrow=TRUE)
a1 <- matrix(c(0, 0, 0), 3, 1)
P1 <- matrix(0, 3, 3)
P1inf <- matrix(c(1,0,0,0,1,0,0,0,1),3,3,byrow=TRUE)
model_gaussian <- SSModel(xx  ~ -1 + SSMcustom(Z = Zt, T = Tt, R = Rt, Q = Qt, a1 = a1, P1 = P1, P1inf = P1inf), H = Ht) 

##############################
# RELEVANT PART OF CODE BELOW
#############################

updatefn <- function(pars, model) {

   model$T[is.na(model$T)] <- pars[1:2]  
   model$Q[is.na(model$Q)] <- pars[3:6] # Matrix Q is the Covariance matrix

   model["Q", etas = "custom"] <- crossprod(model$Q) # Creates covariance matrix from estimated parameters

   model
}


fitSSM(model_gaussian, c(0,0, 0.1,0.1,0.1,0.1), updatefn, method = "L-BFGS-B", lower = c(-0.8, -0.8, 0, -Inf, -Inf,0), upper = c(0.8, 0.8, Inf, Inf, Inf, Inf))

Последняя строка кода (т. Е. fitSSM())) в основном является входом для функции optim () и выводит параметры следующим образом: 

$optim.out
$optim.out$par
[1] -0.2063339 -0.1199854  0.3001893  0.3001893  0.3001893  0.3001893

(значения -0.2063339 -0.1199854 для матрицы T и не имеют отношения к моему вопросу)

Параметры, полученные из функции optim () для ковариационной матрицы, даны как 0.3001893 0.3001893 0.3001893 0.3001893, что не сработает.

Я также пытался изменить мою функцию следующим образом: 

updatefn <- function(pars, model) {

    model$T[is.na(model$T)] <- pars[1:2]
    Q <- matrix(c(NA,NA,0,NA,NA,0,0,0,0),3,3,byrow=TRUE)

    Q[1,1] <- pars[3]
    Q[2,2] <- pars[6]
    Q[1,2] <- pars[4]
    Q[2,1] <- pars[5]
    Q <- crossprod(Q)

    model$Q[is.na(model$Q)] <- c(Q[1,1],Q[1,2],Q[2,1],Q[2,2])

    model
}

Но этот метод все еще не работает - значения Q [1,2] и Q [2,1] должны быть равны, а также Q [1,2], деленное на (sqrt (Q [1,1]) ) * sqrt (Q [2,2])) должен дать корреляцию, которая является значением между -1 и 1 (что, к сожалению, я не получаю. Поэтому, если кто-то может помочь мне с этим, я был бы признателен, спасибо!

...