Dear addresses, I need perform a batch of 10 000 simulations for each of
4 options considered. (The idea is to obtain the parameter estimates in
a heteroskedastic linear regression model - with additive or mixed
heteroskedasticity - via the Kenward-Roger small-sample adjusted
covariance matrix of disturbances). For this purpose I wrote an R
program which would capture all possible options (true
heteroskedasticity = additive or mixed / assumed heteroskedasticity =
additive or mixed): simulate data, compute ML estimates of covariance
components parameters (on the basis of maxLik package), and compute the
resulting beta estimates. Simulations are performed by repeat command
and, after every simulation step, the result is exported to an external
csv file. Since sometimes an error was encountered during the simulation
and the entire program broke down, I made use of the function TRY at
critical evaluation steps to prevent the cessation of the program
evaluation: R evaluates the result and if there should be a problem, it
moves onto another step. Therefore, I distinguish between bad
simulations and good simulations (and have two csv files for the
results). There, however, another problem arises: R seems to be working
(computing), it is not frozen or anything like that, and yet the saving
to the files stops. I should say that it looks as if R was functioning
fully but the simulations stopped. It is possible that R went into an
undesired loop for no reasons that I know of. It requires of me to
service all computers and, after the simulations have frozen, to restart
the repeat command, which is not sound for I cannot be (or perhaps am
not willing to be) with computers non-stop. The computers I use are of
2GB RAM and R is the version 2.10.1 (2009-12-14). I have no idea where
the problem can lie and how to solve it. Can you please give me some
suggestions how to handle this obstacle to speedy simulating? I shall be
much thankful for any suggestion. Later I paste the example of the set
of commands that I use, all though it might be not sensible, because it
all is too complicated. I am convinced that of foremost importance is
the last repeat command. Thank you very much. Faithfully yours, Martin
Boďa (Matej Bel University / Faculty of Natural Sciences / Banská
Bystrica / Slovakia).
# zavádza sériu podporných funkcií
# zavádza trace, sigma, fi, xs, xs, fixs, sxfixs
tr <- function(M) sum(diag(M))
sigma_ADD <- function(alpha) diag(as.vector(alpha %*% t(z)),
nrow=nrow(z), ncol=nrow(z))
sigma_MIX <- function(alpha) diag((as.vector(alpha %*% t(z)))^2,
nrow=nrow(z), ncol=nrow(z))
sigma_EXP <- function(alpha) diag(exp(as.vector(alpha %*% t(z))),
nrow=nrow(z), ncol=nrow(z))
fi <- function(sigma) solve(crossprod(x, crossprod(solve(sigma),x)))
xs <- function(sigma) crossprod(x,solve(sigma))
sx <- function(sigma) crossprod(solve(sigma),x)
fixs <- function(sigma) crossprod(fi(sigma),xs(sigma))
sxfixs <- function(sigma)
crossprod(xs(sigma),crossprod(fi(sigma),xs(sigma)))
# zavádza (ncol(z) x 1) vektor núl, ktorý má na i-tom mieste 1
ei <- function(i) { ei = rep(0,ncol(z))
ei[i] = 1
return(ei) }
# zavádza deriváciu dSIGMA/dalfa_i a zmiešanú deriváciu
(dSIGMA^2)/(ddalfa_i*dalfa_j)
# zavádza maticový súčin (sigma)^(-1)*dSIGMA/dalfa_i
# zavádza maticový súčin (sigma)^(-1)*dSIGMA/dalfa_i*(sigma)^(-1)
dsigma.dalphai_ADD <- function(i) {
diag(as.vector(ei(i) %*% t(z)), nrow=nrow(z), ncol=nrow(z)) }
dsigma.dalphai_MIX <- function(sigma,i) {
2*crossprod(diag(as.vector(ei(i) %*% t(z)), nrow=nrow(z),
ncol=nrow(z)),sigma^(1/2)) }
dsigma.dalphai_EXP <- function(sigma,i) {
crossprod(diag(as.vector(ei(i) %*% t(z)), nrow=nrow(z),
ncol=nrow(z)),sigma) }
d2sigma.dalphaij_MIX <- function(i,j) {
eiz <- diag(as.vector(ei(i) %*% t(z)), nrow=nrow(z), ncol=nrow(z))
ejz <- diag(as.vector(ei(j) %*% t(z)), nrow=nrow(z), ncol=nrow(z))
d2sigma.dalphaij <- 2*crossprod(eiz,ejz)
return(d2sigma.dalphaij) }
d2sigma.dalphaij_EXP <- function(sigma,i,j) {
eiz <- diag(as.vector(ei(i) %*% t(z)), nrow=nrow(z), ncol=nrow(z))
ejz <- diag(as.vector(ei(j) %*% t(z)), nrow=nrow(z), ncol=nrow(z))
d2sigma.dalphaij <- crossprod(crossprod(eiz,ejz),sigma)
return(d2sigma.dalphaij) }
sigmaDERIVi_ADD <- function(sigma,i)
crossprod(solve(sigma),dsigma.dalphai_ADD(i))
sigmaDERIVi_MIX <- function(sigma,i)
crossprod(solve(sigma),dsigma.dalphai_MIX(sigma,i))
sigmaDERIVi_EXP <- function(sigma,i)
crossprod(solve(sigma),dsigma.dalphai_EXP(sigma,i))
sigmaDERIVisigma_ADD <- function(sigma,i) {
crossprod(solve(sigma),crossprod(dsigma.dalphai_ADD(i),solve(sigma))) }
sigmaDERIVisigma_MIX <- function(sigma,i) {
crossprod(solve(sigma),crossprod(dsigma.dalphai_MIX(sigma,i),solve(sigma)))
}
sigmaDERIVisigma_EXP <- function(sigma,i) {
crossprod(solve(sigma),crossprod(dsigma.dalphai_EXP(sigma,i),solve(sigma)))
}
# zavádza pi, qij, rij (iba pre zmiešané a multiplikatívny model
heteroskedasticity)
pi_ADD <- function(sigma,i)
-1*crossprod(x,crossprod(sigmaDERIVisigma_ADD(sigma,i),x))
pi_MIX <- function(sigma,i)
-1*crossprod(x,crossprod(sigmaDERIVisigma_MIX(sigma,i),x))
pi_EXP <- function(sigma,i)
-1*crossprod(x,crossprod(sigmaDERIVisigma_EXP(sigma,i),x))
qij_ADD <- function(sigma,i,j) {
SX <- sx(sigma)
dSIGMA.dalphaj <- dsigma.dalphai_ADD(j)
SIGMAderivi <- sigmaDERIVi_ADD(sigma,i)
qij <- crossprod(SX,crossprod(SIGMAderivi,crossprod(dSIGMA.dalphaj,SX)))
return(qij) }
qij_MIX <- function(sigma,i,j) {
SX <- sx(sigma)
dSIGMA.dalphaj <- dsigma.dalphai_MIX(sigma,j)
SIGMAderivi <- sigmaDERIVi_MIX(sigma,i)
qij <- crossprod(SX,crossprod(SIGMAderivi,crossprod(dSIGMA.dalphaj,SX)))
return(qij) }
qij_EXP <- function(sigma,i,j) {
SX <- sx(sigma)
dSIGMA.dalphaj <- dsigma.dalphai_EXP(sigma,j)
SIGMAderivi <- sigmaDERIVi_EXP(sigma,i)
qij <- crossprod(SX,crossprod(SIGMAderivi,crossprod(dSIGMA.dalphaj,SX)))
return(qij) }
rij_MIX <- function(sigma,i,j)
crossprod(sx(sigma),crossprod(d2sigma.dalphaij_MIX(i,j),sx(sigma)))
rij_EXP <- function(sigma,i,j) {
crossprod(sx(sigma),crossprod(d2sigma.dalphaij_EXP(sigma,i,j),sx(sigma)))
}
# zavádza funkciu pre stanovenie REML likelihood
# je funkciou alpha conditional on x, y, z
#########################################################################################
reml.loglik_ADD <- function(alpha) { SIGMA <- sigma_ADD(alpha)
a <- -0.5*log(det(SIGMA))
b <- -0.5*log(det(crossprod(x,crossprod(solve(SIGMA),x))))
c <-
-0.5*as.numeric(crossprod(y,crossprod(solve(SIGMA)-sxfixs(SIGMA),y)))
logl <- a + b + c
return(logl) }
reml.loglik_MIX <- function(alpha) { SIGMA <- sigma_MIX(alpha)
a <- -0.5*log(det(SIGMA))
b <- -0.5*log(det(crossprod(x,crossprod(solve(SIGMA),x))))
c <-
-0.5*as.numeric(crossprod(y,crossprod(solve(SIGMA)-sxfixs(SIGMA),y)))
logl <- a + b + c
return(logl) }
reml.loglik_EXP <- function(alpha) { SIGMA <- sigma_EXP(alpha)
a <- -0.5*log(det(SIGMA))
b <- -0.5*log(det(crossprod(x,crossprod(solve(SIGMA),x))))
c <-
-0.5*as.numeric(crossprod(y,crossprod(solve(SIGMA)-sxfixs(SIGMA),y)))
logl <- a + b + c
return(logl) }
#########################################################################################
# zavádza funkciu pre stanovenie gradientného vektora pre REML
likelihood
# je funkciou alpha conditional on x, y, z
# obsahuje podfunkciu dl.dalphai, ktorá spočítava deriváciu
likelihood function podľa alpha i
#########################################################################################
######## gradientný vektor aditívneho modelu
############################################
gradvec_ADD = function(alpha) { dl.dalphai <- function(i) {
SIGMA <- sigma_ADD(alpha)
FI <- fi(SIGMA)
FIxs <- fixs(SIGMA)
sxFIxs <- sxfixs(SIGMA)
Pi <- pi_ADD(SIGMA,i)
dSIGMA.dalphai <- dsigma.dalphai_ADD(i)
SIGMAderiviSIGMA <- sigmaDERIVisigma_ADD(SIGMA,i)
dl.dalphai.1 <- -1/2*tr(crossprod(solve(SIGMA),dSIGMA.dalphai ))
dl.dalphai.2 <- -1/2*tr(crossprod(FI,Pi))
dl.dalphai.3 <-
+1/2*as.numeric(crossprod(y,crossprod(SIGMAderiviSIGMA,y)))
INSIDE1 <- crossprod(FIxs,crossprod(Pi,FIxs))
dl.dalphai.4 <- -1/2*as.numeric(crossprod(y,crossprod(INSIDE1,y)))
INSIDE2 <- crossprod(solve(SIGMA),crossprod(dSIGMA.dalphai ,sxFIxs))
dl.dalphai.5 <-
+1/2*(-2)*as.numeric(crossprod(y,crossprod(t(INSIDE2),y)))
dl.dalphai = dl.dalphai.1 + dl.dalphai.2 + dl.dalphai.3 + dl.dalphai.4 +
dl.dalphai.5
return(dl.dalphai)}
grad = rep(dl.dalphai(1),c(1))
for(i in 2:length(alpha)) {
grad = append(grad,rep(dl.dalphai(i),c(1)), after = c(i)) }
return(grad) }
######## gradientný vektor zmiešaného modelu
############################################
gradvec_MIX = function(alpha) { dl.dalphai <- function(i) {
SIGMA <- sigma_MIX(alpha)
FI <- fi(SIGMA)
FIxs <- fixs(SIGMA)
sxFIxs <- sxfixs(SIGMA)
Pi <- pi_MIX(SIGMA,i)
dSIGMA.dalphai <- dsigma.dalphai_MIX(SIGMA,i)
SIGMAderiviSIGMA <- sigmaDERIVisigma_MIX(SIGMA,i)
dl.dalphai.1 <- -1/2*tr(crossprod(solve(SIGMA),dSIGMA.dalphai ))
dl.dalphai.2 <- -1/2*tr(crossprod(FI,Pi))
dl.dalphai.3 <-
+1/2*as.numeric(crossprod(y,crossprod(SIGMAderiviSIGMA,y)))
INSIDE1 <- crossprod(FIxs,crossprod(Pi,FIxs))
dl.dalphai.4 <- -1/2*as.numeric(crossprod(y,crossprod(INSIDE1,y)))
INSIDE2 <- crossprod(solve(SIGMA),crossprod(dSIGMA.dalphai ,sxFIxs))
dl.dalphai.5 <-
+1/2*(-2)*as.numeric(crossprod(y,crossprod(t(INSIDE2),y)))
dl.dalphai = dl.dalphai.1 + dl.dalphai.2 + dl.dalphai.3 + dl.dalphai.4 +
dl.dalphai.5
return(dl.dalphai)}
grad = rep(dl.dalphai(1),c(1))
for(i in 2:length(alpha)) {
grad = append(grad,rep(dl.dalphai(i),c(1)), after = c(i)) }
return(grad) }
######## gradientný vektor multiplikatívneho modelu
#####################################
gradvec_EXP = function(alpha) { dl.dalphai <- function(i) {
SIGMA <- sigma_EXP(alpha)
FI <- fi(SIGMA)
FIxs <- fixs(SIGMA)
sxFIxs <- sxfixs(SIGMA)
Pi <- pi_EXP(SIGMA,i)
dSIGMA.dalphai <- dsigma.dalphai_EXP(SIGMA,i)
SIGMAderiviSIGMA <- sigmaDERIVisigma_EXP(SIGMA,i)
dl.dalphai.1 <- -1/2*tr(crossprod(solve(SIGMA),dSIGMA.dalphai ))
dl.dalphai.2 <- -1/2*tr(crossprod(FI,Pi))
dl.dalphai.3 <-
+1/2*as.numeric(crossprod(y,crossprod(SIGMAderiviSIGMA,y)))
INSIDE1 <- crossprod(FIxs,crossprod(Pi,FIxs))
dl.dalphai.4 <- -1/2*as.numeric(crossprod(y,crossprod(INSIDE1,y)))
INSIDE2 <- crossprod(solve(SIGMA),crossprod(dSIGMA.dalphai ,sxFIxs))
dl.dalphai.5 <-
+1/2*(-2)*as.numeric(crossprod(y,crossprod(t(INSIDE2),y)))
dl.dalphai = dl.dalphai.1 + dl.dalphai.2 + dl.dalphai.3 + dl.dalphai.4 +
dl.dalphai.5
return(dl.dalphai)}
grad = rep(dl.dalphai(1),c(1))
for(i in 2:length(alpha)) {
grad = append(grad,rep(dl.dalphai(i),c(1)), after = c(i)) }
return(grad) }
#########################################################################################
# zavádza funkciu pre stanovenie Hessovej matice pre REML likelihood
# je funkciou alpha conditional on x, y, z
# obsahuje podfunkciu d2l.dalphaij, čo spočítava druhú deriváciu
loglikelihood podľa alpha i a alpha j
#########################################################################################
######## Hessova matica aditívneho modelu
##############################################
hessmat_ADD = function(alpha) {
hess = matrix(nrow = ncol(z), ncol = ncol(z))
d2l.dalphaij <- function(i,j) {
SIGMA <- sigma_ADD(alpha)
FI <- fi(SIGMA)
XS <- xs(SIGMA)
SX <- sx(SIGMA)
FIxs <- fixs(SIGMA)
sxFI <- t(FIxs)
sxFIxs <- sxfixs(SIGMA)
Pi <- pi_ADD(SIGMA,i)
Pj <- pi_ADD(SIGMA,j)
Qij <- qij_ADD(SIGMA,i,j)
dSIGMA.dalphai <- dsigma.dalphai_ADD(i)
dSIGMA.dalphaj <- dsigma.dalphai_ADD(j)
SIGMAderiviSIGMA <- sigmaDERIVisigma_ADD(SIGMA,i)
SIGMAderivjSIGMA <- sigmaDERIVisigma_ADD(SIGMA,j)
SIGMAderivi <- sigmaDERIVi_ADD(SIGMA,i)
SIGMAderivj <- sigmaDERIVi_ADD(SIGMA,j)
d2l.dalphaij.11 = 1/2*tr(crossprod(t(SIGMAderivi),SIGMAderivj))
d2l.dalphaij.12 = 0
d2l.dalpahij.21 = 1/2*tr(crossprod(FI,crossprod(Pi,crossprod(FI,Pj))))
d2l.dalphaij.22 = -1*tr(crossprod(FI,Qij))
d2l.dalphaij.23 = 0
d2l.dalphaij.31 =
-1*crossprod(y,crossprod(t(SIGMAderivj),crossprod(SIGMAderiviSIGMA,y)))
d2l.dalphaij.32 = 0
d2l.dalphaij.41 =
crossprod(y,crossprod(t(SIGMAderivj),crossprod(FIxs,crossprod(Pi,crossprod(sxFI,y)))))
d2l.dalphaij.42 =
crossprod(y,crossprod(FIxs,crossprod(Pi,crossprod(FI,crossprod(Pj,crossprod(sxFI,y))))))
d2l.dalphaij.43 =
-1*crossprod(y,crossprod(FIxs,crossprod(Qij,crossprod(sxFI,y))))
d2l.dalphaij.44 = 0
d2l.dalphaij.51 =
crossprod(y,crossprod(t(SIGMAderivj),crossprod(t(SIGMAderivi),crossprod(sxFIxs,y))))
d2l.dalphaij.52 = 0
d2l.dalphaij.53 =
crossprod(y,crossprod(t(SIGMAderivi),crossprod(t(SIGMAderivj),crossprod(sxFIxs,y))))
d2l.dalphaij.54 =
crossprod(y,crossprod(t(SIGMAderivi),crossprod(FIxs,crossprod(Pj,crossprod(sxFI,y)))))
d2l.dalphaij.55 =
crossprod(y,crossprod(t(SIGMAderivi),crossprod(sxFIxs,crossprod(SIGMAderivj,y))))
d2l.dalphaij = d2l.dalphaij.11 + d2l.dalphaij.12 + d2l.dalpahij.21 +
d2l.dalphaij.22 + d2l.dalphaij.23 +
d2l.dalphaij.31 + d2l.dalphaij.32 + d2l.dalphaij.41 + d2l.dalphaij.42 +
d2l.dalphaij.43 +
d2l.dalphaij.44 + d2l.dalphaij.51 + d2l.dalphaij.52 + d2l.dalphaij.53 +
d2l.dalphaij.54 +
d2l.dalphaij.55
return(d2l.dalphaij)}
for(i in 1:ncol(z)) { for(j in 1:ncol(z)) hess[i,j] = d2l.dalphaij(i,j)
}
return(hess)}
######## Hessova matica zmiešaného modelu
##############################################
hessmat_MIX = function(alpha) {
hess = matrix(nrow = ncol(z), ncol = ncol(z))
d2l.dalphaij <- function(i,j) {
SIGMA <- sigma_MIX(alpha)
FI <- fi(SIGMA)
XS <- xs(SIGMA)
SX <- sx(SIGMA)
FIxs <- fixs(SIGMA)
sxFI <- t(FIxs)
sxFIxs <- sxfixs(SIGMA)
Pi <- pi_MIX(SIGMA,i)
Pj <- pi_MIX(SIGMA,j)
Qij <- qij_MIX(SIGMA,i,j)
Rij<- rij_MIX(SIGMA,i,j)
dSIGMA.dalphai <- dsigma.dalphai_MIX(SIGMA,i)
dSIGMA.dalphaj <- dsigma.dalphai_MIX(SIGMA,j)
d2SIGMA.dalphaij <- d2sigma.dalphaij_MIX(i,j)
SIGMAderiviSIGMA <- sigmaDERIVisigma_MIX(SIGMA,i)
SIGMAderivjSIGMA <- sigmaDERIVisigma_MIX(SIGMA,j)
SIGMAderivi <- sigmaDERIVi_MIX(SIGMA,i)
SIGMAderivj <- sigmaDERIVi_MIX(SIGMA,j)
d2l.dalphaij.11 = 1/2*tr(crossprod(t(SIGMAderivi),SIGMAderivj))
d2l.dalphaij.12 = -1/2*tr(crossprod(solve(SIGMA),d2SIGMA.dalphaij))
d2l.dalpahij.21 = 1/2*tr(crossprod(FI,crossprod(Pi,crossprod(FI,Pj))))
d2l.dalphaij.22 = -1*tr(crossprod(FI,Qij))
d2l.dalphaij.23 = 1/2*tr(crossprod(FI,Rij))
d2l.dalphaij.31 =
-1*crossprod(y,crossprod(t(SIGMAderivj),crossprod(SIGMAderiviSIGMA,y)))
d2l.dalphaij.32 =
1/2*crossprod(y,crossprod(solve(SIGMA),crossprod(d2SIGMA.dalphaij,crossprod(solve(SIGMA),y))))
d2l.dalphaij.41 =
crossprod(y,crossprod(t(SIGMAderivj),crossprod(FIxs,crossprod(Pi,crossprod(sxFI,y)))))
d2l.dalphaij.42 =
crossprod(y,crossprod(FIxs,crossprod(Pi,crossprod(FI,crossprod(Pj,crossprod(sxFI,y))))))
d2l.dalphaij.43 =
-1*crossprod(y,crossprod(FIxs,crossprod(Qij,crossprod(sxFI,y))))
d2l.dalphaij.44 =
1/2*crossprod(y,crossprod(FIxs,crossprod(Rij,crossprod(sxFI,y))))
d2l.dalphaij.51 =
crossprod(y,crossprod(t(SIGMAderivj),crossprod(t(SIGMAderivi),crossprod(sxFIxs,y))))
d2l.dalphaij.52 =
-1*crossprod(y,crossprod(solve(SIGMA),crossprod(d2SIGMA.dalphaij,crossprod(sxFIxs,y))))
d2l.dalphaij.53 =
crossprod(y,crossprod(t(SIGMAderivi),crossprod(t(SIGMAderivj),crossprod(sxFIxs,y))))
d2l.dalphaij.54 =
crossprod(y,crossprod(t(SIGMAderivi),crossprod(FIxs,crossprod(Pj,crossprod(sxFI,y)))))
d2l.dalphaij.55 =
crossprod(y,crossprod(t(SIGMAderivi),crossprod(sxFIxs,crossprod(SIGMAderivj,y))))
d2l.dalphaij = d2l.dalphaij.11 + d2l.dalphaij.12 + d2l.dalpahij.21 +
d2l.dalphaij.22 + d2l.dalphaij.23 +
d2l.dalphaij.31 + d2l.dalphaij.32 + d2l.dalphaij.41 + d2l.dalphaij.42 +
d2l.dalphaij.43 +
d2l.dalphaij.44 + d2l.dalphaij.51 + d2l.dalphaij.52 + d2l.dalphaij.53 +
d2l.dalphaij.54 +
d2l.dalphaij.55
return(d2l.dalphaij)}
for(i in 1:ncol(z)) { for(j in 1:ncol(z)) hess[i,j] = d2l.dalphaij(i,j)
}
return(hess)}
######## Hessova matica multiplikatívneho modelu
#######################################
hessmat_EXP = function(alpha) {
hess = matrix(nrow = ncol(z), ncol = ncol(z))
d2l.dalphaij <- function(i,j) {
SIGMA <- sigma_EXP(alpha)
FI <- fi(SIGMA)
XS <- xs(SIGMA)
SX <- sx(SIGMA)
FIxs <- fixs(SIGMA)
sxFI <- t(FIxs)
sxFIxs <- sxfixs(SIGMA)
Pi <- pi_EXP(SIGMA,i)
Pj <- pi_EXP(SIGMA,j)
Qij <- qij_EXP(SIGMA,i,j)
Rij<- rij_EXP(SIGMA,i,j)
dSIGMA.dalphai <- dsigma.dalphai_EXP(SIGMA,i)
dSIGMA.dalphaj <- dsigma.dalphai_EXP(SIGMA,j)
d2SIGMA.dalphaij <- d2sigma.dalphaij_EXP(SIGMA,i,j)
SIGMAderiviSIGMA <- sigmaDERIVisigma_EXP(SIGMA,i)
SIGMAderivjSIGMA <- sigmaDERIVisigma_EXP(SIGMA,j)
SIGMAderivi <- sigmaDERIVi_EXP(SIGMA,i)
SIGMAderivj <- sigmaDERIVi_EXP(SIGMA,j)
d2l.dalphaij.11 = 1/2*tr(crossprod(t(SIGMAderivi),SIGMAderivj))
d2l.dalphaij.12 = -1/2*tr(crossprod(solve(SIGMA),d2SIGMA.dalphaij))
d2l.dalpahij.21 = 1/2*tr(crossprod(FI,crossprod(Pi,crossprod(FI,Pj))))
d2l.dalphaij.22 = -1*tr(crossprod(FI,Qij))
d2l.dalphaij.23 = 1/2*tr(crossprod(FI,Rij))
d2l.dalphaij.31 =
-1*crossprod(y,crossprod(t(SIGMAderivj),crossprod(SIGMAderiviSIGMA,y)))
d2l.dalphaij.32 =
1/2*crossprod(y,crossprod(solve(SIGMA),crossprod(d2SIGMA.dalphaij,crossprod(solve(SIGMA),y))))
d2l.dalphaij.41 =
crossprod(y,crossprod(t(SIGMAderivj),crossprod(FIxs,crossprod(Pi,crossprod(sxFI,y)))))
d2l.dalphaij.42 =
crossprod(y,crossprod(FIxs,crossprod(Pi,crossprod(FI,crossprod(Pj,crossprod(sxFI,y))))))
d2l.dalphaij.43 =
-1*crossprod(y,crossprod(FIxs,crossprod(Qij,crossprod(sxFI,y))))
d2l.dalphaij.44 =
1/2*crossprod(y,crossprod(FIxs,crossprod(Rij,crossprod(sxFI,y))))
d2l.dalphaij.51 =
crossprod(y,crossprod(t(SIGMAderivj),crossprod(t(SIGMAderivi),crossprod(sxFIxs,y))))
d2l.dalphaij.52 =
-1*crossprod(y,crossprod(solve(SIGMA),crossprod(d2SIGMA.dalphaij,crossprod(sxFIxs,y))))
d2l.dalphaij.53 =
crossprod(y,crossprod(t(SIGMAderivi),crossprod(t(SIGMAderivj),crossprod(sxFIxs,y))))
d2l.dalphaij.54 =
crossprod(y,crossprod(t(SIGMAderivi),crossprod(FIxs,crossprod(Pj,crossprod(sxFI,y)))))
d2l.dalphaij.55 =
crossprod(y,crossprod(t(SIGMAderivi),crossprod(sxFIxs,crossprod(SIGMAderivj,y))))
d2l.dalphaij = d2l.dalphaij.11 + d2l.dalphaij.12 + d2l.dalpahij.21 +
d2l.dalphaij.22 + d2l.dalphaij.23 +
d2l.dalphaij.31 + d2l.dalphaij.32 + d2l.dalphaij.41 + d2l.dalphaij.42 +
d2l.dalphaij.43 +
d2l.dalphaij.44 + d2l.dalphaij.51 + d2l.dalphaij.52 + d2l.dalphaij.53 +
d2l.dalphaij.54 +
d2l.dalphaij.55
return(d2l.dalphaij)}
for(i in 1:ncol(z)) { for(j in 1:ncol(z)) hess[i,j] = d2l.dalphaij(i,j)
}
return(hess)}
#########################################################################################
######## zavádza information matrices a kovariančné matice
aditívneho modelu ############
# expected information matrix
expmat_ADD <- function(alpha) { mat = matrix(nrow = ncol(z), ncol =
ncol(z))
matij <- function(i,j) {
SIGMA <- sigma_ADD(alpha)
FI <- fi(SIGMA)
Pi <- pi_ADD(SIGMA,i)
Pj <- pi_ADD(SIGMA,j)
Qij <- qij_ADD(SIGMA,i,j)
SIGMAderivi <- sigmaDERIVi_ADD(SIGMA,i)
SIGMAderivj <- sigmaDERIVi_ADD(SIGMA,j)
matij.1 = 1/2*tr(crossprod(t(SIGMAderivi),SIGMAderivj))
matij.2 =
-1*tr(crossprod(FI,Qij)-1/2*crossprod(FI,crossprod(Pi,crossprod(FI,Pj))))
matij = matij.1 + matij.2
return(matij)}
for(i in 1:ncol(z)) { for(j in 1:ncol(z)) mat[i,j] = matij(i,j) }
return(mat)}
# observed information matrix & average infromation matrix
obsmat_ADD <- function(alpha) -1*hessmat_ADD(alpha)
avgmat_ADD <- function(alpha) 1/2*(obsmat_ADD(alpha) +
expmat_ADD(alpha))
expcovmat_ADD <- function(alpha) solve(expmat_ADD(alpha))
obscovmat_ADD <- function(alpha) solve(obsmat_ADD(alpha))
avgcovmat_ADD <- function(alpha) solve(avgmat_ADD(alpha))
######## zavádza information matrices a kovariančné matice
zmiešaného modelu ############
# expected information matrix
expmat_MIX = function(alpha) { mat = matrix(nrow = ncol(z), ncol =
ncol(z))
matij <- function(i,j) {
SIGMA <- sigma_MIX(alpha)
FI <- fi(SIGMA)
Pi <- pi_MIX(SIGMA,i)
Pj <- pi_MIX(SIGMA,j)
Qij <- qij_MIX(SIGMA,i,j)
SIGMAderivi <- sigmaDERIVi_MIX(SIGMA,i)
SIGMAderivj <- sigmaDERIVi_MIX(SIGMA,j)
matij.1 = 1/2*tr(crossprod(t(SIGMAderivi),SIGMAderivj))
matij.2 =
-1*tr(crossprod(FI,Qij)-1/2*crossprod(FI,crossprod(Pi,crossprod(FI,Pj))))
matij = matij.1 + matij.2
return(matij)}
for(i in 1:ncol(z)) { for(j in 1:ncol(z)) mat[i,j] = matij(i,j) }
return(mat)}
# observed information matrix & average infromation matrix
obsmat_MIX <- function(alpha) -1*hessmat_MIX(alpha)
avgmat_MIX <- function(alpha) 1/2*(obsmat_MIX(alpha) +
expmat_MIX(alpha))
expcovmat_MIX <- function(alpha) solve(expmat_MIX(alpha))
obscovmat_MIX <- function(alpha) solve(obsmat_MIX(alpha))
avgcovmat_MIX <- function(alpha) solve(avgmat_MIX(alpha))
######## zavádza information matrices a kovariančné matice
mutliplikatívneho modelu #####
# expected information matrix
expmat_EXP = function(alpha) { mat = matrix(nrow = ncol(z), ncol =
ncol(z))
matij <- function(i,j) {
SIGMA <- sigma_EXP(alpha)
FI <- fi(SIGMA)
Pi <- pi_EXP(SIGMA,i)
Pj <- pi_EXP(SIGMA,j)
Qij <- qij_EXP(SIGMA,i,j)
SIGMAderivi <- sigmaDERIVi_EXP(SIGMA,i)
SIGMAderivj <- sigmaDERIVi_EXP(SIGMA,j)
matij.1 = 1/2*tr(crossprod(t(SIGMAderivi),SIGMAderivj))
matij.2 =
-1*tr(crossprod(FI,Qij)-1/2*crossprod(FI,crossprod(Pi,crossprod(FI,Pj))))
matij = matij.1 + matij.2
return(matij)}
for(i in 1:ncol(z)) { for(j in 1:ncol(z)) mat[i,j] = matij(i,j) }
return(mat)}
# observed information matrix & average infromation matrix
obsmat_EXP <- function(alpha) -1*hessmat_EXP(alpha)
avgmat_EXP <- function(alpha) 1/2*(obsmat_EXP(alpha) +
expmat_EXP(alpha))
expcovmat_EXP <- function(alpha) solve(expmat_EXP(alpha))
obscovmat_EXP <- function(alpha) solve(obsmat_EXP(alpha))
avgcovmat_EXP <- function(alpha) solve(avgmat_EXP(alpha))
#########################################################################################
# zavádza knižnicu maxLik
library(maxLik)
# zavedie štartovacie hodnoty jednotky pre alpha
beta_ols <- function(x,y) crossprod(solve(crossprod(x)), crossprod(x,y))
resids_ols <- function(x,y) (y - crossprod(t(x),beta_ols(x,y)))
alpha_ols_ADD <- function(x,y)
crossprod(solve(crossprod(z)),crossprod(z,resids_ols(x,y)^2))
alpha_ols_MIX <- function(x,y)
crossprod(solve(crossprod(z)),crossprod(z,abs(resids_ols(x,y))))
alpha_ols_EXP <- function(x,y)
crossprod(solve(crossprod(z)),crossprod(z,log(resids_ols(x,y)^2)))
coefolsEST_ADD <- function(x,y) t(alpha_ols_ADD(x,y))
coefolsEST_MIX <- function(x,y) t(alpha_ols_MIX(x,y))
coefolsEST_EXP <- function(x,y) t(alpha_ols_EXP(x,y))
kenward_roger_both <- function(model,estimated_alpha,beta_tested,x,y,z)
{ x = x; y = y; z = z
if (model == "ADD") {
inner <- function(i,j,alpha) {
FI <- fi(sigma_ADD(alpha))
Pi <- pi_ADD(sigma_ADD(alpha),i)
Pj <- pi_ADD(sigma_ADD(alpha),j)
Qij <- qij_ADD(sigma_ADD(alpha),i,j)
W <- expcovmat_ADD(alpha)
inner = matrix(nrow = ncol(x), ncol = ncol(x))
inner = W[i,j]*(Qij-crossprod(Pi,crossprod(FI,Pj)))
return(inner) }
INNER = matrix(0, nrow = ncol(x), ncol = ncol(x))
for(i in 1:ncol(z)) { for(j in 1:ncol(z)) INNER = INNER +
inner(i,j,estimated_alpha) }
a1a2ij <- function(i,j,alpha) {
FI <- fi(sigma_ADD(alpha))
Pi <- pi_ADD(sigma_ADD(alpha),i)
Pj <- pi_ADD(sigma_ADD(alpha),j)
W <- expcovmat_ADD(alpha)
a1a2 = vector(mode = "numeric",length = 2)
a1a2[1] <- W[i,j]*tr(crossprod(FI,Pi))*tr(crossprod(FI,Pj))
a1a2[2] <- W[i,j]*tr(crossprod(crossprod(Pi,FI),crossprod(FI,Pj)))
return(a1a2) }
A1A2 = vector(mode = "numeric", length = 2)
for(i in 1:ncol(z)) { for(j in 1:ncol(z)) A1A2 = A1A2 +
a1a2ij(i,j,estimated_alpha) }
A1 = A1A2[1]
A2 = A1A2[2]
FI <- fi(sigma_ADD(estimated_alpha))
p <- ncol(x)
EF <- 1 + 1 / p * A2
DF <- 1 / p^2 * (A1 + 6 * A2 + 2 * p)
m <- 4 + (p + 2) / (p * DF / (2 * EF^2) - 1)
lam <- m / (EF * (m - 2))
beta_egls <-
crossprod(FI,crossprod(x,crossprod(solve(sigma_ADD(estimated_alpha)),y)))
beta_covmat_old <- FI + 2*crossprod(FI,crossprod(INNER,FI))
dif = beta_egls - beta_tested
F.kr = as.numeric(lam / p *
crossprod(dif,crossprod(solve(beta_covmat_old),dif)))
p.v <- pf(F.kr, df1 = p, df2 = m, lower.tail = FALSE, log.p = FALSE)
list <- list(A1 = A1, A2 = A2, EF = EF, DF = DF, df1 = p, df2 = m, F_old
= F.kr,
F_new = F.kr, sig_old = p.v, sig_new = p.v, beta_egls = beta_egls,
beta_covmat_old = beta_covmat_old, beta_covmat_new =
beta_covmat_old,model = model)
list }
else if (model == "MIX") {
inner <- function(i,j,alpha) {
FI <- fi(sigma_MIX(alpha))
Pi <- pi_MIX(sigma_MIX(alpha),i)
Pj <- pi_MIX(sigma_MIX(alpha),j)
Qij <- qij_MIX(sigma_MIX(alpha),i,j)
Rij <- rij_MIX(sigma_MIX(alpha),i,j)
W <- expcovmat_MIX(alpha)
inner = matrix(nrow = ncol(x), ncol = ncol(x))
inner = W[i,j]*(Qij-crossprod(Pi,crossprod(FI,Pj))-1/4*Rij)
return(inner) }
INNER = matrix(0, nrow = ncol(x), ncol = ncol(x))
for(i in 1:ncol(z)) { for(j in 1:ncol(z)) INNER = INNER +
inner(i,j,estimated_alpha) }
a1a2ij <- function(i,j,alpha) {
FI <- fi(sigma_MIX(alpha))
Pi <- pi_MIX(sigma_MIX(alpha),i)
Pj <- pi_MIX(sigma_MIX(alpha),j)
W <- expcovmat_MIX(alpha)
a1a2 = vector(mode = "numeric",length = 2)
a1a2[1] <- W[i,j]*tr(crossprod(FI,Pi))*tr(crossprod(FI,Pj))
a1a2[2] <- W[i,j]*tr(crossprod(crossprod(Pi,FI),crossprod(FI,Pj)))
return(a1a2) }
A1A2 = vector(mode = "numeric", length = 2)
for(i in 1:ncol(z)) { for(j in 1:ncol(z)) A1A2 = A1A2 +
a1a2ij(i,j,estimated_alpha) }
A1 = A1A2[1]
A2 = A1A2[2]
BIAS_COR_MIX <- function(alpha) {
SIGMA <- sigma_MIX(alpha)
FI <- fi(sigma_MIX(alpha))
W <- expcovmat_MIX(alpha)
IN = matrix(0, nrow = nrow(x), ncol = nrow(x))
for(i in 1:ncol(z)) {
for(j in 1:ncol(z)) IN = IN + W[i,j]*d2sigma.dalphaij_MIX(i,j) } #
is symmetric
S <- crossprod(solve(SIGMA),crossprod(IN,solve(SIGMA))) # is
symmetric
V <- function(t) {
ancillary1 <- crossprod(S,dsigma.dalphai_MIX(SIGMA,t)) # is
symmetric
ancillary2 <- crossprod(x,crossprod(S,x)) # is symmetric
v1 <- tr(ancillary1)
v2 <-
-2*tr(crossprod(sx(SIGMA),crossprod(ancillary1,crossprod(t(x),FI))))
v3 <-
-1*tr(crossprod(ancillary2,crossprod(FI,crossprod(pi_MIX(SIGMA,t),FI))))
v <- v1 + v2 + v3
return(v) }
bias = matrix(0, nrow = ncol(x), ncol = ncol(x))
alpha = estimated_alpha
for(t in 1:ncol(z)) {
for(s in 1:ncol(z))
bias = bias +
(-1/4)*W[t,s]*V(t)*crossprod(FI,crossprod(pi_MIX(SIGMA,s),FI)) }
return(bias) }
FI <- fi(sigma_MIX(estimated_alpha))
p <- ncol(x)
EF <- 1 + 1 / p * A2
DF <- 1 / p^2 * (A1 + 6 * A2 + 2 * p)
m <- 4 + (p + 2) / (p * DF / (2 * EF^2) - 1)
lam <- m / (EF * (m - 2))
beta_egls <-
crossprod(FI,crossprod(x,crossprod(solve(sigma_MIX(estimated_alpha)),y)))
beta_covmat_old <- FI + 2*crossprod(FI,crossprod(INNER,FI))
beta_covmat_new <- beta_covmat_old + BIAS_COR_MIX(estimated_alpha)
dif = beta_egls - beta_tested
F.kr_old = as.numeric(lam / p *
crossprod(dif,crossprod(solve(beta_covmat_old),dif)))
p.v_old <- pf(F.kr_old, df1 = p, df2 = m, lower.tail = FALSE, log.p =
FALSE)
F.kr_new = as.numeric(lam / p *
crossprod(dif,crossprod(solve(beta_covmat_new),dif)))
p.v_new <- pf(F.kr_new, df1 = p, df2 = m, lower.tail = FALSE, log.p =
FALSE)
list <- list(A1 = A1, A2 = A2, EF = EF, DF = DF, df1 = p, df2 = m, F_old
= F.kr_old,
F_new = F.kr_new, sig_old = p.v_old, sig_new = p.v_new, beta_egls =
beta_egls,
beta_covmat_old = beta_covmat_old, beta_covmat_new =
beta_covmat_new,model = model)
list }
else if (model == "EXP") {
inner <- function(i,j,alpha) {
FI <- fi(sigma_EXP(alpha))
Pi <- pi_EXP(sigma_EXP(alpha),i)
Pj <- pi_EXP(sigma_EXP(alpha),j)
Qij <- qij_EXP(sigma_EXP(alpha),i,j)
Rij <- rij_EXP(sigma_EXP(alpha),i,j)
W <- expcovmat_EXP(alpha)
inner = matrix(nrow = ncol(x), ncol = ncol(x))
inner = W[i,j]*(Qij-crossprod(Pi,crossprod(FI,Pj))-1/4*Rij)
return(inner) }
INNER = matrix(0, nrow = ncol(x), ncol = ncol(x))
for(i in 1:ncol(z)) { for(j in 1:ncol(z)) INNER = INNER +
inner(i,j,estimated_alpha) }
a1a2ij <- function(i,j,alpha) {
FI <- fi(sigma_EXP(alpha))
Pi <- pi_EXP(sigma_EXP(alpha),i)
Pj <- pi_EXP(sigma_EXP(alpha),j)
W <- expcovmat_EXP(alpha)
a1a2 = vector(mode = "numeric",length = 2)
a1a2[1] <- W[i,j]*tr(crossprod(FI,Pi))*tr(crossprod(FI,Pj))
a1a2[2] <- W[i,j]*tr(crossprod(crossprod(Pi,FI),crossprod(FI,Pj)))
return(a1a2) }
A1A2 = vector(mode = "numeric", length = 2)
for(i in 1:ncol(z)) { for(j in 1:ncol(z)) A1A2 = A1A2 +
a1a2ij(i,j,estimated_alpha) }
A1 = A1A2[1]
A2 = A1A2[2]
BIAS_COR_EXP <- function(alpha) {
SIGMA <- sigma_EXP(alpha)
FI <- fi(sigma_EXP(alpha))
W <- expcovmat_EXP(alpha)
IN = matrix(0, nrow = nrow(x), ncol = nrow(x))
for(i in 1:ncol(z)) {
for(j in 1:ncol(z)) IN = IN + W[i,j]*d2sigma.dalphaij_EXP(SIGMA,i,j) }
# is symmetric
S <- crossprod(solve(SIGMA),crossprod(IN,solve(SIGMA))) # is
symmetric
V <- function(t) {
ancillary1 <- crossprod(S,dsigma.dalphai_EXP(SIGMA,t)) # is
symmetric
ancillary2 <- crossprod(x,crossprod(S,x)) # is symmetric
v1 <- tr(ancillary1)
v2 <-
-2*tr(crossprod(sx(SIGMA),crossprod(ancillary1,crossprod(t(x),FI))))
v3 <-
-1*tr(crossprod(ancillary2,crossprod(FI,crossprod(pi_EXP(SIGMA,t),FI))))
v <- v1 + v2 + v3
return(v) }
bias = matrix(0, nrow = ncol(x), ncol = ncol(x))
for(t in 1:ncol(z)) {
for(s in 1:ncol(z))
bias = bias +
(-1/4)*W[t,s]*V(t)*crossprod(FI,crossprod(pi_EXP(SIGMA,s),FI)) }
return(bias) }
FI <- fi(sigma_EXP(estimated_alpha))
p <- ncol(x)
EF <- 1 + 1 / p * A2
DF <- 1 / p^2 * (A1 + 6 * A2 + 2 * p)
m <- 4 + (p + 2) / (p * DF / (2 * EF^2) - 1)
lam <- m / (EF * (m - 2))
beta_egls <-
crossprod(FI,crossprod(x,crossprod(solve(sigma_EXP(estimated_alpha)),y)))
beta_covmat_old <- FI + 2*crossprod(FI,crossprod(INNER,FI))
beta_covmat_new <- beta_covmat_old + BIAS_COR_EXP(estimated_alpha)
dif = beta_egls - beta_tested
F.kr_old = as.numeric(lam / p *
crossprod(dif,crossprod(solve(beta_covmat_old),dif)))
p.v_old <- pf(F.kr_old, df1 = p, df2 = m, lower.tail = FALSE, log.p =
FALSE)
F.kr_new = as.numeric(lam / p *
crossprod(dif,crossprod(solve(beta_covmat_new),dif)))
p.v_new <- pf(F.kr_new, df1 = p, df2 = m, lower.tail = FALSE, log.p =
FALSE)
list <- list(A1 = A1, A2 = A2, EF = EF, DF = DF, df1 = p, df2 = m, F_old
= F.kr_old,
F_new = F.kr_new, sig_old = p.v_old, sig_new = p.v_new, beta_egls =
beta_egls,
beta_covmat_old = beta_covmat_old, beta_covmat_new =
beta_covmat_new,model = model)
list }
else return(cat("The wrong argument (model) - must be either >ADD<, or
>MIX<, or >EXP<.\n")) }
kenward_roger_both_print <- function(krprocedure) {
A1 = krprocedure$A1; A2 = krprocedure$A2; EF = krprocedure$EF; DF =
krprocedure$DF
p = krprocedure$df1; m = krprocedure$df2; F.kr_old = krprocedure$F_old
F.kr_new = krprocedure$F_new; p.v_new= krprocedure$sig_new; p.v_old =
krprocedure$sig_old
beta_egls = krprocedure$beta_egls; beta_covmat_new =
krprocedure$beta_covmat_new
beta_covmat_old = krprocedure$beta_covmat_old; beta_egls =
krprocedure$beta_egls
if (krprocedure$model == "ADD") {
summary1 <-
cbind(data.frame("Parameter"=c("Estimates")),BETA=t(as.vector(beta_egls)))
summary2 <- cbind(BETA=beta_covmat_old)
cat("**Parameter estimates of the additive heteroskedasticity
model**\n")
cat("-----------------------------------------------------------------------------\n")
print(summary1,row.names=FALSE)
cat("-----------------------------------------------------------------------------\n")
cat("**Both original 1997 and new 2009 Kenward-Roger procedure**\n")
cat(" \n")
cat("**Kenward-Roger small-sample adjusted estimate of the
covmat(beta)**\n")
print(summary2,row.names=FALSE)
cat(" \n")
cat("Small-sample adjusted F-statistics =",F.kr_old,"\n")
cat("df1 (nominator) = ",p,", df2 (denominator) = ",m,", Significance
=",p.v_old,"\n")
cat("-----------------------------------------------------------------------------\n")
}
else if (krprocedure$model == "MIX") {
summary1 <-
cbind(data.frame("Parameter"=c("Estimates")),BETA=t(as.vector(beta_egls)))
summary2_old <- cbind(BETA=beta_covmat_old)
summary2_new <- cbind(BETA=beta_covmat_new)
cat("**Parameter estimates of the mixed heteroskedasticity model**\n")
cat("-----------------------------------------------------------------------------\n")
print(summary1,row.names=FALSE)
cat("-----------------------------------------------------------------------------\n")
cat("**Original 1997 Kenward-Roger procedure**\n")
cat(" \n")
cat("**Kenward-Roger small-sample adjusted estimate of the
covmat(beta)**\n")
print(summary2_old,row.names=FALSE)
cat(" \n")
cat("Small-sample adjusted F-statistics =",F.kr_old,"\n")
cat("df1 (nominator) = ",p,", df2 (denominator) = ",m,", Significance
=",p.v_old,"\n")
cat("-----------------------------------------------------------------------------\n")
cat("**New 2009 Kenward-Roger procedure**\n")
cat(" \n")
cat("**Kenward-Roger small-sample adjusted estimate of the
covmat(beta)**\n")
print(summary2_new,row.names=FALSE)
cat(" \n")
cat("Small-sample adjusted F-statistics =",F.kr_new,"\n")
cat("df1 (nominator) = ",p,", df2 (denominator) = ",m,", Significance
=",p.v_new,"\n")
cat("-----------------------------------------------------------------------------\n")
}
else {
summary1 <-
cbind(data.frame("Parameter"=c("Estimates")),BETA=t(as.vector(beta_egls)))
summary2_old <- cbind(BETA=beta_covmat_old)
summary2_new <- cbind(BETA=beta_covmat_new)
cat("**Parameter estimates of the multiplicative heteroskedasticity
model**\n")
cat("-----------------------------------------------------------------------------\n")
print(summary1,row.names=FALSE)
cat("-----------------------------------------------------------------------------\n")
cat("**Original 1997 Kenward-Roger procedure**\n")
cat(" \n")
cat("**Kenward-Roger small-sample adjusted estimate of the
covmat(beta)**\n")
print(summary2_old,row.names=FALSE)
cat(" \n")
cat("Small-sample adjusted F-statistics =",F.kr_old,"\n")
cat("df1 (nominator) = ",p,", df2 (denominator) = ",m,", Significance
=",p.v_old,"\n")
cat("-----------------------------------------------------------------------------\n")
cat("**New 2009 Kenward-Roger procedure**\n")
cat(" \n")
cat("**Kenward-Roger small-sample adjusted estimate of the
covmat(beta)**\n")
print(summary2_new,row.names=FALSE)
cat(" \n")
cat("Small-sample adjusted F-statistics =",F.kr_new,"\n")
cat("df1 (nominator) = ",p,", df2 (denominator) = ",m,", Significance
=",p.v_new,"\n")
cat("-----------------------------------------------------------------------------\n")
}}
######## generovanie vektora ypsilonov
##################################################
generate_y <- function(x,z,beta,alpha,model) { res =
as.matrix(rnorm(nrow(x)))
if (model == "ADD") {
sigma_add = diag(as.vector(alpha %*% t(z)), nrow=nrow(z), ncol=nrow(z))
y = crossprod(t(x),beta) + crossprod(sigma_add^(1/2),res)
return(y) }
else if (model == "MIX") {
sigma_mix = diag((as.vector(alpha %*% t(z)))^2, nrow=nrow(z),
ncol=nrow(z))
y = crossprod(t(x),beta) + crossprod(sigma_mix^(1/2),res)
return(y) }
else if (model == "EXP") {
sigma_exp = diag(exp(as.vector(alpha %*% t(z))), nrow=nrow(z),
ncol=nrow(z))
y = crossprod(t(x),beta) + crossprod(sigma_exp^(1/2),res)
return(y) }
else return(cat("The wrong argument (type)- must be either >ADD<, or
>MIX<, or >EXP<.\n")) }
######## pomocné funkcie pre simulovanie
################################################
loglk <- function(simtype) { if (simtype == "ADD") reml.loglik_ADD
else if (simtype == "MIX") reml.loglik_MIX
else reml.loglik_EXP }
gradvec <- function(simtype) { if (simtype == "ADD") gradvec_ADD
else if (simtype == "MIX") gradvec_MIX
else gradvec_EXP }
hessmat <- function(simtype) { if (simtype == "ADD") hessmat_ADD
else if (simtype == "MIX") hessmat_MIX
else hessmat_EXP }
starter <- function(simtype,x,y) { if (simtype == "ADD")
abs(coefolsEST_ADD(x,y))
else if (simtype == "MIX") coefolsEST_MIX(x,y)
else coefolsEST_EXP(x,y) }
existencecheck <- function(object)
exists(as.character(substitute(object)))
#######################################################################################
x <- as.matrix(read.table("C:/maddala_70.txt", header = F)) # dané
pevné x
z <- as.matrix(read.table("C:/maddala_70.txt", header = F)) # dané
pevné z
nsimsgood=5000; beta=c(0.85,0.90); beta_tested=c(0.85,0.90)
alpha=c(1, 0.50); realtype="ADD"; simtype="ADD"; prob=0.05
heading = matrix(nrow = 1, ncol = 15 + 2*ncol(z) + 2*ncol(x)); i = 0;
j = 0;
name1 <-
c("SimNo","realHET","simHET",rep("realALPHA",ncol(z)),rep("estALPHA",ncol(z)))
name2 <-
c("estALPHAlikTYPE",rep("realBETA",ncol(x)),rep("estBETA",ncol(x)))
name3 <- c("A1","A2","EF","DF","df1","df2")
name4 <-
c("Fstatisticsold","Fstatisticsnew","Fquantile","inAreaOld","inAreaNew")
heading[1,] <- c(name1,name2,name3,name4)
write.table(heading, file="70addadd__goods.csv", append=T, quote=T,
sep=";", row.names=F, col.names=F)
write.table(c("Unsuccessful simulations"), file="70addadd__bads.csv",
append=T, quote=T, sep=";", row.names=F, col.names=F)
repeat {
y <<- generate_y(x,z,beta,alpha,realtype); x <<- x; z <<- z
row = matrix(nrow = 1, ncol = 15 + 2*ncol(z) + 2*ncol(x))
NR <-
try(maxNR(loglk(simtype),gradvec(simtype),hessmat(simtype),start=starter(simtype,x,y)),silent=TRUE)
if (is.list(NR) == TRUE) kr_nr <-
kenward_roger_both(simtype,NR$estimate,beta_tested,x,y,z)
if (is.list(NR) == TRUE) NRcode <- NR$code else NRcode <- 99
if (existencecheck(kr_nr) == TRUE) kr_nrF_old <- kr_nr$F_old else
kr_nrF_old <- -99
if (existencecheck(kr_nr) == TRUE) kr_nrF_new <- kr_nr$F_new else
kr_nrF_new <- -99
if (existencecheck(kr_nr) == TRUE) qff <-
qf(prob,df1=kr_nr$df1,df2=kr_nr$df2,lower.tail=TRUE) else qff <- NaN
{ if ((NRcode == 0 | NRcode == 1) & (kr_nrF_old >= 0 & kr_nrF_new >= 0 &
is.nan(qff) == FALSE)) { i = i + 1
A1 <- kr_nr$A1; A2 <- kr_nr$A2; EF <- kr_nr$EF; DF <- kr_nr$DF
df1 <- kr_nr$df1; df2 <- kr_nr$df2
Fold <- kr_nr$F_old; Fnew <- kr_nr$F_new
f <- qf(prob,df1=kr_nr$df1,df2=kr_nr$df2,lower.tail=TRUE)
rowi1 = c(i, realtype, simtype, alpha, NR$estimate,
c("Nephton-Raphson"), beta)
rowi2 = c(t(kr_nr$beta_egls), A1, A2, EF, DF, df1, df2)
rowi3 = c(Fold, Fnew, f, Fold < f, Fnew < f)
row[1,] = c(rowi1,rowi2,rowi3)
write.table(row,file="70addadd__goods.csv",append=T,quote=T,sep=";",row.names=F,col.names=F)
rm(y, NR, kr_nr, NRcode, kr_nrF_old, kr_nrF_new, qff) }
else { BFGS <-
try(maxBFGS(loglk(simtype),gradvec(simtype),hessmat(simtype),start=starter(simtype,x,y)),silent=TRUE)
if (is.list(BFGS) == TRUE) kr_bfgs <-
kenward_roger_both(simtype,BFGS$estimate,beta_tested,x,y,z)
if (is.list(BFGS) == TRUE) BFGScode <- BFGS$code else BFGScode <- 99
if (existencecheck(kr_bfgs) == TRUE) kr_bfgsF_old <- kr_bfgs$F_old else
kr_bfgsF_old <- -99
if (existencecheck(kr_bfgs) == TRUE) kr_bfgsF_new <- kr_bfgs$F_new else
kr_bfgsF_new <- -99
if (existencecheck(kr_bfgs) == TRUE) qff <-
qf(prob,df1=kr_bfgs$df1,df2=kr_bfgs$df2,lower.tail=TRUE) else qff <- NaN
{ if ((BFGScode == 0) & (kr_bfgsF_old >= 0 & kr_bfgsF_new >= 0 &
is.nan(qff) == FALSE)) { i = i + 1
A1 <- kr_bfgs$A1; A2 <- kr_bfgs$A2; EF <- kr_bfgs$EF; DF <- kr_bfgs$DF
df1 <- kr_bfgs$df1; df2 <- kr_bfgs$df2
Fold <- kr_bfgs$F_old; Fnew <- kr_bfgs$F_new
f <- qf(prob,df1=kr_bfgs$df1,df2=kr_bfgs$df2,lower.tail=TRUE)
rowi1 = c(i, realtype, simtype, alpha, BFGS$estimate, c("BFGS"), beta)
rowi2 = c(t(kr_bfgs$beta_egls), A1, A2, EF, DF, df1, df2)
rowi3 = c(Fold, Fnew, f, Fold < f, Fnew < f)
row[1,] = c(rowi1,rowi2,rowi3)
write.table(row,file="70addadd__goods.csv",append=T,quote=T,sep=";",row.names=F,col.names=F)
rm(y, BFGS, kr_bfgs, BFGScode, kr_bfgsF_old, kr_bfgsF_new, qff) }
else { NM <-
try(maxNM(loglk(simtype),gradvec(simtype),hessmat(simtype),start=starter(simtype,x,y)),silent=TRUE)
if (is.list(NM) == TRUE) kr_nm <-
kenward_roger_both(simtype,NM$estimate,beta_tested,x,y,z)
if (is.list(NM) == TRUE) NMcode <- NM$code else NMcode <- 99
if (existencecheck(kr_nm) == TRUE) kr_nmF_old <- kr_nm$F_old else
kr_nmF_old <- -99
if (existencecheck(kr_nm) == TRUE) kr_nmF_new <- kr_nm$F_new else
kr_nmF_new <- -99
if (existencecheck(kr_nm) == TRUE) qff <-
qf(prob,df1=kr_nm$df1,df2=kr_nm$df2,lower.tail=TRUE) else qff <- NaN
{ if ((NMcode == 0) & (kr_nmF_old >= 0 & kr_nmF_new >= 0 & is.nan(qff)
== FALSE)) { i = i + 1
A1 <- kr_nm$A1; A2 <- kr_nm$A2; EF <- kr_nm$EF; DF <- kr_nm$DF
df1 <- kr_nm$df1; df2 <- kr_nm$df2
Fold <- kr_nm$F_old; Fnew <- kr_nm$F_new
f <- qf(prob,df1=kr_nm$df1,df2=kr_nm$df2,lower.tail=TRUE)
rowi1 = c(i, realtype, simtype, alpha, NM$estimate, c("Nelder-Mead"),
beta)
rowi2 = c(t(kr_nm$beta_egls), A1, A2, EF, DF, df1, df2)
rowi3 = c(Fold, Fnew, f, Fold < f, Fnew < f)
row[1,] = c(rowi1,rowi2,rowi3)
write.table(row,file="70addadd__goods.csv",append=T,quote=T,sep=";",row.names=F,col.names=F)
rm(y, NN, kr_nm, NMcode, kr_nmF_old, kr_nmF_new, qff) }
else { SANN <-
try(maxSANN(loglk(simtype),gradvec(simtype),hessmat(simtype),start=starter(simtype,x,y)),silent=TRUE)
if (is.list(SANN) == TRUE) kr_sann <-
kenward_roger_both(simtype,SANN$estimate,beta_tested,x,y,z)
if (is.list(SANN) == TRUE) SANNcode <- SANN$code else SANNcode <- 99
if (existencecheck(kr_sann) == TRUE) kr_sannF_old <- kr_sann$F_old else
kr_sannF_old <- -99
if (existencecheck(kr_sann) == TRUE) kr_sannF_new <- kr_sann$F_new else
kr_sannF_new <- -99
if (existencecheck(kr_sann) == TRUE) qff <-
qf(prob,df1=kr_sann$df1,df2=kr_sann$df2,lower.tail=TRUE) else qff <- NaN
{ if ((SANNcode == 0) & (kr_sannF_old >= 0 & kr_sannF_new >= 0 &
is.nan(qff) == FALSE)) { i = i + 1
A1 <- kr_sann$A1; A2 <- kr_sann$A2; EF <- kr_sann$EF; DF <- kr_sann$DF
df1 <- kr_sann$df1; df2 <- kr_sann$df2
Fold <- kr_sann$F_old; Fnew <- kr_sann$F_new
f <- qf(prob,df1=kr_sann$df1,df2=kr_sann$df2,lower.tail=TRUE)
rowi1 = c(i, realtype, simtype, alpha, SANN$estimate, c("SANN"), beta)
rowi2 = c(t(kr_sann$beta_egls), A1, A2, EF, DF, df1, df2)
rowi3 = c(Fold, Fnew, f, Fold < f, Fnew < f)
row[1,] = c(rowi1,rowi2,rowi3)
write.table(row,file="70addadd__goods.csv",append=T,quote=T,sep=";",row.names=F,col.names=F)
rm(y, SANN, kr_sann, SANNcode, kr_sannF_old, kr_sannF_new, qff) }
else { j = j + 1
write.table(j,file="70addadd__bads.csv",append=T,quote=T,sep=";",row.names=F,col.names=F)
rm(y, NR, BFGS, NM, SANN, NRcode, BFGScode, NMcode, SANNcode,
kr_nrF_old, kr_bfgsF_old, kr_nmF_old, kr_sannF_old)
rm(qff, kr_nrF_new, kr_bfgsF_new, kr_nmF_new, kr_sannF_new)
if (existencecheck(kr_nr) == TRUE) rm(kr_nr)
if (existencecheck(kr_bfgs) == TRUE) rm(kr_bfgs)
if (existencecheck(kr_nm) == TRUE) rm(kr_nm)
if (existencecheck(kr_sann) == TRUE) rm(kr_sann) } } } } } } } } }
if (i == nsimsgood) break
cat("-----------------------------------------------------------------------------\n")
cat("Number of successful good simulations = ",i,"\n")
cat("Number of unsuccessful simulations = ",j,"\n")
cat("-----------------------------------------------------------------------------\n")
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
