RMLP <- function(y, S, ml0, mu0, maxitr = 200){
N <- dim(y)[1]
p <- dim(y)[2]
muc <- ml0[2:p] # initial values
mu <- c(mu0, muc)
g1 <- ml0[p+1]
g2 <- ml0[p+2]
Qc0 <- c(muc, g1, g2)
RLL <- function(mur){
mu1 <- c(mu0,mur)
G <- gmat(g1, g2, p)
ll1 <- 0
for(i in 1:N){
yi <- as.vector(y[i, ])
wi <- which(is.na(yi) == FALSE)
pl <- length(wi)
Si <- vmat(S[i, ], p)
yi <- yi[wi]
Si <- pmat(Si, wi)
mui <- mu1[wi]
Gi <- pmat(G, wi)
B1 <- (yi - mui)
B2 <- ginv(Gi + Si)
A1 <- log(det(Gi + Si))
A2 <- t(B1) %*% B2 %*% B1
A3 <- pl * log(2*pi)
ll1 <- ll1 + 0.5*(A1 + A2 + A3)
}
return(ll1)
}
LL1 <- function(g){
#G <- gmat(g, g2, p)
G <- gmat(g, g/2, p)
ll1 <- 0
for(i in 1:N){
yi <- as.vector(y[i, ])
wi <- which(is.na(yi) == FALSE)
pl <- length(wi)
Si <- vmat(S[i, ], p)
yi <- yi[wi]
Si <- pmat(Si, wi)
mui <- mu[wi]
Gi <- pmat(G, wi)
B1 <- (yi - mui)
B2 <- ginv(Gi + Si)
A1 <- log(det(Gi + Si))
A2 <- t(B1) %*% B2 %*% B1
A3 <- pl * log(2*pi)
ll1 <- ll1 + A1 + A2 + A3
}
return(ll1)
}
LL2 <- function(g){
#G <- gmat(g, g2, p)
G <- gmat(g, g/2, p)
ll1 <- 0
for(i in 1:N){
yi <- as.vector(y[i, ])
wi <- which(is.na(yi) == FALSE)
pl <- length(wi)
Si <- vmat(S[i, ], p)
yi <- yi[wi]
Si <- pmat(Si, wi)
mui <- mu[wi]
Gi <- pmat(G, wi)
B1 <- (yi - mui)
B2 <- ginv(Gi + Si)
A1 <- log(det(Gi + Si))
A2 <- t(B1) %*% B2 %*% B1
A3 <- pl * log(2*pi)
ll1 <- ll1 + A1 + A2 + A3
}
return(ll1) # return "minus loglikelihood"
}
for(itr in 1:maxitr){
A1 <- A2 <- A3 <- numeric(p - 1)
A4 <- matrix(numeric((p - 1)*(p - 1)), (p - 1))
G <- gmat(g1, g2, p)
for(i in 1:N){
yi <- as.vector(y[i, ])
wi <- which(is.na(yi) == FALSE)
pl <- length(wi)
Si <- vmat(S[i, ], p)
yi <- yi[wi]
Si <- pmat(Si, wi)
Gi <- pmat(G, wi)
Wi <- ginv(Gi + Si)
Wi <- imat(Wi, wi, p)
yi <- ivec(yi, wi, p)
Wi21 <- Wi[2:p, 1]
Wi22 <- Wi[2:p, 2:p]
A1 <- A1 + Wi21 * yi[1]
A2 <- A2 + Wi22 %*% yi[2:p]
A3 <- A3 + Wi21 * mu0
A4 <- A4 + Wi22
}
muc <- as.vector(A1 + A2 - A3) %*% ginv(A4)
mu <- c(mu0, muc)
g1 <- optimize(LL1, lower = 0, upper = 5)$minimum
g2 <- 0.5*g1
Qc <- c(muc, g1, g2)
rb <- abs(Qc - Qc0)/abs(Qc0); rb[is.nan(rb)] <- 0
if(max(rb) < 10^-4) break
Qc0 <- Qc
}
R7 <- c(mu, g1, g2)
R8 <- as.numeric(-RLL(muc))
return(list("RML" = R7, "Loglikelihood" = R8))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.