Nothing
marginBitsD <- function(params, eta1, eta2, X3, X4, dat, VC, precision, univariate=TRUE, bd=bd) {
sigma <- nu <- df2.sigma <- df2.nu <- dF2.sigma <- dF22.sigma <- dF2.nu <- dF22.nu <- d2f2sigma2 <- NULL
d2f2delta2sigma <- d2f2nu2 <- d2f2delta2nu <- d2f2nusigma <- d2f2sigmanu <- d2F2dsigma2 <- NULL
d2F22dsigma2 <- d2F2ddelta2dsigma <- d2F22ddelta2dsigma <- d2F2dnu2 <- d2F22dnu2 <- d2F2ddelta2dnu <- NULL
d2F22ddelta2dnu <- d2F2dsigmadnu <- d2F22dsigmadnu <- d2f2nusigma <- d2F2dnudsigma <- d2F22dnudsigma <- NULL
F1 <- F2 <- F22 <- dF1 <- d2F1ddelta1delta1 <- dF2 <- dF22 <- dF2.sigma <- dF22.sigma <- dF2.nu <- NULL
dF22.nu <- d2F2ddelta22 <- d2F22ddelta22 <- d2F2dsigma2 <- d2F22dsigma2 <- d2F2ddelta2dsigma <- NULL
d2F22ddelta2dsigma <- d2F2dnu2 <- d2F22dnu2 <- d2F2ddelta2dnu <- d2F22ddelta2dnu <- d2F2dnudsigma <- NULL
d2F22dnudsigma <- d2F2dsigmadnu <- d2F22dsigmadnu <- NULL
etasqv <- etanu <- NULL
eps <- sqrt(.Machine$double.eps)
i1 <- dat[,1]
i0 <- 1-i1
ind <- i1==0
if(univariate==FALSE) {
F1 <- pnorm(-eta1)
F1 <- ifelse(F1>0.00000001,F1,0.00000001)
F1 <- ifelse(F1<0.99999999,F1,0.99999999)
dF1 <- -dnorm(-eta1)
d2F1ddelta1delta1 <- -dF1*eta1
}
#fd.prec <- 10^(-7)
fd.prec <- eps
if(VC$margins[2]=="P") {
i2 <- dat[,2]
mu <- exp(eta2)
f2 <- dpois(i2, mu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- ppois(i2, mu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
# df2, dF2, dF22
gamma <- exp(-exp(eta2))*exp(eta2)^(i2)/(dpois(i2, exp(eta2)))
df2 <- as.vector(1/gamma*(-exp(-exp(eta2))*exp(eta2)^(i2)+exp(-exp(eta2))*i2*exp(eta2)^(i2-1))*exp(eta2))
if(univariate==FALSE) {
df2.int.P <- function(y, mu) {
gamma <- exp(-mu)*mu^(y)/(dpois(y, mu))
as.vector(1/gamma*(-exp(-mu)*mu^(y)+exp(-mu)*y*mu^(y-1))*mu)
}
ly <- length(i2)
cdf.f <- rep(0, ly)
nmu <- rep(exp(eta2), length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
allval <- seq(0, y.y)
pdfall <- df2.int.P(allval, mu = mm)
cdf.f[i] <- sum(pdfall)
}
dF2 <- cdf.f
dF22 <- dF2-df2
}
# Hessian derivative components
gamma.plus <- exp(-exp(eta2+fd.prec))*exp(eta2+fd.prec)^(i2)/(dpois(i2, exp(eta2+fd.prec)))
df2.plus <- as.vector(1/gamma.plus*(-exp(-exp(eta2+fd.prec))*exp(eta2+fd.prec)^(i2)+exp(-exp(eta2+fd.prec))*i2*exp(eta2+fd.prec)^(i2-1))*exp(eta2+fd.prec))
d2f2delta22 <- (df2.plus - df2)/fd.prec
if(univariate==FALSE) {
ly <- length(i2)
cdf.f.plus <- rep(0, ly)
nmu.plus <- rep(exp(eta2 + fd.prec), length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu.plus[i]
allval <- seq(0, y.y)
pdfall <- df2.int.P(allval, mu = mm)
cdf.f.plus[i] <- sum(pdfall)
}
dF2.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
}
} else if (VC$margins[2]=="NB") {
i2 <- dat[,2]
if(univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
}
sigma <- exp(sigma.star)
sigma <- ifelse(sigma>10^-8, sigma, 10^-8)
sigma <- ifelse(sigma<10^7, sigma, 10^7)
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
f2 <- dNBI(i2, sigma = sigma, mu = mu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pNBI(i2, sigma = sigma, mu = mu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
df2 <- f2*(i2*(mu*sigma)^(-1)*sigma*mu-(i2+1/sigma)*(mu*sigma+1)^(-1)*sigma*mu)
df2 <- as.vector(df2)
df2.sigma <- f2*(digamma(i2+1/sigma)*(-1/sigma^2)+i2*(mu*sigma)^(-1)*mu-(digamma(1/sigma)*(-1/sigma^2)+(-1/sigma^2)*log(mu*sigma+1)+(i2+1/sigma)*(1/(mu*sigma+1))*mu))
if(univariate==FALSE) {
df2.mu.int.NB<-function(y, mu, sigma, nu) {
f2 <- dNBI(y, sigma = sigma, mu = mu)
f2 <- ifelse(f2>precision, f2, precision)
df2 <- f2*(y*(mu*sigma)^(-1)*sigma*mu-(y+1/sigma)*(mu*sigma+1)^(-1)*sigma*mu)
df2
}
df2.sigma.int.NB<-function(y, mu, sigma, nu) {
f2 <- dNBI(y, sigma = sigma, mu = mu)
f2 <- ifelse(f2>precision, f2, precision)
df2 <- f2*(digamma(y+1/sigma)*(-1/sigma^2)+y*(mu*sigma)^(-1)*mu-(digamma(1/sigma)*(-1/sigma^2)+(-1/sigma^2)*log(mu*sigma+1)+(y+1/sigma)*(1/(mu*sigma+1))*mu))
df2
}
ly <- length(i2)
cdf.f <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
pdfall <- df2.mu.int.NB(allval, mm, ms)
cdf.f[i] <- sum(pdfall)
}
dF2 <- cdf.f
dF22 <- dF2 - df2
ly <- length(i2)
cdf.f.sigma <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
pdfall <- df2.sigma.int.NB(allval, mm, ms)
cdf.f.sigma[i] <- sum(pdfall)
}
dF2.sigma <- cdf.f.sigma
dF22.sigma <- dF2.sigma - df2.sigma
}
# Hessian derivative components
mu.plus <- ifelse(exp(eta2+fd.prec)>0.0001, exp(eta2+fd.prec), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
f2.plus <- dNBI(i2, sigma = sigma, mu = mu.plus)
f2.plus <- ifelse(f2.plus>precision, f2.plus, precision)
df2.plus <- f2.plus*(i2*(mu.plus*sigma)^(-1)*sigma*mu.plus-(i2+1/sigma)*(mu.plus*sigma+1)^(-1)*sigma*mu.plus)
df2.plus <- as.vector(df2.plus)
d2f2delta22 <- (df2.plus - df2) / fd.prec
sigma.plus <- sigma+fd.prec
f2.plus.sigma <- dNBI(i2, sigma = sigma.plus, mu = mu)
f2.plus.sigma <- ifelse(f2.plus.sigma>precision, f2.plus.sigma, precision)
df2.sigma.plus <- f2.plus.sigma*(digamma(i2+1/sigma.plus)*(-1/sigma.plus^2)+i2*(mu*sigma.plus)^(-1)*mu-(digamma(1/sigma.plus)*(-1/sigma.plus^2)+(-1/sigma.plus^2)*log(mu*sigma.plus+1)+(i2+1/sigma.plus)*(1/(mu*sigma.plus+1))*mu))
d2f2sigma2 <- (df2.sigma.plus - df2.sigma) / fd.prec
df2.dsigma.plus <- f2.plus.sigma*(i2*(mu*sigma.plus)^(-1)*sigma.plus*mu-(i2+1/sigma.plus)*(mu*sigma.plus+1)^(-1)*sigma.plus*mu)
df2.dsigma.plus <- as.vector(df2.dsigma.plus)
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec
if(univariate==FALSE) {
ly <- length(i2)
cdf.f.plus <- rep(0, ly)
nmu <- rep(mu.plus, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
pdfall <- df2.mu.int.NB(allval, mm, ms)
cdf.f.plus[i] <- sum(pdfall)
}
dF2.mu.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.mu.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
ly <- length(i2)
cdf.f.sigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
pdfall <- df2.sigma.int.NB(allval, mm, ms)
cdf.f.sigma.plus[i] <- sum(pdfall)
}
dF2.sigma.plus <- cdf.f.sigma.plus
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
ly <- length(i2)
cdf.f.dsigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
pdfall <- df2.mu.int.NB(allval, mm, ms)
cdf.f.dsigma.plus[i] <- sum(pdfall)
}
dF2.dsigma.plus <- cdf.f.dsigma.plus
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
}
} else if (VC$margins[2]=="D") {
i2 <- dat[,2]
if(univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
if (is.null(VC$X3)) nu.star <- params[(VC$X2.d2 + 2)]
if (!is.null(VC$X3)) nu.star <- etanu <- X4 %*% params[(VC$X2.d2 + VC$X3.d2 + 1):(VC$X2.d2 + VC$X3.d2 + VC$X4.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
if(is.null(VC$X3)) nu.star <- params[(VC$X1.d2+VC$X2.d2+2)]
if(!is.null(VC$X3)) nu.star <- etanu <- X4%*%params[(VC$X1.d2+VC$X2.d2+VC$X3.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2)]
}
nu <- plogis(nu.star)
nu <- ifelse(nu<precision, precision, nu)
nu <- ifelse(nu>(1-precision), 1-precision, nu)
sigma <- exp(sigma.star)
sigma <- ifelse(sigma>precision, sigma, precision)
sigma <- ifelse(sigma<10^7, sigma, 10^7)
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
f2 <- dDEL(i2, mu=mu, sigma=sigma, nu=nu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pDEL(i2, mu=mu, sigma=sigma, nu=nu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
bigS <- getFromNamespace("tofyDEL2", "gamlss.dist")
S_prime_mu <- (exp(bigS(i2, mu+fd.prec, sigma, nu))-exp(bigS(i2, mu, sigma, nu)))/fd.prec
S_prime_eta2 <- mu*S_prime_mu
S <- exp(bigS(i2, mu, sigma, nu))
df2 <- as.vector(f2*(-mu*nu+(-1/sigma)*(1+mu*sigma*(1-nu))^(-1)*mu*sigma*(1-nu)+S_prime_eta2/S))
S_prime_sigma <- (exp(bigS(i2, mu, sigma+fd.prec, nu))-exp(bigS(i2, mu, sigma, nu)))/fd.prec
d_prime <- ((1/sigma^2)*log(1+mu*sigma*(1-nu))-(1/sigma)*(1/(1+mu*sigma*(1-nu)))*mu*(1-nu))*(1+mu*sigma*(1-nu))^(-1/sigma)
df2.sigma <- f2*(d_prime*(1+mu*sigma*(1-nu))^(1/sigma)+S_prime_sigma/S)
S_prime_nu <- (exp(bigS(i2, mu, sigma, nu+fd.prec))-exp(bigS(i2, mu, sigma, nu)))/fd.prec
df2.nu <- f2*(-mu+((-1/sigma)*(1+mu*sigma*(1-nu))^(-1)*(-mu*sigma)+S_prime_nu/S))
if(univariate==FALSE) {
df2.mu.int.D <- function(y, mu, sigma, nu) {
f2 <- dDEL(y, mu=mu, sigma=sigma, nu=nu)
f2 <- ifelse(f2>precision, f2, precision)
S_prime_mu <- (exp(bigS(y, mu+fd.prec, sigma, nu))-exp(bigS(y, mu, sigma, nu)))/fd.prec
S_prime_eta2 <- mu*S_prime_mu
S <- exp(bigS(y, mu, sigma, nu))
df2 <- f2*(-mu*nu+(-1/sigma)*(1+mu*sigma*(1-nu))^(-1)*mu*sigma*(1-nu)+S_prime_eta2/S)
df2
}
df2.sigma.int.D <-function(y, mu, sigma, nu) {
f2 <- dDEL(y, mu=mu, sigma=sigma, nu=nu)
f2 <- ifelse(f2>precision, f2, precision)
S_prime_sigma <- (exp(bigS(y, mu, sigma+fd.prec, nu))-exp(bigS(y, mu, sigma, nu)))/fd.prec
d_prime <- ((1/sigma^2)*log(1+mu*sigma*(1-nu))-(1/sigma)*(1/(1+mu*sigma*(1-nu)))*mu*(1-nu))*(1+mu*sigma*(1-nu))^(-1/sigma)
S <- exp(bigS(y, mu, sigma, nu))
df2 <- f2*(d_prime*(1+mu*sigma*(1-nu))^(1/sigma)+S_prime_sigma/S)
df2
}
df2.nu.int.D <- function(y, mu, sigma, nu) {
f2 <- dDEL(y, mu=mu, sigma=sigma, nu=nu)
f2 <- ifelse(f2>precision, f2, precision)
S_prime_nu <- (exp(bigS(y, mu, sigma, nu+fd.prec))-exp(bigS(y, mu, sigma, nu)))/fd.prec
S <- exp(bigS(y, mu, sigma, nu))
df2 <- f2*(-mu+((-1/sigma)*(1+mu*sigma*(1-nu))^(-1)*(-mu*sigma)+S_prime_nu/S))
df2
}
ly <- length(i2)
cdf.f <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
pdfall <- df2.mu.int.D(allval, mm, ms, mn)
cdf.f[i] <- sum(pdfall)
}
dF2<-cdf.f
dF22<-dF2-df2
ly <- length(i2)
cdf.f.sigma <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
pdfall <- df2.sigma.int.D(allval, mm, ms, mn)
cdf.f.sigma[i] <- sum(pdfall)
}
dF2.sigma <- cdf.f.sigma
dF22.sigma <- dF2.sigma-df2.sigma
ly <- length(i2)
cdf.f.nu <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn <- nnu[i]
allval <- seq(0, y.y)
pdfall <- df2.nu.int.D(allval, mm, ms, mn)
cdf.f.nu[i] <- sum(pdfall)
}
dF2.nu <- cdf.f.nu
dF22.nu <- dF2.nu-df2.nu
}
# Hessian derivative components
#fd.prec2 <- 10^-3
fd.prec2 <- sqrt(eps)
#delta2delta2
mu.plus <- ifelse(exp(eta2+fd.prec2)>0.0001, exp(eta2+fd.prec2), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
f2.plus <- dDEL(i2, sigma = sigma, mu = mu.plus, nu=nu)
f2.plus <- ifelse(f2.plus>precision, f2.plus, precision)
S_prime_mu.plus <- (exp(bigS(i2, mu.plus+fd.prec, sigma, nu))-exp(bigS(i2, mu.plus, sigma, nu)))/fd.prec
S_prime_eta2.plus <- mu.plus*S_prime_mu.plus
S.plus <- exp(bigS(i2, mu.plus, sigma, nu))
df2.plus <- as.vector(f2.plus*(-mu.plus*nu+(-1/sigma)*(1+mu.plus*sigma*(1-nu))^(-1)*mu.plus*sigma*(1-nu)+S_prime_eta2.plus/S.plus))
df2.plus <- as.vector(df2.plus)
d2f2delta22 <- (df2.plus - df2) / fd.prec2
#sigma^2
sigma.plus <- sigma + fd.prec2
f2.plus.sigma <- dDEL(i2, sigma = sigma.plus, mu = mu, nu = nu)
f2.plus.sigma <- ifelse(f2.plus.sigma>precision, f2.plus.sigma, precision)
S_prime_sigma <- (exp(bigS(i2, mu, sigma.plus+fd.prec, nu))-exp(bigS(i2, mu, sigma.plus, nu)))/fd.prec
d_prime <- ((1/sigma.plus^2)*log(1+mu*sigma.plus*(1-nu))-(1/sigma.plus)*(1/(1+mu*sigma.plus*(1-nu)))*mu*(1-nu))*(1+mu*sigma.plus*(1-nu))^(-1/sigma.plus)
S.plus.sigma <- exp(bigS(i2, mu, sigma.plus, nu))
df2.sigma.plus <- f2.plus.sigma*(d_prime*(1+mu*sigma.plus*(1-nu))^(1/sigma.plus)+S_prime_sigma/S.plus.sigma)
d2f2sigma2 <- (df2.sigma.plus - df2.sigma) / fd.prec2
#delta2sigma
S_prime_mu.dsigma <- (exp(bigS(i2, mu+fd.prec, sigma.plus, nu))-exp(bigS(i2, mu, sigma.plus, nu)))/fd.prec
S_prime_eta2.dsigma <- mu*S_prime_mu.dsigma
df2.dsigma.plus <- as.vector(f2.plus.sigma*(-mu*nu+(-1/sigma.plus)*(1+mu*sigma.plus*(1-nu))^(-1)*mu*sigma.plus*(1-nu)+S_prime_eta2.dsigma/S.plus.sigma))
df2.dsigma.plus <- as.vector(df2.dsigma.plus)
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec2
#nu^2
nu.plus <- nu + fd.prec2
nu.plus <- ifelse(nu.plus>=(1-precision), (1-precision), nu.plus)
f2.plus.nu <- dDEL(i2, sigma = sigma, mu = mu, nu = nu.plus)
f2.plus.nu <- ifelse(f2.plus.nu>precision, f2.plus.nu, precision)
S_prime_nu <- (exp(bigS(i2, mu, sigma, nu.plus+fd.prec))-exp(bigS(i2, mu, sigma, nu.plus)))/fd.prec
S.plus.nu <- exp(bigS(i2, mu, sigma, nu.plus))
df2.nu.plus <- f2.plus.nu*(-mu+((-1/sigma)*(1+mu*sigma*(1-nu.plus))^(-1)*(-mu*sigma)+S_prime_nu/S.plus.nu))
d2f2nu2 <- (df2.nu.plus - df2.nu) / fd.prec2
#delta2nu
S_prime_mu.dnu <- (exp(bigS(i2, mu+fd.prec, sigma, nu.plus))-exp(bigS(i2, mu, sigma, nu.plus)))/fd.prec
S_prime_eta2.dnu <- mu*S_prime_mu.dnu
S.dnu <- exp(bigS(i2, mu, sigma, nu.plus))
df2.dnu.plus <- as.vector(f2.plus.nu*(-mu*nu.plus+(-1/sigma)*(1+mu*sigma*(1-nu.plus))^(-1)*mu*sigma*(1-nu.plus)+S_prime_eta2.dnu/S.dnu))
df2.dnu.plus <- as.vector(df2.dnu.plus)
d2f2delta2nu <- (df2.dnu.plus - df2) / fd.prec2
#nusigma
f2.plus.nu.sigma <- dDEL(i2, sigma = sigma.plus, mu = mu, nu = nu)
f2.plus.nu.sigma <- ifelse(f2.plus.nu.sigma>precision, f2.plus.nu.sigma, precision)
S_prime_nu.plus <- (exp(bigS(i2, mu, sigma.plus, nu+fd.prec))-exp(bigS(i2, mu, sigma.plus, nu)))/fd.prec
df2.nu.sigma.plus <- f2.plus.nu.sigma*(-mu+((-1/sigma.plus)*(1+mu*sigma.plus*(1-nu))^(-1)*(-mu*sigma.plus)+S_prime_nu.plus/S.plus.sigma))
d2f2nusigma <- (df2.nu.sigma.plus - df2.nu) / fd.prec2
if(univariate==FALSE) {
ly <- length(i2)
cdf.f.plus <- rep(0, ly)
nmu <- rep(mu.plus, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
pdfall <- df2.mu.int.D(allval, mm, ms, mn)
cdf.f.plus[i] <- sum(pdfall)
}
dF2.mu.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.mu.plus - dF2) / fd.prec2
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
ly <- length(i2)
cdf.f.sigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
pdfall <- df2.sigma.int.D(allval, mm, ms, mn)
cdf.f.sigma.plus[i] <- sum(pdfall)
}
dF2.sigma.plus <- cdf.f.sigma.plus
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec2
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
ly <- length(i2)
cdf.f.dsigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
pdfall <- df2.mu.int.D(allval, mm, ms, mn)
cdf.f.dsigma.plus[i] <- sum(pdfall)
}
dF2.dsigma.plus <- cdf.f.dsigma.plus
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec2
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
ly <- length(i2)
cdf.f.nu.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
pdfall <- df2.nu.int.D(allval, mm, ms, mn)
cdf.f.nu.plus[i] <- sum(pdfall)
}
dF2.nu.plus <- cdf.f.nu.plus
d2F2dnu2 <- (dF2.nu.plus - dF2.nu) / fd.prec2
d2F22dnu2 <- d2F2dnu2 - d2f2nu2
ly <- length(i2)
cdf.f.dnu.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
pdfall <- df2.mu.int.D(allval, mm, ms, mn)
cdf.f.dnu.plus[i] <- sum(pdfall)
}
dF2.dnu.plus <- cdf.f.dnu.plus
d2F2ddelta2dnu <- (dF2.dnu.plus - dF2) / fd.prec2
d2F22ddelta2dnu <- d2F2ddelta2dnu - d2f2delta2nu
ly <- length(i2)
cdf.f.dnu.dsigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
pdfall <- df2.nu.int.D(allval, mm, ms, mn)
cdf.f.dnu.dsigma.plus[i] <- sum(pdfall)
}
dF2.dnu.dsigma.plus <- cdf.f.dnu.dsigma.plus
d2F2dnudsigma <- (dF2.dnu.dsigma.plus - dF2.nu) / fd.prec2
d2F22dnudsigma <- d2F2dnudsigma - d2f2nusigma
}
} else if (VC$margins[2]=="PIG") {
i2 <- dat[,2]
fd.prec <- 10^(-7)
fd.prec2 <- 10^-3
if (univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
}
sigma <- exp(sigma.star)
sigma <- ifelse(sigma>precision, sigma, precision)
sigma <- ifelse(sigma<10^7, sigma, 10^7)
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
f2 <- dPIG(i2, mu=mu, sigma=sigma)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pPIG(i2, mu=mu, sigma=sigma)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
f2.fd.mu <- dPIG(i2, mu = (mu+fd.prec), sigma = sigma)
f2.fd.mu <- ifelse(f2.fd.mu>precision, f2.fd.mu, precision)
df2 <- (f2.fd.mu - f2) / (fd.prec)
df2 <- as.vector(df2)*as.vector(mu)
f2.fd.sigma <- dPIG(i2, mu=mu, sigma=(sigma+fd.prec))
f2.fd.sigma <- ifelse(f2.fd.sigma>precision, f2.fd.sigma, precision)
df2.sigma <- (f2.fd.sigma - f2) / (fd.prec)
if(univariate==FALSE) {
F2.fd.mu <- pPIG(i2, mu = (mu+fd.prec), sigma = sigma)
F2.fd.mu <- ifelse(F2.fd.mu<(1-precision), F2.fd.mu, 1-precision)
dF2 <- (F2.fd.mu-F2) / (fd.prec)
dF2 <- as.vector(dF2)*as.vector(mu)
dF22 <- dF2 - df2
F2.fd.sigma <- pPIG(i2, mu = mu, sigma = (sigma+fd.prec))
F2.fd.sigma <- ifelse(F2.fd.sigma<(1-precision), F2.fd.sigma, 1-precision)
dF2.sigma <- (F2.fd.sigma - F2) / (fd.prec)
dF22.sigma <- dF2.sigma - df2.sigma
}
# Hessian derivative components
#fd.prec2 <- 10^-3
#fd.prec2 <- sqrt(eps)
mu.plus <- ifelse(exp(eta2+fd.prec2)>0.0001, exp(eta2+fd.prec2), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
f2.plus <- dPIG(i2, mu = mu.plus, sigma = sigma)
f2.plus <- ifelse(f2.plus>precision, f2.plus, precision)
f2.fd.mu.plus <- dPIG(i2, mu = (mu.plus+fd.prec), sigma = sigma)
f2.fd.mu.plus <- ifelse(f2.fd.mu.plus>precision, f2.fd.mu.plus, precision)
df2.plus <- (f2.fd.mu.plus - f2.plus) / (fd.prec)
df2.plus <- as.vector(df2.plus)*as.vector(mu.plus)
df2.plus <- as.vector(df2.plus)
d2f2delta22 <- (df2.plus - df2) / fd.prec2
sigma.plus <- sigma + fd.prec2
f2.plus.sigma <- dPIG(i2, mu = mu, sigma = sigma.plus)
f2.plus.sigma <- ifelse(f2.plus.sigma>precision, f2.plus.sigma, precision)
f2.fd.sigma.plus <- dPIG(i2, mu = mu, sigma = (sigma.plus+fd.prec))
f2.fd.sigma.plus <- ifelse(f2.fd.sigma.plus>precision, f2.fd.sigma.plus, precision)
df2.sigma.plus <- (f2.fd.sigma.plus - f2.plus.sigma) / (fd.prec)
d2f2sigma2 <- (df2.sigma.plus - df2.sigma) / fd.prec2
f2.dsigma.plus <- dPIG(i2, mu = mu, sigma = sigma.plus)
f2.dsigma.plus <- ifelse(f2.dsigma.plus>precision, f2.dsigma.plus, precision)
f2.fd.mu.sigma.plus <- dPIG(i2, mu = (mu+fd.prec), sigma = sigma.plus)
f2.fd.mu.sigma.plus <- ifelse(f2.fd.mu.sigma.plus>precision, f2.fd.mu.sigma.plus, precision)
df2.dsigma.plus <- (f2.fd.mu.sigma.plus - f2.dsigma.plus) / (fd.prec)
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*as.vector(mu)
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec2
if(univariate==FALSE) {
F2.plus <- pPIG(i2, mu = mu.plus, sigma = sigma)
F2.fd.mu.plus <- pPIG(i2, mu = (mu.plus+fd.prec), sigma = sigma)
F2.fd.mu.plus <- ifelse(F2.fd.mu.plus<(1-precision), F2.fd.mu.plus, 1-precision)
dF2.plus <- (F2.fd.mu.plus - F2.plus) / (fd.prec)
dF2.plus <- as.vector(dF2.plus)*as.vector(mu.plus)
d2F2ddelta22 <- (dF2.plus - dF2) / fd.prec2
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
F2.sigma.plus <- pPIG(i2, mu = mu, sigma = sigma.plus)
F2.fd.sigma.plus <- pPIG(i2, mu = mu, sigma = (sigma.plus+fd.prec))
F2.fd.sigma.plus <- ifelse(F2.fd.sigma.plus<(1-precision), F2.fd.sigma.plus, 1-precision)
dF2.sigma.plus <- (F2.fd.sigma.plus - F2.sigma.plus) / (fd.prec)
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) /fd.prec2
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
F2.fd.dsigma.plus <- pPIG(i2, mu = (mu+fd.prec), sigma = sigma.plus)
F2.fd.dsigma.plus <- ifelse(F2.fd.dsigma.plus<(1-precision), F2.fd.dsigma.plus, 1-precision)
dF2.dsigma.plus <- (F2.fd.dsigma.plus - F2.sigma.plus) / (fd.prec)
dF2.dsigma.plus <- as.vector(dF2.dsigma.plus)*as.vector(mu)
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec2
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
}
} else if (VC$margins[2]=="S") {
i2 <- dat[,2]
fd.prec <- 10^(-7)
fd.prec2 <- 10^-3
if(univariate==TRUE){
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
if (is.null(VC$X3)) nu <- nu.star <- params[(VC$X2.d2 + 2)]
if (!is.null(VC$X3)) nu <- nu.star <- etanu <- X4 %*% params[(VC$X2.d2 + VC$X3.d2 + 1):(VC$X2.d2 + VC$X3.d2 + VC$X4.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
if(is.null(VC$X3)) nu <- nu.star <- params[(VC$X1.d2+VC$X2.d2+2)]
if(!is.null(VC$X3)) nu <- nu.star <- etanu <- X4%*%params[(VC$X1.d2+VC$X2.d2+VC$X3.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2)]
}
sigma <- exp(sigma.star)
sigma <- ifelse(sigma>precision, sigma, precision)
sigma <- ifelse(sigma<10^7, sigma, 10^7)
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
nu <- ifelse(nu<150, nu, 150)
nu <- ifelse(nu>-150, nu, -150)
f2 <- dSICHEL(i2, mu=mu, sigma=sigma, nu=nu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pSICHEL(i2, mu=mu, sigma=sigma, nu=nu)
F2 <- ifelse(F2>precision, F2, precision)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
f2.fd.mu <- dSICHEL(i2, mu = (mu +fd.prec), sigma = sigma, nu = nu)
f2.fd.mu <- ifelse(f2.fd.mu>precision, f2.fd.mu, precision)
df2 <- (f2.fd.mu - f2) / (fd.prec)
df2 <- as.vector(df2)*as.vector(mu)
f2.fd.sigma <- dSICHEL(i2, mu = mu, sigma = (sigma+fd.prec), nu = nu)
f2.fd.sigma <- ifelse(f2.fd.sigma>precision, f2.fd.sigma, precision)
df2.sigma <- (f2.fd.sigma - f2) / (fd.prec)
f2.fd.nu <- dSICHEL(i2, mu = mu, sigma = sigma, nu =(nu+fd.prec))
f2.fd.nu <- ifelse(f2.fd.nu>precision, f2.fd.nu, precision)
df2.nu <- (f2.fd.nu - f2) / (fd.prec)
if(univariate==FALSE) {
F2.fd.mu <- pSICHEL(i2, mu = (mu+fd.prec), sigma = sigma, nu = nu)
F2.fd.mu <- ifelse(F2.fd.mu<(1-precision), F2.fd.mu, 1-precision)
dF2 <- (F2.fd.mu - F2) / (fd.prec)
dF2 <- as.vector(dF2)*as.vector(mu)
dF22 <- dF2 - df2
F2.fd.sigma <- pSICHEL(i2, mu = mu, sigma = (sigma+fd.prec), nu = nu)
F2.fd.sigma <- ifelse(F2.fd.sigma<(1-precision), F2.fd.sigma, 1-precision)
dF2.sigma <- (F2.fd.sigma - F2) / (fd.prec)
dF22.sigma <- dF2.sigma - df2.sigma
F2.fd.nu <- pSICHEL(i2, mu = mu, sigma = sigma, nu = (nu+fd.prec))
F2.fd.nu <- ifelse(F2.fd.nu<(1-precision), F2.fd.nu, 1-precision)
dF2.nu <- (F2.fd.nu - F2) / (fd.prec)
dF22.nu <- dF2.nu - df2.nu
}
# Hessian derivative components
#second derivative fd precision
#fd.prec2 <- 10^-3
#fd.prec2 <- sqrt(eps)
#pmfs
#delta2^2
mu.plus <- ifelse(exp(eta2+fd.prec2)>0.0001, exp(eta2+fd.prec2), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
f2.plus <- dSICHEL(i2, mu = mu.plus, sigma = sigma, nu = nu)
f2.plus <- ifelse(f2.plus>precision, f2.plus, precision)
f2.fd.mu.plus <- dSICHEL(i2, mu = (mu.plus+fd.prec), sigma = sigma, nu = nu)
f2.fd.mu.plus <- ifelse(f2.fd.mu.plus>precision, f2.fd.mu.plus, precision)
df2.plus <- (f2.fd.mu.plus - f2.plus) / (fd.prec)
df2.plus <- as.vector(df2.plus)*as.vector(mu.plus)
df2.plus <- as.vector(df2.plus)
d2f2delta22 <- (df2.plus - df2) / fd.prec2
#sigma^2
sigma.plus <- sigma+fd.prec2
f2.plus.sigma <- dSICHEL(i2, mu = mu, sigma = sigma.plus, nu = nu)
f2.plus.sigma <- ifelse(f2.plus.sigma>precision, f2.plus.sigma, precision)
f2.fd.sigma.plus <- dSICHEL(i2, mu = mu, sigma = (sigma.plus+fd.prec), nu = nu)
f2.fd.sigma.plus <- ifelse(f2.fd.sigma.plus>precision, f2.fd.sigma.plus, precision)
df2.sigma.plus <- (f2.fd.sigma.plus - f2.plus.sigma) / (fd.prec)
d2f2sigma2 <- (df2.sigma.plus - df2.sigma) / fd.prec2
#delta2sigma
f2.dsigma.plus <- dSICHEL(i2, mu = mu, sigma = sigma.plus, nu = nu)
f2.dsigma.plus <- ifelse(f2.dsigma.plus>precision, f2.dsigma.plus, precision)
f2.fd.mu.sigma.plus <- dSICHEL(i2, mu = (mu+fd.prec), sigma = sigma.plus, nu = nu)
f2.fd.mu.sigma.plus <- ifelse(f2.fd.mu.sigma.plus>precision, f2.fd.mu.sigma.plus, precision)
df2.dsigma.plus <- (f2.fd.mu.sigma.plus - f2.dsigma.plus) / (fd.prec)
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*as.vector(mu)
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec2
#nu^2
nu.plus <- nu+fd.prec2
f2.plus.nu <- dSICHEL(i2, mu = mu, sigma = sigma, nu = nu.plus)
f2.plus.nu <- ifelse(f2.plus.nu>precision, f2.plus.nu, precision)
f2.fd.nu.plus <- dSICHEL(i2, mu = mu, sigma = sigma, nu = nu.plus+fd.prec)
f2.fd.nu.plus <- ifelse(f2.fd.nu.plus>precision, f2.fd.nu.plus, precision)
df2.nu.plus <- (f2.fd.nu.plus - f2.plus.nu) / (fd.prec)
d2f2nu2 <- (df2.nu.plus - df2.nu) / fd.prec2
#delta2nu
f2.dnu.plus <- dSICHEL(i2, mu = mu, sigma = sigma, nu = nu.plus)
f2.dnu.plus <- ifelse(f2.dnu.plus>precision, f2.dnu.plus, precision)
f2.fd.mu.nu.plus <- dSICHEL(i2, mu = (mu+fd.prec), sigma = sigma, nu = nu.plus)
f2.fd.mu.nu.plus <- ifelse(f2.fd.mu.nu.plus>precision, f2.fd.mu.nu.plus, precision)
df2.dnu.plus <- (f2.fd.mu.nu.plus - f2.dnu.plus) / (fd.prec)
df2.dnu.plus <- as.vector(df2.dnu.plus)*as.vector(mu)
d2f2delta2nu <-(df2.dnu.plus - df2) / fd.prec2
#nusigma
f2.fd.sigma.nu.plus <- dSICHEL(i2, mu = mu, sigma = (sigma+fd.prec), nu = nu.plus)
f2.fd.sigma.nu.plus <- ifelse(f2.fd.sigma.nu.plus>precision, f2.fd.sigma.nu.plus, precision)
df2.dnu.sigma.plus <- (f2.fd.sigma.nu.plus - f2.plus.nu) / (fd.prec)
df2.dnu.sigma.plus <- as.vector(df2.dnu.sigma.plus)
d2f2sigmanu <- (df2.dnu.sigma.plus - df2.sigma) / fd.prec2
if(univariate==FALSE) {
#cdfs
#delta2^2
F2.plus <- pSICHEL(i2, mu = mu.plus, sigma = sigma, nu = nu)
F2.fd.mu.plus <- pSICHEL(i2, mu = (mu.plus+fd.prec), sigma = sigma, nu = nu)
F2.fd.mu.plus <- ifelse(F2.fd.mu.plus<(1-precision), F2.fd.mu.plus, 1-precision)
dF2.plus <- (F2.fd.mu.plus - F2.plus) / (fd.prec)
dF2.plus <- as.vector(dF2.plus)*as.vector(mu.plus)
d2F2ddelta22 <-(dF2.plus - dF2) / fd.prec2
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
#sigma^2
F2.sigma.plus <- pSICHEL(i2, mu = mu, sigma = sigma.plus, nu = nu)
F2.fd.sigma.plus <- pSICHEL(i2, mu = mu, sigma = (sigma.plus+fd.prec), nu=nu)
F2.fd.sigma.plus <- ifelse(F2.fd.sigma.plus<(1-precision), F2.fd.sigma.plus, 1-precision)
dF2.sigma.plus <- (F2.fd.sigma.plus - F2.sigma.plus) / (fd.prec)
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec2
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
#delta2sigma
F2.fd.dsigma.plus <- pSICHEL(i2, mu = (mu+fd.prec), sigma = sigma.plus, nu)
F2.fd.dsigma.plus <- ifelse(F2.fd.dsigma.plus<(1-precision), F2.fd.dsigma.plus, 1-precision)
dF2.dsigma.plus <- (F2.fd.dsigma.plus - F2.sigma.plus) / (fd.prec)
dF2.dsigma.plus <- as.vector(dF2.dsigma.plus)*as.vector(mu)
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec2
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
#nu^2
F2.nu.plus<- pSICHEL(i2, mu = mu, sigma = sigma, nu = nu.plus)
F2.fd.nu.plus <- pSICHEL(i2, mu = mu, sigma = sigma, nu = nu.plus+fd.prec)
F2.fd.nu.plus <- ifelse(F2.fd.nu.plus<(1-precision), F2.fd.nu.plus, 1-precision)
dF2.nu.plus <- (F2.fd.nu.plus - F2.nu.plus) / (fd.prec)
d2F2dnu2 <- (dF2.nu.plus - dF2.nu) / fd.prec2
d2F22dnu2 <- d2F2dnu2 - d2f2nu2
#delta2nu
F2.fd.dnu.plus <- pSICHEL(i2, mu = (mu+fd.prec), sigma = sigma, nu.plus)
F2.fd.dnu.plus <- ifelse(F2.fd.dnu.plus<(1-precision), F2.fd.dnu.plus, 1-precision)
dF2.dnu.plus <- (F2.fd.dnu.plus - F2.nu.plus) / (fd.prec)
dF2.dnu.plus <- as.vector(dF2.dnu.plus)*as.vector(mu)
d2F2ddelta2dnu <- (dF2.dnu.plus - dF2) / fd.prec2
d2F22ddelta2dnu <- d2F2ddelta2dnu - d2f2delta2nu
#sigmanu
F2.fd.dnu.sigma.plus <- pSICHEL(i2, mu = mu, sigma = (sigma+fd.prec), nu.plus)
F2.fd.dnu.sigma.plus <- ifelse(F2.fd.dnu.sigma.plus<(1-precision), F2.fd.dnu.sigma.plus, 1-precision)
dF2.dnu.sigma.plus <- (F2.fd.dnu.sigma.plus - F2.nu.plus) / (fd.prec)
dF2.dnu.sigma.plus <- as.vector(dF2.dnu.sigma.plus)
d2F2dsigmadnu <- (dF2.dnu.sigma.plus - dF2.sigma) / fd.prec2
d2F22dsigmadnu <- d2F2dsigmadnu - d2f2sigmanu
}
} else if(VC$margins[2]=="BI") {
i2 <- dat[,2]
mu <- plogis(eta2)
mu <- ifelse(mu<(1-precision), mu, 1-mu)
mu <- ifelse(mu>precision, mu, precision)
f2 <- dBI(i2, mu, bd = bd)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pBI(i2, mu, bd = bd)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
# df2, dF2, dF22
df2 <- -(((1 - mu)^(-1 + bd - i2)*mu^(-1 + i2)*(-i2 + bd*mu)*gamma(1 + bd))/(gamma(1 + bd - i2)*gamma(1 + i2)))*(1-mu)*mu
if(univariate==FALSE) {
df2.int.BI <- function(y, mu, n) {
-(((1 - mu)^(-1 + n - y)*mu^(-1 + y)*(-y + n*mu)*gamma(1 + n))/(gamma(1 + n - y)*gamma(1 + y)))*(1-mu)*mu
}
ly <- length(i2)
cdf.f <- rep(0, ly)
nmu <- rep(mu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
allval <- seq(0, y.y)
pdfall <- df2.int.BI(allval, mu = mm, n=bd)
cdf.f[i] <- sum(pdfall)
}
dF2 <- cdf.f
dF22 <- dF2-df2
}
# Hessian derivative components
mu.plus <- plogis(eta2+fd.prec)
df2.plus <- -(((1 - mu.plus)^(-1 + bd - i2)*mu.plus^(-1 + i2)*(-i2 + bd*mu.plus)*gamma(1 + bd))/(gamma(1 + bd - i2)*gamma(1 + i2)))*(1-mu.plus)*mu.plus
d2f2delta22 <- (df2.plus - df2)/fd.prec
if(univariate==FALSE) {
ly <- length(i2)
cdf.f.plus <- rep(0, ly)
nmu.plus <- rep(mu.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu.plus[i]
allval <- seq(0, y.y)
pdfall <- df2.int.BI(allval, mu = mm, n=bd)
cdf.f.plus[i] <- sum(pdfall)
}
dF2.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
}
} else if (VC$margins[2]=="BB") {
i2 <- dat[,2]
if(univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
}
sigma <- exp(sigma.star)
sigma <- ifelse(sigma>10^-8, sigma, 10^-8)
sigma <- ifelse(sigma<10^7, sigma, 10^7)
mu <- plogis(eta2)
mu <- ifelse(mu<(1-precision), mu, 1-mu)
mu <- ifelse(mu>precision, mu, precision)
f2 <- dBB(i2, sigma = sigma, mu = mu, bd=bd)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pBB(i2, sigma = sigma, mu = mu, bd=bd)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
df2 <- ( (gamma(1 + bd) * gamma(i2 + mu/sigma) * gamma(1/sigma) * gamma(
bd - (-1 + mu + i2 * sigma)/sigma) * (digamma(i2 + mu/sigma) +
digamma((1 - mu)/sigma) - digamma(mu/sigma) -
digamma(bd - (-1 + mu + i2 * sigma)/sigma)))/(sigma * gamma(
1 + bd - i2) * gamma(1 + i2) * gamma(bd + 1/sigma) * gamma((1 - mu)/
sigma) * gamma(mu/sigma)) )*(1-mu)*mu
df2 <- as.vector(df2)
df2.sigma <- (gamma(1 + bd) * gamma(i2 + mu/sigma) * gamma(1/sigma) * gamma(
bd - (-1 + mu + i2 * sigma)/sigma) * (digamma(bd + 1/sigma) -
mu * digamma(i2 + mu/sigma) - digamma(1/sigma) +
digamma((1 - mu)/sigma) - mu * digamma((1 - mu)/sigma) +
mu * digamma(mu/sigma) -
digamma(bd - (-1 + mu + i2 * sigma)/sigma) +
mu * digamma(bd - (-1 + mu + i2 * sigma)/sigma)))/(sigma^2 * gamma(
1 + bd - i2) * gamma(1 + i2) * gamma(bd + 1/sigma) * gamma((1 - mu)/
sigma) * gamma(mu/sigma))
if(univariate==FALSE) {
df2.mu.int.BB<-function(y, mu, sigma, nu, n=bd) {
df2<- ( (gamma(1 + n) * gamma(y + mu/sigma) * gamma(1/sigma) * gamma(
n - (-1 + mu + y * sigma)/sigma) * (digamma(y + mu/sigma) +
digamma((1 - mu)/sigma) - digamma(mu/sigma) -
digamma(n - (-1 + mu + y * sigma)/sigma)))/(sigma * gamma(
1 + n - y) * gamma(1 + y) * gamma(n + 1/sigma) * gamma((1 - mu)/
sigma) * gamma(mu/sigma)) )*(1-mu)*mu
df2
}
df2.sigma.int.BB<-function(y, mu, sigma, nu, n=bd) {
df2 <- (gamma(1 + n) * gamma(y + mu/sigma) * gamma(1/sigma) * gamma(
n - (-1 + mu + y * sigma)/sigma) * (digamma(n + 1/sigma) -
mu * digamma(y + mu/sigma) - digamma(1/sigma) +
digamma((1 - mu)/sigma) - mu * digamma((1 - mu)/sigma) +
mu * digamma(mu/sigma) -
digamma(n - (-1 + mu + y * sigma)/sigma) +
mu * digamma(n - (-1 + mu + y * sigma)/sigma)))/(sigma^2 * gamma(
1 + n - y) * gamma(1 + y) * gamma(n + 1/sigma) * gamma((1 - mu)/
sigma) * gamma(mu/sigma))
df2
}
ly <- length(i2)
cdf.f <- rep(0, ly)
cdf.f.sigma <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f[i] <- sum(df2.mu.int.BB(allval, mm, ms, bd))
cdf.f.sigma[i] <- sum(df2.sigma.int.BB(allval, mm, ms, bd))
}
dF2<-cdf.f
dF22<-dF2-df2
dF2.sigma <- cdf.f.sigma
dF22.sigma <- dF2.sigma-df2.sigma
}
# Hessian derivative components
mu.plus <- ifelse(plogis(eta2+fd.prec)>0.0001, plogis(eta2+fd.prec), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
df2.plus <-( (gamma(1 + bd) * gamma(i2 + mu.plus/sigma) * gamma(1/sigma) * gamma(
bd - (-1 + mu.plus + i2 * sigma)/sigma) * (digamma(i2 + mu.plus/sigma) +
digamma((1 - mu.plus)/sigma) - digamma(mu.plus/sigma) -
digamma(bd - (-1 + mu.plus + i2 * sigma)/sigma)))/(sigma * gamma(
1 + bd - i2) * gamma(1 + i2) * gamma(bd + 1/sigma) * gamma((1 - mu.plus)/
sigma) * gamma(mu.plus/sigma)) )*(1-mu.plus)*mu.plus
df2.plus <- as.vector(df2.plus)
d2f2delta22 <- (df2.plus - df2) / fd.prec
sigma.plus <- sigma+fd.prec
df2.sigma.plus <- (gamma(1 + bd) * gamma(i2 + mu/sigma.plus) * gamma(1/sigma.plus) * gamma(
bd - (-1 + mu + i2 * sigma.plus)/sigma.plus) * (digamma(bd + 1/sigma.plus) -
mu * digamma(i2 + mu/sigma.plus) - digamma(1/sigma.plus) +
digamma((1 - mu)/sigma.plus) - mu * digamma((1 - mu)/sigma.plus) +
mu * digamma(mu/sigma.plus) -
digamma(bd - (-1 + mu + i2 * sigma.plus)/sigma.plus) +
mu * digamma(bd - (-1 + mu + i2 * sigma.plus)/sigma.plus)))/(sigma.plus^2 * gamma(
1 + bd - i2) * gamma(1 + i2) * gamma(bd + 1/sigma.plus) * gamma((1 - mu)/
sigma.plus) * gamma(mu/sigma.plus))
d2f2sigma2 <- (df2.sigma.plus - df2.sigma) / fd.prec
df2.dsigma.plus <-( (gamma(1 + bd) * gamma(i2 + mu/sigma.plus) * gamma(1/sigma.plus) * gamma(
bd - (-1 + mu + i2 * sigma.plus)/sigma.plus) * (digamma(i2 + mu/sigma.plus) +
digamma((1 - mu)/sigma.plus) - digamma(mu/sigma.plus) -
digamma(bd - (-1 + mu + i2 * sigma.plus)/sigma.plus)))/(sigma.plus * gamma(
1 + bd - i2) * gamma(1 + i2) * gamma(bd + 1/sigma.plus) * gamma((1 - mu)/
sigma.plus) * gamma(mu/sigma.plus)) )*(1-mu)*mu
df2.dsigma.plus <- as.vector(df2.dsigma.plus)
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec
if(univariate==FALSE) {
cdf.f.plus <- rep(0, ly)
nmu <- rep(mu.plus, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
pdfall <- df2.mu.int.BB(allval, mm, ms, bd)
cdf.f.plus[i] <- sum(pdfall)
}
dF2.mu.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.mu.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
cdf.f.sigma.plus <- rep(0, ly)
cdf.f.dsigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f.sigma.plus[i] <- sum(df2.sigma.int.BB(allval, mm, ms, mn, bd))
cdf.f.dsigma.plus[i] <- sum(df2.mu.int.BB(allval, mm, ms, mn, bd))
}
dF2.sigma.plus <- cdf.f.sigma.plus
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
dF2.dsigma.plus <- cdf.f.dsigma.plus
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
}
} else if(VC$margins[2]=="GEOM") {
i2 <- dat[,2]
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
f2 <- dGEOM(i2, mu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pGEOM(i2, mu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
# df2, dF2, dF22
p <- 1/(mu+1)
df2 <- ((1 - p)^i2 - (1 - p)^(i2 - 1) * i2 * p)*(-(1/(mu + 1)^2))*mu
if(univariate==FALSE) {
df2.int.GEOM <- function(y, mu) {
p <- 1/(mu+1)
((1 - p)^y - (1 - p)^(y - 1) * y * p)*(-(1/(mu + 1)^2))*mu
}
ly <- length(i2)
cdf.f <- rep(0, ly)
nmu <- rep(mu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
allval <- seq(0, y.y)
pdfall <- df2.int.GEOM(allval, mu = mm)
cdf.f[i] <- sum(pdfall)
}
dF2 <- cdf.f
dF22 <- dF2-df2
}
# Hessian derivative components
mu.plus <- exp(eta2+fd.prec)
p.plus <- 1/(mu.plus+1)
df2.plus <- ((1 - p.plus)^i2 - (1 - p.plus)^(i2 - 1) * i2 * p.plus)*(-(1/(mu.plus + 1)^2))*mu.plus
d2f2delta22 <- (df2.plus - df2)/fd.prec
if(univariate==FALSE) {
ly <- length(i2)
cdf.f.plus <- rep(0, ly)
nmu.plus <- rep(mu.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu.plus[i]
allval <- seq(0, y.y)
pdfall <- df2.int.GEOM(allval, mu = mm)
cdf.f.plus[i] <- sum(pdfall)
}
dF2.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
}
} else if(VC$margins[2]=="LG") {
i2 <- dat[,2]
mu <- plogis(eta2)
mu <- ifelse(mu<(1-precision), mu, 1-mu)
mu <- ifelse(mu>precision, mu, precision)
f2 <- dLG(i2, mu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pLG(i2, mu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
# df2, dF2, dF22
df2 <- ( ((log(1 - mu))^((-1) - 1) * ((-1) * (1/(1 - mu))) * mu^i2 + (-(log(1 - mu))^(-1)) * (mu^(i2 - 1) * i2))/i2 )*(1-mu)*mu
if(univariate==FALSE) {
df2.int.LG <- function(y, mu) {
( ((log(1 - mu))^((-1) - 1) * ((-1) * (1/(1 - mu))) * mu^y + (-(log(1 - mu))^(-1)) * (mu^(y - 1) * y))/y )*(1-mu)*mu
}
ly <- length(i2)
cdf.f <- rep(0, ly)
nmu <- rep(mu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
allval <- seq(1, y.y)
pdfall <- df2.int.LG(allval, mu = mm)
cdf.f[i] <- sum(pdfall)
}
dF2 <- cdf.f
dF22 <- dF2-df2
}
# Hessian derivative components
mu.plus <- plogis(eta2+fd.prec)
df2.plus <- ( ((log(1 - mu.plus))^((-1) - 1) * ((-1) * (1/(1 - mu.plus))) * mu.plus^i2 + (-(log(1 - mu.plus))^(-1)) * (mu.plus^(i2 - 1) * i2))/i2 )*(1-mu.plus)*mu.plus
d2f2delta22 <- (df2.plus - df2)/fd.prec
if(univariate==FALSE) {
ly <- length(i2)
cdf.f.plus <- rep(0, ly)
nmu.plus <- rep(mu.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu.plus[i]
allval <- seq(1, y.y)
pdfall <- df2.int.LG(allval, mu = mm)
cdf.f.plus[i] <- sum(pdfall)
}
dF2.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
}
} else if (VC$margins[2]=="NBII") {
i2 <- dat[,2]
if (univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
}
sigma <- exp(sigma.star)
sigma <- ifelse(sigma>precision, sigma, precision)
sigma <- ifelse(sigma<10^7, sigma, 10^7)
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
f2 <- dNBII(i2, mu=mu, sigma=sigma)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pNBII(i2, mu=mu, sigma=sigma)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
f2.fd.mu <- dNBII(i2, mu = (mu+fd.prec), sigma = sigma)
f2.fd.mu <- ifelse(f2.fd.mu>precision, f2.fd.mu, precision)
df2 <- (f2.fd.mu - f2) / (fd.prec)
df2 <- as.vector(df2)*as.vector(mu)
f2.fd.sigma <- dNBII(i2, mu=mu, sigma=(sigma+fd.prec))
f2.fd.sigma <- ifelse(f2.fd.sigma>precision, f2.fd.sigma, precision)
df2.sigma <- (f2.fd.sigma - f2) / (fd.prec)
if(univariate==FALSE) {
F2.fd.mu <- pNBII(i2, mu = (mu+fd.prec), sigma = sigma)
F2.fd.mu <- ifelse(F2.fd.mu<(1-precision), F2.fd.mu, 1-precision)
dF2 <- (F2.fd.mu-F2) / (fd.prec)
dF2 <- as.vector(dF2)*as.vector(mu)
dF22 <- dF2 - df2
F2.fd.sigma <- pNBII(i2, mu = mu, sigma = (sigma+fd.prec))
F2.fd.sigma <- ifelse(F2.fd.sigma<(1-precision), F2.fd.sigma, 1-precision)
dF2.sigma <- (F2.fd.sigma - F2) / (fd.prec)
dF22.sigma <- dF2.sigma - df2.sigma
}
# Hessian derivative components
#fd.prec2 <- 10^-3
fd.prec2 <- sqrt(eps)
mu.plus <- ifelse(exp(eta2+fd.prec2)>0.0001, exp(eta2+fd.prec2), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
f2.plus <- dNBII(i2, mu = mu.plus, sigma = sigma)
f2.plus <- ifelse(f2.plus>precision, f2.plus, precision)
f2.fd.mu.plus <- dNBII(i2, mu = (mu.plus+fd.prec), sigma = sigma)
f2.fd.mu.plus <- ifelse(f2.fd.mu.plus>precision, f2.fd.mu.plus, precision)
df2.plus <- (f2.fd.mu.plus - f2.plus) / (fd.prec)
df2.plus <- as.vector(df2.plus)*as.vector(mu.plus)
df2.plus <- as.vector(df2.plus)
d2f2delta22 <- (df2.plus - df2) / fd.prec2
sigma.plus <- sigma + fd.prec2
f2.plus.sigma <- dNBII(i2, mu = mu, sigma = sigma.plus)
f2.plus.sigma <- ifelse(f2.plus.sigma>precision, f2.plus.sigma, precision)
f2.fd.sigma.plus <- dNBII(i2, mu = mu, sigma = (sigma.plus+fd.prec))
f2.fd.sigma.plus <- ifelse(f2.fd.sigma.plus>precision, f2.fd.sigma.plus, precision)
df2.sigma.plus <- (f2.fd.sigma.plus - f2.plus.sigma) / (fd.prec)
d2f2sigma2 <- (df2.sigma.plus - df2.sigma) / fd.prec2
f2.dsigma.plus <- dNBII(i2, mu = mu, sigma = sigma.plus)
f2.dsigma.plus <- ifelse(f2.dsigma.plus>precision, f2.dsigma.plus, precision)
f2.fd.mu.sigma.plus <- dNBII(i2, mu = (mu+fd.prec), sigma = sigma.plus)
f2.fd.mu.sigma.plus <- ifelse(f2.fd.mu.sigma.plus>precision, f2.fd.mu.sigma.plus, precision)
df2.dsigma.plus <- (f2.fd.mu.sigma.plus - f2.dsigma.plus) / (fd.prec)
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*as.vector(mu)
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec2
if(univariate==FALSE) {
F2.plus <- pNBII(i2, mu = mu.plus, sigma = sigma)
F2.fd.mu.plus <- pNBII(i2, mu = (mu.plus+fd.prec), sigma = sigma)
F2.fd.mu.plus <- ifelse(F2.fd.mu.plus<(1-precision), F2.fd.mu.plus, 1-precision)
dF2.plus <- (F2.fd.mu.plus - F2.plus) / (fd.prec)
dF2.plus <- as.vector(dF2.plus)*as.vector(mu.plus)
d2F2ddelta22 <- (dF2.plus - dF2) / fd.prec2
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
F2.sigma.plus <- pNBII(i2, mu = mu, sigma = sigma.plus)
F2.fd.sigma.plus <- pNBII(i2, mu = mu, sigma = (sigma.plus+fd.prec))
F2.fd.sigma.plus <- ifelse(F2.fd.sigma.plus<(1-precision), F2.fd.sigma.plus, 1-precision)
dF2.sigma.plus <- (F2.fd.sigma.plus - F2.sigma.plus) / (fd.prec)
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) /fd.prec2
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
F2.fd.dsigma.plus <- pNBII(i2, mu = (mu+fd.prec), sigma = sigma.plus)
F2.fd.dsigma.plus <- ifelse(F2.fd.dsigma.plus<(1-precision), F2.fd.dsigma.plus, 1-precision)
dF2.dsigma.plus <- (F2.fd.dsigma.plus - F2.sigma.plus) / (fd.prec)
dF2.dsigma.plus <- as.vector(dF2.dsigma.plus)*as.vector(mu)
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec2
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
}
} else if (VC$margins[2]=="WARING") {
i2 <- dat[,2]
if (univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
}
sigma <- exp(sigma.star)
sigma <- ifelse(sigma>precision, sigma, precision)
sigma <- ifelse(sigma<10^7, sigma, 10^7)
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
f2 <- dWARING(i2, mu=mu, sigma=sigma)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pWARING(i2, mu=mu, sigma=sigma)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
f2.fd.mu <- dWARING(i2, mu = (mu+fd.prec), sigma = sigma)
f2.fd.mu <- ifelse(f2.fd.mu>precision, f2.fd.mu, precision)
df2 <- (f2.fd.mu - f2) / (fd.prec)
df2 <- as.vector(df2)*as.vector(mu)
f2.fd.sigma <- dWARING(i2, mu=mu, sigma=(sigma+fd.prec))
f2.fd.sigma <- ifelse(f2.fd.sigma>precision, f2.fd.sigma, precision)
df2.sigma <- (f2.fd.sigma - f2) / (fd.prec)
if(univariate==FALSE) {
F2.fd.mu <- pWARING(i2, mu = (mu+fd.prec), sigma = sigma)
F2.fd.mu <- ifelse(F2.fd.mu<(1-precision), F2.fd.mu, 1-precision)
dF2 <- (F2.fd.mu-F2) / (fd.prec)
dF2 <- as.vector(dF2)*as.vector(mu)
dF22 <- dF2 - df2
F2.fd.sigma <- pWARING(i2, mu = mu, sigma = (sigma+fd.prec))
F2.fd.sigma <- ifelse(F2.fd.sigma<(1-precision), F2.fd.sigma, 1-precision)
dF2.sigma <- (F2.fd.sigma - F2) / (fd.prec)
dF22.sigma <- dF2.sigma - df2.sigma
}
# Hessian derivative components
#fd.prec2 <- 10^-3
fd.prec2 <- sqrt(eps)
mu.plus <- ifelse(exp(eta2+fd.prec2)>0.0001, exp(eta2+fd.prec2), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
f2.plus <- dWARING(i2, mu = mu.plus, sigma = sigma)
f2.plus <- ifelse(f2.plus>precision, f2.plus, precision)
f2.fd.mu.plus <- dWARING(i2, mu = (mu.plus+fd.prec), sigma = sigma)
f2.fd.mu.plus <- ifelse(f2.fd.mu.plus>precision, f2.fd.mu.plus, precision)
df2.plus <- (f2.fd.mu.plus - f2.plus) / (fd.prec)
df2.plus <- as.vector(df2.plus)*as.vector(mu.plus)
df2.plus <- as.vector(df2.plus)
d2f2delta22 <- (df2.plus - df2) / fd.prec2
sigma.plus <- sigma + fd.prec2
f2.plus.sigma <- dWARING(i2, mu = mu, sigma = sigma.plus)
f2.plus.sigma <- ifelse(f2.plus.sigma>precision, f2.plus.sigma, precision)
f2.fd.sigma.plus <- dWARING(i2, mu = mu, sigma = (sigma.plus+fd.prec))
f2.fd.sigma.plus <- ifelse(f2.fd.sigma.plus>precision, f2.fd.sigma.plus, precision)
df2.sigma.plus <- (f2.fd.sigma.plus - f2.plus.sigma) / (fd.prec)
d2f2sigma2 <- (df2.sigma.plus - df2.sigma) / fd.prec2
f2.dsigma.plus <- dWARING(i2, mu = mu, sigma = sigma.plus)
f2.dsigma.plus <- ifelse(f2.dsigma.plus>precision, f2.dsigma.plus, precision)
f2.fd.mu.sigma.plus <- dWARING(i2, mu = (mu+fd.prec), sigma = sigma.plus)
f2.fd.mu.sigma.plus <- ifelse(f2.fd.mu.sigma.plus>precision, f2.fd.mu.sigma.plus, precision)
df2.dsigma.plus <- (f2.fd.mu.sigma.plus - f2.dsigma.plus) / (fd.prec)
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*as.vector(mu)
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec2
if(univariate==FALSE) {
F2.plus <- pWARING(i2, mu = mu.plus, sigma = sigma)
F2.fd.mu.plus <- pWARING(i2, mu = (mu.plus+fd.prec), sigma = sigma)
F2.fd.mu.plus <- ifelse(F2.fd.mu.plus<(1-precision), F2.fd.mu.plus, 1-precision)
dF2.plus <- (F2.fd.mu.plus - F2.plus) / (fd.prec)
dF2.plus <- as.vector(dF2.plus)*as.vector(mu.plus)
d2F2ddelta22 <- (dF2.plus - dF2) / fd.prec2
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
F2.sigma.plus <- pWARING(i2, mu = mu, sigma = sigma.plus)
F2.fd.sigma.plus <- pWARING(i2, mu = mu, sigma = (sigma.plus+fd.prec))
F2.fd.sigma.plus <- ifelse(F2.fd.sigma.plus<(1-precision), F2.fd.sigma.plus, 1-precision)
dF2.sigma.plus <- (F2.fd.sigma.plus - F2.sigma.plus) / (fd.prec)
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) /fd.prec2
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
F2.fd.dsigma.plus <- pWARING(i2, mu = (mu+fd.prec), sigma = sigma.plus)
F2.fd.dsigma.plus <- ifelse(F2.fd.dsigma.plus<(1-precision), F2.fd.dsigma.plus, 1-precision)
dF2.dsigma.plus <- (F2.fd.dsigma.plus - F2.sigma.plus) / (fd.prec)
dF2.dsigma.plus <- as.vector(dF2.dsigma.plus)*as.vector(mu)
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec2
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
}
} else if (VC$margins[2]=="YULE") {
i2 <- dat[,2]
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
f2 <- dYULE(i2, mu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pYULE(i2, mu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
f2.fd.mu <- dYULE(i2, mu = (mu+fd.prec))
f2.fd.mu <- ifelse(f2.fd.mu>precision, f2.fd.mu, precision)
df2 <- (f2.fd.mu - f2) / (fd.prec)
df2 <- as.vector(df2)*as.vector(mu)
if(univariate==FALSE) {
F2.fd.mu <- pYULE(i2, mu = (mu+fd.prec))
F2.fd.mu <- ifelse(F2.fd.mu<(1-precision), F2.fd.mu, 1-precision)
dF2 <- (F2.fd.mu-F2) / (fd.prec)
dF2 <- as.vector(dF2)*as.vector(mu)
dF22 <- dF2 - df2
}
# Hessian derivative components
#fd.prec2 <- 10^-3
fd.prec2 <- sqrt(eps)
mu.plus <- ifelse(exp(eta2+fd.prec2)>0.0001, exp(eta2+fd.prec2), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
f2.plus <- dYULE(i2, mu = mu.plus)
f2.plus <- ifelse(f2.plus>precision, f2.plus, precision)
f2.fd.mu.plus <- dYULE(i2, mu = (mu.plus+fd.prec))
f2.fd.mu.plus <- ifelse(f2.fd.mu.plus>precision, f2.fd.mu.plus, precision)
df2.plus <- (f2.fd.mu.plus - f2.plus) / (fd.prec)
df2.plus <- as.vector(df2.plus)*as.vector(mu.plus)
df2.plus <- as.vector(df2.plus)
d2f2delta22 <- (df2.plus - df2) / fd.prec2
if(univariate==FALSE) {
F2.plus <- pYULE(i2, mu = mu.plus)
F2.fd.mu.plus <- pYULE(i2, mu = (mu.plus+fd.prec))
F2.fd.mu.plus <- ifelse(F2.fd.mu.plus<(1-precision), F2.fd.mu.plus, 1-precision)
dF2.plus <- (F2.fd.mu.plus - F2.plus) / (fd.prec)
dF2.plus <- as.vector(dF2.plus)*as.vector(mu.plus)
d2F2ddelta22 <- (dF2.plus - dF2) / fd.prec2
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
}
} else if (VC$margins[2]=="ZABB") {
i2 <- dat[,2]
if(univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
if (is.null(VC$X3)) nu.star <- params[(VC$X2.d2 + 2)]
if (!is.null(VC$X3)) nu.star <- etanu <- X4 %*% params[(VC$X2.d2 + VC$X3.d2 + 1):(VC$X2.d2 + VC$X3.d2 + VC$X4.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
if(is.null(VC$X3)) nu.star <- params[(VC$X1.d2+VC$X2.d2+2)]
if(!is.null(VC$X3)) nu.star <- etanu <- X4%*%params[(VC$X1.d2+VC$X2.d2+VC$X3.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2)]
}
nu <- plogis(nu.star)
nu <- ifelse(nu<precision, precision, nu)
nu <- ifelse(nu>(1-precision), 1-precision, nu)
sigma <- exp(sigma.star)
sigma <- ifelse(sigma>precision, sigma, precision)
sigma <- ifelse(sigma<10^7, sigma, 10^7)
mu <- plogis(eta2)
mu <- ifelse(mu<precision, precision, mu)
mu <- ifelse(mu>(1-precision), 1-precision, mu)
f2 <- dZABB(i2, bd=bd, mu=mu, sigma=sigma, nu=nu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pZABB(i2, bd=bd, mu=mu, sigma=sigma, nu=nu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
f2.fd.mu <- dZABB(i2, mu = (mu +fd.prec), sigma = sigma, nu = nu, bd=bd)
f2.fd.mu <- ifelse(f2.fd.mu>precision, f2.fd.mu, precision)
df2 <- (f2.fd.mu - f2) / (fd.prec)
df2 <- as.vector(df2)*as.vector(mu * (1-mu))
f2.fd.sigma <- dZABB(i2, mu = mu, sigma = (sigma+fd.prec), nu = nu, bd=bd)
f2.fd.sigma <- ifelse(f2.fd.sigma>precision, f2.fd.sigma, precision)
df2.sigma <- (f2.fd.sigma - f2) / (fd.prec)
f2.fd.nu <- dZABB(i2, mu = mu, sigma = sigma, nu =(nu+fd.prec), bd=bd)
f2.fd.nu <- ifelse(f2.fd.nu>precision, f2.fd.nu, precision)
df2.nu <- (f2.fd.nu - f2) / (fd.prec)
if(univariate==FALSE) {
F2.fd.mu <- pZABB(i2, mu = (mu+fd.prec), sigma = sigma, nu = nu, bd=bd)
F2.fd.mu <- ifelse(F2.fd.mu<(1-precision), F2.fd.mu, 1-precision)
dF2 <- (F2.fd.mu - F2) / (fd.prec)
dF2 <- as.vector(dF2)*as.vector(mu * (1-mu))
dF22 <- dF2 - df2
F2.fd.sigma <- pZABB(i2, mu = mu, sigma = (sigma+fd.prec), nu = nu, bd=bd)
F2.fd.sigma <- ifelse(F2.fd.sigma<(1-precision), F2.fd.sigma, 1-precision)
dF2.sigma <- (F2.fd.sigma - F2) / (fd.prec)
dF22.sigma <- dF2.sigma - df2.sigma
F2.fd.nu <- pZABB(i2, mu = mu, sigma = sigma, nu = (nu+fd.prec), bd=bd)
F2.fd.nu <- ifelse(F2.fd.nu<(1-precision), F2.fd.nu, 1-precision)
dF2.nu <- (F2.fd.nu - F2) / (fd.prec)
dF22.nu <- dF2.nu - df2.nu
}
# Hessian derivative components
#second derivative fd precision
#fd.prec2 <- 10^-3
fd.prec2 <- sqrt(eps)
#pmfs
#delta2^2
mu.plus <- ifelse(plogis(eta2+fd.prec2)>0.0001, plogis(eta2+fd.prec2), 0.0001)
mu.plus <- ifelse(mu.plus>(1-precision), 1-precision, mu.plus)
f2.plus <- dZABB(i2, mu = mu.plus, sigma = sigma, nu = nu, bd=bd)
f2.plus <- ifelse(f2.plus>precision, f2.plus, precision)
f2.fd.mu.plus <- dZABB(i2, mu = (mu.plus+fd.prec), sigma = sigma, nu = nu, bd=bd)
f2.fd.mu.plus <- ifelse(f2.fd.mu.plus>precision, f2.fd.mu.plus, precision)
df2.plus <- (f2.fd.mu.plus - f2.plus) / (fd.prec)
df2.plus <- as.vector(df2.plus)*as.vector(mu.plus * (1-mu.plus))
df2.plus <- as.vector(df2.plus)
d2f2delta22 <- (df2.plus - df2) / fd.prec2
#sigma^2
sigma.plus <- sigma+fd.prec2
f2.plus.sigma <- dZABB(i2, mu = mu, sigma = sigma.plus, nu = nu, bd=bd)
f2.plus.sigma <- ifelse(f2.plus.sigma>precision, f2.plus.sigma, precision)
f2.fd.sigma.plus <- dZABB(i2, mu = mu, sigma = (sigma.plus+fd.prec), nu = nu, bd=bd)
f2.fd.sigma.plus <- ifelse(f2.fd.sigma.plus>precision, f2.fd.sigma.plus, precision)
df2.sigma.plus <- (f2.fd.sigma.plus - f2.plus.sigma) / (fd.prec)
d2f2sigma2 <- (df2.sigma.plus - df2.sigma) / fd.prec2
#delta2sigma
f2.dsigma.plus <- dZABB(i2, mu = mu, sigma = sigma.plus, nu = nu, bd=bd)
f2.dsigma.plus <- ifelse(f2.dsigma.plus>precision, f2.dsigma.plus, precision)
f2.fd.mu.sigma.plus <- dZABB(i2, mu = (mu+fd.prec), sigma = sigma.plus, nu = nu, bd=bd)
f2.fd.mu.sigma.plus <- ifelse(f2.fd.mu.sigma.plus>precision, f2.fd.mu.sigma.plus, precision)
df2.dsigma.plus <- (f2.fd.mu.sigma.plus - f2.dsigma.plus) / (fd.prec)
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*as.vector(mu * (1-mu))
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec2
#nu^2
nu.plus <- nu + fd.prec2
nu.plus <- ifelse(nu.plus>=(1-precision), (1-precision), nu.plus)
f2.plus.nu <- dZABB(i2, mu = mu, sigma = sigma, nu = nu.plus, bd=bd)
f2.plus.nu <- ifelse(f2.plus.nu>precision, f2.plus.nu, precision)
f2.fd.nu.plus <- dZABB(i2, mu = mu, sigma = sigma, nu = nu.plus+fd.prec, bd=bd)
f2.fd.nu.plus <- ifelse(f2.fd.nu.plus>precision, f2.fd.nu.plus, precision)
df2.nu.plus <- (f2.fd.nu.plus - f2.plus.nu) / (fd.prec)
d2f2nu2 <- (df2.nu.plus - df2.nu) / fd.prec2
#delta2nu
f2.dnu.plus <- dZABB(i2, mu = mu, sigma = sigma, nu = nu.plus, bd=bd)
f2.dnu.plus <- ifelse(f2.dnu.plus>precision, f2.dnu.plus, precision)
f2.fd.mu.nu.plus <- dZABB(i2, mu = (mu+fd.prec), sigma = sigma, nu = nu.plus, bd=bd)
f2.fd.mu.nu.plus <- ifelse(f2.fd.mu.nu.plus>precision, f2.fd.mu.nu.plus, precision)
df2.dnu.plus <- (f2.fd.mu.nu.plus - f2.dnu.plus) / (fd.prec)
df2.dnu.plus <- as.vector(df2.dnu.plus)*as.vector(mu * (1-mu))
d2f2delta2nu <-(df2.dnu.plus - df2) / fd.prec2
#nusigma
f2.fd.sigma.nu.plus <- dZABB(i2, mu = mu, sigma = (sigma+fd.prec), nu = nu.plus, bd=bd)
f2.fd.sigma.nu.plus <- ifelse(f2.fd.sigma.nu.plus>precision, f2.fd.sigma.nu.plus, precision)
df2.dnu.sigma.plus <- (f2.fd.sigma.nu.plus - f2.plus.nu) / (fd.prec)
df2.dnu.sigma.plus <- as.vector(df2.dnu.sigma.plus)
d2f2sigmanu <- (df2.dnu.sigma.plus - df2.sigma) / fd.prec2
if(univariate==FALSE) {
#cdfs
#delta2^2
F2.plus <- pZABB(i2, mu = mu.plus, sigma = sigma, nu = nu, bd=bd)
F2.fd.mu.plus <- pZABB(i2, mu = (mu.plus+fd.prec), sigma = sigma, nu = nu, bd=bd)
F2.fd.mu.plus <- ifelse(F2.fd.mu.plus<(1-precision), F2.fd.mu.plus, 1-precision)
dF2.plus <- (F2.fd.mu.plus - F2.plus) / (fd.prec)
dF2.plus <- as.vector(dF2.plus)*as.vector(mu.plus * (1-mu.plus))
d2F2ddelta22 <-(dF2.plus - dF2) / fd.prec2
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
#sigma^2
F2.sigma.plus <- pZABB(i2, mu = mu, sigma = sigma.plus, nu = nu, bd=bd)
F2.fd.sigma.plus <- pZABB(i2, mu = mu, sigma = (sigma.plus+fd.prec), nu=nu, bd=bd)
F2.fd.sigma.plus <- ifelse(F2.fd.sigma.plus<(1-precision), F2.fd.sigma.plus, 1-precision)
dF2.sigma.plus <- (F2.fd.sigma.plus - F2.sigma.plus) / (fd.prec)
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec2
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
#delta2sigma
F2.fd.dsigma.plus <- pZABB(i2, mu = (mu+fd.prec), sigma = sigma.plus, nu, bd=bd)
F2.fd.dsigma.plus <- ifelse(F2.fd.dsigma.plus<(1-precision), F2.fd.dsigma.plus, 1-precision)
dF2.dsigma.plus <- (F2.fd.dsigma.plus - F2.sigma.plus) / (fd.prec)
dF2.dsigma.plus <- as.vector(dF2.dsigma.plus)*as.vector(mu * (1-mu))
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec2
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
#nu^2
F2.nu.plus<- pZABB(i2, mu = mu, sigma = sigma, nu = nu.plus, bd=bd)
F2.fd.nu.plus <- pZABB(i2, mu = mu, sigma = sigma, nu = nu.plus+fd.prec, bd=bd)
F2.fd.nu.plus <- ifelse(F2.fd.nu.plus<(1-precision), F2.fd.nu.plus, 1-precision)
dF2.nu.plus <- (F2.fd.nu.plus - F2.nu.plus) / (fd.prec)
d2F2dnu2 <- (dF2.nu.plus - dF2.nu) / fd.prec2
d2F22dnu2 <- d2F2dnu2 - d2f2nu2
#delta2nu
F2.fd.dnu.plus <- pZABB(i2, mu = (mu+fd.prec), sigma = sigma, nu.plus, bd=bd)
F2.fd.dnu.plus <- ifelse(F2.fd.dnu.plus<(1-precision), F2.fd.dnu.plus, 1-precision)
dF2.dnu.plus <- (F2.fd.dnu.plus - F2.nu.plus) / (fd.prec)
dF2.dnu.plus <- as.vector(dF2.dnu.plus)*as.vector(mu * (1-mu))
d2F2ddelta2dnu <- (dF2.dnu.plus - dF2) / fd.prec2
d2F22ddelta2dnu <- d2F2ddelta2dnu - d2f2delta2nu
#sigmanu
F2.fd.dnu.sigma.plus <- pZABB(i2, mu = mu, sigma = (sigma+fd.prec), nu.plus, bd=bd)
F2.fd.dnu.sigma.plus <- ifelse(F2.fd.dnu.sigma.plus<(1-precision), F2.fd.dnu.sigma.plus, 1-precision)
dF2.dnu.sigma.plus <- (F2.fd.dnu.sigma.plus - F2.nu.plus) / (fd.prec)
dF2.dnu.sigma.plus <- as.vector(dF2.dnu.sigma.plus)
d2F2dsigmadnu <- (dF2.dnu.sigma.plus - dF2.sigma) / fd.prec2
d2F22dsigmadnu <- d2F2dsigmadnu - d2f2sigmanu
}
} else if (VC$margins[2]=="ZABI") {
i2 <- dat[,2]
if(univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
}
sigma <- plogis(sigma.star)
sigma <- ifelse(sigma<(1-precision), sigma, 1-sigma)
sigma <- ifelse(sigma>precision, sigma, precision)
mu <- plogis(eta2)
mu <- ifelse(mu<(1-precision), mu, 1-mu)
mu <- ifelse(mu>precision, mu, precision)
f2 <- dZABI(i2, sigma = sigma, mu = mu, bd=bd)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pZABI(i2, sigma = sigma, mu = mu, bd=bd)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
df2<- ifelse(i2>0, ((1 - mu)^(-1 + bd - i2) * mu^(-1 + i2) * (i2 * (-1 + (1 - mu)^bd) + bd * mu) *
(-1 + sigma) * gamma(1 + bd))/((-1 + (1 - mu)^bd)^2 * gamma(1 + bd - i2) * gamma(1 + i2)), 0)
df2 <- as.vector(df2)*(1-mu)*mu
df2.sigma <- ifelse(i2>0, -(gamma(bd + 1) * mu^(i2) * (1 - mu)^(bd - i2)/((1 - (1 - mu)^(bd)) * gamma(i2 + 1) * gamma(bd - i2 + 1))), 1)
if(univariate==FALSE) {
df2.mu.int.ZABI<-function(y, mu, sigma, n) {
df2<- ifelse(y>0, ((1 - mu)^(-1 + n - y) * mu^(-1 + y) * (y * (-1 + (1 - mu)^n) + n * mu) *
(-1 + sigma) * gamma(1 + n))/((-1 + (1 - mu)^n)^2 *
gamma(1 + n - y) * gamma(1 + y)), 0)
df2 <- as.vector(df2)*(1-mu)*mu
df2
}
df2.sigma.int.ZABI<-function(y, mu, sigma, n) {
df2<- ifelse(y>0, -(gamma(n + 1) * mu^(y) * (1 - mu)^(n - y)/((1 - (1 - mu)^(n)) * gamma(y + 1) * gamma(n - y + 1))), 1)
df2
}
ly <- length(i2)
cdf.f <- rep(0, ly)
cdf.f.sigma <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f[i] <- sum(df2.mu.int.ZABI(allval, mm, ms, bd))
cdf.f.sigma[i] <- sum(df2.sigma.int.ZABI(allval, mm, ms, bd))
}
dF2<-cdf.f
dF22<-dF2-df2
dF2.sigma <- cdf.f.sigma
dF22.sigma <- dF2.sigma-df2.sigma
}
# Hessian derivative components
mu.plus <- ifelse(plogis(eta2+fd.prec)>0.0001, plogis(eta2+fd.prec), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
df2.plus <- ifelse(i2>0, ((1 - mu.plus)^(-1 + bd - i2) * mu.plus^(-1 + i2) * (i2 * (-1 +
(1 - mu.plus)^bd) + bd * mu.plus) * (-1 + sigma) * gamma(1 + bd))/((-1 + (1 - mu.plus)^bd)^2 *
gamma(1 + bd - i2) * gamma(1 + i2)), 0)
df2.plus <- as.vector(df2.plus)*(1-mu.plus)*mu.plus
d2f2delta22 <- (df2.plus - df2) / fd.prec
sigma.plus <- sigma + fd.prec
d2f2sigma2 <- 0
df2.dsigma.plus <- ifelse(i2>0, ((1 - mu)^(-1 + bd - i2) * mu^(-1 + i2) * (i2 * (-1 + (1 - mu)^bd) + bd * mu) *
(-1 + sigma.plus) * gamma(1 + bd))/((-1 + (1 - mu)^bd)^2 *
gamma(1 + bd - i2) * gamma(1 + i2)), 0)
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*(1-mu)*mu
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec
if(univariate==FALSE) {
cdf.f.plus <- rep(0, ly)
nmu <- rep(mu.plus, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f.plus[i] <- sum(df2.mu.int.ZABI(allval, mm, ms, bd))
}
dF2.mu.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.mu.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
cdf.f.sigma.plus <- rep(0, ly)
cdf.f.dsigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f.sigma.plus[i] <- sum(df2.sigma.int.ZABI(allval, mm, ms, bd))
cdf.f.dsigma.plus[i] <- sum(df2.mu.int.ZABI(allval, mm, ms, bd))
}
dF2.sigma.plus <- cdf.f.sigma.plus
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
dF2.dsigma.plus <- cdf.f.dsigma.plus
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
}
} else if (VC$margins[2]=="ZALG") {
i2 <- dat[,2]
if(univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
}
sigma <- plogis(sigma.star)
sigma <- ifelse(sigma<(1-precision), sigma, 1-sigma)
sigma <- ifelse(sigma>precision, sigma, precision)
mu <- plogis(eta2)
mu <- ifelse(mu<(1-precision), mu, 1-mu)
mu <- ifelse(mu>precision, mu, precision)
f2 <- dZALG(i2, sigma = sigma, mu = mu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pZALG(i2, sigma = sigma, mu = mu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
df2<- ifelse(i2>0, (mu^(-1 + i2) * (-1 + sigma) * (-mu + i2 * (-1 + mu) * log(1 - mu)))/(i2 * (-1 + mu) * log(1 - mu)^2), 0)
df2 <- as.vector(df2)*(1-mu)*mu
df2.sigma <- ifelse(i2>0, -((-(log(1 - mu))^(-1)) * mu^(i2)/i2), 1)
if(univariate==FALSE) {
df2.mu.int.ZALG<-function(y, mu, sigma) {
df2<- ifelse(y>0, (mu^(-1 + y) * (-1 + sigma) * (-mu + y * (-1 + mu) * log(1 - mu)))/(y * (-1 + mu) * log(1 - mu)^2), 0)
df2 <- as.vector(df2)*(1-mu)*mu
df2
}
df2.sigma.int.ZALG<-function(y, mu, sigma) {
df2<- ifelse(y>0, -((-(log(1 - mu))^(-1)) * mu^(y)/y), 1)
df2
}
ly <- length(i2)
cdf.f <- rep(0, ly)
cdf.f.sigma <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f[i] <- sum(df2.mu.int.ZALG(allval, mm, ms))
cdf.f.sigma[i] <- sum(df2.sigma.int.ZALG(allval, mm, ms))
}
dF2<-cdf.f
dF22<-dF2-df2
dF2.sigma <- cdf.f.sigma
dF22.sigma <- dF2.sigma-df2.sigma
}
# Hessian derivative components
mu.plus <- ifelse(plogis(eta2+fd.prec)>0.0001, plogis(eta2+fd.prec), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
df2.plus <- ifelse(i2>0, (mu.plus^(-1 + i2) * (-1 + sigma) * (-mu.plus + i2 * (-1 + mu.plus) * log(1 - mu.plus)))/(i2 *
(-1 + mu.plus) * log(1 - mu.plus)^2), 0)
df2.plus <- as.vector(df2.plus)*(1-mu.plus)*mu.plus
d2f2delta22 <- (df2.plus - df2) / fd.prec
sigma.plus <- sigma + fd.prec
d2f2sigma2 <- 0
df2.dsigma.plus <- ifelse(i2>0, (mu^(-1 + i2) * (-1 + sigma.plus) * (-mu + i2 * (-1 + mu) * log(1 - mu)))/(i2 *
(-1 + mu) * log(1 - mu)^2), 0)
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*(1-mu)*mu
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec
if(univariate==FALSE) {
cdf.f.plus <- rep(0, ly)
nmu <- rep(mu.plus, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f.plus[i] <- sum(df2.mu.int.ZALG(allval, mm, ms))
}
dF2.mu.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.mu.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
cdf.f.sigma.plus <- rep(0, ly)
cdf.f.dsigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f.sigma.plus[i] <- sum(df2.sigma.int.ZALG(allval, mm, ms))
cdf.f.dsigma.plus[i] <- sum(df2.mu.int.ZALG(allval, mm, ms))
}
dF2.sigma.plus <- cdf.f.sigma.plus
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
dF2.dsigma.plus <- cdf.f.dsigma.plus
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
}
} else if (VC$margins[2]=="ZANBI") {
i2 <- dat[,2]
if(univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
if (is.null(VC$X3)) nu.star <- params[(VC$X2.d2 + 2)]
if (!is.null(VC$X3)) nu.star <- etanu <- X4 %*% params[(VC$X2.d2 + VC$X3.d2 + 1):(VC$X2.d2 + VC$X3.d2 + VC$X4.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
if(is.null(VC$X3)) nu.star <- params[(VC$X1.d2+VC$X2.d2+2)]
if(!is.null(VC$X3)) nu.star <- etanu <- X4%*%params[(VC$X1.d2+VC$X2.d2+VC$X3.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2)]
}
nu <- plogis(nu.star)
nu <- ifelse(nu<precision, precision, nu)
nu <- ifelse(nu>(1-precision), 1-precision, nu)
sigma <- exp(sigma.star)
sigma <- ifelse(sigma>precision, sigma, precision)
sigma <- ifelse(sigma<10^7, sigma, 10^7)
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
f2 <- dZANBI(i2, mu=mu, sigma=sigma, nu=nu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pZANBI(i2, mu=mu, sigma=sigma, nu=nu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
df2 <- ifelse(i2>0, ((-1 + nu) * (1/(1 + mu * sigma))^(1 + 1/sigma) * ((mu * sigma)/(1 + mu * sigma))^i2 * (mu +
i2 * (-1 + (1/(1 + mu * sigma))^(1/sigma))) * gamma(i2 + 1/sigma))/(mu *
(-1 + (1/(1 + mu * sigma))^(1/sigma))^2 * gamma(1 + i2) * gamma(1/sigma)), 0)
df2 <- as.vector(df2) * mu
df2.sigma <- ifelse(i2>0, -(((-1 + nu) * (1/(1 + mu * sigma))^(1 + 1/sigma) * ((mu * sigma)/(1 + mu * sigma))^i2 *
gamma(i2 + 1/sigma) * (i2 * sigma - mu * sigma -
i2 * sigma * (1/(1 + mu * sigma))^(1/sigma) - log(1/(1 + mu * sigma)) -
mu * sigma * log(1/(1 + mu * sigma)) + (1 + mu * sigma) * (-1 + (1/(1 + mu * sigma))^(1/
sigma)) * digamma(i2 + 1/sigma) - (1 + mu * sigma) * (-1 + (1/(1 + mu * sigma))^(1/
sigma)) * digamma(1/sigma)))/(sigma^2 * (-1 + (1/(1 + mu * sigma))^(1/sigma))^2 * gamma(1 + i2) *
gamma(1/sigma))), 0)
df2.nu <- ifelse(i2>0, (1/(1 + mu * sigma))^(1/sigma) * (mu * sigma/(1 + mu * sigma))^i2 *
gamma(1/sigma + i2)/((-1 + (1/(1 + sigma * mu))^(1/sigma)) *
gamma(1/sigma) * gamma(1 + i2)), 1)
if(univariate==FALSE) {
df2.mu.int.ZANBI <- function(y, mu, sigma, nu) {
df2 <- ifelse(y>0, ((-1 + nu) * (1/(1 + mu * sigma))^(1 + 1/sigma) * ((mu * sigma)/(1 + mu * sigma))^y * (mu +
y * (-1 + (1/(1 + mu * sigma))^(1/sigma))) * gamma(y + 1/sigma))/(mu *
(-1 + (1/(1 + mu * sigma))^(1/sigma))^2 * gamma(1 + y) * gamma(1/sigma)), 0)
df2 <- as.vector(df2) * mu
}
df2.sigma.int.ZANBI <-function(y, mu, sigma, nu) {
df2 <- ifelse(y>0, -(((-1 + nu) * (1/(1 + mu * sigma))^(1 + 1/sigma) * ((mu * sigma)/(1 + mu * sigma))^y *
gamma(y + 1/sigma) * (y * sigma - mu * sigma -
y * sigma * (1/(1 + mu * sigma))^(1/sigma) - log(1/(1 + mu * sigma)) -
mu * sigma * log(1/(1 + mu * sigma)) + (1 + mu * sigma) * (-1 + (1/(1 + mu * sigma))^(1/
sigma)) * digamma(y + 1/sigma) - (1 + mu * sigma) * (-1 + (1/(1 + mu * sigma))^(1/
sigma)) * digamma(1/sigma)))/(sigma^2 * (-1 + (1/(1 + mu * sigma))^(1/sigma))^2 * gamma(1 + y) *
gamma(1/sigma))), 0)
df2
}
df2.nu.int.ZANBI <- function(y, mu, sigma, nu) {
df2 <- ifelse(y>0, (1/(1 + mu * sigma))^(1/sigma) * (mu * sigma/(1 + mu * sigma))^y *
gamma(1/sigma + y)/((-1 + (1/(1 + sigma * mu))^(1/sigma)) *
gamma(1/sigma) * gamma(1 + y)), 1)
df2
}
ly <- length(i2)
cdf.f <- rep(0, ly)
cdf.f.sigma <- rep(0, ly)
cdf.f.nu <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
cdf.f[i] <- sum(df2.mu.int.ZANBI(allval, mm, ms, mn))
cdf.f.sigma[i] <- sum(df2.sigma.int.ZANBI(allval, mm, ms, mn))
cdf.f.nu[i] <- sum(df2.nu.int.ZANBI(allval, mm, ms, mn))
}
dF2<-cdf.f
dF22<-dF2-df2
dF2.sigma <- cdf.f.sigma
dF22.sigma <- dF2.sigma-df2.sigma
dF2.nu <- cdf.f.nu
dF22.nu <- dF2.nu-df2.nu
}
# Hessian derivative components
#fd.prec2 <- 10^-3
fd.prec2 <- sqrt(eps)
#delta2delta2
mu.plus <- ifelse(exp(eta2+fd.prec)>0.0001, exp(eta2+fd.prec), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
df2.plus <- ifelse(i2>0, ((-1 + nu) * (1/(1 + mu.plus * sigma))^(1 + 1/sigma) * ((mu.plus * sigma)/(1 + mu.plus * sigma))^i2 *
(mu.plus + i2 * (-1 + (1/(1 + mu.plus * sigma))^(1/sigma))) * gamma(i2 + 1/sigma))/(mu.plus *
(-1 + (1/(1 + mu.plus * sigma))^(1/sigma))^2 * gamma(1 + i2) * gamma(1/sigma)), 0)
df2.plus <- as.vector(df2.plus) * mu.plus
d2f2delta22 <- (df2.plus - df2) / fd.prec
#sigma^2
sigma.plus <- sigma + fd.prec
df2.sigma.plus <- ifelse(i2>0, -(((-1 + nu) * (1/(1 + mu * sigma.plus))^(1 + 1/sigma.plus) * ((mu * sigma.plus)/(1 + mu *
sigma.plus))^i2 * gamma(i2 + 1/sigma.plus) * (i2 * sigma.plus - mu * sigma.plus -
i2 * sigma.plus * (1/(1 + mu * sigma.plus))^(1/sigma.plus) - log(1/(1 + mu * sigma.plus)) -
mu * sigma.plus * log(1/(1 + mu * sigma.plus)) + (1 + mu * sigma.plus) * (-1 + (1/(1 + mu * sigma.plus))^(1/
sigma.plus)) * digamma(i2 + 1/sigma.plus) - (1 + mu * sigma.plus) * (-1 + (1/(1 + mu * sigma.plus))^(1/
sigma.plus)) * digamma(1/sigma.plus)))/(sigma.plus^2 * (-1 + (1/(1 + mu * sigma.plus))^(1/sigma.plus))^2 *
gamma(1 + i2) * gamma(1/sigma.plus))), 0)
d2f2sigma2 <- (df2.sigma.plus - df2.sigma) / fd.prec
#delta2sigma
df2.dsigma.plus <- ifelse(i2>0, ((-1 + nu) * (1/(1 + mu * sigma.plus))^(1 + 1/sigma.plus) *
((mu * sigma.plus)/(1 + mu * sigma.plus))^i2 * (mu +
i2 * (-1 + (1/(1 + mu * sigma.plus))^(1/sigma.plus))) * gamma(i2 + 1/sigma.plus))/(mu *
(-1 + (1/(1 + mu * sigma.plus))^(1/sigma.plus))^2 * gamma(1 + i2) * gamma(1/sigma.plus)), 0)
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*mu
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec
#nu^2
nu.plus <- nu + fd.prec
nu.plus <- ifelse(nu.plus>=(1-precision), (1-precision), nu.plus)
d2f2nu2 <- 0
#delta2nu
df2.dnu.plus <- ifelse(i2>0, ((-1 + nu.plus) * (1/(1 + mu * sigma))^(1 + 1/sigma) * ((mu * sigma)/(1 + mu * sigma))^i2 *
(mu + i2 * (-1 + (1/(1 + mu * sigma))^(1/sigma))) * gamma(i2 + 1/sigma))/(mu *
(-1 + (1/(1 + mu * sigma))^(1/sigma))^2 * gamma(1 + i2) * gamma(1/sigma)), 0)
df2.dnu.plus <- as.vector(df2.dnu.plus)*mu
d2f2delta2nu <- (df2.dnu.plus - df2) / fd.prec
#nusigma
df2.nu.sigma.plus <- ifelse(i2>0, (1/(1 + mu * sigma.plus))^(1/sigma.plus) * (mu * sigma.plus/(1 + mu * sigma.plus))^i2 *
gamma(1/sigma.plus + i2)/((-1 + (1/(1 + sigma.plus * mu))^(1/sigma.plus)) *
gamma(1/sigma.plus) * gamma(1 + i2)), 1)
d2f2sigmanu <- (df2.nu.sigma.plus - df2.nu) / fd.prec
if(univariate==FALSE) {
cdf.f.plus <- rep(0, ly)
nmu <- rep(mu.plus, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
pdfall <- df2.mu.int.ZANBI(allval, mm, ms, mn)
cdf.f.plus[i] <- sum(pdfall)
}
dF2.mu.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.mu.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
cdf.f.sigma.plus <- rep(0, ly)
cdf.f.dsigma.plus <- rep(0, ly)
cdf.f.dnu.dsigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
cdf.f.sigma.plus[i] <- sum(df2.sigma.int.ZANBI(allval, mm, ms, mn))
cdf.f.dsigma.plus[i] <- sum(df2.mu.int.ZANBI(allval, mm, ms, mn))
cdf.f.dnu.dsigma.plus[i] <- sum(df2.nu.int.ZANBI(allval, mm, ms, mn))
}
dF2.sigma.plus <- cdf.f.sigma.plus
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
dF2.dsigma.plus <- cdf.f.dsigma.plus
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
dF2.dnu.dsigma.plus <- cdf.f.dnu.dsigma.plus
d2F2dsigmadnu <- (dF2.dnu.dsigma.plus - dF2.nu) / fd.prec
d2F22dsigmadnu <- d2F2dsigmadnu - d2f2sigmanu
cdf.f.nu.plus <- rep(0, ly)
cdf.f.dnu.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
cdf.f.nu.plus[i] <- sum(df2.nu.int.ZANBI(allval, mm, ms, mn))
cdf.f.dnu.plus[i] <- sum(df2.mu.int.ZANBI(allval, mm, ms, mn))
}
dF2.nu.plus <- cdf.f.nu.plus
d2F2dnu2 <- (dF2.nu.plus - dF2.nu) / fd.prec
d2F22dnu2 <- d2F2dnu2 - d2f2nu2
dF2.dnu.plus <- cdf.f.dnu.plus
d2F2ddelta2dnu <- (dF2.dnu.plus - dF2) / fd.prec
d2F22ddelta2dnu <- d2F2ddelta2dnu - d2f2delta2nu
}
} else if (VC$margins[2]=="ZAP") {
i2 <- dat[,2]
if(univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
}
sigma <- plogis(sigma.star)
sigma <- ifelse(sigma<(1-precision), sigma, 1-sigma)
sigma <- ifelse(sigma>precision, sigma, precision)
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
f2 <- dZAP(i2, sigma = sigma, mu = mu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pZAP(i2, sigma = sigma, mu = mu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
df2 <- ifelse(i2>0, (mu^(-1 + i2) * (i2 - exp(mu) * i2 + exp(mu)* mu) * (-1 + sigma))/((-1 + exp(mu))^2 * gamma(1 + i2)), 0)
df2 <- as.vector(df2)*mu
df2.sigma <- ifelse(i2>0, -(exp(-mu) * mu^i2/(gamma(i2 + 1) * (1 - exp(-mu)))), 1)
if(univariate==FALSE) {
df2.mu.int.ZAP<-function(y, mu, sigma) {
df2 <- ifelse(y>0, (mu^(-1 + y) * (y - exp(mu) * y + exp(mu)* mu) * (-1 + sigma))/((-1 + exp(mu))^2 * gamma(1 + y)), 0)
df2 <- as.vector(df2)*mu
df2
}
df2.sigma.int.ZAP<-function(y, mu, sigma) {
df2 <- ifelse(y>0, -(exp(-mu) * mu^y/(gamma(y + 1) * (1 - exp(-mu)))), 1)
df2
}
ly <- length(i2)
cdf.f <- rep(0, ly)
cdf.f.sigma <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f[i] <- sum(df2.mu.int.ZAP(allval, mm, ms))
cdf.f.sigma[i] <- sum(df2.sigma.int.ZAP(allval, mm, ms))
}
dF2<-cdf.f
dF22<-dF2-df2
dF2.sigma <- cdf.f.sigma
dF22.sigma <- dF2.sigma-df2.sigma
}
# Hessian derivative components
mu.plus <- ifelse(exp(eta2+fd.prec)>0.0001, exp(eta2+fd.prec), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
df2.plus <- ifelse(i2>0, (mu.plus^(-1 + i2) * (i2 - exp(mu.plus) * i2 + exp(mu.plus)* mu.plus) *
(-1 + sigma))/((-1 + exp(mu.plus))^2 * gamma(1 + i2)), 0)
df2.plus <- as.vector(df2.plus)*mu.plus
d2f2delta22 <- (df2.plus - df2) / fd.prec
sigma.plus <- sigma+fd.prec
d2f2sigma2 <- 0
df2.dsigma.plus <- ifelse(i2>0, (mu^(-1 + i2) * (i2 - exp(mu) * i2 + exp(mu)* mu) *
(-1 + sigma.plus))/((-1 + exp(mu))^2 * gamma(1 + i2)), 0)
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*mu
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec
if(univariate==FALSE) {
cdf.f.plus <- rep(0, ly)
nmu <- rep(mu.plus, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f.plus[i] <- sum(df2.mu.int.ZAP(allval, mm, ms))
}
dF2.mu.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.mu.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
cdf.f.sigma.plus <- rep(0, ly)
cdf.f.dsigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f.sigma.plus[i] <- sum(df2.sigma.int.ZAP(allval, mm, ms))
cdf.f.dsigma.plus[i] <- sum(df2.mu.int.ZAP(allval, mm, ms))
}
dF2.sigma.plus <- cdf.f.sigma.plus
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
dF2.dsigma.plus <- cdf.f.dsigma.plus
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
}
} else if (VC$margins[2]=="ZIBB") {
i2 <- dat[,2]
if(univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
if (is.null(VC$X3)) nu.star <- params[(VC$X2.d2 + 2)]
if (!is.null(VC$X3)) nu.star <- etanu <- X4 %*% params[(VC$X2.d2 + VC$X3.d2 + 1):(VC$X2.d2 + VC$X3.d2 + VC$X4.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
if(is.null(VC$X3)) nu.star <- params[(VC$X1.d2+VC$X2.d2+2)]
if(!is.null(VC$X3)) nu.star <- etanu <- X4%*%params[(VC$X1.d2+VC$X2.d2+VC$X3.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2)]
}
nu <- plogis(nu.star)
nu <- ifelse(nu<precision, precision, nu)
nu <- ifelse(nu>(1-precision), 1-precision, nu)
sigma <- exp(sigma.star)
sigma <- ifelse(sigma>precision, sigma, precision)
sigma <- ifelse(sigma<10^7, sigma, 10^7)
mu <- plogis(eta2)
mu <- ifelse(mu<precision, precision, mu)
mu <- ifelse(mu>(1-precision), 1-precision, mu)
f2 <- dZIBB(i2, bd=bd, mu=mu, sigma=sigma, nu=nu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pZIBB(i2, bd=bd, mu=mu, sigma=sigma, nu=nu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
df2 <- ifelse(i2>0, ((-1 + nu)*gamma(1 + bd)*gamma(i2 + mu/sigma)*gamma(1/sigma)*
gamma(bd - (-1 + mu + i2*sigma)/sigma)*(-digamma(i2 + mu/sigma) -
digamma((1 - mu)/sigma) + digamma(mu/sigma) + digamma(bd - (-1 + mu + i2*sigma)/sigma)))/
(sigma*gamma(1 + bd - i2)*gamma(1 + i2)*gamma(bd + 1/sigma)*gamma((1 - mu)/sigma)*gamma(mu/sigma)),
((-1 + nu)*gamma(bd + (1 - mu)/sigma)*gamma(1/sigma)*(digamma(bd + (1 - mu)/sigma) -
digamma((1 - mu)/sigma)))/(sigma*gamma(bd + 1/sigma)*gamma((1 - mu)/sigma)))
df2 <- as.vector(df2) * mu * (1-mu)
df2.sigma <- ifelse(i2>0, -(((-1 + nu)*gamma(1 + bd)*gamma(i2 + mu/sigma)*gamma(1/sigma)*
gamma(bd - (-1 + mu + i2*sigma)/sigma)*(digamma(bd + 1/sigma) -
mu*digamma(i2 + mu/sigma) - digamma(1/sigma) +
digamma((1 - mu)/sigma) - mu*digamma((1 - mu)/sigma) +
mu*digamma(mu/sigma) - digamma(bd - (-1 + mu + i2*sigma)/sigma) +
mu*digamma(bd - (-1 + mu + i2*sigma)/sigma)))/
(sigma^2*gamma(1 + bd - i2)*gamma(1 + i2)*gamma(bd + 1/sigma)*
gamma((1 - mu)/sigma)*gamma(mu/sigma))),
-(((-1 + nu)*gamma(bd + (1 - mu)/sigma)*gamma(1/sigma)*
(digamma(bd + 1/sigma) + (-1 + mu)*digamma(bd + (1 - mu)/sigma) -
digamma(1/sigma) + digamma((1 - mu)/sigma) -
mu*digamma((1 - mu)/sigma)))/(sigma^2*gamma(bd + 1/sigma)*
gamma((1 - mu)/sigma))))
df2.nu <- ifelse(i2>0, -(gamma(bd + 1)/(gamma(i2 + 1) * gamma(bd - i2 + 1)) * (gamma(1/sigma) *
gamma(i2 + mu/sigma) * gamma(bd + (1 - mu)/sigma - i2))/(gamma(bd +
1/sigma) * gamma(mu/sigma) * gamma((1 - mu)/sigma))),
1 - (gamma(bd + 1)/(gamma(1) * gamma(bd + 1)) * (gamma(1/sigma) *
gamma(mu/sigma) * gamma(bd + (1 - mu)/sigma))/(gamma(bd + 1/sigma) *
gamma(mu/sigma) * gamma((1 - mu)/sigma))))
if(univariate==FALSE) {
df2.mu.int.ZIBB <- function(y, mu, sigma, nu, n) {
df2 <- ifelse(y>0, ((-1 + nu)*gamma(1 + n)*gamma(y + mu/sigma)*gamma(1/sigma)*
gamma(n - (-1 + mu + y*sigma)/sigma)*(-digamma(y + mu/sigma) -
digamma((1 - mu)/sigma) + digamma(mu/sigma) + digamma(n - (-1 + mu + y*sigma)/sigma)))/
(sigma*gamma(1 + n - y)*gamma(1 + y)*gamma(n + 1/sigma)*gamma((1 - mu)/sigma)*gamma(mu/sigma)),
((-1 + nu)*gamma(n + (1 - mu)/sigma)*gamma(1/sigma)*(digamma(n + (1 - mu)/sigma) -
digamma((1 - mu)/sigma)))/(sigma*gamma(n + 1/sigma)*gamma((1 - mu)/sigma)))
df2 <- as.vector(df2) * mu*(1-mu)
}
df2.sigma.int.ZIBB<-function(y, mu, sigma, nu, n) {
df2 <- ifelse(y>0, -(((-1 + nu)*gamma(1 + n)*gamma(y + mu/sigma)*gamma(1/sigma)*
gamma(n - (-1 + mu + y*sigma)/sigma)*(digamma(n + 1/sigma) -
mu*digamma(y + mu/sigma) - digamma(1/sigma) +
digamma((1 - mu)/sigma) - mu*digamma((1 - mu)/sigma) +
mu*digamma(mu/sigma) - digamma(n - (-1 + mu + y*sigma)/sigma) +
mu*digamma(n - (-1 + mu + y*sigma)/sigma)))/
(sigma^2*gamma(1 + n - y)*gamma(1 + y)*gamma(n + 1/sigma)*
gamma((1 - mu)/sigma)*gamma(mu/sigma))),
-(((-1 + nu)*gamma(n + (1 - mu)/sigma)*gamma(1/sigma)*
(digamma(n + 1/sigma) + (-1 + mu)*digamma(n + (1 - mu)/sigma) -
digamma(1/sigma) + digamma((1 - mu)/sigma) -
mu*digamma((1 - mu)/sigma)))/(sigma^2*gamma(n + 1/sigma)*
gamma((1 - mu)/sigma))))
df2
}
df2.nu.int.ZIBB <- function(y, mu, sigma, nu, n) {
df2 <- ifelse(y>0, -(gamma(n + 1)/(gamma(y + 1) * gamma(n - y + 1)) * (gamma(1/sigma) *
gamma(y + mu/sigma) * gamma(n + (1 - mu)/sigma - y))/(gamma(n +
1/sigma) * gamma(mu/sigma) * gamma((1 - mu)/sigma))), 1 - (gamma(n + 1)/(gamma(1) * gamma(n + 1)) * (gamma(1/sigma) *
gamma(mu/sigma) * gamma(n + (1 - mu)/sigma))/(gamma(n + 1/sigma) *
gamma(mu/sigma) * gamma((1 - mu)/sigma))))
df2
}
ly <- length(i2)
cdf.f <- rep(0, ly)
cdf.f.sigma <- rep(0, ly)
cdf.f.nu <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
cdf.f[i] <- sum(df2.mu.int.ZIBB(allval, mm, ms, mn, bd))
cdf.f.sigma[i] <- sum(df2.sigma.int.ZIBB(allval, mm, ms, mn, bd))
cdf.f.nu[i] <- sum(df2.nu.int.ZIBB(allval, mm, ms, mn, bd))
}
dF2<-cdf.f
dF22<-dF2-df2
dF2.sigma <- cdf.f.sigma
dF22.sigma <- dF2.sigma-df2.sigma
dF2.nu <- cdf.f.nu
dF22.nu <- dF2.nu-df2.nu
}
# Hessian derivative components
#delta2delta2
mu.plus <- plogis(eta2+fd.prec)
mu.plus <- ifelse(mu.plus<precision, precision, mu.plus)
mu.plus <- ifelse(mu.plus>(1-precision), 1-precision, mu.plus)
df2.plus <- ifelse(i2>0, ((-1 + nu)*gamma(1 + bd)*gamma(i2 + mu.plus/sigma)*gamma(1/sigma)*
gamma(bd - (-1 + mu.plus + i2*sigma)/sigma)*(-digamma(i2 + mu.plus/sigma) -
digamma((1 - mu.plus)/sigma) + digamma(mu.plus/sigma) +
digamma(bd - (-1 + mu.plus + i2*sigma)/sigma)))/(sigma*gamma(1 + bd - i2)*gamma(1 + i2)*gamma(bd + 1/sigma)*
gamma((1 - mu.plus)/sigma)*gamma(mu.plus/sigma)),
((-1 + nu)*gamma(bd + (1 - mu.plus)/sigma)*gamma(1/sigma)*(digamma(bd + (1 - mu.plus)/sigma) -
digamma((1 - mu.plus)/sigma)))/(sigma*gamma(bd + 1/sigma)*gamma((1 - mu.plus)/sigma)))
df2.plus <- as.vector(df2.plus) * mu.plus*(1-mu.plus)
d2f2delta22 <- (df2.plus - df2) / fd.prec
#sigma^2
sigma.plus <- sigma + fd.prec
df2.sigma.plus <- ifelse(i2>0, -(((-1 + nu)*gamma(1 + bd)*gamma(i2 + mu/sigma.plus)*gamma(1/sigma.plus)*
gamma(bd - (-1 + mu + i2*sigma.plus)/sigma.plus)*(digamma(bd + 1/sigma.plus) -
mu*digamma(i2 + mu/sigma.plus) - digamma(1/sigma.plus) +
digamma((1 - mu)/sigma.plus) - mu*digamma((1 - mu)/sigma.plus) +
mu*digamma(mu/sigma.plus) - digamma(bd - (-1 + mu + i2*sigma.plus)/sigma.plus) +
mu*digamma(bd - (-1 + mu + i2*sigma.plus)/sigma.plus)))/
(sigma.plus^2*gamma(1 + bd - i2)*gamma(1 + i2)*gamma(bd + 1/sigma.plus)*
gamma((1 - mu)/sigma.plus)*gamma(mu/sigma.plus))),
-(((-1 + nu)*gamma(bd + (1 - mu)/sigma.plus)*gamma(1/sigma.plus)*
(digamma(bd + 1/sigma.plus) + (-1 + mu)*digamma(bd + (1 - mu)/sigma.plus) -
digamma(1/sigma.plus) + digamma((1 - mu)/sigma.plus) -
mu*digamma((1 - mu)/sigma.plus)))/(sigma.plus^2*gamma(bd + 1/sigma.plus)*
gamma((1 - mu)/sigma.plus))))
d2f2sigma2 <- (df2.sigma.plus - df2.sigma) / fd.prec
#delta2sigma
df2.dsigma.plus <- ifelse(i2>0, ((-1 + nu)*gamma(1 + bd)*gamma(i2 + mu/sigma.plus)*gamma(1/sigma.plus)*
gamma(bd - (-1 + mu + i2*sigma.plus)/sigma.plus)*(-digamma(i2 + mu/sigma.plus) -
digamma((1 - mu)/sigma.plus) + digamma(mu/sigma.plus) +
digamma(bd - (-1 + mu + i2*sigma.plus)/sigma.plus)))/
(sigma.plus*gamma(1 + bd - i2)*gamma(1 + i2)*gamma(bd + 1/sigma.plus)*
gamma((1 - mu)/sigma.plus)*gamma(mu/sigma.plus)),
((-1 + nu)*gamma(bd + (1 - mu)/sigma.plus)*gamma(1/sigma.plus)*(digamma(bd + (1 - mu)/sigma.plus) -
digamma((1 - mu)/sigma.plus)))/(sigma.plus*gamma(bd + 1/sigma.plus)*gamma((1 - mu)/sigma.plus)))
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*(mu*(1-mu))
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec
#nu^2
nu.plus <- nu + fd.prec
nu.plus <- ifelse(nu.plus>=(1-precision), (1-precision), nu.plus)
df2.nu.plus <- ifelse(i2>0, -(gamma(bd + 1)/(gamma(i2 + 1) * gamma(bd - i2 + 1)) * (gamma(1/sigma) *
gamma(i2 + mu/sigma) * gamma(bd + (1 - mu)/sigma - i2))/(gamma(bd +
1/sigma) * gamma(mu/sigma) * gamma((1 - mu)/sigma))),
1 - (gamma(bd + 1)/(gamma(1) * gamma(bd + 1)) * (gamma(1/sigma) *
gamma(mu/sigma) * gamma(bd + (1 - mu)/sigma))/(gamma(bd + 1/sigma) *
gamma(mu/sigma) * gamma((1 - mu)/sigma))))
d2f2nu2 <- (df2.nu.plus - df2.nu) / fd.prec
#delta2nu
df2.dnu.plus <- ifelse(i2>0, ((-1 + nu.plus)*gamma(1 + bd)*gamma(i2 + mu/sigma)*gamma(1/sigma)*
gamma(bd - (-1 + mu + i2*sigma)/sigma)*(-digamma(i2 + mu/sigma) -
digamma((1 - mu)/sigma) + digamma(mu/sigma) +
digamma(bd - (-1 + mu + i2*sigma)/sigma)))/
(sigma*gamma(1 + bd - i2)*gamma(1 + i2)*gamma(bd + 1/sigma)*
gamma((1 - mu)/sigma)*gamma(mu/sigma)),
((-1 + nu.plus)*gamma(bd + (1 - mu)/sigma)*gamma(1/sigma)*(digamma(bd + (1 - mu)/sigma) -
digamma((1 - mu)/sigma)))/(sigma*gamma(bd + 1/sigma)*gamma((1 - mu)/sigma)))
df2.dnu.plus <- as.vector(df2.dnu.plus)*(mu*(1-mu))
d2f2delta2nu <- (df2.dnu.plus - df2) / fd.prec
#nusigma
df2.nu.sigma.plus <- ifelse(i2>0, -(gamma(bd + 1)/(gamma(i2 + 1) * gamma(bd - i2 + 1)) * (gamma(1/sigma.plus) *
gamma(i2 + mu/sigma.plus) * gamma(bd + (1 - mu)/sigma.plus - i2))/(gamma(bd +
1/sigma.plus) * gamma(mu/sigma.plus) * gamma((1 - mu)/sigma.plus))) ,
1 - (gamma(bd + 1)/(gamma(1) * gamma(bd + 1)) * (gamma(1/sigma.plus) *
gamma(mu/sigma.plus) * gamma(bd + (1 - mu)/sigma.plus))/(gamma(bd + 1/sigma.plus) *
gamma(mu/sigma.plus) * gamma((1 - mu)/sigma.plus))) )
d2f2sigmanu <- (df2.nu.sigma.plus - df2.nu) / fd.prec
if(univariate==FALSE) {
cdf.f.plus <- rep(0, ly)
nmu <- rep(mu.plus, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
pdfall <- df2.mu.int.ZIBB(allval, mm, ms, mn, bd)
cdf.f.plus[i] <- sum(pdfall)
}
dF2.mu.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.mu.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
cdf.f.sigma.plus <- rep(0, ly)
cdf.f.dsigma.plus <- rep(0, ly)
cdf.f.dnu.dsigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
cdf.f.sigma.plus[i] <- sum(df2.sigma.int.ZIBB(allval, mm, ms, mn, bd))
cdf.f.dsigma.plus[i] <- sum(df2.mu.int.ZIBB(allval, mm, ms, mn, bd))
cdf.f.dnu.dsigma.plus[i] <- sum(df2.nu.int.ZIBB(allval, mm, ms, mn, bd))
}
dF2.sigma.plus <- cdf.f.sigma.plus
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
dF2.dsigma.plus <- cdf.f.dsigma.plus
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
dF2.dnu.dsigma.plus <- cdf.f.dnu.dsigma.plus
d2F2dsigmadnu <- (dF2.dnu.dsigma.plus - dF2.nu) / fd.prec
d2F22dsigmadnu <- d2F2dsigmadnu - d2f2sigmanu
cdf.f.nu.plus <- rep(0, ly)
cdf.f.dnu.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
cdf.f.nu.plus[i] <- sum(df2.nu.int.ZIBB(allval, mm, ms, mn, bd))
cdf.f.dnu.plus[i] <- sum(df2.mu.int.ZIBB(allval, mm, ms, mn, bd))
}
dF2.nu.plus <- cdf.f.nu.plus
d2F2dnu2 <- (dF2.nu.plus - dF2.nu) / fd.prec
d2F22dnu2 <- d2F2dnu2 - d2f2nu2
dF2.dnu.plus <- cdf.f.dnu.plus
d2F2ddelta2dnu <- (dF2.dnu.plus - dF2) / fd.prec
d2F22ddelta2dnu <- d2F2ddelta2dnu - d2f2delta2nu
}
} else if (VC$margins[2]=="ZIBI") {
i2 <- dat[,2]
if(univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
}
sigma <- plogis(sigma.star)
sigma <- ifelse(sigma<(1-precision), sigma, 1-sigma)
sigma <- ifelse(sigma>precision, sigma, precision)
mu <- plogis(eta2)
mu <- ifelse(mu<(1-precision), mu, 1-mu)
mu <- ifelse(mu>precision, mu, precision)
f2 <- dZIBI(i2, sigma = sigma, mu = mu, bd=bd)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pZIBI(i2, sigma = sigma, mu = mu, bd=bd)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
df2<- ifelse(i2>0, ((1 - mu)^(-1 + bd - i2) * mu^(-1 + i2) * (-i2 + bd * mu) * (-1 + sigma) * gamma(1 + bd))/(gamma(1 + bd - i2)
* gamma(1 + i2)),
-((1 - sigma) * ((1 - mu)^(bd - 1) * bd)))
df2 <- as.vector(df2)*(1-mu)*mu
df2.sigma <- ifelse(i2>0, -(gamma(bd + 1) * mu^i2 * (1 - mu)^(bd - i2)/(gamma(i2 + 1) * gamma(bd - i2 + 1))), 1 - (1 - mu)^bd)
if(univariate==FALSE) {
df2.mu.int.ZIBI<-function(y, mu, sigma, n) {
df2<- ifelse(y>0, ((1 - mu)^(-1 + n - y) * mu^(-1 + y) * (-y + n * mu) * (-1 + sigma) * gamma(1 + n))/(gamma(1 + n - y) *
gamma(1 + y)),
-((1 - sigma) * ((1 - mu)^(n - 1) * n)))
df2 <- as.vector(df2)*(1-mu)*mu
df2
}
df2.sigma.int.ZIBI<-function(y, mu, sigma, n) {
df2<- ifelse(y>0, -(gamma(n + 1) * mu^y * (1 - mu)^(n - y)/(gamma(y + 1) * gamma(n -
y + 1))), 1 - (1 - mu)^n)
df2
}
ly <- length(i2)
cdf.f <- rep(0, ly)
cdf.f.sigma <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f[i] <- sum(df2.mu.int.ZIBI(allval, mm, ms, bd))
cdf.f.sigma[i] <- sum(df2.sigma.int.ZIBI(allval, mm, ms, bd))
}
dF2<-cdf.f
dF22<-dF2-df2
dF2.sigma <- cdf.f.sigma
dF22.sigma <- dF2.sigma-df2.sigma
}
# Hessian derivative components
mu.plus <- ifelse(plogis(eta2+fd.prec)>0.0001, plogis(eta2+fd.prec), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
df2.plus <- ifelse(i2>0, ((1 - mu.plus)^(-1 + bd - i2) * mu.plus^(-1 + i2) * (-i2 + bd * mu.plus) * (-1 + sigma) *
gamma(1 + bd))/(gamma(1 + bd - i2) * gamma(1 + i2)),
-((1 - sigma) * ((1 - mu.plus)^(bd - 1) * bd)))
df2.plus <- as.vector(df2.plus)*(1-mu.plus)*mu.plus
d2f2delta22 <- (df2.plus - df2) / fd.prec
sigma.plus <- sigma+fd.prec
d2f2sigma2 <- 0
df2.dsigma.plus <- ifelse(i2>0, ((1 - mu)^(-1 + bd - i2) * mu^(-1 + i2) * (-i2 + bd * mu) * (-1 + sigma.plus) *
gamma(1 + bd))/(gamma(1 + bd - i2) * gamma(1 + i2)),
-((1 - sigma.plus) * ((1 - mu)^(bd - 1) * bd)))
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*(1-mu)*mu
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec
if(univariate==FALSE) {
cdf.f.plus <- rep(0, ly)
nmu <- rep(mu.plus, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f.plus[i] <- sum(df2.mu.int.ZIBI(allval, mm, ms, bd))
}
dF2.mu.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.mu.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
cdf.f.sigma.plus <- rep(0, ly)
cdf.f.dsigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f.sigma.plus[i] <- sum(df2.sigma.int.ZIBI(allval, mm, ms, bd))
cdf.f.dsigma.plus[i] <- sum(df2.mu.int.ZIBI(allval, mm, ms, bd))
}
dF2.sigma.plus <- cdf.f.sigma.plus
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
dF2.dsigma.plus <- cdf.f.dsigma.plus
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
}
} else if (VC$margins[2]=="ZINBI") {
i2 <- dat[,2]
if(univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
if (is.null(VC$X3)) nu.star <- params[(VC$X2.d2 + 2)]
if (!is.null(VC$X3)) nu.star <- etanu <- X4 %*% params[(VC$X2.d2 + VC$X3.d2 + 1):(VC$X2.d2 + VC$X3.d2 + VC$X4.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
if(is.null(VC$X3)) nu.star <- params[(VC$X1.d2+VC$X2.d2+2)]
if(!is.null(VC$X3)) nu.star <- etanu <- X4%*%params[(VC$X1.d2+VC$X2.d2+VC$X3.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2)]
}
nu <- plogis(nu.star)
nu <- ifelse(nu<precision, precision, nu)
nu <- ifelse(nu>(1-precision), 1-precision, nu)
sigma <- exp(sigma.star)
sigma <- ifelse(sigma>precision, sigma, precision)
sigma <- ifelse(sigma<10^7, sigma, 10^7)
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
f2 <- dZINBI(i2, mu=mu, sigma=sigma, nu=nu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pZINBI(i2, mu=mu, sigma=sigma, nu=nu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
df2 <- ifelse(i2>0, ((-i2 + mu) * (-1 + nu) * (1/(1 + mu * sigma))^(1 + 1/sigma) * ((mu * sigma)/(1 + mu * sigma))^i2 *
gamma(i2 + 1/sigma))/(mu * gamma(1 + i2) * gamma(1/sigma)),
(-1 + nu) * ((1/(1 + mu * sigma))^((1/sigma) - 1) * ((1/sigma) *
(sigma/(1 + mu * sigma)^2))) )
df2 <- as.vector(df2) * mu
df2.sigma <- ifelse(i2>0, ((-1 + nu) * (1/(1 + mu * sigma))^(1 + 1/sigma) * ((mu * sigma)/(1 + mu * sigma))^
i2 * gamma(i2 + 1/sigma) * (-i2 * sigma + mu * sigma + log(1/(1 + mu * sigma)) +
mu * sigma * log(1/(1 + mu * sigma)) + (1 + mu * sigma) * digamma(
i2 + 1/sigma) - (1 + mu * sigma) * digamma(1/
sigma)))/(sigma^2 * gamma(1 + i2) * gamma(1/sigma)),
((-1 + nu) * (1/(1 + mu * sigma))^(1 + 1/sigma) * (mu * sigma + (1 + mu * sigma) * log(1/(1 + mu * sigma))))/sigma^2)
df2.nu <- ifelse(i2>0, -((1/(1 + sigma * mu))^(1/sigma) * (mu * sigma/(1 + mu * sigma))^i2 *
gamma(1/sigma + i2)/(gamma(1/sigma) * gamma(1 + i2))),
1 - (1/(1 + mu * sigma))^(1/sigma))
if(univariate==FALSE) {
df2.mu.int.ZINBI <- function(y, mu, sigma, nu) {
df2 <- ifelse(y>0, -(((-1 + nu) * (1/(1 + sigma * mu))^(1/sigma) * ((mu * sigma/(1 +
mu * sigma))^(y - 1) * (y * (sigma/(1 + mu * sigma) - mu *
sigma * sigma/(1 + mu * sigma)^2))) - (-1 + nu) * ((1/(1 +
sigma * mu))^((1/sigma) - 1) * ((1/sigma) * (sigma/(1 + sigma *
mu)^2))) * (mu * sigma/(1 + mu * sigma))^y) * gamma(1/sigma +
y)/(gamma(1/sigma) * gamma(1 + y))),
(-1 + nu) * ((1/(1 + mu * sigma))^((1/sigma) - 1) * ((1/sigma) *
(sigma/(1 + mu * sigma)^2))) )
df2 <- as.vector(df2) * mu
}
df2.sigma.int.ZINBI<-function(y, mu, sigma, nu) {
df2 <- ifelse(y>0, ((-1 + nu) * (1/(1 + mu * sigma))^(1 + 1/sigma) * ((mu * sigma)/(1 + mu * sigma))^
y * gamma(y + 1/sigma) * (-y * sigma + mu * sigma + log(1/(1 + mu * sigma)) +
mu * sigma * log(1/(1 + mu * sigma)) + (1 + mu * sigma) * digamma(
y + 1/sigma) - (1 + mu * sigma) * digamma(1/
sigma)))/(sigma^2 * gamma(1 + y) * gamma(1/sigma)),
((-1 + nu) * (1/(1 + mu * sigma))^(1 + 1/sigma) * (mu * sigma + (1 + mu * sigma) *
log(1/(1 + mu * sigma))))/sigma^2)
df2
}
df2.nu.int.ZINBI <- function(y, mu, sigma, nu) {
df2 <- ifelse(y>0, -((1/(1 + sigma * mu))^(1/sigma) * (mu * sigma/(1 + mu * sigma))^y *
gamma(1/sigma + y)/(gamma(1/sigma) * gamma(1 + y))),
1 - (1/(1 + mu * sigma))^(1/sigma))
df2
}
ly <- length(i2)
cdf.f <- rep(0, ly)
cdf.f.sigma <- rep(0, ly)
cdf.f.nu <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
cdf.f[i] <- sum(df2.mu.int.ZINBI(allval, mm, ms, mn))
cdf.f.sigma[i] <- sum(df2.sigma.int.ZINBI(allval, mm, ms, mn))
cdf.f.nu[i] <- sum(df2.nu.int.ZINBI(allval, mm, ms, mn))
}
dF2<-cdf.f
dF22<-dF2-df2
dF2.sigma <- cdf.f.sigma
dF22.sigma <- dF2.sigma-df2.sigma
dF2.nu <- cdf.f.nu
dF22.nu <- dF2.nu-df2.nu
}
# Hessian derivative components
#fd.prec2 <- 10^-3
fd.prec2 <- sqrt(eps)
#delta2delta2
mu.plus <- ifelse(exp(eta2+fd.prec)>0.0001, exp(eta2+fd.prec), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
df2.plus <- ifelse(i2>0, ((-i2 + mu.plus) * (-1 + nu) * (1/(1 + mu.plus * sigma))^(1 + 1/sigma) *
((mu.plus * sigma)/(1 + mu.plus * sigma))^i2 * gamma(i2 + 1/sigma))/(mu.plus * gamma(1 + i2) * gamma(1/sigma)),
(-1 + nu) * ((1/(1 + mu.plus * sigma))^((1/sigma) - 1) * ((1/sigma) *
(sigma/(1 + mu.plus * sigma)^2))) )
df2.plus <- as.vector(df2.plus) * mu.plus
d2f2delta22 <- (df2.plus - df2) / fd.prec
#sigma^2
sigma.plus <- sigma + fd.prec
df2.sigma.plus <- ifelse(i2>0, ((-1 + nu) * (1/(1 + mu * sigma.plus))^(1 + 1/sigma.plus) * ((mu * sigma.plus)/(1 + mu *
sigma.plus))^i2 * gamma(i2 + 1/sigma.plus) * (-i2 * sigma.plus + mu * sigma.plus + log(1/(1 + mu * sigma.plus)) +
mu * sigma.plus * log(1/(1 + mu * sigma.plus)) + (1 + mu * sigma.plus) * digamma(i2 + 1/sigma.plus) -
(1 + mu * sigma.plus) * digamma(1/sigma.plus)))/(sigma.plus^2 * gamma(1 + i2) * gamma(1/sigma.plus)),
((-1 + nu) * (1/(1 + mu * sigma.plus))^(1 + 1/sigma.plus) * (mu * sigma.plus + (1 + mu * sigma.plus) *
log(1/(1 + mu * sigma.plus))))/sigma.plus^2)
d2f2sigma2 <- (df2.sigma.plus - df2.sigma) / fd.prec
#delta2sigma
df2.dsigma.plus <- ifelse(i2>0, ((-i2 + mu) * (-1 + nu) * (1/(1 + mu * sigma.plus))^(1 + 1/sigma.plus) *
((mu * sigma.plus)/(1 + mu * sigma.plus))^i2 * gamma(i2 + 1/sigma.plus))/(mu * gamma(1 + i2) * gamma(1/sigma.plus)),
(-1 + nu) * ((1/(1 + mu * sigma.plus))^((1/sigma.plus) - 1) * ((1/sigma.plus) *
(sigma.plus/(1 + mu * sigma.plus)^2))) )
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*mu
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec
#nu^2
nu.plus <- nu + fd.prec
nu.plus <- ifelse(nu.plus>=(1-precision), (1-precision), nu.plus)
d2f2nu2 <- 0
#delta2nu
df2.dnu.plus <- ifelse(i2>0, ((-i2 + mu) * (-1 + nu.plus) * (1/(1 + mu * sigma))^(1 + 1/sigma) *
((mu * sigma)/(1 + mu * sigma))^i2 * gamma(i2 + 1/sigma))/(mu * gamma(1 + i2) * gamma(1/sigma)),
(-1 + nu.plus) * ((1/(1 + mu * sigma))^((1/sigma) - 1) * ((1/sigma) *
(sigma/(1 + mu * sigma)^2))) )
df2.dnu.plus <- as.vector(df2.dnu.plus)*mu
d2f2delta2nu <- (df2.dnu.plus - df2) / fd.prec
#nusigma
df2.nu.sigma.plus <- ifelse(i2>0, -((1/(1 + sigma.plus * mu))^(1/sigma.plus) * (mu * sigma.plus/(1 + mu * sigma.plus))^i2 *
gamma(1/sigma.plus + i2)/(gamma(1/sigma.plus) * gamma(1 + i2))),
1 - (1/(1 + mu * sigma.plus))^(1/sigma.plus))
d2f2sigmanu <- (df2.nu.sigma.plus - df2.nu) / fd.prec
if(univariate==FALSE) {
cdf.f.plus <- rep(0, ly)
nmu <- rep(mu.plus, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
pdfall <- df2.mu.int.ZINBI(allval, mm, ms, mn)
cdf.f.plus[i] <- sum(pdfall)
}
dF2.mu.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.mu.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
cdf.f.sigma.plus <- rep(0, ly)
cdf.f.dsigma.plus <- rep(0, ly)
cdf.f.dnu.dsigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
nnu <- rep(nu, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
cdf.f.sigma.plus[i] <- sum(df2.sigma.int.ZINBI(allval, mm, ms, mn))
cdf.f.dsigma.plus[i] <- sum(df2.mu.int.ZINBI(allval, mm, ms, mn))
cdf.f.dnu.dsigma.plus[i] <- sum(df2.nu.int.ZINBI(allval, mm, ms, mn))
}
dF2.sigma.plus <- cdf.f.sigma.plus
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
dF2.dsigma.plus <- cdf.f.dsigma.plus
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
dF2.dnu.dsigma.plus <- cdf.f.dnu.dsigma.plus
d2F2dsigmadnu <- (dF2.dnu.dsigma.plus - dF2.nu) / fd.prec
d2F22dsigmadnu <- d2F2dsigmadnu - d2f2sigmanu
cdf.f.nu.plus <- rep(0, ly)
cdf.f.dnu.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
nnu <- rep(nu.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
mn<- nnu[i]
allval <- seq(0, y.y)
cdf.f.nu.plus[i] <- sum(df2.nu.int.ZINBI(allval, mm, ms, mn))
cdf.f.dnu.plus[i] <- sum(df2.mu.int.ZINBI(allval, mm, ms, mn))
}
dF2.nu.plus <- cdf.f.nu.plus
d2F2dnu2 <- (dF2.nu.plus - dF2.nu) / fd.prec
d2F22dnu2 <- d2F2dnu2 - d2f2nu2
dF2.dnu.plus <- cdf.f.dnu.plus
d2F2ddelta2dnu <- (dF2.dnu.plus - dF2) / fd.prec
d2F22ddelta2dnu <- d2F2ddelta2dnu - d2f2delta2nu
}
} else if (VC$margins[2]=="ZIP") {
i2 <- dat[,2]
if(univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
}
sigma <- plogis(sigma.star)
sigma <- ifelse(sigma<(1-precision), sigma, 1-sigma)
sigma <- ifelse(sigma>precision, sigma, precision)
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
f2 <- dZIP(i2, sigma = sigma, mu = mu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pZIP(i2, sigma = sigma, mu = mu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
df2 <- ifelse(i2>0, (1 - sigma) * (mu^((i2) - 1) * (i2))/gamma(i2 + 1) * exp(-mu) -
(1 - sigma) * mu^(i2)/gamma(i2 + 1) * exp(-mu), -((1 - sigma) * exp(-mu)))
df2 <- as.vector(df2)*mu
df2.sigma <- ifelse(i2>0, -(mu^(i2)/gamma(i2 + 1) * exp(-mu)), 1 - exp(-mu))
if(univariate==FALSE) {
df2.mu.int.ZIP<-function(y, mu, sigma) {
df2 <- ifelse(y>0, (1 - sigma) * (mu^((y) - 1) * (y))/gamma(y + 1) * exp(-mu) -
(1 - sigma) * mu^(y)/gamma(y + 1) * exp(-mu), -((1 - sigma) * exp(-mu)))
df2 <- as.vector(df2)*mu
df2
}
df2.sigma.int.ZIP<-function(y, mu, sigma) {
df2 <- ifelse(y>0, -(mu^(y)/gamma(y + 1) * exp(-mu)), 1 - exp(-mu))
df2
}
ly <- length(i2)
cdf.f <- rep(0, ly)
cdf.f.sigma <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f[i] <- sum(df2.mu.int.ZIP(allval, mm, ms))
cdf.f.sigma[i] <- sum(df2.sigma.int.ZIP(allval, mm, ms))
}
dF2<-cdf.f
dF22<-dF2-df2
dF2.sigma <- cdf.f.sigma
dF22.sigma <- dF2.sigma-df2.sigma
}
# Hessian derivative components
mu.plus <- ifelse(exp(eta2+fd.prec)>0.0001, exp(eta2+fd.prec), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
df2.plus <- ifelse(i2>0, (1 - sigma) * (mu.plus^((i2) - 1) * (i2))/gamma(i2 + 1) * exp(-mu.plus) -
(1 - sigma) * mu.plus^(i2)/gamma(i2 + 1) * exp(-mu.plus), -((1 - sigma) * exp(-mu.plus)))
df2.plus <- as.vector(df2.plus)*mu.plus
d2f2delta22 <- (df2.plus - df2) / fd.prec
sigma.plus <- sigma + fd.prec
d2f2sigma2 <- 0
df2.dsigma.plus <- ifelse(i2>0, (1 - sigma.plus) * (mu^((i2) - 1) * (i2))/gamma(i2 + 1) * exp(-mu) -
(1 - sigma.plus) * mu^(i2)/gamma(i2 + 1) * exp(-mu), -((1 - sigma.plus) * exp(-mu)))
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*mu
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec
if(univariate==FALSE) {
cdf.f.plus <- rep(0, ly)
nmu <- rep(mu.plus, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f.plus[i] <- sum(df2.mu.int.ZIP(allval, mm, ms))
}
dF2.mu.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.mu.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
cdf.f.sigma.plus <- rep(0, ly)
cdf.f.dsigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f.sigma.plus[i] <- sum(df2.sigma.int.ZIP(allval, mm, ms))
cdf.f.dsigma.plus[i] <- sum(df2.mu.int.ZIP(allval, mm, ms))
}
dF2.sigma.plus <- cdf.f.sigma.plus
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
dF2.dsigma.plus <- cdf.f.dsigma.plus
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
}
} else if (VC$margins[2]=="ZIP2") {
i2 <- dat[,2]
if(univariate==TRUE) {
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
}
sigma <- plogis(sigma.star)
sigma <- ifelse(sigma<(1-precision), sigma, 1-sigma)
sigma <- ifelse(sigma>precision, sigma, precision)
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
f2 <- dZIP2(i2, sigma = sigma, mu = mu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pZIP2(i2, sigma = sigma, mu = mu)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F2 <- ifelse(F2>precision, F2, precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
df2 <- ifelse(i2>0, -(exp(mu/(-1 + sigma))*mu^(-1 + i2)*(mu + i2*(-1 + sigma)))/((1 - sigma)^i2*gamma(1 + i2)),
-((1 - sigma) * (exp(-(mu/(1 - sigma))) * (1/(1 - sigma)))) )
df2 <- as.vector(df2)*mu
df2.sigma <- ifelse(i2>0, -(exp(mu/(-1 + sigma))*mu^i2*(1 - sigma)^(-1 - i2)*(1 + mu + i2*(-1 + sigma) - sigma))/gamma(1 + i2),
1 - ((1 - sigma) * (exp(-(mu/(1 - sigma))) * (mu/(1 - sigma)^2)) +
exp(-(mu/(1 - sigma)))))
if(univariate==FALSE) {
df2.mu.int.ZIP2<-function(y, mu, sigma) {
df2 <- ifelse(y>0, -(exp(mu/(-1 + sigma))*mu^(-1 + y)*(mu + y*(-1 + sigma)))/((1 - sigma)^y*gamma(1 + y)),
-((1 - sigma) * (exp(-(mu/(1 - sigma))) * (1/(1 - sigma)))) )
df2 <- as.vector(df2)*mu
df2
}
df2.sigma.int.ZIP2<-function(y, mu, sigma) {
df2 <- ifelse(y>0, -(exp(mu/(-1 + sigma))*mu^y*(1 - sigma)^(-1 - y)*(1 + mu + y*(-1 + sigma) - sigma))/gamma(1 + y),
1 - ((1 - sigma) * (exp(-(mu/(1 - sigma))) * (mu/(1 - sigma)^2)) +
exp(-(mu/(1 - sigma)))))
df2
}
ly <- length(i2)
cdf.f <- rep(0, ly)
cdf.f.sigma <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f[i] <- sum(df2.mu.int.ZIP2(allval, mm, ms))
cdf.f.sigma[i] <- sum(df2.sigma.int.ZIP2(allval, mm, ms))
}
dF2<-cdf.f
dF22<-dF2-df2
dF2.sigma <- cdf.f.sigma
dF22.sigma <- dF2.sigma-df2.sigma
}
# Hessian derivative components
mu.plus <- ifelse(exp(eta2+fd.prec)>0.0001, exp(eta2+fd.prec), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
df2.plus <- ifelse(i2>0, -(exp(mu.plus/(-1 + sigma))*mu.plus^(-1 + i2)*(mu.plus + i2*(-1 + sigma)))/((1 - sigma)^i2*gamma(1 + i2)),
-((1 - sigma) * (exp(-(mu.plus/(1 - sigma))) * (1/(1 - sigma)))) )
df2.plus <- as.vector(df2.plus)*mu.plus
d2f2delta22 <- (df2.plus - df2) / fd.prec
sigma.plus <- sigma + fd.prec
df2.sigma.plus <- ifelse(i2>0, -(exp(mu/(-1 + sigma.plus))*mu^i2*(1 - sigma.plus)^(-1 - i2)*(1 + mu + i2*(-1 + sigma.plus) - sigma.plus))/gamma(1 + i2),
1 - ((1 - sigma.plus) * (exp(-(mu/(1 - sigma.plus))) * (mu/(1 - sigma.plus)^2)) + exp(-(mu/(1 - sigma.plus)))))
d2f2sigma2 <- (df2.sigma.plus - df2.sigma) / fd.prec
df2.dsigma.plus <- ifelse(i2>0, -(exp(mu/(-1 + sigma.plus))*mu^(-1 + i2)*(mu + i2*(-1 + sigma.plus)))/((1 - sigma.plus)^i2*gamma(1 + i2)),
-((1 - sigma.plus) * (exp(-(mu/(1 - sigma.plus))) * (1/(1 - sigma.plus)))) )
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*mu
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec
if(univariate==FALSE) {
cdf.f.plus <- rep(0, ly)
nmu <- rep(mu.plus, length = ly)
nsigma <- rep(sigma, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f.plus[i] <- sum(df2.mu.int.ZIP2(allval, mm, ms))
}
dF2.mu.plus <- cdf.f.plus
d2F2ddelta22 <- (dF2.mu.plus - dF2) / fd.prec
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
cdf.f.sigma.plus <- rep(0, ly)
cdf.f.dsigma.plus <- rep(0, ly)
nmu <- rep(mu, length = ly)
nsigma <- rep(sigma.plus, length = ly)
for (i in 1:ly) {
y.y <- i2[i]
mm <- nmu[i]
ms <- nsigma[i]
allval <- seq(0, y.y)
cdf.f.sigma.plus[i] <- sum(df2.sigma.int.ZIP2(allval, mm, ms))
cdf.f.dsigma.plus[i] <- sum(df2.mu.int.ZIP2(allval, mm, ms))
}
dF2.sigma.plus <- cdf.f.sigma.plus
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
dF2.dsigma.plus <- cdf.f.dsigma.plus
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
}
} else if (VC$margins[2]=="ZIPIG") {
i2 <- dat[,2]
if(univariate==TRUE){
if (is.null(VC$X3)) sigma.star <- params[(VC$X2.d2 + 1)]
if (!is.null(VC$X3)) sigma.star <- etasqv <- X3 %*% params[(VC$X2.d2 + 1):(VC$X2.d2 + VC$X3.d2)]
if (is.null(VC$X3)) nu.star <- params[(VC$X2.d2 + 2)]
if (!is.null(VC$X3)) nu.star <- etanu <- X4 %*% params[(VC$X2.d2 + VC$X3.d2 + 1):(VC$X2.d2 + VC$X3.d2 + VC$X4.d2)]
} else {
if(is.null(VC$X3)) sigma.star <- params[(VC$X1.d2+VC$X2.d2+1)]
if(!is.null(VC$X3)) sigma.star <- etasqv <- X3%*%params[(VC$X1.d2+VC$X2.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2)]
if(is.null(VC$X3)) nu.star <- params[(VC$X1.d2+VC$X2.d2+2)]
if(!is.null(VC$X3)) nu.star <- etanu <- X4%*%params[(VC$X1.d2+VC$X2.d2+VC$X3.d2+1):(VC$X1.d2+VC$X2.d2+VC$X3.d2+VC$X4.d2)]
}
sigma <- exp(sigma.star)
sigma <- ifelse(sigma>precision, sigma, precision)
sigma <- ifelse(sigma<10^7, sigma, 10^7)
mu <- ifelse(exp(eta2)>0.0001, exp(eta2), 0.0001)
mu <- ifelse(mu<10^7, mu, 10^7)
nu <- plogis(nu.star)
nu <- ifelse(nu<precision, precision, nu)
nu <- ifelse(nu>(1-precision), 1-precision, nu)
f2 <- dZIPIG(i2, mu=mu, sigma=sigma, nu=nu)
f2 <- ifelse(f2>precision, f2, precision)
if(univariate==FALSE) {
F2 <- pZIPIG(i2, mu=mu, sigma=sigma, nu=nu)
F2 <- ifelse(F2>precision, F2, precision)
F2 <- ifelse(F2<(1-precision), F2, 1-precision)
F22 <- F2 - f2
F22 <- ifelse(F22<(1-precision), F22, 1-precision)
F22 <- ifelse(F22>precision, F22, precision)
}
f2.fd.mu <- dZIPIG(i2, mu = (mu +fd.prec), sigma = sigma, nu = nu)
f2.fd.mu <- ifelse(f2.fd.mu>precision, f2.fd.mu, precision)
df2 <- (f2.fd.mu - f2) / (fd.prec)
df2 <- as.vector(df2)*as.vector(mu)
f2.fd.sigma <- dZIPIG(i2, mu = mu, sigma = (sigma+fd.prec), nu = nu)
f2.fd.sigma <- ifelse(f2.fd.sigma>precision, f2.fd.sigma, precision)
df2.sigma <- (f2.fd.sigma - f2) / (fd.prec)
f2.fd.nu <- dZIPIG(i2, mu = mu, sigma = sigma, nu =(nu+fd.prec))
f2.fd.nu <- ifelse(f2.fd.nu>precision, f2.fd.nu, precision)
df2.nu <- (f2.fd.nu - f2) / (fd.prec)
if(univariate==FALSE) {
F2.fd.mu <- pZIPIG(i2, mu = (mu+fd.prec), sigma = sigma, nu = nu)
F2.fd.mu <- ifelse(F2.fd.mu<(1-precision), F2.fd.mu, 1-precision)
dF2 <- (F2.fd.mu - F2) / (fd.prec)
dF2 <- as.vector(dF2)*as.vector(mu)
dF22 <- dF2 - df2
F2.fd.sigma <- pZIPIG(i2, mu = mu, sigma = (sigma+fd.prec), nu = nu)
F2.fd.sigma <- ifelse(F2.fd.sigma<(1-precision), F2.fd.sigma, 1-precision)
dF2.sigma <- (F2.fd.sigma - F2) / (fd.prec)
dF22.sigma <- dF2.sigma - df2.sigma
F2.fd.nu <- pZIPIG(i2, mu = mu, sigma = sigma, nu = (nu+fd.prec))
F2.fd.nu <- ifelse(F2.fd.nu<(1-precision), F2.fd.nu, 1-precision)
dF2.nu <- (F2.fd.nu - F2) / (fd.prec)
dF22.nu <- dF2.nu - df2.nu
}
# Hessian derivative components
#second derivative fd precision
#fd.prec2 <- 10^-3
fd.prec2 <- sqrt(eps)
#pmfs
#delta2^2
mu.plus <- ifelse(exp(eta2+fd.prec2)>0.0001, exp(eta2+fd.prec2), 0.0001)
mu.plus <- ifelse(mu.plus<10^7, mu.plus, 10^7)
f2.plus <- dZIPIG(i2, mu = mu.plus, sigma = sigma, nu = nu)
f2.plus <- ifelse(f2.plus>precision, f2.plus, precision)
f2.fd.mu.plus <- dZIPIG(i2, mu = (mu.plus+fd.prec), sigma = sigma, nu = nu)
f2.fd.mu.plus <- ifelse(f2.fd.mu.plus>precision, f2.fd.mu.plus, precision)
df2.plus <- (f2.fd.mu.plus - f2.plus) / (fd.prec)
df2.plus <- as.vector(df2.plus)*as.vector(mu.plus)
df2.plus <- as.vector(df2.plus)
d2f2delta22 <- (df2.plus - df2) / fd.prec2
#sigma^2
sigma.plus <- sigma+fd.prec2
f2.plus.sigma <- dZIPIG(i2, mu = mu, sigma = sigma.plus, nu = nu)
f2.plus.sigma <- ifelse(f2.plus.sigma>precision, f2.plus.sigma, precision)
f2.fd.sigma.plus <- dZIPIG(i2, mu = mu, sigma = (sigma.plus+fd.prec), nu = nu)
f2.fd.sigma.plus <- ifelse(f2.fd.sigma.plus>precision, f2.fd.sigma.plus, precision)
df2.sigma.plus <- (f2.fd.sigma.plus - f2.plus.sigma) / (fd.prec)
d2f2sigma2 <- (df2.sigma.plus - df2.sigma) / fd.prec2
#delta2sigma
f2.dsigma.plus <- dZIPIG(i2, mu = mu, sigma = sigma.plus, nu = nu)
f2.dsigma.plus <- ifelse(f2.dsigma.plus>precision, f2.dsigma.plus, precision)
f2.fd.mu.sigma.plus <- dZIPIG(i2, mu = (mu+fd.prec), sigma = sigma.plus, nu = nu)
f2.fd.mu.sigma.plus <- ifelse(f2.fd.mu.sigma.plus>precision, f2.fd.mu.sigma.plus, precision)
df2.dsigma.plus <- (f2.fd.mu.sigma.plus - f2.dsigma.plus) / (fd.prec)
df2.dsigma.plus <- as.vector(df2.dsigma.plus)*as.vector(mu)
d2f2delta2sigma <- (df2.dsigma.plus - df2) / fd.prec2
#nu^2
nu.plus <- nu + fd.prec2
nu.plus <- ifelse(nu.plus>=(1-precision), (1-precision), nu.plus)
f2.plus.nu <- dZIPIG(i2, mu = mu, sigma = sigma, nu = nu.plus)
f2.plus.nu <- ifelse(f2.plus.nu>precision, f2.plus.nu, precision)
nu.plus2 <- nu.plus + fd.prec
nu.plus2 <- ifelse(nu.plus2>(1-precision), 1-precision, nu.plus2)
f2.fd.nu.plus <- dZIPIG(i2, mu = mu, sigma = sigma, nu = nu.plus2)
f2.fd.nu.plus <- ifelse(f2.fd.nu.plus>precision, f2.fd.nu.plus, precision)
df2.nu.plus <- (f2.fd.nu.plus - f2.plus.nu) / (fd.prec)
d2f2nu2 <- (df2.nu.plus - df2.nu) / fd.prec2
#delta2nu
f2.dnu.plus <- dZIPIG(i2, mu = mu, sigma = sigma, nu = nu.plus)
f2.dnu.plus <- ifelse(f2.dnu.plus>precision, f2.dnu.plus, precision)
f2.fd.mu.nu.plus <- dZIPIG(i2, mu = (mu+fd.prec), sigma = sigma, nu = nu.plus)
f2.fd.mu.nu.plus <- ifelse(f2.fd.mu.nu.plus>precision, f2.fd.mu.nu.plus, precision)
df2.dnu.plus <- (f2.fd.mu.nu.plus - f2.dnu.plus) / (fd.prec)
df2.dnu.plus <- as.vector(df2.dnu.plus)*as.vector(mu)
d2f2delta2nu <-(df2.dnu.plus - df2) / fd.prec2
#nusigma
f2.fd.sigma.nu.plus <- dZIPIG(i2, mu = mu, sigma = (sigma+fd.prec), nu = nu.plus)
f2.fd.sigma.nu.plus <- ifelse(f2.fd.sigma.nu.plus>precision, f2.fd.sigma.nu.plus, precision)
df2.dnu.sigma.plus <- (f2.fd.sigma.nu.plus - f2.plus.nu) / (fd.prec)
df2.dnu.sigma.plus <- as.vector(df2.dnu.sigma.plus)
d2f2sigmanu <- (df2.dnu.sigma.plus - df2.sigma) / fd.prec2
if(univariate==FALSE) {
#cdfs
#delta2^2
F2.plus <- pZIPIG(i2, mu = mu.plus, sigma = sigma, nu = nu)
F2.fd.mu.plus <- pZIPIG(i2, mu = (mu.plus+fd.prec), sigma = sigma, nu = nu)
F2.fd.mu.plus <- ifelse(F2.fd.mu.plus<(1-precision), F2.fd.mu.plus, 1-precision)
dF2.plus <- (F2.fd.mu.plus - F2.plus) / (fd.prec)
dF2.plus <- as.vector(dF2.plus)*as.vector(mu.plus)
d2F2ddelta22 <-(dF2.plus - dF2) / fd.prec2
d2F22ddelta22 <- d2F2ddelta22 - d2f2delta22
#sigma^2
F2.sigma.plus <- pZIPIG(i2, mu = mu, sigma = sigma.plus, nu = nu)
F2.fd.sigma.plus <- pZIPIG(i2, mu = mu, sigma = (sigma.plus+fd.prec), nu=nu)
F2.fd.sigma.plus <- ifelse(F2.fd.sigma.plus<(1-precision), F2.fd.sigma.plus, 1-precision)
dF2.sigma.plus <- (F2.fd.sigma.plus - F2.sigma.plus) / (fd.prec)
d2F2dsigma2 <- (dF2.sigma.plus - dF2.sigma) / fd.prec2
d2F22dsigma2 <- d2F2dsigma2 - d2f2sigma2
#delta2sigma
F2.fd.dsigma.plus <- pZIPIG(i2, mu = (mu+fd.prec), sigma = sigma.plus, nu)
F2.fd.dsigma.plus <- ifelse(F2.fd.dsigma.plus<(1-precision), F2.fd.dsigma.plus, 1-precision)
dF2.dsigma.plus <- (F2.fd.dsigma.plus - F2.sigma.plus) / (fd.prec)
dF2.dsigma.plus <- as.vector(dF2.dsigma.plus)*as.vector(mu)
d2F2ddelta2dsigma <- (dF2.dsigma.plus - dF2) / fd.prec2
d2F22ddelta2dsigma <- d2F2ddelta2dsigma - d2f2delta2sigma
#nu^2
F2.nu.plus<- pZIPIG(i2, mu = mu, sigma = sigma, nu = nu.plus)
F2.fd.nu.plus <- pZIPIG(i2, mu = mu, sigma = sigma, nu = nu.plus2)
F2.fd.nu.plus <- ifelse(F2.fd.nu.plus<(1-precision), F2.fd.nu.plus, 1-precision)
dF2.nu.plus <- (F2.fd.nu.plus - F2.nu.plus) / (fd.prec)
d2F2dnu2 <- (dF2.nu.plus - dF2.nu) / fd.prec2
d2F22dnu2 <- d2F2dnu2 - d2f2nu2
#delta2nu
F2.fd.dnu.plus <- pZIPIG(i2, mu = (mu+fd.prec), sigma = sigma, nu.plus)
F2.fd.dnu.plus <- ifelse(F2.fd.dnu.plus<(1-precision), F2.fd.dnu.plus, 1-precision)
dF2.dnu.plus <- (F2.fd.dnu.plus - F2.nu.plus) / (fd.prec)
dF2.dnu.plus <- as.vector(dF2.dnu.plus)*as.vector(mu)
d2F2ddelta2dnu <- (dF2.dnu.plus - dF2) / fd.prec2
d2F22ddelta2dnu <- d2F2ddelta2dnu - d2f2delta2nu
#sigmanu
F2.fd.dnu.sigma.plus <- pZIPIG(i2, mu = mu, sigma = (sigma+fd.prec), nu.plus)
F2.fd.dnu.sigma.plus <- ifelse(F2.fd.dnu.sigma.plus<(1-precision), F2.fd.dnu.sigma.plus, 1-precision)
dF2.dnu.sigma.plus <- (F2.fd.dnu.sigma.plus - F2.nu.plus) / (fd.prec)
dF2.dnu.sigma.plus <- as.vector(dF2.dnu.sigma.plus)
d2F2dsigmadnu <- (dF2.dnu.sigma.plus - dF2.sigma) / fd.prec2
d2F22dsigmadnu <- d2F2dsigmadnu - d2f2sigmanu
}
}
list(i0=i0, i1=i1, i2=i2, ind=ind, F1=F1, dF1=dF1, d2F1ddelta1delta1=d2F1ddelta1delta1, mu=mu, sigma=sigma, nu=nu,
F2=F2, f2=f2, F22=F22, df2=df2, df2.sigma=df2.sigma, df2.nu=df2.nu, dF2=dF2, dF22=dF22, dF2.sigma=dF2.sigma,
dF22.sigma=dF22.sigma, dF2.nu=dF2.nu, dF22.nu=dF22.nu, d2f2delta22=d2f2delta22, d2f2sigma2=d2f2sigma2,
d2f2delta2sigma=d2f2delta2sigma, d2f2nu2=d2f2nu2, d2f2delta2nu=d2f2delta2nu, d2f2nusigma=d2f2nusigma,
d2f2sigmanu=d2f2sigmanu, d2F2ddelta22=d2F2ddelta22, d2F22ddelta22=d2F22ddelta22, d2F2dsigma2=d2F2dsigma2,
d2F22dsigma2=d2F22dsigma2, d2F2ddelta2dsigma=d2F2ddelta2dsigma, d2F22ddelta2dsigma=d2F22ddelta2dsigma, d2F2dnu2=d2F2dnu2,
d2F22dnu2=d2F22dnu2, d2F2ddelta2dnu=d2F2ddelta2dnu, d2F22ddelta2dnu=d2F22ddelta2dnu, d2F2dnudsigma=d2F2dnudsigma,
d2F22dnudsigma=d2F22dnudsigma, d2F2dsigmadnu=d2F2dsigmadnu, d2F22dsigmadnu=d2F22dsigmadnu, etasqv=etasqv, etanu=etanu)
}
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.