## File Name: SRM_OPTIMIZE.R
## File Version: 0.293
SRM_OPTIMIZE <- function(x, parm.table, parm_list, symm_matrices, npar,
data_list, NOP, parm_table_free, use_rcpp=TRUE, shortcut=TRUE,
compute_gradient=TRUE, compute_hessian=FALSE, method="ml", inv_type="pinv",
use_woodbury=FALSE)
{
if (!compute_gradient){
compute_hessian <- FALSE
}
#* include parameters in parameter table and parameter list
parm.table <- SRM_INCLUDE_PARAMETERS_PARM_TABLE( parm=x, parm.table=parm.table)
parm_list <- SRM_INCLUDE_PARAMETERS_PARM_LIST( parm.table=parm.table,
parm_list=parm_list, symm_matrices=symm_matrices, include_fixed=FALSE )
ll_new <- 0
grad <- rep(0, npar)
hessian <- matrix(0, nrow=npar, ncol=npar)
G <- length(data_list)
for (gg in seq_len(G) ){
data_list_gg <- data_list[[gg]]
calculate_gg <- (! shortcut ) | data_list_gg$calculate_gg
# compute covariance matrix at individual level
SIGMA_U <- SRM_COMPUTE_SEM_SIGMA( parm_list=parm_list, level="U", group=gg )
# compute covariance matrix at dyad level
SIGMA_D <- SRM_COMPUTE_SEM_SIGMA( parm_list=parm_list, level="D", group=gg )
nrr <- data_list_gg$nrr
NI <- data_list_gg$NI
ND <- data_list_gg$ND
null_list <- as.list( seq_len(NOP) )
SIGMA_der <- null_list
if (compute_gradient){
for (nn in 1:NOP){
parm_type <- parm_table_free[nn, "mat"]
parm_pos <- c( parm_table_free[nn,"row"], parm_table_free[nn,"col"])
parm_level <- parm_table_free[nn, "level"]
SIGMA_der[[nn]] <- SRM_COMPUTE_SEM_SIGMA_GRADIENT( parm_list=parm_list,
parm_type=parm_type, parm_pos=parm_pos, parm_level=parm_level, group=gg )
}
}
#* loop over round robin groups rr
hess_list <- as.list( seq_len(npar) )
mu_y_der_list <- as.list( seq_len(npar) )
mu_y_der_bool_list <- as.list( seq_len(npar) )
sigma_y_der_list <- as.list( seq_len(npar) )
Zis_rcpp <- Wds_rcpp <- NULL
res_ll <- NULL
SIGMA_Y_saved <- as.list( seq_len(nrr) )
SIGMA_Y_der_saved <- null_list
res_log_like_grad_saved <- null_list
der_bool_saved <- null_list
for (rr in seq_len(nrr) ){
calculate_gg_rr <- calculate_gg[rr]
y <- data_list_gg$y_rr[[rr]]
Xs <- data_list_gg$Xs_rr[[rr]]
Zis <- data_list_gg$Zis_rr[[rr]]
Wds <- data_list_gg$Wds_rr[[rr]]
Z_ind <- data_list_gg$Z_ind_rr_rcpp[[rr]]
if (use_rcpp){
Zis_rcpp <- data_list_gg$Zis_rr_rcpp[[rr]]
Wds_rcpp <- data_list_gg$Wds_rr_rcpp[[rr]]
}
NI_rr <- NI[rr,2]
ND_rr <- ND[rr,2]
NY <- length(y)
hessian_rr <- matrix(0, nrow=npar, ncol=npar)
M0 <- matrix( 0, nrow=NY, ncol=NY)
mu0 <- rep(0, NY)
for (pp in 1:npar){
hess_list[[pp]] <- M0
mu_y_der_list[[pp]] <- mu0
mu_y_der_bool_list[[pp]] <- FALSE
sigma_y_der_list[[pp]] <- M0
}
# compute SIGMA_Y
if (calculate_gg_rr){
SIGMA_Y <- SRM_COMPUTE_SIGMA_Y( Zis=Zis, Wds=Wds, NI=NI_rr, ND=ND_rr,
SIGMA_U=SIGMA_U, SIGMA_D=SIGMA_D, NY=NY, use_rcpp=use_rcpp,
Zis_rcpp=Zis_rcpp, Wds_rcpp=Wds_rcpp, Z_ind=Z_ind )
SIGMA_Y_saved[[rr]] <- SIGMA_Y
} else {
SIGMA_Y <- SIGMA_Y_saved[[rr]]
}
# compute muy and ey
res <- SRM_COMPUTE_MU_Y(parm_list=parm_list, y=y, Xs=Xs, group=gg)
muy <- res$muy
ey <- res$ey
# compute log likelihood
args_eval <- list( y=y, muy=muy, ey=ey, SIGMA_Y=SIGMA_Y, use_rcpp=use_rcpp,
calculate_full=calculate_gg_rr, res_ll=res_ll,
inv_type=inv_type, Zis_rcpp=Zis_rcpp, Wds_rcpp=Wds_rcpp,
SIGMA_U=SIGMA_U, SIGMA_D=SIGMA_D, Z_ind=Z_ind,
NI=NI_rr, ND=ND_rr, use_woodbury=use_woodbury)
res_ll <- SRM_COMPUTE_OPTIM_FUNCTION_EVAL(args_eval=args_eval, method=method)
SIGMA_Y <- res_ll$SIGMA_Y
SIGMA_Y_inv <- res_ll$SIGMA_Y_inv
cov_resid <- res_ll$cov_resid
ll_new <- ll_new + res_ll$ll
#-- compute gradient with respect to all parameters
if (compute_gradient){
t1b <- NULL
for (nn in 1:NOP){
index_nn <- parm_table_free[nn,"index"]
#- derivative with respect to SIGMA_Y
if (calculate_gg_rr){
SIGMA_Y_der <- SRM_COMPUTE_SIGMA_Y_GRADIENT( SIGMA_der=SIGMA_der,
parm_table_free=parm_table_free, nn=nn, rr=rr, Zis=Zis, Wds=Wds, NY=NY, NI=NI, ND=ND,
use_rcpp=use_rcpp, Zis_rcpp=Zis_rcpp, Wds_rcpp=Wds_rcpp )
SIGMA_Y_der_saved[[nn]] <- SIGMA_Y_der
} else {
SIGMA_Y_der <- SIGMA_Y_der_saved[[nn]]
}
sigma_y_der_list[[index_nn]] <- sigma_y_der_list[[index_nn]] + SIGMA_Y_der
#- derivative with respect to mu
res <- SRM_COMPUTE_MU_Y_GRADIENT(y=y, Xs=Xs, parm_list=parm_list,
parm_table_free=parm_table_free, nn=nn, group=gg)
der_in_mu <- res$der_in_mu
MU_Y_der <- res$MU_Y_der
#-- non-zero gradient in SIGMA_Y_der
if (calculate_gg_rr & use_rcpp){
der_bool <- SRM_RCPP_SRM_COMPUTE_NONZERO_GRADIENT_INDICES(
sigma_y_der=SIGMA_Y_der, eps=1e-12)
der_bool_saved[[nn]] <- der_bool
} else {
der_bool <- der_bool_saved[[nn]]
}
#-- derivative log-likelihood and Hessian
args_grad <- list( y=y, muy=muy, ey=ey,
SIGMA_Y_inv=SIGMA_Y_inv, SIGMA_Y_der=SIGMA_Y_der,
MU_Y_der=MU_Y_der, der_in_mu=der_in_mu, use_rcpp=use_rcpp, t1b=t1b,
calculate_full=calculate_gg_rr, res_saved=res_log_like_grad_saved[[nn]],
der_bool=der_bool, cov_resid=cov_resid)
res <- SRM_COMPUTE_OPTIM_FUNCTION_GRAD(args_grad=args_grad, method=method)
res_log_like_grad_saved[[nn]] <- res
t1b <- res$t1b
grad[index_nn] <- grad[index_nn] + res$ll_grad_pos
if (compute_hessian){
hess_list[[index_nn]] <- hess_list[[index_nn]] + res$ll_matrix
}
mu_y_der_list[[index_nn]] <- mu_y_der_list[[index_nn]] + MU_Y_der
mu_y_der_bool_list[[index_nn]] <- mu_y_der_bool_list[[index_nn]] | der_in_mu
} # end nn ... free parameters
} # end compute_gradient
#-- compute Hessian
if (compute_hessian){
calc_hess_rr <- calculate_gg_rr
if (calc_hess_rr){
# compute Hessian for round robin group rr
hessian_rr <- SRM_COMPUTE_HESSIAN_RR(hess_list=hess_list, mu_y_der_list=mu_y_der_list,
mu_y_der_bool_list=mu_y_der_bool_list, SIGMA_Y_inv=SIGMA_Y_inv,
npar=npar, use_rcpp=use_rcpp)
hessian_rr_temp <- hessian_rr
} else {
hessian_rr <- hessian_rr_temp
}
# add to total Hessian
hessian <- hessian + hessian_rr
} # end compute_hessian
} # end rr
} ## end gg
#--- output
if (! compute_hessian){
hessian <- NULL
}
if (! compute_gradient){
grad <- NULL
}
res <- list(ll_new=ll_new, grad=grad, expected_infomat=hessian, parm.table=parm.table,
parm_list=parm_list, G=G )
return(res)
}
# cat(" *** hessian") ; z1 <- Sys.time(); print(z1-z0) ; z0 <- z1
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.