Nothing
hessianLogistic<-function(beta,y,w,X,variant){
n<-NROW(w)
J<-NCOL(w)
K<-NCOL(X)
y<-sapply(1:J,function(j) y)
eta<-sapply(1:J,function(j){
Xj<-cbind(X[,,j])
betaj<-beta[,j,drop=FALSE]
Xj%*%betaj
})
p<-1/(1+exp(-eta))
h<-w*exp(y*eta)/(1+exp(eta))
g<-rowSums(h)
dlogP<-function(k){
if (variant[k]) colSums(dh(k)/g) else sum(rowSums(dh(k))/g)
}
dh<-function(k) h*(y-p)*X[,k,]
S<-NULL
for (k in 1:K){
Sk<-dlogP(k)
names(Sk)<-rep(paste("var",k,sep=""),length(Sk))
S<-c(S,Sk)
}
d2h<-function(k,kprime,j,jprime) X[,k,j]*dh(kprime)[,jprime]*(y-p)[,j]-X[,k,j]*X[,kprime,jprime]*h[,j]*(p-p^2)[,jprime]
d2logP<-function(k,kprime,variantk,variantkprime){
if (variantk & variantkprime){
res<-matrix(0,J,J)
for (j in 1:J){
for (jprime in 1:J){
dhk<-dh(k)[,j]
dhkprime<-dh(kprime)[,jprime]
d2hkkprime<-d2h(k,kprime,j,jprime)
if (j!=jprime) res[j,jprime]<- sum(-dhk*dhkprime/g^2)
if (j==jprime) res[j,jprime]<- sum((d2hkkprime*g-dhk*dhkprime)/g^2)
}
}
}
if (variantk & !variantkprime){
res<-matrix(NA,J,1)
for (j in 1:J){
dhk<-dh(k)[,j]
dhkprime<-rowSums(dh(kprime))
d2hkkprime<-d2h(k,kprime,j,j)
res[j,1]<-sum((d2hkkprime*g-dhk*dhkprime)/g^2)
}
res<-cbind(res)
}
if (!variantk & variantkprime){
res<-matrix(NA,1,J)
for (j in 1:J){
dhk<-rowSums(dh(k))
dhkprime<-dh(kprime)[,j]
d2hkkprime<-d2h(k,kprime,j,j)
res[1,j]<-sum((d2hkkprime*g-dhk*dhkprime)/g^2)
}
res<-rbind(res)
}
if (!variantk & !variantkprime){
dhk<-rowSums(dh(k))
dhkprime<-rowSums(dh(kprime))
d2hkkprime<-rowSums(d2h(k,kprime,1:J,1:J))
res<-cbind(sum((d2hkkprime*g-dhk*dhkprime)/g^2))
}
res
}
H<-NULL
for (k in 1:K){
Hfilak<-NULL
for (kprime in 1:K){
res<-d2logP(k,kprime,variant[k],variant[kprime])
colnames(res)<-rep(paste("var",kprime,sep=""),NCOL(res))
rownames(res)<-rep(paste("var",k,sep=""),NROW(res))
Hfilak<-cbind(Hfilak,res)
}
H<-rbind(H,Hfilak)
}
return(list(S=S,H=H))
}
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.