Nothing
FI.MSTEC <-
function(P,y,X,nu)
{
if (missing(X)) {
X <- array(1, c(ncol(y), 1, nrow(y)))
}
y <- as.matrix(y)
if (!is.matrix(y))
stop("y must have at least one element")
if (is.null(X)) {
X <- array(c(diag(ncol(y))), c(ncol(y), ncol(y), nrow(y)))
}
if (is.array(X) == FALSE & is.list(X) == FALSE)
stop("X must be an array or a list")
if (is.array(X)) {
Xs <- list()
for (i in 1:nrow(y)) {
Xs[[i]] <- matrix(t(X[, , i]), nrow = ncol(y))
}
X <- Xs
}
if (ncol(y) != nrow(X[[1]]))
stop("y does not have the same number of columns than X")
if (nrow(y) != length(X))
stop("y does not have the same number of observations than X")
if (nu < 0)
stop("nu must be positive")
FI.MSMSNC<- function(P,y,X,nu,dist=c("STEC","SSLEC","SCEC"))
{
dtgamma1<-function(x,shape,rate,b=1,log=TRUE)
{
ll<-dgamma(x,shape=shape, rate=rate,log=TRUE)-pgamma(b,shape=shape,rate=rate,log.p=TRUE)
if(log==FALSE)
{ll<-exp(ll)}
ll
}
Iphi.MSMSNC<-function(a,d,p,w,dist=c("STEC","SSLEC","SCEC"),nu=1)
{
aux.Iphi.MSMSNC<-function(x,a,d,p,w,dist=c("STEC","SSLEC","SCEC"),nu=1){
if(dist=="STEC" || dist=="SSLEC")
logf<-switch(dist,
STEC=dgamma(x,shape=(nu+p)/2,rate=(nu+d)/2),
SSLEC=dtgamma1(x,shape=(2*nu+p)/2,rate=d/2,log=FALSE))
log.int<-x^w*exp(-0.5*x*d)*dcauchy(sqrt(x)*a)*logf
return(log.int)}
if(dist=="STEC" | dist=="SSLEC"){
max.int<-switch(dist,STEC=Inf,SSLEC=1)}
integrate(aux.Iphi.MSMSNC,lower=0,upper=max.int,a,d,p,w,dist=dist,nu=nu,abs.tol=1e-12)$value
}
IPhi.MSMSNC<-function(a,d,p,w,dist=c("STEC","SSLEC","SCEC"),nu=1)
{
aux.IPhi.MSMSNC<-function(x,a,d,p,w,dist=c("STEC","SSLEC","SCEC"),nu=1){
if(dist=="STEC" || dist=="SSLEC")
logf<-switch(dist,
STEC=dgamma(x,shape=(nu+p)/2,rate=(nu+d)/2),
SSLEC=dtgamma1(x,shape=(2*nu+p)/2,rate=d/2,log=FALSE))
log.int<-x^w*exp(-0.5*x*d)*pcauchy(sqrt(x)*a)*logf
return(log.int)}
if(dist=="STEC" | dist=="SSLEC") max.int<-switch(dist,STEC=Inf,SSLEC=1)
integrate(aux.IPhi.MSMSNC,lower=0,upper=max.int,a,d,p,w,dist=dist,nu=nu,abs.tol=1e-12)$value
}
if(is.array(X))
{Xs<-list()
if(ncol(y)>1 | !is.matrix(X)){
for (i in 1:nrow(y)){
Xs[[i]]<- matrix(t(X[,,i]),nrow=ncol(y))}}
if(ncol(y)==1 & is.matrix(X)){
for (i in 1:nrow(y)){
Xs[[i]]<- matrix(t(X[i,]),nrow=1)}}
X<-Xs}
# theta=list(beta=beta.new,Sigma=Sigma.new,eta=eta.new,nu=nu,dif=dif,iter=i)
p=ncol(y)
n=nrow(y)
m=ncol(X[[1]])
beta=P[1:m]
q0=length(beta)
Sigma<-matrix(0,ncol=p, nrow=p)
Sigma[lower.tri(Sigma, diag = TRUE)]<-P[1:(p*(p+1)/2)+m]
Sigma[upper.tri(Sigma, diag = TRUE)]<-P[1:(p*(p+1)/2)+m]
Sinv=solve2(Sigma) # Sigma^{-1}
C=matrix.sqrt(Sigma)
Cinv=solve2(C)# Sigma^{-1/2}
eta=P[m+(p*(p+1)/2)+1:p]
lambda=C%*%eta
#Delta=Sminv%*%lambda # eta
q1=p*(p+1)/2
pth=q0+p+q1
dlogS=matrix(0,pth,1)
# dlog Sigma
for (r in 1:q1){
indr=rep(0,q1)
indr[r]=1
Cr=xpnd(indr)
auxi=Sinv%*%(Cr%*%C+C%*%Cr)
dlogS[q0+r]=sum(diag(auxi))
}
dAi=matrix(0,pth,1)
ddi=matrix(0,pth,1)
MI=matrix(0,pth,pth)
for (i in 1:n){
resi=y[i,]-X[[i]]%*%beta # p x 1
ai= as.numeric(t(eta)%*%resi)
di= as.numeric(t(resi)%*%Sinv%*%resi)
if(dist!="SCEC")
{
Iphi.MSMSNC.i<-Iphi.MSMSNC(a=ai,d=di,p=p,w=(p+1)/2,dist=dist,nu=nu)
IPhi.MSMSNC.i<-IPhi.MSMSNC(a=ai,d=di,p=p,w=(p+2)/2,dist=dist,nu=nu)
Ki<-IPhi.MSMSNC(a=ai,d=di,p=p,w=p/2,dist=dist,nu=nu)
}
else
{ #print(nu)
nu1=nu[1]
nu2=nu[2]
Iphi.MSMSNC.i<-nu1*(nu2^((p+1)/2))*exp(-0.5*nu2^0.5*di)*dcauchy(nu2^0.25*ai)+(1-nu1)*exp(-0.5*di)*dcauchy(ai)
IPhi.MSMSNC.i<-nu1*(nu2^((p+2)/2))*exp(-0.5*nu2^0.5*di)*pcauchy(nu2^0.25*ai)+(1-nu1)*exp(-0.5*di)*pcauchy(ai)
Ki<-nu1*(nu2^(p/2))*exp(-0.5*nu2^0.5*di)*pcauchy(nu2^0.25*ai)+(1-nu1)*exp(-0.5*di)*pcauchy(ai)
}
#
dAi[1:q0]=t(X[[i]])%*%Cinv%*%lambda
dAi[(q0+q1+1):pth]=Cinv%*%resi
ddi[1:q0]=-2*t(X[[i]])%*%Cinv%*%Cinv%*%resi
for (r in 1:q1){
indr=rep(0,q1)
indr[r]=1
Cr=xpnd(indr)
ddi[q0+r]=-t(resi)%*%Cinv%*%(Cr%*%C+C%*%Cr)%*%Cinv%*%resi
dAi[q0+r]=-t(lambda)%*%Cinv%*%Cr%*%Cinv%*%resi
}
scoretheta=-0.5*dlogS +(Iphi.MSMSNC.i*dAi-0.5*IPhi.MSMSNC.i*ddi)/Ki
MI=MI+scoretheta%*%t(scoretheta)
}
return(-MI)
}
FI.MSMSNC(P,y,X,nu, dist="STEC")
}
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.