Package 'HelpersMG'

Description

Contains miscellaneous functions useful for managing
'NetCDF' files (see http://en.wikipedia.org/wiki/NetCDF),
get tide levels on any point of the globe,
get moon phase and time for sun rise and fall,
analyse and reconstruct daily time series of temperature
with irregular sinusoidal pattern,
show scales and wind rose in plot with change of color of text,
Metropolis-Hastings algorithm for Bayesian MCMC analysis,
plot graphs or boxplot with error bars,
search files in disk by there names or their content,
read the contents of all files from a folder at one time,
calculate IC50 for ecotoxicological studies,
calculate the probability mass function of the sum of negative binomial
distributions, calculate distribution of unobserved values in censored or truncated distributions.
The latest version of this package can always be installed using:
install.packages("https://hebergement.universite-paris-saclay.fr/marcgirondot/CRAN/HelpersMG.tar.gz", repos=NULL, type="source")

Details

Helpers functions for several packages

Author(s)

Marc Girondot [email protected]

Examples

## Not run: 
library(HelpersMG)
print('----------------------------------')
print('Examples for mcmcComposite objects')
print('----------------------------------')
require(coda)
x <- rnorm(30, 10, 2)
dnormx <- function(x, par) return(-sum(dnorm(x, mean=par['mean'], sd=par['sd'], log=TRUE)))
parameters_mcmc <- data.frame(Density=c('dnorm', 'dlnorm'), 
Prior1=c(10, 0.5), Prior2=c(2, 0.5), SDProp=c(0.35, 0.2), 
Min=c(-3, 0), Max=c(100, 10), Init=c(10, 2), stringsAsFactors = FALSE, 
row.names=c('mean', 'sd'))
mcmc_run <- MHalgoGen(n.iter=100000, parameters=parameters_mcmc, data=x, 
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=1)
plot(mcmc_run, xlim=c(0, 20))
plot(mcmc_run, xlim=c(0, 10), parameters="sd")
mcmcforcoda <- as.mcmc(mcmc_run)
# Optimal rejection rate should be 0.234
rejectionRate(mcmcforcoda)
heidel.diag(mcmcforcoda)
raftery.diag(mcmcforcoda)
autocorr.diag(mcmcforcoda)
acf(mcmcforcoda[[1]][,"mean"], lag.max=20, bty="n", las=1)
acf(mcmcforcoda[[1]][,"sd"], lag.max=20, bty="n", las=1)
batchSE(mcmcforcoda, batchSize=100)
# The batch standard error procedure is usually thought to 
# be not as accurate as the time series methods used in summary
summary(mcmcforcoda)$statistics[,"Time-series SE"]
summary(mcmc_run)
as.parameters(mcmc_run)
lastp <- as.parameters(mcmc_run, index="last")
parameters_mcmc[,"Init"] <- lastp
# The n.adapt set to 1 is used to not record the first set of parameters
# then it is not duplicated (as it is also the last one for 
# the object mcmc_run)
mcmc_run2 <- MHalgoGen(n.iter=10000, parameters=parameters_mcmc, data=x, 
likelihood=dnormx, n.chains=1, n.adapt=1, thin=1, trace=1)
mcmc_run3 <- merge(mcmc_run, mcmc_run2)
####### no adaptation, n.adapt must be 0
parameters_mcmc[,"Init"] <- c(mean(x), sd(x))
mcmc_run3 <- MHalgoGen(n.iter=10000, parameters=parameters_mcmc, data=x, 
likelihood=dnormx, n.chains=1, n.adapt=0, thin=1, trace=1)
print('------------------------------------------')
print('Examples for Daily patterns of temperature')
print('------------------------------------------')
# Generate a timeserie of time
time.obs <- NULL
for (i in 0:9) time.obs <- c(time.obs, c(0, 6, 12, 18)+i*24)
# For these time, generate a timeseries of temperatures
temp.obs <- rep(NA, length(time.obs))
temp.obs[3+(0:9)*4] <- rnorm(10, 25, 3)
temp.obs[1+(0:9)*4] <- rnorm(10, 10, 3)
for (i in 1:(length(time.obs)-1)) 
  if (is.na(temp.obs[i])) 
  temp.obs[i] <- mean(c(temp.obs[i-1], temp.obs[i+1]))
  if (is.na(temp.obs[length(time.obs)])) 
  temp.obs[length(time.obs)] <- temp.obs[length(time.obs)-1]/2
observed <- data.frame(time=time.obs, temperature=temp.obs)
# Search for the minimum and maximum values
r <- minmax.periodic(time.minmax.daily=c(Min=2, Max=15), 
observed=observed, period=24)

# Estimate all the temperatures for these values
t <- temperature.periodic(minmax=r)

plot_errbar(x=t[,"time"], y=t[,"temperature"],
errbar.y=ifelse(is.na(t[,"sd"]), 0, 2*t[,"sd"]),
type="l", las=1, bty="n", errbar.y.polygon = TRUE, 
xlab="hours", ylab="Temperatures", ylim=c(0, 35), 
errbar.y.polygon.list = list(col="grey"))

plot_add(x=t[,"time"], y=t[,"temperature"], type="l")

# How many times this package has been download
library(cranlogs)
HelpersMG <- cran_downloads("HelpersMG", from = "2015-04-07", 
                            to = Sys.Date() - 1) 
sum(HelpersMG$count)
plot(HelpersMG$date, HelpersMG$count, type="l", bty="n")

## End(Not run)

Description

Add a S3 class as first class to an object.

Usage

addS3Class(x, class = NULL)

Arguments

Details

addS3Class add a S3 class to an object

Value

The same object with the new class as first class

Author(s)

Marc Girondot [email protected]

Examples

print.CF <- function(x) {cat("print.CF ", x)}
result <- "Je suis donc je pense"
result <- addS3Class(result, class="CF")
class(result)
print(result)
result <- addS3Class(result, class=c("ECF", "OCF"))
class(result)
print(result)

Description

Take a mcmcComposite object and create a mcmc.list object to be used with coda package.

Usage

## S3 method for class 'mcmcComposite'
as.mcmc(x, ...)

Arguments

Details

as.mcmc Extract mcmc object from the result of phenology_MHmcmc to be used with coda package

Value

A mcmc.list object

Author(s)

Marc Girondot [email protected]

See Also

Other mcmcComposite functions: MHalgoGen(), as.parameters(), as.quantiles(), merge.mcmcComposite(), plot.PriorsmcmcComposite(), plot.mcmcComposite(), setPriors(), summary.mcmcComposite()

Examples

## Not run: 
library(HelpersMG)
require(coda)
x <- rnorm(30, 10, 2)
dnormx <- function(data, x) {
 data <- unlist(data)
 return(-sum(dnorm(data, mean=x['mean'], sd=x['sd'], log=TRUE)))
}
parameters_mcmc <- data.frame(Density=c('dnorm', 'dlnorm'), 
Prior1=c(10, 0.5), Prior2=c(2, 0.5), SDProp=c(1, 1), 
Min=c(-3, 0), Max=c(100, 10), Init=c(10, 2), stringsAsFactors = FALSE, 
row.names=c('mean', 'sd'))
mcmc_run <- MHalgoGen(n.iter=1000, parameters=parameters_mcmc, data=x, 
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=1)
plot(mcmc_run, xlim=c(0, 20))
plot(mcmc_run, xlim=c(0, 10), parameters="sd")
mcmcforcoda <- as.mcmc(mcmc_run)
#' heidel.diag(mcmcforcoda)
raftery.diag(mcmcforcoda)
autocorr.diag(mcmcforcoda)
acf(mcmcforcoda[[1]][,"mean"], lag.max=20, bty="n", las=1)
acf(mcmcforcoda[[1]][,"sd"], lag.max=20, bty="n", las=1)
batchSE(mcmcforcoda, batchSize=100)
# The batch standard error procedure is usually thought to 
# be not as accurate as the time series methods used in summary
summary(mcmcforcoda)$statistics[,"Time-series SE"]
summary(mcmc_run)
as.parameters(mcmc_run)
lastp <- as.parameters(mcmc_run, index="last")
parameters_mcmc[,"Init"] <- lastp
# The n.adapt set to 1 is used to not record the first set of parameters
# then it is not duplicated (as it is also the last one for 
# the object mcmc_run)
mcmc_run2 <- MHalgoGen(n.iter=1000, parameters=parameters_mcmc, data=x, 
likelihood=dnormx, n.chains=1, n.adapt=1, thin=1, trace=1)
mcmc_run3 <- merge(mcmc_run, mcmc_run2)
####### no adaptation, n.adapt must be 0
parameters_mcmc[,"Init"] <- c(mean(x), sd(x))
mcmc_run3 <- MHalgoGen(n.iter=1000, parameters=parameters_mcmc, data=x, 
likelihood=dnormx, n.chains=1, n.adapt=0, thin=1, trace=1)

## End(Not run)

Description

Take a mcmcComposite object and create a vector object with parameter value at specified iteration.
If index="best", the function will return the parameters for the highest likelihood. It also indicates at which iteration the maximum lihelihood has been observed.
If index="last", the function will return the parameters for the last likelihood.
If index="median", the function will return the median value of the parameter.
if index="quantile", the function will return the probs defined by quantiles parameter.
If index="mode", the function will return the mode value of the parameter based on Asselin de Beauville (1978) method.
index can also be a numeric value.
This function uses the complete iterations available except the adaptation part, even if thin parameter is not equal to 1.

Usage

as.parameters(x, index = "best", chain = 1, probs = 0.025, silent = FALSE)

Arguments

Value

A vector with parameters at maximum likelihood or index position

Author(s)

Marc Girondot [email protected]

References

Asselin de Beauville J.-P. (1978). Estimation non paramétrique de la densité et du mode, exemple de la distribution Gamma. Revue de Statistique Appliquée, 26(3):47-70.

See Also

Other mcmcComposite functions: MHalgoGen(), as.mcmc.mcmcComposite(), as.quantiles(), merge.mcmcComposite(), plot.PriorsmcmcComposite(), plot.mcmcComposite(), setPriors(), summary.mcmcComposite()

Examples

## Not run: 
library(HelpersMG)
require(coda)
x <- rnorm(30, 10, 2)
dnormx <- function(data, x) {
 data <- unlist(data)
 return(-sum(dnorm(data, mean=x['mean'], sd=x['sd'], log=TRUE)))
}
parameters_mcmc <- data.frame(Density=c('dnorm', 'dlnorm'), 
                              Prior1=c(10, 0.5), 
                              Prior2=c(2, 0.5), 
                              SDProp=c(1, 1), 
                              Min=c(-3, 0), 
                              Max=c(100, 10), 
                              Init=c(10, 2), 
                              stringsAsFactors = FALSE, 
                              row.names=c('mean', 'sd'))
mcmc_run <- MHalgoGen(n.iter=1000, parameters=parameters_mcmc, data=x, 
                      likelihood=dnormx, n.chains=1, n.adapt=100, 
                      thin=1, trace=1)
plot(mcmc_run, xlim=c(0, 20))
plot(mcmc_run, xlim=c(0, 10), parameters="sd")
mcmcforcoda <- as.mcmc(mcmc_run)
#' heidel.diag(mcmcforcoda)
raftery.diag(mcmcforcoda)
autocorr.diag(mcmcforcoda)
acf(mcmcforcoda[[1]][,"mean"], lag.max=20, bty="n", las=1)
acf(mcmcforcoda[[1]][,"sd"], lag.max=20, bty="n", las=1)
batchSE(mcmcforcoda, batchSize=100)

# The batch standard error procedure is usually thought to 
# be not as accurate as the time series methods used in summary
summary(mcmcforcoda)$statistics[,"Time-series SE"]
summary(mcmc_run)
as.parameters(mcmc_run)
lastp <- as.parameters(mcmc_run, index="last")
parameters_mcmc[,"Init"] <- lastp

# The n.adapt set to 1 is used to not record the first set of parameters
# then it is not duplicated (as it is also the last one for 
# the object mcmc_run)
mcmc_run2 <- MHalgoGen(n.iter=1000, parameters=parameters_mcmc, data=x, 
                       likelihood=dnormx, n.chains=1, n.adapt=1, 
                       thin=1, trace=1)
mcmc_run3 <- merge(mcmc_run, mcmc_run2)

####### no adaptation, n.adapt must be 0
parameters_mcmc[,"Init"] <- c(mean(x), sd(x))
mcmc_run3 <- MHalgoGen(n.iter=1000, parameters=parameters_mcmc, data=x, 
                       likelihood=dnormx, n.chains=1, n.adapt=0, 
                       thin=1, trace=1)
# With index being median, it returns the median value for each parameter
as.parameters(mcmc_run3, index="median")
as.parameters(mcmc_run3, index="mode")
as.parameters(mcmc_run3, index="best")
as.parameters(mcmc_run3, index="quantile", probs=0.025)
as.parameters(mcmc_run3, index="quantile", probs=0.975)
as.parameters(mcmc_run3, index="quantile", probs=c(0.025, 0.975))

## End(Not run)

Description

Extract quantile distribution from mcmcComposite object

Usage

as.quantiles(
  x,
  chain = 1,
  fun = function(...) return(as.numeric(list(...))),
  probs = c(0.025, 0.975),
  xlim = NULL,
  nameparxlim = NULL,
  namepar = NULL
)

Arguments

Value

A data.frame with quantiles

Author(s)

Marc Girondot [email protected]

See Also

Other mcmcComposite functions: MHalgoGen(), as.mcmc.mcmcComposite(), as.parameters(), merge.mcmcComposite(), plot.PriorsmcmcComposite(), plot.mcmcComposite(), setPriors(), summary.mcmcComposite()

Examples

## Not run: 
library(HelpersMG)
require(coda)
x <- rnorm(30, 10, 2)
dnormx <- function(data, x) {
 data <- unlist(data)
 return(-sum(dnorm(data, mean=x['mean'], sd=x['sd'], log=TRUE)))
}
parameters_mcmc <- data.frame(Density=c('dnorm', 'dlnorm'), 
Prior1=c(10, 0.5), Prior2=c(2, 0.5), SDProp=c(1, 1), 
Min=c(-3, 0), Max=c(100, 10), Init=c(10, 2), stringsAsFactors = FALSE, 
row.names=c('mean', 'sd'))
mcmc_run <- MHalgoGen(n.iter=10000, parameters=parameters_mcmc, data=x, 
likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=1)
k <- as.quantiles(x=mcmc_run, namepar="mean")
k <- as.quantiles(x=mcmc_run, namepar="mean", 
                 xlim=c(1:5), nameparxlim="sd", 
                 fun=function(...) return(sum(as.numeric(list(...)))))

## End(Not run)

Description

Return the codes (in UTF-8) of a string.

Usage

asc(x)

Arguments

Details

asc returns the codes (in UTF-8) of a string

Value

A vector with ITF-8 codes of a string

Author(s)

Based on this blog: http://datadebrief.blogspot.com/2011/03/ascii-code-table-in-r.html

See Also

Other Characters: char(), d(), tnirp()

Examples

asc("abcd")
asc("ABCD")

Description

To plot data, just use it as a normal barplot but add the errbar.y values or errbar.y.minus, errbar.y.plus if bars for y axis are asymetric. Use y.plus and y.minus to set absolut limits for error bars. Note that y.plus and y.minus have priority over errbar.y, errbar.y.minus and errbar.y.plus.

Usage

barplot_errbar(
  ...,
  errbar.y = NULL,
  errbar.y.plus = NULL,
  errbar.y.minus = NULL,
  y.plus = NULL,
  y.minus = NULL,
  errbar.tick = 1/50,
  errbar.lwd = par("lwd"),
  errbar.lty = par("lty"),
  errbar.col = par("fg"),
  add = FALSE
)

Arguments

Details

barplot_errbar plot a barplot with error bar on y

Value

A numeric vector (or matrix, when beside = TRUE), say mp, giving the coordinates of all the bar midpoints drawn, useful for adding to the graph.
If beside is true, use colMeans(mp) for the midpoints of each group of bars, see example.

Author(s)

Marc Girondot [email protected]

See Also

plot_errorbar

Other plot and barplot functions: ScalePreviousPlot(), plot_add(), plot_errbar(), show_name()

Examples

## Not run: 
barplot_errbar(rnorm(10, 10, 3), 
	xlab="axe x", ylab="axe y", bty="n", 
		errbar.y.plus=rnorm(10, 1, 0.1), col=rainbow(10), 
		names.arg=paste("Group",1:10), cex.names=0.6)
y <- rnorm(10, 10, 3)
barplot_errbar(y, 
               	xlab="axe x", ylab="axe y", bty="n", 
            		y.plus=y+2)

## End(Not run)

Description

Draw a curved line with arrowhead.

Usage

cArrows(
  x1,
  y1,
  x2,
  y2,
  code = 2,
  size = 1,
  width = 1.2/4/cin,
  open = TRUE,
  sh.adj = 0.1,
  sh.lwd = 1,
  sh.col = par("fg"),
  sh.lty = 1,
  h.col = sh.col,
  h.col.bo = sh.col,
  h.lwd = sh.lwd,
  h.lty = sh.lty,
  curved = FALSE,
  beautiful.arrow = 2/3
)

Arguments

Details

cArrows draws curved lines with arrowhead

Value

A list wit lab.x and lab.y being the position where to draw label

Author(s)

Modified from iGraph

Examples

plot(c(1, 10), c(1, 10), type="n", bty="n")
cArrows(x1=2, y1=2, x2=6, y2=6, curved=1)
cArrows(x1=2, y1=2, x2=6, y2=6, curved=0)
cArrows(x1=2, y1=2, x2=6, y2=6, curved=1, sh.adj=1)
cArrows(x1=2, y1=2, x2=6, y2=6, curved=-1, open=FALSE)
cArrows(x1=9, y1=2, x2=6, y2=6, curved=-1, open=FALSE, sh.col="red")
cArrows(x1=9, y1=9, x2=6, y2=6, curved=-1, open=FALSE, h.col="red")
cArrows(x1=2, y1=9, x2=6, y2=6, curved=1, open=FALSE, h.col="red", h.col.bo="red")

Description

Return a value in a changed coordinate system.

Usage

ChangeCoordinate(
  x = stop("At least one value to convert must be provided"),
  initial = stop("Set of two values must be provided as references"),
  transformed = stop("Set of two transformed values must be provided")
)

Arguments

Details

ChangeCoordinate returns a value in a changed coordinate

Value

A value in the new system

Author(s)

Marc Girondot [email protected]

Examples

ChangeCoordinate(x=c(10, 20), initial=c(1, 100), transformed=c(0, 1))

Description

Return a string with characters defined by the codes.

Usage

char(n)

Arguments

Details

char returns the characters defined by the codes

Value

A string with characters defined by the codes

Author(s)

Based on this blog: http://datadebrief.blogspot.com/2011/03/ascii-code-table-in-r.html

See Also

Other Characters: asc(), d(), tnirp()

Examples

char(65:75)
char(unlist(tapply(144:175, 144:175, function(x) {c(208, x)})))

Description

Run a shiny application for basic functions of comparison.

Usage

compare()

Details

compare runs a shiny application for basic functions of comparison

Value

Nothing

Author(s)

Marc Girondot [email protected]

References

Girondot, M., Guillon, J.-M., 2018. The w-value: An alternative to t- and X2 tests. Journal of Biostatistics & Biometrics 1, 1-3.

See Also

Other w-value functions: contingencyTable.compare(), series.compare()

Examples

## Not run: 
library(HelpersMG)
compare()

## End(Not run)

Description

This function is used to compare the AIC of several outputs obtained with the same data but with different set of parameters.
The parameters must be lists with $aic or $AIC or $value and $par elements or if AIC(element) is defined.
if $value and $par are present in the object, the AIC is calculated as 2*factor.value*value+2*length(par). If $value is -log(likeihood), then factor.value must be 1 and if $value is log(likeihood), then factor.value must be -1.
If several objects are within the same list, their AIC are summed.
For example, compare_AIC(g1=list(group), g2=list(separe1, separe2)) can be used to compare a single model onto two different sets of data against each set of data fited with its own set of parameters.
Take a look at ICtab in package bbmle which is similar.

Usage

compare_AIC(
  ...,
  factor.value = 1,
  silent = FALSE,
  FUN = function(x) specify_decimal(x, decimals = 2)
)

Arguments

Details

compare_AIC compares the AIC of several outputs obtained with the same data.

Value

A list with DeltaAIC and Akaike weight for the models.

Author(s)

Marc Girondot [email protected]

See Also

Other AIC: ExtractAIC.glm(), FormatCompareAIC(), compare_AICc(), compare_BIC()

Examples

## Not run: 
library("HelpersMG")
# Here two different models are fitted
x <- 1:30
y <- rnorm(30, 10, 2)+log(x)
plot(x, y)
d <- data.frame(x=x, y=y)
m1 <- lm(y ~ x, data=d)
m2 <- lm(y ~ log(x), data=d)
compare_AIC(linear=m1, log=m2)
# Here test if two datasets can be modeled with a single model
x2 <- 1:30
y2 <- rnorm(30, 15, 2)+log(x2)
plot(x, y, ylim=c(5, 25))
plot_add(x2, y2, col="red")
d2 <- data.frame(x=x2, y=y2)
m1_2 <- lm(y ~ x, data=d2)
x_grouped <- c(x, x2)
y_grouped <- c(y, y2)
d_grouped <- data.frame(x=x_grouped, y=y_grouped)
m1_grouped <- lm(y ~ x, data=d_grouped)
compare_AIC(separate=list(m1, m1_2), grouped=m1_grouped)

## End(Not run)

Description

This function is used to compare the AICc of several outputs obtained with the same data but with different set of parameters.
Each object must have associated logLik() method with df and nobs attributes.
AICc for object x will be calculated as 2*factor.value*logLik(x)+(2*attributes(logLik(x))$df*(attributes(logLik(x))$df+1)/(attributes(logLik(x))$nobs-attributes(logLik(x))$df-1).

Usage

compare_AICc(
  ...,
  factor.value = -1,
  silent = FALSE,
  FUN = function(x) specify_decimal(x, decimals = 2)
)

Arguments

Details

compare_AICc compares the AICc of several outputs obtained with the same data.

Value

A list with DeltaAICc and Akaike weight for the models.

Author(s)

Marc Girondot [email protected]

See Also

Other AIC: ExtractAIC.glm(), FormatCompareAIC(), compare_AIC(), compare_BIC()

Examples

## Not run: 
library("HelpersMG")
# Here two different models are fitted
x <- 1:30
y <- rnorm(30, 10, 2)+log(x)
plot(x, y)
d <- data.frame(x=x, y=y)
m1 <- lm(y ~ x, data=d)
m2 <- lm(y ~ log(x), data=d)
compare_BIC(linear=m1, log=m2, factor.value=-1)
# Here test if two datasets can be modeled with a single model
x2 <- 1:30
y2 <- rnorm(30, 15, 2)+log(x2)
plot(x, y, ylim=c(5, 25))
plot_add(x2, y2, col="red")
d2 <- data.frame(x=x2, y=y2)
m1_2 <- lm(y ~ x, data=d2)
x_grouped <- c(x, x2)
y_grouped <- c(y, y2)
d_grouped <- data.frame(x=x_grouped, y=y_grouped)
m1_grouped <- lm(y ~ x, data=d_grouped)
compare_AICc(separate=list(m1, m1_2), grouped=m1_grouped, factor.value=-1)
# Or simply
compare_AICc(m1=list(AICc=100), m2=list(AICc=102))

## End(Not run)

Description

This function is used to compare the BIC of several outputs obtained with the same data but with different set of parameters.
Each object must have associated logLik() method with df and nobs attributes.
BIC for object x will be calculated as 2*factor.value*sum(logLik(x))+sum(attributes(logLik(x))$df)*log(attributes(logLik(x))$nobs)).
When several data (i..n) are included, the global BIC is calculated as:
2*factor.value*sum(logLik(x)) for i..n+sum(attributes(logLik(x))$df) for i..n*log(attributes(logLik(x))$nobs for i..n))

Usage

compare_BIC(
  ...,
  factor.value = -1,
  silent = FALSE,
  FUN = function(x) specify_decimal(x, decimals = 2)
)

Arguments

Details

compare_BIC compares the BIC of several outputs obtained with the same data.

Value

A list with DeltaBIC and Akaike weight for the models.

Author(s)

Marc Girondot [email protected]

See Also

Other AIC: ExtractAIC.glm(), FormatCompareAIC(), compare_AIC(), compare_AICc()

Examples

## Not run: 
library("HelpersMG")
# Here two different models are fitted
x <- 1:30
y <- rnorm(30, 10, 2)+log(x)
plot(x, y)
d <- data.frame(x=x, y=y)
m1 <- lm(y ~ x, data=d)
m2 <- lm(y ~ log(x), data=d)
compare_BIC(linear=m1, log=m2, factor.value=-1)
# Here test if two datasets can be modeled with a single model
x2 <- 1:30
y2 <- rnorm(30, 15, 2)+log(x2)
plot(x, y, ylim=c(5, 25))
plot_add(x2, y2, col="red")
d2 <- data.frame(x=x2, y=y2)
m1_2 <- lm(y ~ x, data=d2)
x_grouped <- c(x, x2)
y_grouped <- c(y, y2)
d_grouped <- data.frame(x=x_grouped, y=y_grouped)
m1_grouped <- lm(y ~ x, data=d_grouped)
compare_BIC(separate=list(m1, m1_2), grouped=m1_grouped, factor.value=-1)

## End(Not run)

Description

This function is used as a replacement of chisq.test() to not use p-value.

Usage

contingencyTable.compare(
  table,
  criterion = c("AIC", "AICc", "BIC"),
  probs = NULL
)

Arguments

Details

contingencyTable.compare compares contingency table using Akaike weight.

Value

The probability that a single proportion model is sufficient to explain the data

Author(s)

Marc Girondot [email protected]

References

Girondot, M., Guillon, J.-M., 2018. The w-value: An alternative to t- and X2 tests. Journal of Biostatistics & Biometrics 1, 1-4.

See Also

Other w-value functions: compare(), series.compare()

Examples

## Not run: 
library("HelpersMG")

# Symmetry of Lepidochelys olivacea scutes
table <- t(data.frame(SriLanka=c(200, 157), AfricaAtl=c(19, 12), 
                      Guyana=c(8, 6), Suriname=c(162, 88), 
                      MexicoPac1984=c(42, 34), MexicoPac2014Dead=c(8, 9),
                      MexicoPac2014Alive=c(13, 12), 
                      row.names =c("Symmetric", "Asymmetric")))
table
contingencyTable.compare(table)

table <- t(data.frame(SriLanka=c(200, 157), AfricaAtl=c(19, 12), Guyana=c(8, 6),
                      Suriname=c(162, 88), MexicoPac1984=c(42, 34), 
                      MexicoPac2014Dead=c(8, 9),
                      MexicoPac2014Alive=c(13, 12), Lepidochelys.kempii=c(99, 1), 
                      row.names =c("Symmetric", "Asymmetric")))
table
contingencyTable.compare(table)

# Conformity to a model
table <- matrix(c(33, 12, 25, 75), ncol = 2, byrow = TRUE)
probs <- c(0.5, 0.5)
contingencyTable.compare(table, probs=probs)

# Conformity to a model
table <- matrix(c(33, 12), ncol = 2, byrow = TRUE)
probs <- c(0.5, 0.5)
contingencyTable.compare(table, probs=probs)

# Conformity to a model
table <- matrix(c(33, 12, 8, 25, 75, 9), ncol = 3, byrow = TRUE)
probs <- c(0.8, 0.1, 0.1)
contingencyTable.compare(table, probs=probs)

# Comparison of chisq.test() and this function
table <- matrix(c(NA, NA, 25, 75), ncol = 2, byrow = TRUE)

pv <- NULL
aw <- NULL
par(new=FALSE)
n <- 100

for (GroupA in 0:n) {
  table[1, 1] <- GroupA
  table[1, 2] <- n-GroupA
  pv <- c(pv, chisq.test(table)$p.value)
  aw <- c(aw, contingencyTable.compare(table, criterion="BIC")[1])
}

x <- 0:n
y <- pv
y2 <- aw
plot(x=x, y=y, type="l", bty="n", las=1, xlab="Number of type P in Group B", ylab="Probability", 
     main="", lwd=2)
lines(x=x, y=y2, type="l", col="red", lwd=2)

# w-value
(l1 <- x[which(aw>0.05)[1]])
(l2 <- rev(x)[which(rev(aw)>0.05)[1]])

aw[l1]
pv[l1]

aw[l2+2]
pv[l2+2]

# p-value
l1 <- which(pv>0.05)[1]
l2 <- max(which(pv>0.05))

aw[l1]
pv[l1]

aw[l2]
pv[l2]

y[which(y2>0.05)[1]]
y[which(rev(y2)>0.05)[1]]

par(xpd=TRUE)
text(x=25, y=1.15, labels="Group A: 25 type P / 100", pos=1)

segments(x0=25, y0=0, x1=25, y1=1, lty=3)

# plot(1, 1)

v1 <- c(expression(italic("p")*"-value"), expression("after "*chi^2*"-test"))
v2 <- c(expression(italic("w")*"-value for A"), expression("and B identical models"))
legend("topright", legend=c(v1, v2), 
       y.intersp = 1, 
       col=c("black", "black", "red", "red"), bty="n", lty=c(1, 0, 1, 0))

segments(x0=0, x1=n, y0=0.05, y1=0.05, lty=2)
text(x=101, y=0.05, labels = "0.05", pos=4)

## End(Not run)

Description

Convert one Date-Time from one timezone to another.
Available timezones can be shown using OlsonNames().

Usage

convert.tz(x, tz = Sys.timezone())

Arguments

Details

convert.tz Convert one Date-Time from one timezone to another

Value

A POSIXlt or POSIXct date converted

Author(s)

Marc Girondot [email protected]

See Also

Function with_tz() from lubridate package does the same. I keep it here only for compatibility with old scripts.

Examples

d <- as.POSIXlt("2010-01-01 17:34:20", tz="UTC")
convert.tz(d, tz="America/Guatemala")

Description

If observations is a data.frame, it can have 4 columns:
A column for the measurements;
A column for the lower detection limit;
A column for the upper detection limit;
A column for the truncated of censored nature of the data.
The names of the different columns are in the observations.colname, lower_detection_limit.colname, upper_detection_limit.colname and cut_method.colname.
If lower_detection_limit.colname is NULL or if the column does not exist, the data are supposed to not be left-cut and if upper_detection_limit.colname is NULL or if the column does not exist, the data are supposed to not be right-cut.
If observations is a vector, then the parameters lower_detection_limit and/or upper_detection_limit must be given. Then cut_method must be also provided.
In abservations, -Inf must be used to indicate a value below the lower detection limit and +Inf must be used for a value above the upper detection limit.
Be careful: NA is used to represent a missing data and not a value below of above the detection limit.
If lower_detection_limit, upper_detection_limit or cut_method are only one value, they are supposed to be used for all the observations.
Definitions for censored or truncated distribution vary, and the two terms are sometimes used interchangeably. Let the following data set:
1 1.25 2 4 5

Censoring: some observations will be censored, meaning that we only know that they are below (or above) some bound. This can for instance occur if we measure the concentration of a chemical in a water sample. If the concentration is too low, the laboratory equipment cannot detect the presence of the chemical. It may still be present though, so we only know that the concentration is below the laboratory's detection limit.

If the detection limit is 1.5, so that observations that fall below this limit is censored, our example data set would become:
<1.5 <1.5 2 4 5; that is, we don't know the actual values of the first two observations, but only that they are smaller than 1.5.

Truncation: the process generating the data is such that it only is possible to observe outcomes above (or below) the truncation limit. This can for instance occur if measurements are taken using a detector which only is activated if the signals it detects are above a certain limit. There may be lots of weak incoming signals, but we can never tell using this detector.

If the truncation limit is 1.5, our example data set would become: 2 4 5; and we would not know that there in fact were two signals which were not recorded.
If n.iter is NULL, no Bayesian MCMC is performed but credible interval will not be available.

Usage

cutter(
  observations = stop("Observations must be provided"),
  observations.colname = "Observations",
  lower_detection_limit.colname = "LDL",
  upper_detection_limit.colname = "UDL",
  cut_method.colname = "Cut",
  par = NULL,
  lower_detection_limit = NULL,
  upper_detection_limit = NULL,
  cut_method = "censored",
  distribution = "gamma",
  n.mixture = 1,
  n.iter = 5000,
  n.adapt = 100,
  debug = FALSE,
  progress.bar = TRUE,
  session = NULL
)

Arguments

Details

cutter returns the fitted distribution without cut

Value

The parameters of distribution of values below or above the detection limit.

Author(s)

Marc Girondot [email protected]

See Also

Other Distributions: dSnbinom(), dbeta_new(), dcutter(), dggamma(), logLik.cutter(), plot.cutter(), print.cutter(), r2norm(), rcutter(), rmnorm(), rnbinom_new()

Examples

## Not run: 
library(HelpersMG)
# _______________________________________________________________
# right censored distribution with gamma distribution
# _______________________________________________________________
# Detection limit
DL <- 100
# Generate 100 random data from a gamma distribution
obc <- rgamma(100, scale=20, shape=2)
# remove the data below the detection limit
obc[obc>DL] <- +Inf
# search for the parameters the best fit these censored data
result <- cutter(observations=obc, upper_detection_limit=DL, 
                           cut_method="censored")
result
plot(result, xlim=c(0, 150), breaks=seq(from=0, to=150, by=10), col.mcmc=NULL)
plot(result, xlim=c(0, 150), breaks=seq(from=0, to=150, by=10))
# _______________________________________________________________
# The same data seen as truncated data with gamma distribution
# _______________________________________________________________
obc <- obc[is.finite(obc)]
# search for the parameters the best fit these truncated data
result <- cutter(observations=obc, upper_detection_limit=DL, 
                           cut_method="truncated")
result
plot(result, xlim=c(0, 150), breaks=seq(from=0, to=150, by=10))
# _______________________________________________________________
# left censored distribution with gamma distribution
# _______________________________________________________________
# Detection limit
DL <- 10
# Generate 100 random data from a gamma distribution
obc <- rgamma(100, scale=20, shape=2)
# remove the data below the detection limit
obc[obc<DL] <- -Inf
# search for the parameters the best fit these truncated data
result <- cutter(observations=obc, lower_detection_limit=DL, 
                          cut_method="censored")
result
plot(result, xlim=c(0, 200), breaks=seq(from=0, to=300, by=10))
# _______________________________________________________________
# left censored distribution with mixture of gamma distribution
# _______________________________________________________________
#' # Detection limit
library(HelpersMG)
# Generate 200 random data from a gamma distribution
set.seed(1234)
obc <- c(rgamma(100, scale=10, shape=5), rgamma(100, scale=20, shape=10))
LDL <- 20
l <- seq(from=0, to=LDL, length.out=1001)
p <- pgamma(l, scale=10, shape=5)*0.5+pgamma(l, scale=20, shape=10)
deltal <- l[2]-l[1]
expected_LDL <- sum((l[-1]-deltal/2)*(p[-1]-p[-length(p)]))/sum((p[-1]-p[-length(p)]))
# remove the data below the detection limit
obc[obc<LDL] <- -Inf

UDL <- 300
l <- seq(from=UDL, to=1000, length.out=1001)
p <- pgamma(l, scale=10, shape=5)*0.5+pgamma(l, scale=20, shape=10)
deltal <- l[2]-l[1]
expected_UDL <- sum((l[-1]-deltal/2)*(p[-1]-p[-length(p)]))/sum((p[-1]-p[-length(p)]))
obc[obc>UDL] <- +Inf

# search for the parameters the best fit these truncated data
result1_gamma <- cutter(observations=obc, lower_detection_limit=LDL, 
                          upper_detection_limit = UDL, 
                          distribution="gamma", 
                          cut_method="censored", n.iter=5000, debug=0)
result1_normal <- cutter(observations=obc, lower_detection_limit=LDL, 
                          upper_detection_limit = UDL, 
                          distribution="normal", 
                          cut_method="censored", n.iter=5000)
result1_lognormal <- cutter(observations=obc, lower_detection_limit=LDL, 
                          upper_detection_limit = UDL, 
                          distribution="lognormal", 
                          cut_method="censored", n.iter=5000)
result1_Weibull <- cutter(observations=obc, lower_detection_limit=LDL, 
                          upper_detection_limit = UDL,  
                          distribution="Weibull", 
                          cut_method="censored", n.iter=5000)
result1_generalized.gamma <- cutter(observations=obc, lower_detection_limit=LDL, 
                          upper_detection_limit = UDL, 
                          distribution="generalized.gamma", 
                          cut_method="censored", n.iter=5000)
result2_gamma <- cutter(observations=obc, lower_detection_limit=LDL, 
                          upper_detection_limit = UDL, 
                          distribution="gamma", 
                          n.mixture=2, 
                          cut_method="censored", n.iter=5000)
result2_normal <- cutter(observations=obc, lower_detection_limit=LDL, 
                          upper_detection_limit = UDL, 
                          distribution="normal", 
                          n.mixture=2, 
                          cut_method="censored", n.iter=5000)
result2_lognormal <- cutter(observations=obc, lower_detection_limit=LDL, 
                          upper_detection_limit = UDL, 
                          distribution="lognormal", 
                          n.mixture=2, 
                          cut_method="censored", n.iter=5000)
result2_Weibull <- cutter(observations=obc, lower_detection_limit=LDL, 
                          upper_detection_limit = UDL, 
                          distribution="Weibull", 
                          n.mixture=2, 
                          cut_method="censored", n.iter=5000)
result2_generalized.gamma <- cutter(observations=obc, lower_detection_limit=LDL, 
                          upper_detection_limit = UDL, 
                          distribution="generalized.gamma", 
                          n.mixture=2, 
                          cut_method="censored", n.iter=5000)
                          
compare_AIC(nomixture.gamma=result1_gamma, 
             nomixture.normal=result1_normal, 
             nomixture.lognormal=result1_lognormal, 
             nomixture.Weibull=result1_Weibull, 
             nomixture.generalized.gamma=result1_generalized.gamma, 
             mixture.gamma=result2_gamma, 
             mixture.normal=result2_normal, 
             mixture.lognormal=result2_lognormal, 
             mixture.Weibull=result2_Weibull, 
             mixture.generalized.gamma=result2_generalized.gamma)
             
plot(result2_gamma, xlim=c(0, 600), breaks=seq(from=0, to=600, by=10))
plot(result2_generalized.gamma, xlim=c(0, 600), breaks=seq(from=0, to=600, by=10))

# _______________________________________________________________
# left and right censored distribution
# _______________________________________________________________
# Generate 100 random data from a gamma distribution
obc <- rgamma(100, scale=20, shape=2)
# Detection limit
LDL <- 10
# remove the data below the detection limit
obc[obc<LDL] <- -Inf
# Detection limit
UDL <- 100
# remove the data below the detection limit
obc[obc>UDL] <- +Inf
# search for the parameters the best fit these censored data
result <- cutter(observations=obc, lower_detection_limit=LDL, 
                           upper_detection_limit=UDL, 
                          cut_method="censored")
result
plot(result, xlim=c(0, 150), col.DL=c("black", "grey"), 
                             col.unobserved=c("green", "blue"), 
     breaks=seq(from=0, to=150, by=10))
# _______________________________________________________________
# Example with two values for lower detection limits
# corresponding at two different methods of detection for example
# with gamma distribution
# _______________________________________________________________
obc <- rgamma(50, scale=20, shape=2)
# Detection limit for sample 1 to 50
LDL1 <- 10
# remove the data below the detection limit
obc[obc<LDL1] <- -Inf
obc2 <- rgamma(50, scale=20, shape=2)
# Detection limit for sample 1 to 50
LDL2 <- 20
# remove the data below the detection limit
obc2[obc2<LDL2] <- -Inf
obc <- c(obc, obc2)
# search for the parameters the best fit these censored data
result <- cutter(observations=obc, 
                           lower_detection_limit=c(rep(LDL1, 50), rep(LDL2, 50)), 
                          cut_method="censored")
result
# It is difficult to choose the best set of colors
plot(result, xlim=c(0, 150), col.dist="red", 
     col.unobserved=c(rgb(red=1, green=0, blue=0, alpha=0.1), 
                      rgb(red=1, green=0, blue=0, alpha=0.2)), 
     col.DL=c(rgb(red=0, green=0, blue=1, alpha=0.5), 
                      rgb(red=0, green=0, blue=1, alpha=0.9)), 
     breaks=seq(from=0, to=200, by=10))
     
# ________________________________________________________________________
# left censored distribution comparison of normal, lognormal, 
# weibull, generalized gamma, and gamma without Bayesian MCMC
# Comparison with Akaike Information Criterion
# ________________________________________________________________________
# Detection limit
DL <- 10
# Generate 100 random data from a gamma distribution
obc <- rgamma(100, scale=20, shape=2)
# remove the data below the detection limit
obc[obc<DL] <- -Inf

result_gamma <- cutter(observations=obc, lower_detection_limit=DL, 
                          cut_method="censored", distribution="gamma", 
                          n.iter=NULL)
plot(result_gamma)

result_lognormal <- cutter(observations=obc, lower_detection_limit=DL, 
                          cut_method="censored", distribution="lognormal", 
                          n.iter=NULL)
plot(result_lognormal)

result_weibull <- cutter(observations=obc, lower_detection_limit=DL, 
                          cut_method="censored", distribution="weibull", 
                          n.iter=NULL)
plot(result_weibull)

result_normal <- cutter(observations=obc, lower_detection_limit=DL, 
                          cut_method="censored", distribution="normal", 
                          n.iter=NULL)
plot(result_normal)

result_generalized.gamma <- cutter(observations=obc, lower_detection_limit=DL, 
                          cut_method="censored", distribution="generalized.gamma", 
                          n.iter=NULL)
plot(result_generalized.gamma)

compare_AIC(gamma=result_gamma, 
            lognormal=result_lognormal, 
            normal=result_normal, 
            Weibull=result_weibull, 
            Generalized.gamma=result_generalized.gamma)

# ______________________________________________________________________________
# left censored distribution comparison of normal, lognormal, 
# weibull, generalized gamma, and gamma
# ______________________________________________________________________________
# Detection limit
DL <- 10
# Generate 100 random data from a gamma distribution
obc <- rgamma(100, scale=20, shape=2)
# remove the data below the detection limit
obc[obc<DL] <- -Inf
# search for the parameters the best fit these truncated data
result_gamma <- cutter(observations=obc, lower_detection_limit=DL, 
                          cut_method="censored", distribution="gamma")
result_gamma
plot(result_gamma, xlim=c(0, 250), breaks=seq(from=0, to=250, by=10))

result_lognormal <- cutter(observations=obc, lower_detection_limit=DL, 
                          cut_method="censored", distribution="lognormal")
result_lognormal
plot(result_lognormal, xlim=c(0, 250), breaks=seq(from=0, to=250, by=10))

result_weibull <- cutter(observations=obc, lower_detection_limit=DL, 
                          cut_method="censored", distribution="weibull")
result_weibull
plot(result_weibull, xlim=c(0, 250), breaks=seq(from=0, to=250, by=10))

result_normal <- cutter(observations=obc, lower_detection_limit=DL, 
                          cut_method="censored", distribution="normal")
result_normal
plot(result_normal, xlim=c(0, 250), breaks=seq(from=0, to=250, by=10))

result_generalized.gamma <- cutter(observations=obc, lower_detection_limit=DL, 
                          cut_method="censored", distribution="generalized.gamma")
result_generalized.gamma
plot(result_generalized.gamma, xlim=c(0, 250), breaks=seq(from=0, to=250, by=10))

# ___________________________________________________________________
# Test for similarity in gamma left censored distribution between two
# datasets
# ___________________________________________________________________
obc1 <- rgamma(100, scale=20, shape=2)
# Detection limit for sample 1 to 50
LDL <- 10
# remove the data below the detection limit
obc1[obc1<LDL] <- -Inf
obc2 <- rgamma(100, scale=10, shape=2)
# remove the data below the detection limit
obc2[obc2<LDL] <- -Inf
# search for the parameters the best fit these censored data
result1 <- cutter(observations=obc1, 
                  distribution="gamma", 
                  lower_detection_limit=LDL, 
                  cut_method="censored", n.iter=NULL)
logLik(result1)
plot(result1, xlim=c(0, 200), 
     breaks=seq(from=0, to=200, by=10))
result2 <- cutter(observations=obc2, 
                  distribution="gamma", 
                  lower_detection_limit=LDL, 
                  cut_method="censored", n.iter=NULL)
logLik(result2)
plot(result2, xlim=c(0, 200), 
     breaks=seq(from=0, to=200, by=10))
result_totl <- cutter(observations=c(obc1, obc2), 
                      distribution="gamma", 
                      lower_detection_limit=LDL, 
                      cut_method="censored", n.iter=NULL)
logLik(result_totl)
plot(result_totl, xlim=c(0, 200), 
     breaks=seq(from=0, to=200, by=10))
     
compare_AIC(Separate=list(result1, result2), 
            Common=result_totl, factor.value=1)
compare_BIC(Separate=list(result1, result2), 
            Common=result_totl, factor.value=1)           

## End(Not run)

Description

Writes an ASCII text representation of an R object.
It can be used as a replacement of dput() for named vectors.
The controls "keepNA", "keepInteger" and "showAttributes" are utilized for named vectors.

Usage

d(
  x,
  file = "",
  control = c("keepNA", "keepInteger", "showAttributes"),
  collapse = ", \n  "
)

Arguments

Details

d Write an ASCII Representation of a vector object

Value

A string

Author(s)

Marc Girondot [email protected]

See Also

Other Characters: asc(), char(), tnirp()

Examples

d(c(A=10, B=20))
dput(c(A=10, B=20))

Description

Density for the Beta distribution with parameters mu and v or shape1 and shape2 (and optional non-centrality parameter ncp).
The returned object has three attributes:
shape1, shape2, and ncp
Note that if x has other attributes, they are preserved.

Usage

dbeta_new(
  x,
  mu = NULL,
  v = NULL,
  shape1,
  shape2,
  ncp = 0,
  log = FALSE,
  silent = FALSE
)

Arguments

Details

dbeta_new returns the density for the Beta distributions

The Beta distribution with parameters shape1 = a and shape2 = b has density
gamma(a+b)/(gamma(a)gamma(b))x^(a-1)(1-x)^(b-1)
for a > 0, b > 0 and 0 <= x <= 1 where the boundary values at x=0 or x=1 are defined as by continuity (as limits).
The mean is a/(a+b) and the variance is ab/((a+b)^2 (a+b+1)). These moments and all distributional properties can be defined as limits.

Value

dbeta_new gives the density for the Beta distributions

Author(s)

Marc Girondot [email protected]

See Also

Other Distributions: cutter(), dSnbinom(), dcutter(), dggamma(), logLik.cutter(), plot.cutter(), print.cutter(), r2norm(), rcutter(), rmnorm(), rnbinom_new()

Examples

pi <- rbeta(100, shape1=0.48, shape2=0.12)
hist(pi, freq=FALSE, breaks=seq(from=0, to=1, by=0.1), ylim=c(0, 8), las=1)
library("HelpersMG")
mx <- ScalePreviousPlot()$ylim["end"]/
      max(dbeta_new(seq(from=0.01, to=0.99, by=0.01), mu = 0.8, v=0.1))
curve(dbeta_new(x, mu = 0.8, v=0.1)*mx, add=TRUE, col="red")
out <- dbeta_new(0.1, mu = 0.8, v=0.1)
out
attributes(out)$shape1; attributes(out)$shape2; attributes(out)$ncp
dbeta(0.1, shape1=attributes(out)$shape1, shape2=attributes(out)$shape2, 
      ncp=attributes(out)$ncp)

# It can be used to generate random numbers using mu and v
out <- dbeta_new(0.1, mu = 0.8, v=0.1, silent=TRUE)
pi <- rbeta(100, shape1=attributes(out)$shape1, shape2=attributes(out)$shape2, 
      ncp=attributes(out)$ncp)
hist(pi, freq=FALSE, breaks=seq(from=0, to=1, by=0.1), ylim=c(0, 8), las=1)

Description

If observations must be a data.frame with 4 columns:
observations: A column for the measurements;
LDL: A column for the lower detection limit;
UDL: A column for the upper detection limit;
Cut: A column for the truncated of censored nature of the data.

Usage

dcutter(
  par,
  observations = NULL,
  distribution = "gamma",
  n.mixture = NULL,
  debug = FALSE,
  limits.lower = NULL,
  limits.upper = NULL,
  log = TRUE
)

Arguments

Details

dcutter returns the density of the cutter function

Value

The density of the cutter function according to observations.

Author(s)

Marc Girondot [email protected]

See Also

Other Distributions: cutter(), dSnbinom(), dbeta_new(), dggamma(), logLik.cutter(), plot.cutter(), print.cutter(), r2norm(), rcutter(), rmnorm(), rnbinom_new()

Examples

## Not run: 
library(HelpersMG)
par <- c('shape1' = 0.42265849507444225, 
         'scale1' = 14.139457094879594, 
         'shape2' = 1.667131542489706, 
         'scale2' = 0.10763344388223803, 
         'p1' = 0.12283307526788023)
obs <- data.frame(Observations=c(0.755, 1.013, 2.098, 6.265, 4.708, 0.078, 2.169, 0.403, 1.251, 
                                 0.008, 1.419, 1.078, 2.744, 81.534, 1.426, 13.486, 7.813, 0.165, 
                                 0.118, 0.864, 0.369, 7.159, 2.605, 1.579, 1.646, 0.484, 4.492, 
                                 0.139, 0.28, 0.154, 0.106, 0.104, 4.185, 0.735, 0.149, 0.183, 
                                 0.062, 8.246, 0.165, 0.121, 0.109, 0.092, 0.162, 0.108, 0.139, 
                                 0.141, 0.124, 0.124, 0.151, 0.141, 0.364, 0.295, 0.09, 0.135, 
                                 0.154, 0.218, 0.167, -Inf, 0.203, 0.228, 0.107, 0.162, 0.194, 
                                 0.322, 0.351, 0.17, 0.236, 0.176, 0.107, 0.12, 0.095, 0.27, 0.194, 
                                 0.125, 0.123, 0.085, 0.164, 0.106, 0.079, 0.162), 
                 LDL=0.001, UDL=NA, Cut="censored")
dcutter(par=par, observations=obs, distribution="gamma", 
        n.mixture=NULL, debug=FALSE, limits.lower=NULL, 
        limits.upper=NULL,log=FALSE)
dcutter(par=par, observations=obs, distribution="gamma", 
        n.mixture=NULL, debug=FALSE, limits.lower=NULL, 
        limits.upper=NULL, log=TRUE)
        

## End(Not run)

Description

Generalized gamma distribution

Usage

dggamma(x, theta, kappa, delta, log = FALSE)

pggamma(q, theta, kappa, delta, lower.tail = TRUE, log.p = FALSE)

qggamma(p, theta, kappa, delta, lower.tail = TRUE, log.p = FALSE)

rggamma(n, theta, kappa, delta)

Arguments

Details

pggamma, qggamma, dggamma, and rggamma are used to model the generalized gamma distribution.

The code is modified from https://rpubs.com/FJRubio/GG.

Value

dggamma gives the density, pggamma gives the distribution function, qggamma gives the quantile function, and rggamma generates random deviates.

Functions

• dggamma(): Density of the generalized gamma.

• pggamma(): Distribution function of the generalized gamma.

• qggamma(): Quantile of the generalized gamma.

• rggamma(): Random of the generalized gamma.

More details here

The generalized gamma is described here https://en.wikipedia.org/wiki/Generalized_gamma_distribution.
With a being theta, b being kappa, and p being delta.
theta, kappa and delta must be all > 0.

Author(s)

Marc Girondot [email protected]

See Also

Other Distributions: cutter(), dSnbinom(), dbeta_new(), dcutter(), logLik.cutter(), plot.cutter(), print.cutter(), r2norm(), rcutter(), rmnorm(), rnbinom_new()

Examples

# To reproduce the wikipedia page graphic
x <- seq(from=0, to=8, by=0.1)
plot(x, dggamma(x, theta=2, kappa=0.5, delta=0.5), lty=1, col="blue", 
     type="l", lwd=2, xlab="x", ylab="PDF")
lines(x, dggamma(x, theta=1, kappa=1, delta=0.5), lty=1, col="green", lwd=2)
lines(x, dggamma(x, theta=2, kappa=1, delta=2), lty=1, col="red", lwd=2)
lines(x, dggamma(x, theta=5, kappa=1, delta=5), lty=1, col="yellow", lwd=2)
lines(x, dggamma(x, theta=7, kappa=1, delta=7), lty=1, col="grey", lwd=2)
legend("topright", legend=c("a=2, d=0.5, p=0.5", "a=1, d=1, p=0.5", 
                            "a=2, d=1, p=2", "a=5, d=1, p=5", "a=7, d=1, p=7"), 
                            col=c("blue", "green", "red", "yellow", "grey"), 
                            lty=1, lwd=2, bty="n")
par <- c(theta=2, kappa=0.5, delta=0.5)
# Mean, var and sd
mean.ggamma <- function(theta, kappa, delta) 
       return(theta*(gamma((kappa+1)/delta))/gamma(kappa/delta))
var.ggamma <- function(theta, kappa, delta) 
       return(theta^2* ( ( (gamma((kappa+2)/delta))/gamma(kappa/delta) ) - 
                ( (gamma((kappa+1)/delta))/gamma(kappa/delta) )^2 ) )
sd.ggamma <- function(theta, kappa, delta) 
       return(sqrt(theta^2* ( ( (gamma((kappa+2)/delta))/gamma(kappa/delta) ) - 
                ( (gamma((kappa+1)/delta))/gamma(kappa/delta) )^2 ) ))

Description

Return an index of quantitative asymmetry and complexity.
Higher is the value, higher is the complexity (number of objects) and diversity (difference between them).
The indice is based on the product of the average angular distance of Edwards (1971) for all permutations of measures for both sides with the geometric mean of the inverse of Shannon entropy H for both sides. Let p1 and p2 two vectors of relative measures of objects with sum(p1) = 1 and sum(p2)=1 and n1 being the number of objects in p1 and n2 being the number of objects in p2.
Edwards distance for all permutations of p1 and p2 objects are computed and the average value E is calculated.
The maximun possible Shannon index for identical n1 is max1 = sum((1/n1) * log(1/n1)).
Shannon index is v1 = sum(p1 * log(p1)).
If version == 2, the complementary of Shannon index for these n1 objects is used: c1 = 2 * max1 - v1
If version == 1, the Shannon index is used directly.
The geometry mean between both sides defined the measure of diversity within each side: S=sqrt(c1 * c2)
The Developmental Instability Index is then S * E

Usage

DIx(l1, l2, details = FALSE, version = 1)

Arguments

Details

DIx returns an index of quantitative asymmetry and complexity

Value

A numeric value

Author(s)

Marc Girondot [email protected]

References

Edwards, A.W.F., 1971. Distances between populations on the basis of gene frequencies. Biometrics 27, 873–881.
Shannon C.E. 1948 A mathematical theory of communication. Bell System Technical Journal 27(3), 379-423.

Examples

## Not run: 
l1 <- c(0.1, 0.1, 0.05, 0.2, 0.3, 0.25)
l2 <- c(0.2, 0.3, 0.5)
DIx(l1, l2)

l1 <- c(0.1, 0.1, 0.05, 0.2, 0.3, 0.25)
l2 <- c(0.1, 0.1, 0.05, 0.2, 0.3, 0.25)
DIx(l1, l2)

l1 <- c(0.2, 0.3, 0.5)
l2 <- c(0.2, 0.3, 0.5)
DIx(l1, l2)

l1 <- c(0.2, 0.2, 0.2, 0.2, 0.2)
l2 <- c(0.2, 0.3, 0.5)
DIx(l1, l2)

l1 <- c(0.2, 0.2, 0.2, 0.2, 0.2)
l2 <- c(0.3333, 0.3333, 0.3333)
DIx(l1, l2)

l1 <- c(0.2, 0.2, 0.2, 0.2, 0.2)
l2 <- c(0.2, 0.2, 0.2, 0.2, 0.2)
DIx(l1, l2)

l1 <- c(0.3333, 0.3333, 0.3333)
l2 <- c(0.3333, 0.3333, 0.3333)
DIx(l1, l2)

## End(Not run)

Description

Density for the negative binomial distribution with parameters mu, sd, var, size or prob. See dnbinom.

Usage

dnbinom_new(
  x,
  size = NULL,
  prob = NULL,
  mu = NULL,
  sd = NULL,
  var = NULL,
  log = FALSE
)

Arguments

Details

dnbinom_new returns density for the negative binomial distribution

Value

Random numbers for the negative binomial distribution

Author(s)

Marc Girondot [email protected]

Examples

## Not run: 
library("HelpersMG")
set.seed(1)
x <- rnbinom_new(n=100, mu=2, sd=3)
LnL <- NULL
df <- data.frame(mu=seq(from=0.1, to=8, by=0.1), "-LnL"=NA)
for (mu in df[, "mu"])
LnL <- c(LnL, -sum(dnbinom_new(x=x, mu=mu, sd=3, log=TRUE)))
df[, "-LnL"] <- LnL
ggplot(data = df, aes(x = .data[["mu"]], y = .data[["-LnL"]])) + geom_line()
# Examples of wrong parametrization
dnbinom_new(x=x, mu=c(1, 2), sd=3, log=TRUE)

## End(Not run)

Description

Distribution of the sum of random variable with negative binomial distributions.
Technically the sum of random variable with negative binomial distributions is a convolution of negative binomial random variables.
dSnbinom returns the density for the sum of random variable with negative binomial distributions.
pSnbinom returns the distribution function for the sum of random variable with negative binomial distributions.
qSnbinom returns the quantile function for the sum of random variable with negative binomial distributions.
rSnbinom returns random numbers for the sum of random variable with negative binomial distributions.
If all prob values are the same, exact probabilities are estimated.
Estimate using Vellaisamy&Upadhye method uses parallel computing depending on value of parallel. The number of cores in usage can be defined using options(mc.cores = c) with c being the number of cores to be used. By default it will use all the available cores. Forking will be used in Unix system and no forking on Windows systems.
When Furman method is in use, it will return the progress of Pr(S = x) during recursion in an attribute if verbose is TRUE (see examples).

Usage

dSnbinom(
  x = stop("You must provide at least one x value"),
  size = NULL,
  prob = NULL,
  mu = NULL,
  log = FALSE,
  tol = NULL,
  method = "Furman",
  normalize = TRUE,
  max.iter = NULL,
  mean = NULL,
  sd = NULL,
  n.random = 1e+06,
  parallel = FALSE,
  verbose = FALSE
)

pSnbinom(
  q = stop("At least one quantile must be provided"),
  size = NULL,
  prob = NULL,
  mu = NULL,
  lower.tail = TRUE,
  log.p = FALSE,
  tol = NULL,
  method = "Furman",
  normalize = TRUE
)

qSnbinom(
  p = stop("At least one probability must be provided"),
  size = stop("size parameter is mandatory"),
  prob = NULL,
  mu = NULL,
  lower.tail = TRUE,
  log.p = FALSE,
  tol = NULL,
  method = "Furman"
)

rSnbinom(n = 1, size = NULL, prob = NULL, mu = NULL)

Arguments

Details

Distribution of the Sum of Independent Negative Binomial Random Variables.

Value

dSnbinom gives the density, pSnbinom gives the distribution function, qSnbinom gives the quantile function, and rSnbinom generates random deviates.

Functions

• dSnbinom(): Density for the sum of random variable with negative binomial distributions.

• pSnbinom(): Distribution function for the sum of random variable with negative binomial distributions.

• qSnbinom(): Quantile function for the sum of random variable with negative binomial distributions.

• rSnbinom(): Random numbers for the sum of random variable with negative binomial distributions.

Author(s)

Marc Girondot [email protected] and Jon Barry [email protected]

References

Furman, E., 2007. On the convolution of the negative binomial random variables. Statistics & Probability Letters 77, 169-172.

Vellaisamy, P. & Upadhye, N.S. 2009. On the sums of compound negative binomial and gamma random variables. Journal of Applied Probability, 46, 272-283.

Girondot M, Barry J. 2023. Computation of the distribution of the sum of independent negative binomial random variables. Mathematical and Computational Applications 2023, 28, 63, doi:10.3390/mca28030063

See Also

Other Distributions: cutter(), dbeta_new(), dcutter(), dggamma(), logLik.cutter(), plot.cutter(), print.cutter(), r2norm(), rcutter(), rmnorm(), rnbinom_new()

Examples

## Not run: 
library(HelpersMG)
alpha <- c(1, 2, 5, 1, 2)
p <- c(0.1, 0.12, 0.13, 0.14, 0.14)
# By default, the Furman method with tol=1E-40 is used
dSnbinom(20, size=alpha, prob=p)
# Note the attribute is the dynamics of convergence of Pr(X=x)
attributes(dSnbinom(20, size=alpha, prob=p, verbose=TRUE))$Pk[, 1]

mutest <- c(0.01, 0.02, 0.03)
sizetest <- 2
x <- 20
# Exact probability
dSnbinom(x, size=sizetest, mu=mutest, method="vellaisamy&upadhye")
dSnbinom(x, size=sizetest, mu=mutest, method="vellaisamy&upadhye", log=TRUE)
# With Furman method and tol=1E-12, when probability 
# is very low, it will be biased
dSnbinom(x, size=sizetest, mu=mutest, method="Furman", tol=1E-12)
# The solution is to use a tolerance lower than the estimate
dSnbinom(x, size=sizetest, mu=mutest, method="Furman", tol=1E-45)
# Here the estimate used a first estimation by saddlepoint approximation
dSnbinom(x, size=sizetest, mu=mutest, method="Furman", tol=NULL)
# Or a huge number of iterations; but it is not the best solution
dSnbinom(x, size=sizetest, mu=mutest, method="Furman", 
         tol=1E-12, max.iter=10000)
# With the saddle point approximation method
dSnbinom(x, size=sizetest, mu=mutest, method="saddlepoint", log=FALSE)
dSnbinom(x, size=sizetest, mu=mutest, method="saddlepoint", log=TRUE)

# Another example
sizetest <- c(1, 1, 0.1)
mutest <- c(2, 1, 10)
x <- 5
(exact <- dSnbinom(x=x, size=sizetest, mu=mutest, method="Vellaisamy&Upadhye"))
(sp <- dSnbinom(x=x, size=sizetest, mu=mutest, method="saddlepoint"))
paste0("Saddlepoint approximation: Error of ", specify_decimal(100*abs(sp-exact)/exact, 2), "%")
(furman <- dSnbinom(x=x, size=sizetest, mu=mutest, method="Furman"))
paste0("Inversion of mgf: Error of ", specify_decimal(100*abs(furman-exact)/exact, 2), "%")
(na <- dSnbinom(x=x, size=sizetest, mu=mutest, method="approximate.normal")) 
paste0("Gaussian approximation: Error of ", specify_decimal(100*abs(na-exact)/exact, 2), "%")
(nb <- dSnbinom(x=x, size=sizetest, mu=mutest, method="approximate.negativebinomial"))
paste0("NB approximation: Error of ", specify_decimal(100*abs(nb-exact)/exact, 2), "%")

plot(0:20, dSnbinom(0:20, size=sizetest, mu=mutest, method="furman"), bty="n", type="h", 
     xlab="x", ylab="Density", ylim=c(0, 0.2), las=1)
points(x=0:20, y=dSnbinom(0:20, size=sizetest, mu=mutest, 
                           method="saddlepoint"), pch=1, col="blue")
points(x=0:20, y=dSnbinom(0:20, size=sizetest, mu=mutest, 
                           method="approximate.negativebinomial"), 
                           col="red")
points(x=0:20, y=dSnbinom(0:20, size=sizetest, mu=mutest, 
                           method="approximate.normal"), 
                           col="green")

# Test with a single distribution
dSnbinom(20, size=1, mu=20)
# when only one distribution is available, it is the same as dnbinom()
dnbinom(20, size=1, mu=20)

# If a parameter is supplied as only one value, it is supposed to be constant
dSnbinom(20, size=1, mu=c(14, 15, 10))
dSnbinom(20, size=c(1, 1, 1), mu=c(14, 15, 10))

# The functions are vectorized:
plot(0:200, dSnbinom(0:200, size=alpha, prob=p, method="furman"), bty="n", type="h", 
     xlab="x", ylab="Density")
points(0:200, dSnbinom(0:200, size=alpha, prob=p, method="saddlepoint"), 
     col="red", pch=3)
     
# Comparison with simulated distribution using rep replicates
alpha <- c(2.1, 2.05, 2)
mu <- c(10, 30, 20)
rep <- 100000
distEmpirique <- rSnbinom(rep, size=alpha, mu=mu)
tabledistEmpirique <- rep(0, 301)
names(tabledistEmpirique) <- as.character(0:300)
tabledistEmpirique[names(table(distEmpirique))] <- table(distEmpirique)/rep

plot(0:300, dSnbinom(0:300, size=alpha, mu=mu, method="furman"), type="h", bty="n", 
   xlab="x", ylab="Density", ylim=c(0,0.02))
plot_add(0:(length(tabledistEmpirique)-1), tabledistEmpirique, type="l", col="red")
legend(x=200, y=0.02, legend=c("Empirical", "Theoretical"), 
   text.col=c("red", "black"), bty="n")


# Example from Vellaisamy, P. & Upadhye, N.S. (2009) - Table 1
# Note that computing time for k = 7 using exact method is very long
k <- 2:7
x <- c(3, 5, 8, 10, 15)
table1_Vellaisamy <- matrix(NA, ncol=length(x), nrow=length(k))
rownames(table1_Vellaisamy) <- paste0("n = ", as.character(k))
colnames(table1_Vellaisamy) <- paste0("x = ", as.character(x))
table1_approximateObservations <- table1_Vellaisamy
table1_Furman3 <- table1_Vellaisamy
table1_Furman6 <- table1_Vellaisamy
table1_Furman9 <- table1_Vellaisamy
table1_Furman12 <- table1_Vellaisamy
table1_Furman40 <- table1_Vellaisamy
table1_Furman40 <- table1_Vellaisamy
table1_FurmanAuto <- table1_Vellaisamy
table1_FurmanAuto_iter <- table1_Vellaisamy
table1_Vellaisamy_parallel <- table1_Vellaisamy
table1_Approximate_Normal <- table1_Vellaisamy
table1_saddlepoint <- table1_Vellaisamy

st_Furman3 <- rep(NA, length(k))
st_Furman6 <- rep(NA, length(k))
st_Furman9 <- rep(NA, length(k))
st_Furman12 <- rep(NA, length(k))
st_Furman40 <- rep(NA, length(k))
st_FurmanAuto <- rep(NA, length(k))
st_approximateObservations <- rep(NA, length(k))
st_Vellaisamy <- rep(NA, length(k))
st_Vellaisamy_parallel <- rep(NA, length(k))
st_Approximate_Normal <- rep(NA, length(k))
st_saddlepoint <- rep(NA, length(k))

for (n in k) {
    print(n)
    alpha <- 1:n
    p <- (1:n)/10
    st_Vellaisamy[which(n == k)] <- 
        system.time({
        table1_Vellaisamy[which(n == k), ] <- dSnbinom(x=x, prob=p, size=alpha, 
                                     method="Vellaisamy&Upadhye", log=FALSE, verbose=FALSE)
        })[1]
    st_Vellaisamy_parallel[which(n == k)] <- 
        system.time({
        table1_Vellaisamy_parallel[which(n == k), ] <- dSnbinom(x=x, prob=p, size=alpha, 
                                     parallel=TRUE, 
                                     method="Vellaisamy&Upadhye", log=FALSE, verbose=FALSE)
        })[1]
    st_approximateObservations[which(n == k)] <- 
        system.time({
        table1_approximateObservations[which(n == k), ] <- dSnbinom(x=x, prob=p, size=alpha, 
                                     method="approximate.RandomObservations", log=FALSE, 
                                     verbose=FALSE)
        })[1]
    st_Furman3[which(n == k)] <- 
        system.time({
            table1_Furman3[which(n == k), ] <- dSnbinom(x=x, prob=p, size=alpha, 
                                     method="Furman", tol=1E-3, log=FALSE, 
                                     verbose=FALSE)
        })[1]
    st_Furman6[which(n == k)] <- 
        system.time({
            table1_Furman6[which(n == k), ] <- dSnbinom(x=x, prob=p, size=alpha, 
                                     method="Furman", tol=1E-6, log=FALSE, 
                                     verbose=FALSE)
        })[1]
    st_Furman9[which(n == k)] <- 
        system.time({
            table1_Furman9[which(n == k), ] <- dSnbinom(x=x, prob=p, size=alpha, 
                                     method="Furman", tol=1E-9, log=FALSE, 
                                     verbose=FALSE)
        })[1]
    st_Furman12[which(n == k)] <- 
        system.time({
            table1_Furman12[which(n == k), ] <- dSnbinom(x=x, prob=p, size=alpha, 
                                     method="Furman", tol=1E-12, log=FALSE, 
                                     verbose=FALSE)
        })[1]
    st_Furman40[which(n == k)] <- 
        system.time({
            table1_Furman40[which(n == k), ] <- dSnbinom(x=x, prob=p, size=alpha, 
                                     method="Furman", tol=1E-40, log=FALSE, 
                                     verbose=FALSE)
        })[1]

    st_FurmanAuto[which(n == k)] <- 
        system.time({
            table1_FurmanAuto[which(n == k), ] <- dSnbinom(x=x, prob=p, size=alpha, 
                                     method="Furman", tol=NULL, log=FALSE, 
                                     verbose=FALSE)
        })[1]

    st_Approximate_Normal[which(n == k)] <- 
        system.time({
            table1_Approximate_Normal[which(n == k), ] <- dSnbinom(x=x, prob=p, size=alpha, 
                                     method="approximate.normal", tol=1E-12, log=FALSE, 
                                     verbose=FALSE)
        })[1]
        st_saddlepoint[which(n == k)] <- 
        system.time({
            table1_saddlepoint[which(n == k), ] <- dSnbinom(x=x, prob=p, size=alpha, 
                                     method="saddlepoint", tol=1E-12, log=FALSE, 
                                     verbose=FALSE)
        })[1]

for (xc in x) {
 essai <- dSnbinom(x=xc, prob=p, size=alpha, method="Furman", tol=NULL, log=FALSE, verbose=TRUE)
 table1_FurmanAuto_iter[which(n == k), which(xc == x)] <- nrow(attributes(essai)[[1]])
}
}

cbind(table1_Vellaisamy, st_Vellaisamy)
cbind(table1_Vellaisamy_parallel, st_Vellaisamy_parallel)
cbind(table1_Furman3, st_Furman3)
cbind(table1_Furman6, st_Furman6)
cbind(table1_Furman9, st_Furman9)
cbind(table1_Furman12, st_Furman12)
cbind(table1_Furman40, st_Furman40)
cbind(table1_FurmanAuto, st_FurmanAuto)
cbind(table1_approximateObservations, st_approximateObservations)
cbind(table1_Approximate_Normal, st_Approximate_Normal)
cbind(table1_saddlepoint, st_saddlepoint)


# Test of different methods
n <- 9
x <- 17
alpha <- 1:n
p <- (1:n)/10

# Parallel computing is not always performant
# Here it is very performant
system.time({print(dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", log=FALSE, 
         verbose=TRUE, parallel=TRUE))})
system.time({print(dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", log=FALSE, 
         verbose=TRUE, parallel=FALSE))})
         
# Test of different methods
n <- 7
x <- 8
alpha <- 1:n
p <- (1:n)/10

# Parallel computing is not always performant
# Here it is approximately the same time of execution
system.time({print(dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", log=FALSE, 
         verbose=TRUE, parallel=TRUE))})
system.time({print(dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", log=FALSE, 
         verbose=TRUE, parallel=FALSE))})
         
# Test of different methods
n <- 7
x <- 15
alpha <- 1:n
p <- (1:n)/10

# Parallel computing is sometimes very performant
# Here parallel computing is 7 times faster (with a 8 cores computer) 
#             for vellaisamy&upadhye method
system.time(dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", log=FALSE, 
         verbose=TRUE, parallel=TRUE))
system.time(dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", log=FALSE, 
         verbose=TRUE, parallel=FALSE))
         
# Test of different methods
n <- 2
x <- 3
alpha <- 1:n
p <- (1:n)/10

# Parallel computing is sometimes very performant
# Here parallel computing is 7 times faster (with a 8 cores computer) 
#             for vellaisamy&upadhye method
system.time(dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", log=FALSE, 
         verbose=TRUE, parallel=TRUE))
system.time(dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", log=FALSE, 
         verbose=TRUE, parallel=FALSE))
         
# Test for different tolerant values
n <- 7
x <- 8
alpha <- 1:n
p <- (1:n)/10
dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", log=FALSE, verbose=TRUE)
as.numeric(dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, tol=1E-3, verbose=TRUE))
as.numeric(dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, tol=1E-6, verbose=TRUE))
as.numeric(dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, tol=1E-9, verbose=TRUE))
as.numeric(dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, verbose=TRUE))
as.numeric(dSnbinom(x=x, prob=p, size=alpha, method="Saddlepoint", log=FALSE, verbose=TRUE))
as.numeric(dSnbinom(x=x, prob=p, size=alpha, method="approximate.RandomObservations", 
         log=FALSE, verbose=TRUE))

# Test for criteria of convergence
Pr_exact <- dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", 
                   log=FALSE, verbose=TRUE)
Pr_Furman <- dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, 
                 verbose=TRUE)
Pr_exact;as.numeric(Pr_Furman)
plot(1:length(attributes(Pr_Furman)$Pk), 
     log10(abs(attributes(Pr_Furman)$Pk-Pr_exact)), type="l", xlab="Iterations", 
     ylab="Abs log10", bty="n")
lines(1:(length(attributes(Pr_Furman)$Pk)-1), 
      log10(abs(diff(attributes(Pr_Furman)$Pk))), col="red")
legend("bottomleft", legend=c("Log10 Convergence to true value", "Log10 Rate of change"), 
       col=c("black", "red"), 
       lty=1)
       
n <- 7
x <- 6
alpha <- 1:n
p <- (1:n)/10
Pr_exact <- dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", 
                   log=FALSE, verbose=TRUE)
Pr_Furman <- dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, 
                 verbose=TRUE)
Pr_saddlepoint <- dSnbinom(x=x, prob=p, size=alpha, method="saddlepoint", log=FALSE,  
                 verbose=TRUE)                  
                 
# pdf("figure.pdf", width=7, height=7, pointsize=14)                 
ylab <- as.expression(bquote(.("P(S=6) x 10")^"6"))
layout(1:2)
par(mar=c(3, 4, 1, 1))
plot(1:length(attributes(Pr_Furman)$Pk), 
     attributes(Pr_Furman)$Pk*1E6, type="l", xlab="", 
     ylab="", bty="n", las=1, xlim=c(0, 80))
mtext(text=ylab, side=2, line=2.5)
segments(x0=25, x1=80, y0=Pr_exact*1E6, y1=Pr_exact*1E6, lty=3)
par(xpd=TRUE)
text(x=0, y=6, labels="Exact probability", pos=4)
txt <- "        Approximate probability\n    based on Furman (2007)\nrecursive iterations"
text(x=20, y=2, labels=txt, pos=4)
text(x=75, y=5, labels="A", cex=2)
par(mar=c(4, 4, 1, 1))
ylab <- as.expression(bquote("log"["10"]*""*"(P"["k+1"]*""*" - P"["k"]*""*")"))
plot(1:(length(attributes(Pr_Furman)$Pk)-1), 
      log10(diff(attributes(Pr_Furman)$Pk)), col="black", xlim=c(0, 80), type ="l", 
      bty="n", las=1, xlab="Iterations", ylab="")
mtext(text=ylab, side=2, line=2.5)
 peak <- (1:(length(attributes(Pr_Furman)$Pk)-1))[which.max(
            log10(abs(diff(attributes(Pr_Furman)$Pk))))]
segments(x0=peak, x1=peak, y0=-12, y1=-5, lty=2)
text(x=0, y=-6, labels="Positive trend", pos=4)
text(x=30, y=-6, labels="Negative trend", pos=4)
segments(x0=0, x1=43, y0=-12, y1=-12, lty=4)
segments(x0=62, x1=80, y0=-12, y1=-12, lty=4)
text(x=45, y=-12, labels="Tolerance", pos=4)
text(x=75, y=-7, labels="B", cex=2)
# dev.off()

# pdf("figure 2.pdf", width=7, height=7, pointsize=14)                 
ylab <- as.expression(bquote(.("P(S=6) x 10")^"6"))
layout(1:2)
par(mar=c(3, 4, 1, 1))
plot(1:length(attributes(Pr_Furman)$Pk), 
     attributes(Pr_Furman)$Pk*1E6, type="l", xlab="", 
     ylab="", bty="n", las=1, xlim=c(0, 80))
mtext(text=ylab, side=2, line=2.5)
segments(x0=25, x1=80, y0=Pr_exact*1E6, y1=Pr_exact*1E6, lty=3)
par(xpd=TRUE)
text(x=0, y=6, labels="Exact probability", pos=4)
txt <- "        Approximate probability\n    based on Furman (2007)\nrecursive iterations"
text(x=20, y=2, labels=txt, pos=4)
text(x=75, y=5, labels="A", cex=2)
par(mar=c(4, 4, 1, 1))
ylab <- as.expression(bquote("(P"["k+1"]*""*" - P"["k"]*""*") x 10"^"7"))
plot(1:(length(attributes(Pr_Furman)$Pk)-1), 
      diff(attributes(Pr_Furman)$Pk)*1E7, 
      col="black", xlim=c(0, 80), type ="l", 
      bty="n", las=1, xlab="Iterations", ylab="")
mtext(text=ylab, side=2, line=2.5)
 peak <- (1:(length(attributes(Pr_Furman)$Pk)-1))[which.max(diff(attributes(Pr_Furman)$Pk))]
segments(x0=peak, x1=peak, y0=0, y1=3.5, lty=2)
text(x=-2, y=3.5, labels="Positive trend", pos=4)
text(x=30, y=3.5, labels="Negative trend", pos=4)
segments(x0=0, x1=22, y0=1E-12, y1=1E-12, lty=4)
segments(x0=40, x1=80, y0=1E-12, y1=1E-12, lty=4)
text(x=22, y=1E-12+0.2, labels="Tolerance", pos=4)
text(x=75, y=3, labels="B", cex=2)
# dev.off()

# Test of different methods
n <- 2
x <- 15
alpha <- 1:n
p <- (1:n)/10
dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", log=FALSE, verbose=TRUE)
as.numeric(dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, tol=1E-3, verbose=TRUE))
as.numeric(dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, tol=1E-6, verbose=TRUE))
as.numeric(dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, tol=1E-9, verbose=TRUE))
as.numeric(dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, verbose=TRUE))
dSnbinom(x=x, prob=p, size=alpha, method="approximate.RandomObservations", 
         log=FALSE, verbose=TRUE)


n <- 50
x <- 300
alpha <- (1:n)/100
p <- (1:n)/1000
# Produce an error
Pr_exact <- dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", 
                   log=FALSE, verbose=TRUE)
Pr_Furman <- dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, 
                 verbose=FALSE)
Pr_ApproximateNormal <- dSnbinom(x=x, prob=p, size=alpha, method="approximate.normal", 
                        log=FALSE, 
                        verbose=TRUE)
Pr_ApproximateRandom <- dSnbinom(x=x, prob=p, size=alpha, method="approximate.RandomObservations", 
                        log=FALSE, n.random=1E6, 
                        verbose=TRUE)
                        
n <- 500
x <- 3000
alpha <- (1:n)/100
p <- (1:n)/1000
# Produce an error
Pr_exact <- dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", 
                   log=FALSE, verbose=TRUE)
Pr_Furman <- dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, 
                 verbose=FALSE)
Pr_ApproximateNormal <- dSnbinom(x=x, prob=p, size=alpha, method="approximate.normal", 
                        log=FALSE, 
                        verbose=TRUE)
Pr_ApproximateNegativeBinomial <- dSnbinom(x=x, prob=p, size=alpha, 
                        method="approximate.negativebinomial", 
                        log=FALSE, 
                        verbose=TRUE)
Pr_ApproximateRandom <- dSnbinom(x=x, prob=p, size=alpha, 
                        method="approximate.RandomObservations", 
                        log=FALSE, n.random=1E6, 
                        verbose=TRUE)
Pr_ApproximateSaddlepoint <- dSnbinom(x=x, prob=p, size=alpha, 
                        method="saddlepoint", 
                        log=FALSE, 
                        verbose=TRUE)
             
             
             
layout(matrix(1:4, ncol=2, byrow=TRUE))
par(mar=c(3, 4.5, 1, 1))
alpha <- seq(from=10, to=100, length.out=3)
p <- seq(from=0.5, to=0.9, length.out=3)

p_nb <- dSnbinom(0:100, prob=p, size=alpha, method="vellaisamy&upadhye", verbose=TRUE)
p_Furman <- dSnbinom(0:100, prob=p, size=alpha, method="Furman", verbose=FALSE)
p_normal <- dSnbinom(0:100, prob=p, size=alpha, method="approximate.normal", verbose=TRUE)
p_aNB <- dSnbinom(0:100, prob=p, size=alpha, method="approximate.negativebinomial", verbose=TRUE)
p_SA <- dSnbinom(0:100, prob=p, size=alpha, method="saddlepoint", verbose=TRUE)

lab_PSnx <- bquote(italic("P(S"* "" [n] * "=x)"))

plot(1, 1, las=1, bty="n", col="grey", xlab="", 
       xlim=c(10, 70), ylim=c(0, 0.05), 
       ylab=lab_PSnx, type="n")
par(xpd=FALSE)
segments(x0=(0:100), x1=(0:100), 
         y0=0, y1=as.numeric(p_nb), col="black")
         
plot(x=p_nb, y=p_normal, pch=4, cex=0.5, las=1, bty="n", 
     xlab=bquote(italic("P" * "" [exact] * "(S"* "" [n] * "=x)")), 
     ylab=bquote(italic("P" * "" [approximate] * "(S"* "" [n] * "=x)")),  
     xlim=c(0, 0.05), ylim=c(0, 0.05))
points(x=p_nb, y=p_aNB, pch=5, cex=0.5)
points(x=p_nb, y=p_SA, pch=6, cex=0.5)
points(x=p_Furman, y=p_SA, pch=19, cex=0.5)

n <- 2
x <- 15
alpha <- 1:n
p <- (1:n)/10

p_nb <- dSnbinom(0:80, prob=p, size=alpha, method="vellaisamy&upadhye", verbose=TRUE)
p_Furman <- dSnbinom(0:80, prob=p, size=alpha, method="Furman", verbose=FALSE)
p_normal <- dSnbinom(0:80, prob=p, size=alpha, method="approximate.normal", verbose=TRUE)
p_aNB <- dSnbinom(0:80, prob=p, size=alpha, method="approximate.negativebinomial", verbose=TRUE)
p_SA <- dSnbinom(0:80, prob=p, size=alpha, method="saddlepoint", verbose=TRUE)


par(mar=c(4, 4.5, 1, 1))
plot(1, 1, las=1, bty="n", col="grey", xlab="x", 
       xlim=c(0, 60), ylim=c(0, 0.05), 
       ylab=lab_PSnx, type="n")
par(xpd=FALSE)
segments(x0=(0:80), x1=(0:80), 
         y0=0, y1=as.numeric(p_nb), col="black")
         
plot(x=p_nb, y=p_normal, pch=4, cex=0.5, las=1, bty="n", 
     xlab=bquote(italic("P" * "" [exact] * "(S"* "" [n] * "=x)")), 
     ylab=bquote(italic("P" * "" [approximate] * "(S"* "" [n] * "=x)")), 
     xlim=c(0, 0.05), ylim=c(0, 0.05))
points(x=p_nb, y=p_aNB, pch=5, cex=0.5)
points(x=p_nb, y=p_SA, pch=6, cex=0.5)
points(x=p_Furman, y=p_SA, pch=19, cex=0.5)

# pdf("figure 1.pdf", width=7, height=7, pointsize=14)    

layout(1:2)
par(mar=c(3, 4, 1, 1))
alpha <- seq(from=10, to=100, length.out=3)
p <- seq(from=0.5, to=0.9, length.out=3)

p_nb <- dSnbinom(0:100, prob=p, size=alpha, method="vellaisamy&upadhye", verbose=TRUE)
p_Furman <- dSnbinom(0:100, prob=p, size=alpha, method="Furman", verbose=FALSE)
p_normal <- dSnbinom(0:100, prob=p, size=alpha, method="approximate.normal", verbose=TRUE)
p_aNB <- dSnbinom(0:100, prob=p, size=alpha, method="approximate.negativebinomial", verbose=TRUE)
p_SA <- dSnbinom(0:100, prob=p, size=alpha, method="saddlepoint", verbose=TRUE)

lab_PSnx <- bquote(italic("P(S"* "" [n] * "=x)"))

plot(1, 1, las=1, bty="n", col="grey", xlab="", 
       xlim=c(10, 70), ylim=c(0, 0.09), 
       ylab="", type="n", yaxt="n")
axis(2, at=seq(from=0, to=0.05, by=0.01), las=1)
mtext(lab_PSnx, side = 2, adj=0.3, line=3)
par(xpd=FALSE)
segments(x0=(0:100), x1=(0:100), 
         y0=0, y1=as.numeric(p_nb), col="black")
errr <- (abs((100*(p_Furman-p_nb)/p_nb)))/1000+0.055
errr <- ifelse(is.infinite(errr), NA, errr)
lines(x=(0:100), y=errr, lty=5, col="red", lwd=2)
errr <- (abs((100*(p_normal-p_nb)/p_nb)))/1000+0.055
lines(x=(0:100), y=errr, lty=2, col="blue", lwd=2)
errr <- (abs((100*(p_aNB-p_nb)/p_nb)))/1000+0.055
lines(x=(0:100), y=errr, lty=3, col="purple", lwd=2)
errr <- (abs((100*(p_SA-p_nb)/p_nb)))/1000+0.055
lines(x=(0:100), y=errr, lty=4, col="green", lwd=2)
axis(2, at=seq(from=0, to=40, by=10)/1000+0.055, las=1, 
     labels=as.character(seq(from=0, to=40, by=10)))
mtext("|% error|", side = 2, adj=0.9, line=3)
par(xpd=TRUE)
legend(x=30, y=0.1, legend=c("Inversion of mgf", "Saddlepoint", "Normal", "Negative binomial"), 
       lty=c(5, 4, 2, 3), bty="n", cex=0.8, col=c("red", "green", "blue", "purple"), lwd=2)
legend(x=10, y=0.05, legend=c("Exact"), lty=c(1), bty="n", cex=0.8)
par(xpd=TRUE)
text(x=ScalePreviousPlot(x = 0.95, y = 0.1)$x, 
     y=ScalePreviousPlot(x = 0.95, y = 0.1)$y, labels="A", cex=2)
     
     
# When normal approximation will fail
n <- 2
x <- 15
alpha <- 1:n
p <- (1:n)/10

p_nb <- dSnbinom(0:80, prob=p, size=alpha, method="vellaisamy&upadhye", verbose=TRUE)
p_Furman <- dSnbinom(0:80, prob=p, size=alpha, method="Furman", verbose=FALSE)
p_normal <- dSnbinom(0:80, prob=p, size=alpha, method="approximate.normal", verbose=TRUE)
p_aNB <- dSnbinom(0:80, prob=p, size=alpha, method="approximate.negativebinomial", verbose=TRUE)
p_SA <- dSnbinom(0:80, prob=p, size=alpha, method="saddlepoint", verbose=TRUE)

par(mar=c(4, 4, 1, 1))
plot(1, 1, las=1, bty="n", col="grey", xlab="x", 
       xlim=c(0, 60), ylim=c(0, 0.09), 
       ylab="", type="n", yaxt="n")
axis(2, at=seq(from=0, to=0.05, by=0.01), las=1)
mtext(lab_PSnx, side = 2, adj=0.3, line=3)
par(xpd=FALSE)
segments(x0=(0:80), x1=(0:80), 
         y0=0, y1=as.numeric(p_nb), col="black")
errr <- (abs((100*(p_Furman-p_nb)/p_nb)))/1000+0.055
errr <- ifelse(is.infinite(errr), NA, errr)
lines(x=(0:80), y=errr, lty=5, col="red", lwd=2)
errr <- (abs((100*(p_normal-p_nb)/p_nb)))/1000+0.055
lines(x=(0:80), y=errr, lty=2, col="blue", lwd=2)
errr <- (abs((100*(p_aNB-p_nb)/p_nb)))/1000+0.055
lines(x=(0:80), y=errr, lty=3, col="purple", lwd=2)
errr <- (abs((100*(p_SA-p_nb)/p_nb)))/1000+0.055
lines(x=(0:80), y=errr, lty=4, col="green", lwd=2)
axis(2, at=seq(from=0, to=40, by=10)/1000+0.055, las=1, 
     labels=as.character(seq(from=0, to=40, by=10)))
mtext("|% error|", side = 2, adj=0.9, line=3)
legend(x=30, y=0.055, 
       legend=c("Exact", "Inversion of mgf", "Saddlepoint", "Normal", "Negative binomial"), 
       lty=c(1, 5, 4, 2, 3), bty="n", cex=0.8, col=c("black", "red", "green", "blue", "purple"), 
       lwd=c(1, 2, 2, 2, 2))
par(xpd=TRUE)
text(x=ScalePreviousPlot(x = 0.95, y = 0.1)$x, 
     y=ScalePreviousPlot(x = 0.95, y = 0.1)$y, labels="B", cex=2)
     

# dev.off()


# Test for criteria of convergence
Pr_exact <- dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", 
                   log=FALSE, verbose=TRUE)
Pr_Furman <- dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, 
                 verbose=TRUE)
Pr_exact;as.numeric(Pr_Furman)
plot(1:length(attributes(Pr_Furman)$Pk), 
     log10(abs(attributes(Pr_Furman)$Pk-Pr_exact)), type="l", xlab="Iterations", 
     ylab="Abs log10", bty="n")
lines(1:(length(attributes(Pr_Furman)$Pk)-1), 
      log10(abs(diff(attributes(Pr_Furman)$Pk))), col="red")
legend("bottomleft", legend=c("Log10 Convergence to true value", "Log10 Rate of change"), 
       col=c("black", "red"), 
       lty=1)


# Test of different methods
alpha <- c(2.05, 2)
mu <- c(10, 30)
test <- rSnbinom(n=100000, size=alpha, mu=mu)
plot(0:200, table(test)[as.character(0:200)]/sum(table(test), na.rm=TRUE), 
     bty="n", type="h", xlab="x", ylab="Density")
lines(x=0:200, dSnbinom(0:200, size=alpha, mu=mu, log=FALSE, method="Furman"), col="blue")
lines(x=0:200, y=dSnbinom(0:200, size=alpha, mu=mu, log=FALSE, 
                          method="vellaisamy&upadhye"), col="red")
lines(x=0:200, y=dSnbinom(0:200, size=alpha, mu=mu, log=FALSE, 
                          method="approximate.randomobservations"), col="green")

# Test for criteria of convergence for x = 50
x <- 50
# Test for criteria of convergence
Pr_exact <- dSnbinom(x=x, prob=p, size=alpha, method="vellaisamy&upadhye", 
                   log=FALSE, verbose=TRUE)
Pr_Furman <- dSnbinom(x=x, prob=p, size=alpha, method="Furman", log=FALSE, tol=1E-12, 
                 verbose=TRUE)
Pr_exact;as.numeric(Pr_Furman)
plot(1:length(attributes(Pr_Furman)$Pk), 
     log10(abs(attributes(Pr_Furman)$Pk-Pr_exact)), type="l", xlab="Iterations", 
     ylab="Abs log10", bty="n")
lines(1:(length(attributes(Pr_Furman)$Pk)-1), 
      log10(abs(diff(attributes(Pr_Furman)$Pk))), col="red")
legend("bottomleft", legend=c("Log10 Convergence to true value", "Log10 Rate of change"), 
       col=c("black", "red"), 
       lty=1)

# Another example more complicated
set.seed(2)
mutest <- c(56, 6.75, 1)
ktest <- c(50, 50, 50)
nr <- 100000
test <- rSnbinom(nr, size=ktest, mu=mutest)
system.time({pr_vellaisamy <- dSnbinom(x=0:150, size=ktest, mu=mutest, 
                          method = "vellaisamy&upadhye", verbose=FALSE, parallel=FALSE)})
# Parallel computing is not efficient
system.time({pr_vellaisamy <- dSnbinom(x=0:150, size=ktest, mu=mutest, 
                          method = "vellaisamy&upadhye", verbose=FALSE, parallel=TRUE)})
system.time({pr_furman <- dSnbinom(x=0:150, size=ktest, mu=mutest, prob=NULL, 
                      method = "furman", verbose=FALSE, log=FALSE)})
pr_approximateObservations <- dSnbinom(0:150, size=ktest, mu=mutest, 
                                       method = "approximate.randomobservations")

plot(table(test), xlab="N", ylab="Density", las=1, bty="n", ylim=c(0, 4000), xlim=c(0, 150))
lines(0:150, pr_vellaisamy*nr, col="red")
lines(0:150, pr_furman*nr, col="blue")
lines(0:150, pr_approximateObservations*nr, col="green")

dSnbinom(x=42, size=ktest, mu=mutest, prob=NULL, 
         method = "vellaisamy&upadhye", verbose=TRUE)
as.numeric(dSnbinom(x=42, size=ktest, mu=mutest, prob=NULL, 
           method = "Furman", verbose=TRUE))
dSnbinom(x=42, size=ktest, mu=mutest, prob=NULL, 
         method = "approximate.randomobservations", verbose=TRUE)

x <- 100
# Test for criteria of convergence
Pr_exact <- dSnbinom(x=x, size=ktest, mu=mutest, method="vellaisamy&upadhye", 
                   log=FALSE, verbose=TRUE)
Pr_Furman <- dSnbinom(x=x, size=ktest, mu=mutest, method="Furman", log=FALSE, 
                 verbose=TRUE)
Pr_exact;as.numeric(Pr_Furman)
plot(1:length(attributes(Pr_Furman)$Pk), 
     log10(abs(attributes(Pr_Furman)$Pk-Pr_exact)), type="l", xlab="Iterations", 
     ylab="Abs log10", bty="n", ylim=c(-100, 0))
lines(1:(length(attributes(Pr_Furman)$Pk)-1), 
      log10(abs(diff(attributes(Pr_Furman)$Pk))), col="red")
legend("bottomright", legend=c("Log10 Convergence to true value", "Log10 Rate of change"), 
       col=c("black", "red"), 
       lty=1)
       
     
# example to fit a distribution
data <- rnbinom(1000, size=1, mu=10)
hist(data)
ag <- rep(1:100, 10)
r <- aggregate(data, by=list(ag), FUN=sum)
hist(r[,2])

parx <- c(size=1, mu=10)

dSnbinomx <- function(x, par) {
  -sum(dSnbinom(x=x[,2], mu=rep(par["mu"], 10), size=par["size"], log=TRUE))
}

fit_mu_size <- optim(par = parx, fn=dSnbinomx, x=r, method="BFGS", control=c(trace=TRUE))
fit_mu_size$par

alpha <- c(2.1, 2.05, 2)
mu <- c(10, 30, 20)
p <- pSnbinom(q=10, size=alpha, mu=mu, lower.tail = TRUE)

alpha <- c(2.1, 2.05, 2)
mu <- c(10, 30, 20)
q <- qSnbinom(p=0.1, size=alpha, mu=mu, lower.tail = TRUE)

alpha <- c(2.1, 2.05, 2)
mu <- c(10, 30, 20)
rep <- 100000
distEmpirique <- rSnbinom(n=rep, size=alpha, mu=mu)
tabledistEmpirique <- rep(0, 301)
names(tabledistEmpirique) <- as.character(0:300)
tabledistEmpirique[names(table(distEmpirique))] <- table(distEmpirique)/rep

plot(0:300, dSnbinom(0:300, size=alpha, mu=mu), type="h", bty="n", 
   xlab="x", ylab="Density", ylim=c(0,0.02))
plot_add(0:300, tabledistEmpirique, type="l", col="red")
legend(x=200, y=0.02, legend=c("Empirical", "Theoretical"), 
   text.col=c("red", "black"), bty="n")
   
   
# Test if saddlepoint approximation must be normalized
# Yes, it must be
n <- 7
alpha <- 1:n
p <- (1:n)/10
dSnbinom(x=10, prob=p, size=alpha, method="saddlepoint", log=FALSE,  
                 verbose=TRUE)
dSnbinom(x=10, prob=p, size=alpha, method="saddlepoint", log=FALSE,  
                 verbose=TRUE, normalize=FALSE)
                 
# Test for saddlepoint when x=0
n <- 7
alpha <- 1:n
p <- (1:n)/10
dSnbinom(x=0, prob=p, size=alpha, method="saddlepoint", log=FALSE,  
                 verbose=TRUE)
dSnbinom(x=1, prob=p, size=alpha, method="saddlepoint", log=FALSE,  
                 verbose=TRUE)
dSnbinom(x=c(0, 1), prob=p, size=alpha, method="saddlepoint", log=FALSE,  
                 verbose=TRUE)
dSnbinom(x=c(0, 1), prob=p, size=alpha, method="saddlepoint", log=FALSE,  
                 verbose=FALSE)
                 
# Test when prob are all the same
p <- rep(0.2, 7)
n <- 7
alpha <- 1:n
dSnbinom(x=0:10, prob=p, size=alpha, method="saddlepoint", log=FALSE,  
                 verbose=TRUE)
dSnbinom(x=0:10, prob=p, size=alpha, method="furman", log=FALSE,  
                 verbose=TRUE)
dSnbinom(x=0:10, prob=p, size=alpha, method="exact", log=FALSE,  
                 verbose=TRUE)
                 
# Test when n=1
p <- 0.2
n <- 1
alpha <- 1:n
dSnbinom(x=0:10, prob=p, size=alpha, method="saddlepoint", log=FALSE,  
                 verbose=TRUE)
dSnbinom(x=0:10, prob=p, size=alpha, method="furman", log=FALSE,  
                 verbose=TRUE)
dSnbinom(x=0:10, prob=p, size=alpha, method="exact", log=FALSE,  
                 verbose=TRUE)
                 

## End(Not run)