Nothing
#library(mvtnorm); library(lqr); library(quantreg); library(psych); library(ald); library(progress)
QSAEM_NL = function(y,x,nj,initial,exprNL,covar=NA,p=0.5,
precision = 0.0001,M=20,pc=0.25,MaxIter=500,
beta,sigmae,D,nlmodel,d,q){
start.time <- Sys.time()
n = length(nj)
N = sum(nj)
beta = c(1.24332278, -0.03787132, 1.42456439, 0.02953380, 0.40366525,
2.67806346, 0.03564155)
#NEW CODE
if(all(is.na(covar)==TRUE)){n.covar = 0}else{n.covar = dim(covar)[2]}
delta1 = 0.001
delta2 = precision
#assymetry
vp = (1-2*p)/(p*(1-p))
tp = sqrt(2/(p*(1-p)))
MDel = MElim(q)
ndiag = (q*(1+q)/2)
npar = d+1+ndiag
###################################
exprFX = gsub("|]","",gsub("[|[]","",as.character(exprNL)))
paste1 = paste("\"fixed",1:d,"\"",sep = "",collapse = ",")
paste2 = paste("fixed",1:d,sep = "",collapse = ",")
pasteA = paste("covar",1:n.covar,sep = "",collapse = ",")
paste3 = ifelse(n.covar==0,"",pasteA)
paste4 = ifelse(n.covar==0,paste("x",paste2,sep =","),paste("x",paste2,paste3,sep =","))
paste5 = paste("random",1:q,sep = "",collapse = ",")
paste6 = paste("deriv( ~ ",exprFX,",c(",paste1,"),function(",paste(paste4,",",paste5,sep = ""),"){})",sep = "")
dd = eval(parse(text = paste6))
critval = 1
count = 0
teta = c(beta,sigmae,D[upper.tri(D, diag = T)])
tetam = matrix(data=NA,nrow=npar,ncol=MaxIter)
EPV = matrix(0,nrow = npar,ncol = MaxIter)
if(pc==1){
seqq=rep(1,pc*MaxIter)
}else{
seqq = c(rep(1,pc*MaxIter),(1/((((pc*MaxIter)+1):MaxIter)-(pc*MaxIter))))
seqq = c(rep(1,MaxIter-length(seqq)),seqq)
}
SAEM_bb = array(data=0,dim=c(MaxIter+1,n,q,q))
SAEM_bb2 = array(data=0,dim=c(MaxIter+1,n,q,q))
SAEM_bi = array(data=0,dim=c(MaxIter+1,n,q))
SAEM_VGi = array(data=0,dim=c(MaxIter+1,n,npar))
SAEM_sig = array(data=0,dim=c(MaxIter+1,n))
SAEM_Hbeta = array(data=0,dim=c(MaxIter+1,n,d,d))
for(j in 1:n){SAEM_bb[count+1,j,,] = diag(q)}
# Progress bar ------------------------------------------------------------
# pb = tkProgressBar(title = "QRNLMM via SAEM", min = 0,max = MaxIter, width = 300)
# setTkProgressBar(pb, 0, label=paste("Iter ",0,"/",MaxIter," - ",0,"% done",sep = ""))
#
pb <- progress_bar$new(
format = ":what [:bar] :percent eta: :eta",
total = MaxIter, clear = FALSE,
width= 60, show_after = 0)
pb$tick(len = -1,tokens = list(what = "Working..."))
while(critval < 3)
{
count = count + 1
sumb1 = matrix(data=0,nrow=d,ncol=1)
sumb2 = matrix(data=0,nrow=d,ncol=d)
sumD = matrix(data=0,nrow=q,ncol=q)
sumsig = 0
IE = 0
x = as.matrix(x)
for (j in 1:n)
{
print(j)
y1 = y[(sum(nj[1:j-1])+1):(sum(nj[1:j]))]
x1 = matrix(x[(sum(nj[1:j-1])+1):(sum(nj[1:j])),],ncol=1)
if(n.covar>0){cov1 = matrix(covar[(sum(nj[1:j-1])+1):(sum(nj[1:j])),],ncol=n.covar)}
v1s = matrix(data=1,nrow=nj[j],ncol=1)
##########################################################################
#PASSO E
##########################################################################
ui = matrix(0,nrow = nj[j],ncol = 1)
Dui = matrix(0,nrow = nj[j],ncol = nj[j])
VGi = matrix(data = 0,nrow = npar,ncol = M)
Msig = array(NA,M)
Mbb = array(NA,dim = c(q,q,M))
Mbb2 = array(NA,dim = c(q,q,M))
Hbeta = array(NA,dim = c(d,d,M))
bmetro = matrix(MHbi2(j=j,M=M,y1,x1,cov1,
bi=as.matrix(SAEM_bi[count,j,]),
bibi=as.matrix(SAEM_bb[count,j,,]),
d=d,q=q,p=p,nj=nj,
beta=beta,sigmae=sigmae,D=D,
nlmodel=nlmodel),q,M)
paste6 = paste("beta[",1:d,"]",sep = "",collapse = ",")
paste7 = ifelse(n.covar==0,"",paste(",",paste("cov1[,",1:n.covar,"]",sep = "",collapse = ","),sep=""))
paste8 = paste("attr(dd(x1,",paste6,paste7,",",paste(rep(0,q),sep="",collapse = ","),"),\"gradient\")",sep = "")
grad = eval(parse(text = paste8))
for(l in 1:M)
{
for(k in 1:nj[j])
{
chi = ((y1[k] - nlmodel(x1[k],beta,bmetro[,l],cov1[k,]))^2)/(sigmae*tp^2)
if(chi==0){chi=10^-100}
psi = (tp^2)/(4*sigmae)
ui[k] = Egig(lambda = 0.5,chi = chi,psi = psi,func = "x")
Dui[k,k] = Egig(lambda = 0.5,chi = chi,psi = psi,func = "1/x")
}
Msig[l] = t(y1)%*%Dui%*%y1 - 2*vp*t(y1)%*%v1s + ((tp^4)/4)*t(ui)%*%v1s -
2*t(y1)%*%Dui%*%nlmodel(x1,beta,bmetro[,l],cov1) +
2*vp*t(v1s)%*%nlmodel(x1,beta,bmetro[,l],cov1) +
t(nlmodel(x1,beta,bmetro[,l],cov1))%*%Dui%*%nlmodel(x1,beta,bmetro[,l],cov1)
Mbb[,,l] = bmetro[,l]%*%t(bmetro[,l])
GG1 = -(1/(sigmae*tp^2))*t(grad)%*%(2*vp*v1s + 2*Dui%*%nlmodel(x1,beta,bmetro[,l],cov1) -2*Dui%*%y1)
GG2 = -(3/2)*(nj[j])*(1/sigmae) + (1/(2*(tp^2)*(sigmae^2)))*Msig[l]
GG3 = (1/2)*MElim(q)%*%(kronecker(X = solve(D),Y = solve(D)))%*%as.vector(Mbb[,,l]-D)
VGi[,l] = rbind(GG1,GG2,GG3)
Hbeta[,,l] = -(1/(sigmae*tp^2))*t(grad)%*%Dui%*%grad
}
E_bi = apply(bmetro,1,mean)
E_bb = apply(Mbb,c(1,2),mean)
E_sig = mean(Msig)
E_VGi = apply(VGi,1,mean)
E_Hbeta = apply(Hbeta,c(1,2),mean)
SAEM_bb[count+1,j,,] = SAEM_bb[count,j,,] + seqq[count]*(E_bb - SAEM_bb[count,j,,])
SAEM_sig[count+1,j] = SAEM_sig[count,j] + seqq[count]*(E_sig - SAEM_sig[count,j])
SAEM_bi[count+1,j,] = SAEM_bi[count,j,] + seqq[count]*(E_bi - SAEM_bi[count,j,])
SAEM_VGi[count+1,j,] = SAEM_VGi[count,j,] + seqq[count]*(E_VGi - SAEM_VGi[count,j,])
SAEM_Hbeta[count+1,j,,] = SAEM_Hbeta[count,j,,] + seqq[count]*(E_Hbeta - SAEM_Hbeta[count,j,,])
##########################################################################
#PASSO M
##########################################################################
#PASSO M sigmae
sumsig = sumsig + SAEM_sig[count+1,j]
#PASSO M matriz D
sumD = sumD + SAEM_bb[count+1,j,,]
#PASSO M betas
sumb1 = sumb1 + SAEM_VGi[count+1,j,1:d]#Gradiente
sumb2 = sumb2 + SAEM_Hbeta[count+1,j,,]#Hessiana
#SUM do prod vector gradiente
IE = IE + SAEM_VGi[count+1,j,]%*%t(SAEM_VGi[count+1,j,])
}
ssumb2 = tryCatch(solve(-sumb2), error=function(e) stop("Error trying to compute an inverse of a matrix involved in beta. Frequently, for Nonlinear models this happens when the initial values are not properly."))
beta = beta + ssumb2%*%sumb1
D = sumD/n
sigmae = sumsig/(3*N*tp^2)
EP = sqrt(diag(tryCatch(solve(IE), error=function(e) diag(npar))))
EPV[,count] = EP
param = teta
teta = c(beta,sigmae,D[upper.tri(D, diag = T)])
criterio = abs(teta-param)/(abs(param)+delta1)
print(teta)
print(criterio)
if(max(criterio) < delta2){critval=critval+1}else{critval=0}
#PRUEBAS
#############################################################################
tetam[,count] = teta
#setTkProgressBar(pb, count, label=paste("Iter ",count,"/",MaxIter," - ",round(count/MaxIter*100,0),"% done",sep = ""))
pb$tick(tokens = list(what = paste0("Quantile ",p,": Working...")))
if (count == MaxIter){critval=10}
}
#loglik = logveroMC(beta,sigmae,D,y,x,covar,nj,MC=1000,n=n,d=d,q=q,p=p,nlmodel=nlmodel)
loglik = logveroIS(beta,sigmae,D,y,x,covar,nj,bi=SAEM_bi[count+1,,],bibi=SAEM_bb[count+1,,,],MIS=500,n=n,d=d,q=q,p=p,n.covar,nlmodel=nlmodel)
AIC = -2*loglik +2*npar
BIC = -2*loglik +log(N)*npar
HQ = -2*loglik +2*log(log(N))*npar
table = data.frame(beta,EP[1:d],beta/EP[1:d],2*pnorm(abs(beta/EP[1:d]),lower.tail = F))
rownames(table) = paste("beta",1:d)
colnames(table) = c("Estimate","Std. Error","z value","Pr(>|z|)")
end.time <- Sys.time()
time.taken <- end.time - start.time
res = list(quantile = p,nlmodel = nlmodel,iter = count,
criterio = max(criterio),beta = beta,
weights=SAEM_bi[count+1,,],
sigmae= sigmae,D = D,EP=EP,
table = table,
loglik=loglik,AIC=AIC,BIC=BIC,HQ=HQ,
time = time.taken)
conv = list(teta = tetam[,1:count],EPV = EPV[,1:count])
obj.out = list(conv=conv,res = res)
if (count == MaxIter){
#setTkProgressBar(pb, MaxIter, label=paste("MaxIter reached ",count,"/",MaxIter," - 100 % done",sep = ""))
#Sys.sleep(2)
#close(pb)
pb$tick(MaxIter,tokens = list(what =
paste0("Quantile ",p,": MaxIter reached!")))
Sys.sleep(1/20)
cat("\n")
#cat("\n")
}else{
#setTkProgressBar(pb, MaxIter, label=paste("Convergence at Iter ",count,"/",MaxIter," - 100 % done",sep = ""))
#Sys.sleep(2)
#close(pb)
pb$tick(MaxIter,tokens = list(what =
paste0("Quantile ",p,": Convergence reached!")))
Sys.sleep(1/20)
cat("\n")
#cat("\n")
}
class(obj.out) = "QSAEM_NL"
return(obj.out)
}
####################################################################
####################################################################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.