Nothing
#' Negative log Likelihood functions for Poisson, negative binomial,
#' Delaporte, Poisson-inverse Gaussian and Poisson-beta distributions
#'
#' The negative log Likelihood functions for Poisson, negative binomial,
#' Delaporte, Poisson-inverse Gaussian and Poisson-beta distributions.
#' Mixing two distributions of the same kind and/or adding zero-inflation
#' allows to take characteristics of real data into account.
#' Additionally, one population and two population mixtures - with and
#' without zero-inflation - allow distribution fitting of the Poisson,
#' negative binomial, Delaporte, Poisson-inverse Gaussian and the
#' Poisson-beta distribution.
#'
#'
#' @details
#' Functions nlogL_pois, nlogL_nb, nlogL_del, nlogL_pig, nlogL_pb compute the negative
#' log-likelihood of Poisson, negative binomial, Poisson-inverse Gaussian and
#' the Poisson-beta distributions given the data.
#' Functions nlogL_pois2, nlogL_nb2, nlogL_del2, nlogL_pig2 and nlogL_pb2 compute the negative
#' log-likelihood values for a two population mixture of distributions whereas
#' nlogL_zipois, nlogL_zinb, nlogL_zidel, nlogL_zipig, nlogL_zipb compute the same for the
#' zero-inflated distributions. Furthermore, nlogL_zipois2, nlogL_zinb2, nlogL_zidel2,
#' nlogL_zipig2 and nlogL_zipb2 are for two population mixtures with zero-inflation.
#' @param data Vector containing the discrete observations
#' @param par.pois Scalar containing the lambda parameter
#' of the Poisson distribution
#' @param par.nb Vector of length 2, containing the size and the mu
#' parameter of the negative binomial distribution
#' @param par.del Vector of length 3, containing the mu, sigma and the nu
#' parameter of the Delaporte distribution
#' @param par.pig Vector of length 2, containing the mu and the sigma
#' parameter of the Poisson-inverse Gaussian distribution
#' @param par.pb Vector of length 3, containing the alpha, beta
#' and c parameter of the Poisson-beta distribution
#' @param par.pois2,par.nb2,par.del2,par.pig2,par.pb2 Vector containing the parameters
#' of the two mixing distributions. First entry represents the
#' fraction of the first distribution, followed by all parameters
#' of the first, then all of the second distribution.
#' @param par.zipois,par.zinb,par.zidel,par.zipig,par.zipb Vector containing the respective
#' zero-inflated distribution parameters. The additional first
#' entry is the inflation parameter for all cases.
#' @param par.zipois2,par.zinb2,par.zidel2,par.zipig2,par.zipb2 Parameters for the zero-inflated
#' two population model.
#'
#' @keywords likelihood negative binomial Poisson-beta
#'
#' @name nlogL
#' @importFrom stats dpois dnbinom rnorm rpois rnbinom
#' @importFrom gamlss.dist dPIG dZIPIG rPIG rZIPIG dDEL rDEL
#' @export
#' @examples
#' x <- rpois(100, 11)
#' nl1 <- nlogL_pois(x, 11)
#' nl2 <- nlogL_pois(x, 13)
nlogL_pois <- function(data, par.pois) {
if (par.pois <= 0) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
} else {
nl <- -sum(dpois(x = data, lambda = par.pois, log = TRUE))
if (is.infinite(nl))
return (nl_inf + (rnorm(1, 10000, 20) ^ 2))
else
return(nl)
}
}
#' @rdname nlogL
#' @export
#' @examples
#' x <- rnbinom(100, size = 13, mu = 9)
#' nl <- nlogL_nb(x, c(13, 9))
nlogL_nb <- function(data, par.nb) {
if (par.nb[1] <= 0 || par.nb[2] < 0) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
} else {
nl <- -sum(dnbinom(x = data, size = par.nb[1], mu = par.nb[2], log = TRUE))
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else
return(nl)
}
}
#' @rdname nlogL
#' @export
#' @examples
#' x <- gamlss.dist::rDEL(100, mu = 5, sigma = 0.2, nu= 0.5)
#' nl <- nlogL_del(x, c(5, 0.2, 0.5))
nlogL_del <- function(data, par.del) {
if (par.del[1] <= 0 || par.del[2] <= 0 || par.del[3] <= 0 || par.del[3] >= 1) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
} else {
nl <- -sum(dDEL(x = data, mu = par.del[1], sigma = par.del[2], nu = par.del[3], log = TRUE))
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else
return(nl)
}
}
#' @rdname nlogL
#' @export
#' @examples
#' x <- gamlss.dist::rPIG(100, mu = 5, sigma = 0.2)
#' nl <- nlogL_pig(x, c(5, 0.2))
nlogL_pig <- function(data, par.pig) {
if (par.pig[1] <= 0 || par.pig[2] < 0) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
} else {
nl <- -sum(dPIG(x = data, mu = par.pig[1], sigma = par.pig[2], log = TRUE))
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else
return(nl)
}
}
#' @rdname nlogL
#' @export
#' @examples
#' x <- rpb(n = 1000, alpha=5, beta= 3, c=20)
#' nl <- nlogL_pb(x, c(5, 3, 20))
nlogL_pb <- function(data, par.pb) {
if (par.pb[1] < 0 ||
par.pb[2] < 0 || par.pb[3] <= 0) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
} else {
nl <- -sum(dpb(x = data, alpha = par.pb[1], beta = par.pb[2], c = par.pb[3], log = TRUE))
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else
return(nl)
}
}
sum_2pop_terms <- function(t1, t2) {
b1 <- (t1)/log(10) -300
b3 <- (t2)/log(10) -300
b13 <- pmax(b1,b3)
b2 <- (t1)/log(10) +300
b4 <- (t2)/log(10) +300
b24 <- pmin(b2,b4)
b <- rowMeans(matrix(c(b1*is.finite(b1),b2*is.finite(b2),b3*is.finite(b3),b4*is.finite(b4)),ncol = 4, byrow= FALSE), na.rm = TRUE )
t1_b <- (t1 / log(10)-b)
t2_b <- (t2 / log(10)-b)
t_b_check_1<- (t2_b - t1_b > 600)
t_b_check_2<- (t1_b - t2_b > 600)
nl <- sum( b[!t_b_check_1 & !t_b_check_2]*log(10) + log( 10 ^(t1[!t_b_check_1 & ! t_b_check_2] / log(10)-b[!t_b_check_1 & !t_b_check_2]) + 10 ^(t2[ !t_b_check_1 & !t_b_check_2] / log(10)-b[!t_b_check_1 & !t_b_check_2])))+
sum( (t2[t_b_check_1] )) + sum( (t1[t_b_check_2] ))
return(-nl)
}
#' @rdname nlogL
#' @export
#' @examples
#' s <- sample(x = c(0,1), size = 100, replace = TRUE, prob = c(0.3,0.7))
#' x <- s*rpois(100, 7) + (1-s)*rpois(100, 13)
#' nl <- nlogL_pois2(x, c(0.3, 13, 7))
nlogL_pois2 <- function(data, par.pois2) {
if (par.pois2[2] <= 0 ||
par.pois2[3] <= 0 ||
par.pois2[1] < 0 ||
par.pois2[1] > 1) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else if (par.pois2[1] == 0) {
new.par.pois2 <- c(1, par.pois2[c(3, 2)])
return(nlogL_pois2(data, new.par.pois2))
}
else {
t1 <- log(par.pois2[1])+dpois(x = data, lambda = par.pois2[2], log = TRUE)
t2 <- log(1-par.pois2[1])+dpois(x = data, lambda = par.pois2[3], log = TRUE)
nl <- sum_2pop_terms(t1, t2)
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else{
return(nl)
}
}
}
#' @rdname nlogL
#' @export
#' @examples
#' s <- sample(x = c(0,1), size = 100, replace = TRUE, prob = c(0.3,0.7))
#' x <- s*rnbinom(100, size = 13, mu = 9) + (1-s)*rnbinom(100, size = 17, mu = 29)
#' nl <- nlogL_nb2(x, c(0.3, 17, 29, 13, 9))
nlogL_nb2 <- function(data, par.nb2) {
if (par.nb2[2] <= 0 ||
par.nb2[3] < 0 ||
par.nb2[4] <= 0 ||
par.nb2[5] < 0 ||
par.nb2[1] < 0 ||
par.nb2[1] > 1) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else if (par.nb2[1] == 0) {
new.par.nb2 <- c(1, par.nb2[c(4, 5, 2, 3)])
return(nlogL_nb2(data, new.par.nb2))
}
else {
t1 <- log(par.nb2[1]) + dnbinom(x = data,size = par.nb2[2],mu = par.nb2[3],log = TRUE)
t2 <- log(1-par.nb2[1]) + dnbinom(x = data,size = par.nb2[4],mu = par.nb2[5],log = TRUE)
nl <- sum_2pop_terms(t1, t2)
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else
return(nl)
}
}
#' @rdname nlogL
#' @export
#' @examples
#' s <- sample(x = c(0,1), size = 100, replace = TRUE, prob = c(0.3,0.7))
#' x <- s * gamlss.dist::rDEL(100, mu = 5, sigma = 0.2, nu = 0.5) +
#' (1 - s) * gamlss.dist::rDEL(100, mu = 20, sigma = 2, nu = 0.1)
#' nl <- nlogL_del2(x, c(0.7,5, 0.2, 20, 2))
nlogL_del2 <- function(data, par.del2) {
if (par.del2[2] <= 0 ||
par.del2[3] <= 0 ||
par.del2[4] <= 0 ||
par.del2[4] >= 1 ||
par.del2[5] <= 0 ||
par.del2[6] <= 0 ||
par.del2[7] <= 0 ||
par.del2[7] >= 1 ||
par.del2[1] < 0 ||
par.del2[1] > 1) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else if (par.del2[1] == 0) {
new.par.del2 <- c(1, par.del2[c(5,6,7, 2, 3, 4)])
return(nlogL_del2(data, new.par.del2))
}
else {
t1 <- log(par.del2[1]) + dDEL(x = data, mu = par.del2[2], sigma = par.del2[3], nu = par.del2[4], log = TRUE)
t2 <- log(1 - par.del2[1]) + dDEL(x = data, mu = par.del2[5], sigma = par.del2[6], nu = par.del2[7], log = TRUE)
nl <- sum_2pop_terms(t1, t2)
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else
return(nl)
}
}
#' @rdname nlogL
#' @export
#' @examples
#' s <- sample(x = c(0,1), size = 100, replace = TRUE, prob = c(0.3,0.7))
#' x <- s * gamlss.dist::rPIG(100, mu = 5, sigma = 0.2) +
#' (1 - s) * gamlss.dist::rPIG(100, mu = 20, sigma = 2)
#' nl <- nlogL_pig2(x, c(0.7, 5, 0.2, 20, 2))
nlogL_pig2 <- function(data, par.pig2) {
if (par.pig2[2] < 0 ||
par.pig2[3] < 0 ||
par.pig2[4] < 0 ||
par.pig2[5] < 0 ||
par.pig2[1] < 0 ||
par.pig2[1] > 1) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else if (par.pig2[1] == 0) {
new.par.pig2 <- c(1, par.pig2[c(4, 5, 2, 3)])
return(nlogL_pig2(data, new.par.pig2))
}
else {
t1 <- log(par.pig2[1]) + dPIG(x = data,mu = par.pig2[2],sigma = par.pig2[3],log = TRUE)
t2 <- log(1-par.pig2[1]) + dPIG(x = data,mu = par.pig2[4],sigma = par.pig2[5],log = TRUE)
nl <- sum_2pop_terms(t1, t2)
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else
return(nl)
}
}
#' @rdname nlogL
#' @export
#' @examples
#' s <- sample(x = c(0,1), size = 100, replace = TRUE, prob = c(0.3,0.7))
#' x <- s*rpb(100, 5, 3, 20) + (1-s)*rpb(100, 7, 13, 53)
#' nl <- nlogL_pb2(x, c(0.7, 7, 13, 53, 5, 3, 20))
nlogL_pb2 <- function(data, par.pb2) {
if (par.pb2[2] < 0 ||
par.pb2[3] < 0 ||
par.pb2[4] <= 0 ||
par.pb2[5] < 0 ||
par.pb2[6] < 0 ||
par.pb2[7] <= 0 ||
par.pb2[1] < 0 ||
par.pb2[1] > 1) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else if (par.pb2[1] == 0) {
new.par.pb2 <- c(1, par.pb2[c(5:7, 2:4)])
return(nlogL_pb2(data, new.par.pb2))
}
else {
t1 <- log(par.pb2[1]) + dpb(x = data, alpha = par.pb2[2], beta = par.pb2[3], c = par.pb2[4], log = TRUE)
t2 <- log(1-par.pb2[1]) + dpb(x = data, alpha = par.pb2[5], beta = par.pb2[6], c = par.pb2[7], log = TRUE)
nl <- sum_2pop_terms(t1, t2)
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else{
return(nl)
}
}
}
#' @rdname nlogL
#' @export
#' @examples
#' x <- c(rep(0, 10), rpois(90, 7))
#' nl <- nlogL_zipois(x, c(0.1, 7))
nlogL_zipois <- function(data, par.zipois) {
if (par.zipois[2] <= 0 ||
par.zipois[1] < 0 ||
par.zipois[1] > 1) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else {
n <- length(data)
n0 <- length(which(data == 0))
non_zero <- data[which(data != 0)]
if(n0 == 0) {
t1 = 0
}
else {
l <- log(par.zipois[1] + (1 - par.zipois[1]) * exp(-par.zipois[2]))
# This can happen when the inflation parameter is 0 and the probability
# of a 0 is very low. For an explanation, see this:
#>>>> > exp(-1259)
#>>>> [1] 0
if(is.nan(l)){
t1 = -par.zipois[2]
}
else {
t1 <- l
}
}
if(n != n0){
nl <- n0 * t1 + (n - n0) * log(1 - par.zipois[1]) + sum(dpois(x = non_zero, lambda = par.zipois[2], log = TRUE))
} else {nl <- n0 * t1}
nl <- -nl
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else{
return(nl)
}
}
}
#' @rdname nlogL
#' @export
#' @examples
#' x <- c(rep(0,10), rnbinom(90, size = 13, mu = 9))
#' nl <- nlogL_zinb(x, c(0.1, 13, 9))
nlogL_zinb <- function(data, par.zinb) {
if (par.zinb[2] <= 0 ||
par.zinb[3] < 0 ||
par.zinb[1] < 0 ||
par.zinb[1] > 1) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else {
n <- length(data)
n0 <- length(which(data == 0))
non_zero <- data[which(data != 0)]
if(n0 == 0) {
t1 = 0
}
else {
l <- log(par.zinb[1] + (1 - par.zinb[1])*dnbinom(0, size = par.zinb[2], mu = par.zinb[3]))
# Motivated by conclusions from zipois, but the same scenario
# occurs in less cases here. The nestorowa dataset had convergent fits
# for all cases without this condition.
if(is.nan(l)){
t1 = dnbinom(0, size = par.zinb[2], mu = par.zinb[3], log = TRUE)
}
else {
t1 <- l
}
}
if(n != n0){
nl <- n0 * t1 + (n-n0)*log(1-par.zinb[1])+sum(dnbinom(x = non_zero, size = par.zinb[2], mu = par.zinb[3], log = TRUE))
} else {nl <- n0 * t1}
nl <- -nl
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else{
return(nl)
}
}
}
#' @rdname nlogL
#' @export
#' @examples
#' x <- c(rep(0,10), gamlss.dist::rDEL(90, mu = 13, sigma = 2, nu = 0.5))
#' nl <- nlogL_zidel(x, c(0.1, 13, 2, 0.5))
nlogL_zidel <- function(data, par.zidel) {
if (par.zidel[2] <= 0 ||
par.zidel[3] <= 0 ||
par.zidel[4] <= 0 ||
par.zidel[4] >= 1 ||
par.zidel[1] < 0 ||
par.zidel[1] > 1) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else {
n <- length(data)
n0 <- length(which(data == 0))
non_zero <- data[which(data != 0)]
if(n0 == 0) {
t1 = 0
}
else {
l <- log(par.zidel[1] + (1 - par.zidel[1]) * dDEL(0, mu = par.zidel[2], sigma = par.zidel[3], nu = par.zidel[4]))
# Motivated by conclusions from zipois, but the same scenario
# occurs in less cases here. The nestorowa dataset had convergent fits
# for all cases without this condition.
if(is.nan(l)){
t1 = dDEL(x = 0, mu = par.zidel[2], sigma = par.zidel[3], nu = par.zidel[4], log = TRUE)
}
else {
t1 <- l
}
}
if(n != n0){
nl <- n0 * t1 + (n - n0) * log(1 - par.zidel[1]) + sum(dDEL(x = non_zero, mu = par.zidel[2], sigma = par.zidel[3], nu = par.zidel[4], log = TRUE))
} else {nl<- n0 * t1}
nl <- -nl
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else{
return(nl)
}
}
}
#' @rdname nlogL
#' @export
#' @examples
#' x <- c(rep(0,10), gamlss.dist::rPIG(90, mu = 13, sigma = 2))
#' nl <- nlogL_zipig(x, c(0.1, 13, 2))
nlogL_zipig <- function(data, par.zipig) {
if (par.zipig[2] < 0 ||
par.zipig[3] < 0 ||
par.zipig[1] < 0 ||
par.zipig[1] > 1) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
} else {
n0 <- length(which(data == 0))
if(par.zipig[1]==0) {
nl= -sum(dPIG(x = data, mu = par.zipig[2], sigma = par.zipig[3], log = TRUE))
}
else {
nl <- -sum(dZIPIG(data, mu = par.zipig[2], sigma = par.zipig[3], nu = par.zipig[1], log = TRUE))
}
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else
return(nl)
}
}
#' @rdname nlogL
#' @export
#' @examples
#' x <- c(rep(0, 10), rpb(n = 90, alpha=5, beta= 3, c=20))
#' nl <- nlogL_zipb(x, c(0.1, 5, 3, 20))
nlogL_zipb <- function(data, par.zipb) {
if (par.zipb[2] < 0 ||
par.zipb[3] < 0 ||
par.zipb[4] <= 0 ||
par.zipb[1] < 0 ||
par.zipb[1] > 1) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else {
n <- length(data)
n0 <- length(which(data == 0))
non_zero <- data[which(data != 0)]
if(n0 == 0) {
t1 = 0
}
else{
l <- log(par.zipb[1] + (1 - par.zipb[1])*dpb(0, par.zipb[2], par.zipb[3], par.zipb[4]))
if(is.nan(l)){
# Going by the change from zipois to zinb [more parameters => better fit?],
# maybe even lesser cases happen here. Inference not tested.
t1 = dpb(0, par.zipb[2], par.zipb[3], par.zipb[4], log = TRUE)
}
else {
t1 <- l
}
}
if(n != n0){
nl <- n0 * t1 + (n-n0)*log(1-par.zipb[1])+sum(dpb(x = non_zero, par.zipb[2], par.zipb[3], par.zipb[4], log = TRUE))
} else { nl <- n0 * t1}
nl <- -nl
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else{
return(nl)
}
}
}
#' @rdname nlogL
#' @export
#' @examples
#' s <- sample(x = c(0, 1), size = 90, replace = TRUE, prob = c(0.3, 0.7))
#' x <- c(rep(0, 10), s * rpois(90, 7) + (1 - s) * rpois(90, 13))
#' nl1 <- nlogL_zipois2(x, c(0.1, 0.63, 7, 13))
nlogL_zipois2 <- function(data, par.zipois2) {
if (par.zipois2[1] < 0 ||
par.zipois2[1] > 1 ||
par.zipois2[2] < 0 ||
par.zipois2[2] > 1 ||
par.zipois2[1] + par.zipois2[2] > 1 ||
par.zipois2[3] <= 0 ||
par.zipois2[4] <= 0) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else if (par.zipois2[2] == 0) {
new.par.zipois2 <- c(par.zipois2[1], 1- par.zipois2[1], par.zipois2[c(4, 3)])
return(nlogL_zipois2(data, new.par.zipois2))
}
else {
n <- length(data)
n0 <- length(which(data == 0))
non_zero <- data[which(data != 0)]
t1 <- log(par.zipois2[2]) + dpois(x = 0, lambda = par.zipois2[3], log = TRUE)
t2 <- log(1-(par.zipois2[1] + par.zipois2[2]))+dpois(x = 0, lambda = par.zipois2[4], log = TRUE)
nl_zero <- -sum_2pop_terms(t1, t2)
if(n != n0){
t1 <- log(par.zipois2[2]) + dpois(x = non_zero, lambda = par.zipois2[3], log = TRUE)
t2 <- log(1-(par.zipois2[1] + par.zipois2[2]))+dpois(x = non_zero, lambda = par.zipois2[4], log = TRUE)
nl_non_zero <- -sum_2pop_terms(t1, t2)
} else {nl_non_zero <- 0}
# reduce expression when no zero-inflation
if(par.zipois2[1] == 0)
nl <- n0 * nl_zero + nl_non_zero
else
nl <- n0 * log(par.zipois2[1] + exp(nl_zero ) ) + nl_non_zero
nl <- -nl
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else{
return(nl)
}
}
}
#' @rdname nlogL
#' @export
#' @examples
#' s <- sample(x = c(0, 1), size = 90, replace = TRUE, prob = c(0.3, 0.7))
#' x <- c(rep(0, 10), s * rnbinom(90, size = 13, mu = 9) + (1 - s) * rnbinom(90, size = 17, mu = 29))
#' nl <- nlogL_zinb2(x, c(0.1, 0.63, 13, 9, 17, 29))
nlogL_zinb2 <- function(data, par.zinb2) {
if (par.zinb2[1] < 0 ||
par.zinb2[1] > 1 ||
par.zinb2[2] < 0 ||
par.zinb2[2] > 1 ||
par.zinb2[1] + par.zinb2[2] > 1 ||
par.zinb2[3] <= 0 ||
par.zinb2[4] < 0 ||
par.zinb2[5] <= 0 ||
par.zinb2[6] < 0) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else if (par.zinb2[2] == 0) {
new.par.zinb2 <- c(par.zinb2[1], 1 - par.zinb2[1], par.zinb2[c(5, 6, 3, 4)])
return(nlogL_zinb2(data, new.par.zinb2))
}
else {
n <- length(data)
n0 <- length(which(data == 0))
non_zero <- data[which(data != 0)]
t1 <- log(par.zinb2[2]) + dnbinom(x = 0, size = par.zinb2[3], mu = par.zinb2[4], log = TRUE)
t2 <- log(1 - (par.zinb2[1] + par.zinb2[2])) + dnbinom(x = 0, size = par.zinb2[5], mu = par.zinb2[6], log = TRUE)
nl_zero <- -sum_2pop_terms(t1, t2)
if(n != n0){
t1 <- log(par.zinb2[2]) + dnbinom(x = non_zero, size = par.zinb2[3], mu = par.zinb2[4], log = TRUE)
t2 <- log(1 - (par.zinb2[1] + par.zinb2[2])) + dnbinom(x = non_zero, size = par.zinb2[5], mu = par.zinb2[6], log = TRUE)
nl_non_zero <- -sum_2pop_terms(t1, t2)
} else {nl_non_zero <- 0}
if(par.zinb2[1] == 0)
nl <- n0 * nl_zero + nl_non_zero
else
nl <- n0 * log(par.zinb2[1] + exp(nl_zero ) ) + nl_non_zero
nl <- -nl
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else{
return(nl)
}
}
}
#' @rdname nlogL
#' @export
#' @examples
#' s <- sample(x = c(0, 1), size = 90, replace = TRUE, prob = c(0.3, 0.7))
#' x <- c(rep(0, 10), s * gamlss.dist::rDEL(90, mu = 13, sigma = 9, nu = 0.5) +
#' (1 - s) * gamlss.dist::rDEL(90, mu = 17, sigma = 29, nu = 0.1))
#' nl <- nlogL_zidel2(x, c(0.1, 0.63, 13, 9, 0.5, 17, 29, 0.1))
nlogL_zidel2 <- function(data, par.zidel2) {
if (par.zidel2[1] < 0 ||
par.zidel2[1] > 1 ||
par.zidel2[2] < 0 ||
par.zidel2[2] > 1 ||
par.zidel2[1] + par.zidel2[2] > 1 ||
par.zidel2[3] <= 0 ||
par.zidel2[4] <= 0 ||
par.zidel2[5] <= 0 ||
par.zidel2[5] >= 1 ||
par.zidel2[6] <= 0 ||
par.zidel2[7] <= 0 ||
par.zidel2[8] <= 0 ||
par.zidel2[8] >= 1 ) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else if (par.zidel2[2] == 0) {
new.par.zidel2 <- c(par.zidel2[1], 1 - par.zidel2[1], par.zidel2[c(6, 7, 8, 3, 4, 5)])
return(nlogL_zidel2(data, new.par.zidel2))
}
else {
n <- length(data)
n0 <- length(which(data == 0))
non_zero <- data[which(data != 0)]
t1 <- log(par.zidel2[2]) + dDEL(x = 0, mu = par.zidel2[3], sigma = par.zidel2[4], nu = par.zidel2[5], log = TRUE) #dnbinom(x = 0, size = par.zinb2[3], mu = par.zinb2[4], log = TRUE)
t2 <- log(1 - (par.zidel2[1] + par.zidel2[2])) + dDEL(x = 0, mu = par.zidel2[6], sigma = par.zidel2[7], nu = par.zidel2[8], log = TRUE) #dnbinom(x = 0, size = par.zinb2[5], mu = par.zinb2[6], log = TRUE)
nl_zero <- -sum_2pop_terms(t1, t2)
if(n != n0){
t1 <- log(par.zidel2[2]) + dDEL(x = non_zero, mu = par.zidel2[3], sigma = par.zidel2[4], nu = par.zidel2[5], log = TRUE) #dnbinom(x = non_zero, size = par.zinb2[3], mu = par.zinb2[4], log = TRUE)
t2 <- log(1 - (par.zidel2[1] + par.zidel2[2])) + dDEL(x = non_zero, mu = par.zidel2[6], sigma = par.zidel2[7], nu = par.zidel2[8], log = TRUE) #dnbinom(x = non_zero, size = par.zinb2[5], mu = par.zinb2[6], log = TRUE)
nl_non_zero <- -sum_2pop_terms(t1, t2)
} else {nl_non_zero <- 0}
if(par.zidel2[1] == 0)
nl <- n0 * nl_zero + nl_non_zero
else
nl <- n0 * log(par.zidel2[1] + exp(nl_zero ) ) + nl_non_zero
nl <- -nl
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else{
return(nl)
}
}
}
#' @rdname nlogL
#' @export
#' @examples
#' s <- sample(x = c(0,1), size = 90, replace = TRUE, prob = c(0.3, 0.7))
#' x <- c(rep(0, 10), s * gamlss.dist::rPIG(90, mu = 13, sigma = 0.2) +
#' (1-s) * gamlss.dist::rPIG(90, mu = 17, sigma = 2))
#' nl <- nlogL_zipig2(x, c(0.1, 0.63, 13, 0.2, 17, 2))
nlogL_zipig2 <- function(data, par.zipig2) {
if (par.zipig2[1] < 0 ||
par.zipig2[1] > 1 ||
par.zipig2[2] < 0 ||
par.zipig2[2] > 1 ||
par.zipig2[1] + par.zipig2[2] > 1 ||
par.zipig2[3] < 0 ||
par.zipig2[4] < 0 ||
par.zipig2[5] < 0 ||
par.zipig2[6] < 0) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else if (par.zipig2[2] == 0) {
new.par.zipig2 <- c(par.zipig2[1], 1 - par.zipig2[1], par.zipig2[c(5, 6, 3, 4)])
return(nlogL_zipig2(data, new.par.zipig2))
}
else {
n <- length(data)
n0 <- length(which(data == 0))
non_zero <- data[which(data != 0)]
t1 <- log(par.zipig2[2]) + dPIG(x = 0, mu = par.zipig2[3], sigma = par.zipig2[4], log = TRUE)
t2 <- log(1 - (par.zipig2[1] + par.zipig2[2])) + dPIG(x = 0, mu = par.zipig2[5], sigma = par.zipig2[6], log = TRUE)
nl_zero <- -sum_2pop_terms(t1, t2)
if(n != n0){
t1 <- log(par.zipig2[2]) + dPIG(x = non_zero, mu = par.zipig2[3], sigma = par.zipig2[4], log = TRUE)
t2 <- log(1 - (par.zipig2[1] + par.zipig2[2])) + dPIG(x = non_zero, mu = par.zipig2[5], sigma = par.zipig2[6], log = TRUE)
nl_non_zero <- -sum_2pop_terms(t1, t2)
} else {nl_non_zero <- 0}
if(par.zipig2[1] == 0)
nl <- n0 * nl_zero + nl_non_zero
else
nl <- n0 * log(par.zipig2[1] + exp(nl_zero ) ) + nl_non_zero
nl <- -nl
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else{
return(nl)
}
}
}
#' @rdname nlogL
#' @export
#' @examples
#' s <- sample(x = c(0,1), size = 90, replace = TRUE, prob = c(0.3,0.7))
#' x <- c(rep(0,10), s*rpb(90, 5, 3, 20) + (1-s)*rpb(90, 7, 13, 53))
#' nl <- nlogL_zipb2(x, c(0.1, 0.63, 7, 13, 53, 5, 3, 20))
nlogL_zipb2 <- function(data, par.zipb2) {
if (par.zipb2[1] < 0 ||
par.zipb2[1] > 1 ||
par.zipb2[2] < 0 ||
par.zipb2[2] > 1 ||
par.zipb2[1] + par.zipb2[2] > 1 ||
par.zipb2[3] < 0 ||
par.zipb2[4] < 0 ||
par.zipb2[5] <= 0 ||
par.zipb2[6] < 0 ||
par.zipb2[7] < 0 ||
par.zipb2[8] <= 0) {
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
}
else if (par.zipb2[2] == 0) {
new.par.zipb2 <- c(par.zipb2[1], 1 - par.zipb2[1], par.zipb2[c(6:8, 3:5)])
return(nlogL_zipb2(data, new.par.zipb2))
}
else {
n <- length(data)
n0 <- length(which(data == 0))
non_zero <- data[which(data != 0)]
t1 <- log(par.zipb2[2]) + dpb(x = 0, alpha = par.zipb2[3], beta = par.zipb2[4], c = par.zipb2[5], log = TRUE)
t2 <- log(1 - (par.zipb2[1] + par.zipb2[2])) + dpb(x = 0, alpha = par.zipb2[6], beta = par.zipb2[7], c = par.zipb2[8], log = TRUE)
nl_zero <- -sum_2pop_terms(t1, t2)
if(n != n0){
t1 <- log(par.zipb2[2]) + dpb(x = non_zero, alpha = par.zipb2[3], beta = par.zipb2[4], c = par.zipb2[5], log = TRUE)
t2 <- log(1 - (par.zipb2[1] + par.zipb2[2])) + dpb(x = non_zero, alpha = par.zipb2[6], beta = par.zipb2[7], c = par.zipb2[8], log = TRUE)
nl_non_zero <- -sum_2pop_terms(t1, t2)
} else {nl_non_zero <- 0}
if(par.zipb2[1] == 0)
nl <- n0 * nl_zero + nl_non_zero
else
nl <- n0 * log(par.zipb2[1] + exp(nl_zero ) ) + nl_non_zero
nl <- -nl
if (is.infinite(nl))
return(nl_inf + (rnorm(1, 10000, 20) ^ 2))
else{
return(nl)
}
}
}
nl_inf <- 1e+100
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.