#' @importFrom stats optim var coef cov
update_alpha.lasso <- function(X, Y, Z, alpha.old, sigma.square, theta, maxstep_inner,
tolerance_inner) {
## initial
alpha.inner.old = alpha.old
k_inner = 1
n = nrow(X)
p = ncol(X)
while (k_inner < maxstep_inner) {
# given alpha update delta
gamma = 2 * exp(-2 * Z %*% alpha.inner.old)
sd_y <- sqrt(var(Y) * (n - 1)/n)
C = sum(1/gamma)/p * sd_y * sigma.square/n
delta.est = coef(glmnet(X, Y, alpha = 0, penalty.factor = 1/gamma, lambda = C,
standardize = F, intercept = FALSE))[-1]
## given delta update alpha
alpha.inner.new <- optim(alpha.old, likelihood.alpha.theta.lasso, likelihood.alpha.theta.gradient.lasso,
Z = Z, theta = theta, delta = delta.est, method = "BFGS")$par
if (sum(abs(alpha.inner.new - alpha.inner.old)) < tolerance_inner) {
break
}
k_inner = k_inner + 1
alpha.inner.old <- alpha.inner.new
}
return(list(alpha.est = alpha.inner.old, inner_iter = k_inner))
}
likelihood.alpha.theta.lasso <- function(Z, alpha, theta, delta) {
gamma = 2 * exp(-2 * Z %*% alpha)
return(as.numeric(t(theta) %*% gamma + delta^2 %*% (1/gamma)))
}
likelihood.alpha.theta.gradient.lasso <- function(Z, alpha, theta, delta) {
gamma = 2 * exp(-2 * Z %*% alpha)
dev_gamma = (theta - delta^2/(gamma^2))
return(crossprod(dev_gamma, as.vector(gamma) * Z) * (-2))
}
approx_likelihood.lasso <- function(to_estimate, X, Y, Z, sigma.square.est) {
n = nrow(X)
gamma = 2/exp(2 * Z %*% to_estimate) ## to_estimate:alpha estimates
K = sigma.square.est * diag(n) + X %*% diag(c(gamma)) %*% t(X)
logdetK = determinant(K)$modulus[1]
part1 = t(Y) %*% solve(K, Y)
normapprox = 1/2 * (part1 + logdetK)
return(-as.numeric(normapprox))
}
update_alpha.ridge <- function(X, Y, Z, alpha.old, sigma.square, theta, maxstep_inner,
tolerance_inner) {
## initial
alpha.inner.old = alpha.old
k_inner = 1
n = nrow(X)
p = ncol(X)
while (k_inner < maxstep_inner) {
# given alpha update delta
gamma = exp(-Z %*% alpha.inner.old)
sd_y <- sqrt(var(Y) * (n - 1)/n)
C = sum(1/gamma)/p * sd_y * sigma.square/n
delta.est = coef(glmnet(X, Y, alpha = 0, penalty.factor = 1/gamma, lambda = C,
standardize = F, intercept = FALSE))[-1]
## given delta update alpha
alpha.inner.new <- optim(alpha.old, likelihood.alpha.theta.ridge, likelihood.alpha.theta.gradient.ridge,
Z = Z, theta = theta, delta = delta.est, method = "BFGS")$par
if (sum(abs(alpha.inner.new - alpha.inner.old)) < tolerance_inner) {
break
}
k_inner = k_inner + 1
alpha.inner.old <- alpha.inner.new
}
return(list(alpha.est = alpha.inner.old, inner_iter = k_inner))
}
likelihood.alpha.theta.ridge <- function(Z, alpha, theta, delta) {
gamma = exp(-Z %*% alpha)
return(as.numeric(t(theta) %*% gamma + delta^2 %*% (1/gamma)))
}
likelihood.alpha.theta.gradient.ridge <- function(Z, alpha, theta, delta) {
gamma = exp(-Z %*% alpha)
dev_gamma = (theta - delta^2/(gamma^2))
return(-crossprod(dev_gamma, as.vector(gamma) * Z))
}
approx_likelihood.ridge <- function(to_estimate, X, Y, Z, sigma.square.est) {
n = nrow(X)
gamma = exp(-Z %*% to_estimate) ## to_estimate:alpha estimates
K = sigma.square.est * diag(n) + X %*% diag(c(gamma)) %*% t(X)
logdetK = determinant(K)$modulus[1]
part1 = t(Y) %*% solve(K, Y)
normapprox = 1/2 * (part1 + logdetK)
return(-as.numeric(normapprox))
}
get_mse <- function(estimation, true) {
return(mean((as.vector(estimation) - as.vector(true))^2))
}
score_function <- function(to_estimate, input_X, input_Y, input_Z, sigma2_est) {
X = input_X
Y = input_Y
Z = input_Z
sigma_square = sigma2_est
n = nrow(X)
p = ncol(X)
gamma = 2/exp(2 * Z %*% to_estimate)
Rinv <- backsolve(chol(crossprod(X)/sigma_square + diag(c(1/gamma))), diag(1,
p))
diagSigma <- rowSums(Rinv^2)
mu_vec <- (Rinv %*% (crossprod(Rinv, crossprod(X, Y))))/sigma_square
dev_gamma = (gamma - diagSigma - mu_vec^2)/(gamma^2)
return(-crossprod(dev_gamma, as.vector(gamma) * Z))
}
approx_likelihood <- function(to_estimate, input_X, input_Y, input_Z, sigma2_est) {
X = input_X
Y = input_Y
Z = input_Z
sigma_s = sigma2_est
n = nrow(X)
gamma = 2/exp(2 * Z %*% to_estimate) ## to_estimate:alpha estimates
K = sigma_s * diag(n) + X %*% diag(c(gamma)) %*% t(X)
logdetK = determinant(K)$modulus[1]
part1 = t(Y) %*% solve(K, Y)
normapprox = 1/2 * (part1 + logdetK)
return(as.numeric(normapprox))
}
beta.lda <- function(beta,X,Y){
n<-length(Y)
n1<-n-sum(Y)
n2<-sum(Y)
beta1<-beta[-1]
sel<-which(abs(beta1)!=0)
mu1<-as.matrix(apply(X[Y==0,sel,drop=F],2,mean))
mu2<-as.matrix(apply(X[Y==1,sel,drop=F],2,mean))
sigma.hat<-(cov(X[Y==0,sel,drop=F])*(n1-1)+cov(X[Y==1,sel,drop=F])*(n2-1))/(n-2)
beta1[sel]<-sweep(as.matrix(beta1[sel,drop=F]),2,sign(t(mu2-mu1)%*%beta1[sel,drop=F]),"*")
beta0<--t((mu1+mu2))%*%beta1[sel,drop=F]/2-log(n1/n2)*diag(t(beta1[sel,drop=F])%*%sigma.hat%*%beta1[sel,drop=F])/(t(mu2-mu1)%*%beta1[sel,drop=F])
beta0[is.na(beta0)]<-n2-n1
beta[1,]<-beta0
return(beta)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.