#' Hessian Matrix for Truncated Normal Distribution with Censoring with Linear Equation Mean
#'
#' @inheritParams tcensReg_llike
#'
#' @importFrom stats dnorm pnorm
#'
#' @return Matrix of Hessian values for the nth iterate
#' @keywords internal
tcensReg_hess <- function(theta, y, X, a = -Inf, v = NULL){
H <- matrix(nrow = length(theta), ncol = length(theta)) #creating an empty vector to store the gradient
#assume that there are a total of p parameters, then the first p-1 are assumed to be beta
#and the last parameter is log_sigma
p <- length(theta)
log_sigma <- theta[p]
n <- length(y) #total number of observations
if(is.null(v) == FALSE){
uncens <- which(y > v) #identifying which values are uncensored
cens <- which(y == v) #idenfitying which values are censored
n_0 <- length(cens) #number of censored observations
n_1 <- length(uncens) #number of uncensored observations
a_stand <- (a - X %*% theta[seq_len(p - 1)]) / exp(log_sigma) #standardized value with respect to truncated value
v_stand <- (v - X %*% theta[seq_len(p - 1)]) / exp(log_sigma) #standardized value with respect to censored value
y_stand <- (y - X %*% theta[seq_len(p - 1)]) / exp(log_sigma)
if(a == -Inf){
#calculating the hessian for each beta parameter
for(j in seq_len(p - 1)){
for(k in j:(p - 1)){
H1 <- 0
H2 <- sum((X[cens, j] * X[cens, k] * exp(-2 * log_sigma)) * (((pnorm(v_stand[cens])) * (dnorm(v_stand[cens]) * v_stand[cens]) + (dnorm(v_stand[cens])) ^ 2) / (pnorm(v_stand[cens])) ^ 2))
H3 <- (X[uncens, j] %*% X[uncens, k]) * exp(-2 * log_sigma)
H[j, k] <- H1 - H2 - H3
}
}
#calculating the hessian for beta with respect to log_sigma
for(j in seq_len(p - 1)){
H1 <- 0
H2 <- sum((-X[cens, j] * exp(-log_sigma)) * ((pnorm(v_stand[cens]))*((dnorm(v_stand[cens]) * (1 - v_stand[cens] ^ 2))) - ((dnorm(v_stand[cens])) * ((dnorm(v_stand[cens]) * v_stand[cens])))) / ((pnorm(v_stand[cens])) ^ 2))
H3 <- -2 * exp(-log_sigma) * (X[uncens, j] %*% y_stand[uncens])
H[j, p] <- H1 - H2 + H3
}
#this completes the upper triangle of the hessian matrix
#calculating the hessian with respect to log_sigma
#calculating the second derivative with respect to log_sigma squared
H1 <- 0
H2 <- sum(((pnorm(v_stand[cens])) * (((dnorm(v_stand[cens]) * (v_stand[cens] ^ 3 - v_stand[cens])))) + ((dnorm(v_stand[cens]) * v_stand[cens])) ^ 2) / ((pnorm(v_stand[cens])) ^ 2))
H3 <- 2 * sum((y_stand[uncens] ^ 2))
H[p, p] <- H1 - H2 - H3
#completing the matrix using the symmetric property
H[lower.tri(H)] <- H[upper.tri(H)]
}else{
#calculating the hessian with respect to the beta parameters
for(j in seq_len(p - 1)){
for(k in j:(p - 1)){
H1 <- sum((-X[, j] * X[, k] * exp(-2 * log_sigma)) * ((pnorm(-a_stand) * dnorm(a_stand) * a_stand - dnorm(a_stand) ^ 2) / pnorm(-a_stand) ^ 2))
H2 <- sum((X[cens, j] * X[cens, k] * exp(-2 * log_sigma)) * (((pnorm(v_stand[cens]) - pnorm(a_stand[cens])) * (dnorm(v_stand[cens]) * v_stand[cens] - dnorm(a_stand[cens]) * a_stand[cens]) + (dnorm(v_stand[cens]) - dnorm(a_stand[cens])) ^ 2) / (pnorm(v_stand[cens]) - pnorm(a_stand[cens])) ^ 2))
H3 <- (X[uncens, j] %*% X[uncens, k]) * exp(-2 * log_sigma)
H[j, k] <- H1 - H2 - H3
}
}
#calculating the hessian for beta with respect to log_sigma
for(j in seq_len(p - 1)){
H1 <- sum((X[,j]*exp(-log_sigma))*((pnorm(-a_stand)*(dnorm(a_stand)*(1-a_stand^2))+a_stand*dnorm(a_stand)^2)/pnorm(-a_stand)^2))
H2 <- sum((-X[cens,j]*exp(-log_sigma))*((pnorm(v_stand[cens])-pnorm(a_stand[cens]))*((dnorm(v_stand[cens])*(1-v_stand[cens]^2))-(dnorm(a_stand[cens])*(1-a_stand[cens]^2)))-((dnorm(v_stand[cens])-dnorm(a_stand[cens]))*((dnorm(v_stand[cens])*v_stand[cens])-(dnorm(a_stand[cens])*a_stand[cens]))))/((pnorm(v_stand[cens])-pnorm(a_stand[cens]))^2))
H3 <- -2*exp(-log_sigma)*(X[uncens,j]%*%y_stand[uncens])
H[j, p] <- H1 - H2 + H3
}
#this completes the upper triangle of the hessian matrix
#calculating the hessian with respect to log_sigma
#calculating the second derivative with respect to log_sigma squared
H1 <- sum((a_stand*(pnorm(-a_stand)*(dnorm(a_stand)*(1-a_stand^2))+a_stand*dnorm(a_stand)^2))/(pnorm(-a_stand)^2))
H2 <- sum(((pnorm(v_stand[cens])-pnorm(a_stand[cens]))*(((dnorm(v_stand[cens])*(v_stand[cens]^3-v_stand[cens]))-((dnorm(a_stand[cens])*(a_stand[cens]^3-a_stand[cens])))))+((dnorm(v_stand[cens])*v_stand[cens])-(dnorm(a_stand[cens])*a_stand[cens]))^2)/((pnorm(v_stand[cens])-pnorm(a_stand[cens]))^2))
H3 <- 2*sum((y_stand[uncens]^2))
H[p,p] <- H1 - H2 - H3
#completing the matrix using the symmetric property
H[lower.tri(H)]<- H[upper.tri(H)]
}
}else if(is.null(v)==TRUE){
a_stand <- (a-X%*%theta[seq_len(p - 1)])/exp(log_sigma) #standardized value with respect to truncated value
y_stand <- (y-X%*%theta[seq_len(p - 1)])/exp(log_sigma)
if(a == -Inf){
#calculating the hessian for each beta parameter
for(j in seq_len(p - 1)){
for(k in j:(p-1)){
H1 <- 0
H3 <- (X[,j]%*%X[,k])*exp(-2*log_sigma)
H[j, k] <- H1 - H3
}
}
#calculating the hessian for beta with respect to log_sigma
for(j in seq_len(p - 1)){
H1 <- 0
H3 <- -2*exp(-log_sigma)*(X[,j]%*%y_stand)
H[j, p] <- H1 + H3
}
#this completes the upper triangle of the hessian matrix
#calculating the hessian with respect to log_sigma
#calculating the second derivative with respect to log_sigma squared
H1 <- 0
H3 <- 2*sum((y_stand^2))
H[p,p] <- H1 - H3
#completing the matrix using the symmetric property
H[lower.tri(H)]<- H[upper.tri(H)]
}else{
#calculating the hessian for each beta parameter
for(j in seq_len(p - 1)){
for(k in j:(p-1)){
H1 <- sum((-X[,j]*X[,k]*exp(-2*log_sigma))*((pnorm(-a_stand)*dnorm(a_stand)*a_stand-dnorm(a_stand)^2)/pnorm(-a_stand)^2))
H3 <- (X[,j]%*%X[,k])*exp(-2*log_sigma)
H[j, k] <- H1 - H3
}
}
#calculating the hessian for beta with respect to log_sigma
for(j in seq_len(p - 1)){
H1 <- sum((X[,j]*exp(-log_sigma))*((pnorm(-a_stand)*(dnorm(a_stand)*(1-a_stand^2))+a_stand*dnorm(a_stand)^2)/pnorm(-a_stand)^2))
H3 <- -2*exp(-log_sigma)*(X[,j]%*%y_stand)
H[j, p] <- H1 + H3
}
#this completes the upper triangle of the hessian matrix
#calculating the hessian with respect to log_sigma
#calculating the second derivative with respect to log_sigma squared
H1 <- sum((a_stand*(pnorm(-a_stand)*(dnorm(a_stand)*(1-a_stand^2))+a_stand*dnorm(a_stand)^2))/(pnorm(-a_stand)^2))
H3 <- 2*sum((y_stand^2))
H[p,p] <- H1 - H3
#completing the matrix using the symmetric property
H[lower.tri(H)]<- H[upper.tri(H)]
}
}
return(H)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.