R help - Mar 2010 - R runs in a usual way, but simulations are not performed

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")

The script you sent seems to be too large for me to spend time debugging,
plus it seems that there are files being read that are not part of the
data.  I would suggest that if you have a long running program that you put
messages (cat('reached part 1\n')) throughout the program so you can
trace
what it is doing.  R typically does not stop on its own; it is only
following the directions of the program you have written and there is a bug
in it somewhere.  So you will have to apply some basic debugging skills to
help you track it down.

There are a lot of things to look for.  Does it take longer doing through
each loop (e.g., you have a vector that is growing and taking longer to
access).  Is the system paging because you have run out of physical memory?
Do you have some condition in a 'while' loop that is never satisfied? 
Only
you will be able to determine this after you have put some trace statements
in your program and localized where the problem is.

On Sat, Mar 27, 2010 at 7:19 AM, Nomen Omen <hcd@azet.sk> wrote:
> 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<http://www.r-project.org/posting-guide.html>
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Jim Holtman
Cincinnati, OH
+1 513 646 9390

What is the problem that you are trying to solve?

	[[alternative HTML version deleted]]