Nothing
#
# rmutil : A Library of Special Functions for Repeated Measurements
# Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public Licence as published by
# the Free Software Foundation; either version 2 of the Licence, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public Licence for more details.
#
# You should have received a copy of the GNU General Public Licence
# along with this program; if not, write to the Free Software
# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
# SYNOPSIS
#
# pinvgauss(q, m, s)
# dinvgauss(y, m, s, log=FALSE)
# qinvgauss(p, m, s)
# rinvgauss(n, m, s)
# plaplace(q, m=0, s=1)
# dlaplace(y, m=0, s=1)
# qlaplace(p, m=0, s=1)
# rlaplace(n, m=0, s=1)
# plevy(q, m=0, s=1)
# dlevy(y, m=0, s=1)
# qlevy(p, m=0, s=1)
# rlevy(n, m=0, s=1)
# ppareto(q, m, s)
# dpareto(y, m, s)
# qpareto(p, m, s)
# rpareto(n, m, s)
# psimplex(q, m, s)
# dsimplex(y, m, s)
# qsimplex(p, m, s)
# rsimplex(n, m, s)
# ptwosidedpower(q, m, s=2)
# dtwosidedpower(y, m, s=2)
# qtwosidedpower(p, m, s=2)
# rtwosidedpower(n, m, s=2)
#
# pboxcox(q, m=0, s=1, f=1)
# dboxcox(y, m=0, s=1, f=1)
# qboxcox(p, m=0, s=1, f=1)
# rboxcox(n, m=0, s=1, f=1)
# pburr(q, m, s, f)
# dburr(y, m, s, f)
# qburr(p, m, s, f)
# rburr(n, m, s, f)
# pgextval(q, s, m, f)
# dgextval(y, s, m, f)
# qgextval(p, s, m, f)
# rgextval(n, s, m, f)
# pggamma(q, s, m, f)
# dggamma(y, s, m, f)
# qggamma(p, s, m, f)
# rggamma(n, s, m, f)
# pginvgauss(q, m, s, f)
# dginvgauss(y, m, s, f)
# qginvgauss(p, m, s, f)
# rginvgauss(n, m, s, f)
# pglogis(q, m=0, s=1, f=1)
# dglogis(y, m=0, s=1, f=1)
# qglogis(p, m=0, s=1, f=1)
# rglogis(n, m=0, s=1, f=1)
# pgweibull(q, s, m, f)
# dgweibull(y, s, m, f)
# qgweibull(p, s, m, f)
# rgweibull(n, s, m, f)
# phjorth(q, m, s, f)
# dhjorth(y, m, s, f)
# qhjorth(p, m, s, f)
# rhjorth(n, m, s, f)
# ppowexp(q, m=0, s=1, f=1)
# dpowexp(y, m=0, s=1, f=1)
# qpowexp(p, m=0, s=1, f=1)
# rpowexp(n, m=0, s=1, f=1)
# pskewlaplace(q, m=0, s=1, f=1)
# dskewlaplace(y, m=0, s=1, f=1)
# qskewlaplace(p, m=0, s=1, f=1)
# rskewlaplace(n, m=0, s=1, f=1)
#
# pbetabinom(q, size, m, s)
# dbetabinom(y, size, m, s)
# qbetabinom(p, size, m, s)
# rbetabinom(n, size, m, s)
# pdoublebinom(q, size, m, s)
# ddoublebinom(y, size, m, s)
# qdoublebinom(p, size, m, s)
# rdoublebinom(n, size, m, s)
# pmultbinom(q, size, m, s)
# dmultbinom(y, size, m, s)
# qmultbinom(p, size, m, s)
# rmultbinom(n, size, m, s)
# pdoublepois(q, m, s)
# ddoublepois(y, m, s)
# qdoublepois(p, m, s)
# rdoublepois(n, m, s)
# pmultpois(q, m, s)
# dmultpois(y, m, s)
# qmultpois(p, m, s)
# rmultpois(n, m, s)
# ppvfpois(q, m, s, f)
# dpvfpois(y, m, s, f)
# qpvfpois(p, m, s, f)
# rpvfpois(n, m, s, f)
# pgammacount(q, m, s)
# dgammacount(y, m, s)
# qgammacount(p, m, s)
# rgammacount(n, m, s)
# pconsul(q, m, s)
# dconsul(y, m, s)
# qconsul(p, m, s)
# rconsul(n, m, s)
#
# DESCRIPTION
#
# Functions to compute the probability and cumulative probability
# functions for
# continuous two parameter distributions:
# inverse Gaussian, Laplace, Levy, Pareto, simplex
# continuous three parameter distributions:
# Box-Cox, Burr, generalized extreme value, generalized gamma, generalized
# inverse Gaussian, generalized logistic, generalized Weibull, Hjorth,
# power exponential, skew Laplace
# discrete two parameter distributions:
# double Poisson, multiplicative Poisson, double binomial, Consul,
# multiplicative binomial, beta binomial
### continuous two-parameter distributions
###
### inverse Gaussian distribution
###
pinvgauss <- function(q, m, s){
if(any(q<=0))stop("q must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
t <- q/m
v <- sqrt(q*s)
pnorm((t-1)/v)+exp(2/(m*s))*pnorm(-(t+1)/v)}
dinvgauss <- function(y, m, s, log=FALSE){
if(any(y<=0))stop("y must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
tmp <- -(y-m)^2/(2*y*s*m^2)-(log(2*pi*s)+3*log(y))/2
if(!log)tmp <- exp(tmp)
tmp}
qinvgauss <- function(p, m, s){
h <- function(y){
t <- y/m[i]
v <- sqrt(y*s[i])
pnorm((t-1)/v)+exp(2/(m[i]*s[i]))*pnorm(-(t+1)/v)-p[i]}
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0))stop("m must be positive")
if(any(s<0))stop("s must be positive")
len <- max(length(p),length(m),length(s))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
interval <- c(.Machine$double.xmin,20)
while(h(interval[1])*h(interval[2])>0)interval <- 2*interval
tmp[i] <- uniroot(h,interval)$root}
tmp}
rinvgauss <- function(n=1, m, s) qinvgauss(runif(n),m=m,s=s)
### Laplace distribution
###
plaplace <- function(q, m=0, s=1){
if(any(s<=0))stop("s must be positive")
u <- (q-m)/s
t <- exp(-abs(u))/2
ifelse(u<0,t,1-t)}
dlaplace <- function(y, m=0, s=1, log=FALSE){
if(any(s<=0))stop("s must be positive")
tmp <- -abs(y-m)/s-log(2*s)
if(!log)tmp <- exp(tmp)
tmp}
qlaplace <- function(p, m=0, s=1){
h <- function(y){
#u <- (y-m[i])/s[i]
#t <- exp(-abs(u))/2
#ifelse(u<0,t,1-t)-p[i]
## bruce edit:
u <- (y-m)/s
t <- exp(-abs(u))/2
ifelse(u<0,t,1-t)-p
}
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(s<=0))stop("s must be positive")
ifelse(p<0.5,s*log(2*p)+m,-s*log(2*(1-p))+m)}
rlaplace <- function(n = 1, m = 0, s = 1){
if(any(s<=0))stop("s must be positive")
q <- runif(n)
ifelse(q<0.5,s*log(2*q)+m,-s*log(2*(1-q))+m)}
### Levy distribution
###
plevy <- function(q, m=0, s=1){
if(any(q<=m))stop("some y <= m")
if(any(s<=0))stop("s must be positive")
2*(1-pnorm(1/sqrt((q-m)/s)))}
dlevy <- function(y, m=0, s=1, log=FALSE){
if(any(y<=m))stop("some y <= m")
if(any(s<=0))stop("s must be positive")
#sqrt(s/(2*pi*(y-m)^3))*exp(-s/(2*(y-m)))}
tmp <- log(s/(2*pi))/2-3*log(y-m)/2-s/(2*(y-m))
if(!log)tmp <- exp(tmp)
tmp}
qlevy <- function(p, m=0, s=1){
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(s<0))stop("s must be positive")
s/qnorm(1-p/2)^2+m}
rlevy <- function(n=1, m=0, s=1) {
if(any(s<0))stop("s must be positive")
s/qnorm(1-runif(n)/2)^2+m}
### Pareto distribution
###
ppareto <- function(q, m, s){
if(any(q<0))stop("q must be >= 0")
if(any(m<=0))stop("m must be positive")
if(any(s<=1))stop("s must be > 1")
1-(1+q/(m*(s-1)))^-s}
dpareto <- function(y, m, s, log=FALSE){
if(any(y<0))stop("y must be >= 0")
if(any(m<=0))stop("m must be positive")
if(any(s<=1))stop("s must be > 1")
m <- m*(s-1)
tmp <- log(s)-(s+1)*log(1+y/m)-log(m)
if(!log)tmp <- exp(tmp)
tmp}
qpareto <- function(p, m, s){
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0))stop("m must be positive")
if(any(s<=1))stop("s must be >1")
((1-p)^(-1/s)-1)*m*(s-1)}
rpareto <- function(n=1, m, s){
if(any(m<=0))stop("m must be positive")
if(any(s<=1))stop("s must be >1")
((1-runif(n))^(-1/s)-1)*m*(s-1)}
### simplex distribution
###
psimplex <- function(q, m, s){
if(any(q<=0)||any(q>=1))stop("q must contain values between 0 and 1")
if(any(m<=0)||any(m>=1))stop("m must contain values between 0 and 1")
if(any(s<=0))stop("s must be positive")
len <- max(length(q),length(m),length(s))
if(length(q)!=len){
if(length(q)==1)q <- rep(q,len)
else stop("length of q incorrect")}
if(length(m)!=len){
if(length(m)!=1)stop("m has incorrect length")
else m <- rep(m,len)}
if(length(s)!=len){
if(length(s)!=1)stop("s has incorrect length")
else s <- rep(s,len)}
z <- .C("psimplex_c",
as.double(q),
as.double(m),
as.double(s),
as.double(1),
len=as.integer(len),
eps=as.double(1.0e-6),
pts=as.integer(5),
max=as.integer(16),
err=integer(1),
res=double(len),
## DUP=FALSE,
PACKAGE="rmutil")
if(z$err==1)warning("Unable to allocate memory for integration")
if(z$err==2)warning("Division by zero in integration")
else if(z$err==3)warning("No convergence in integration")
z$res}
dsimplex <- function(y, m, s, log=FALSE){
if(any(y<=0)||any(y>=1))stop("y must contain values between 0 and 1")
if(any(m<=0)||any(m>=1))stop("m must contain values between 0 and 1")
if(any(s<=0))stop("s must be positive")
tmp <- -((y-m)/(m*(1-m)))^2/(2*y*(1-y)*s)-(log(2*pi*s)+3*(log(y)+log(1-y)))/2
if(!log)tmp <- exp(tmp)
tmp}
qsimplex <- function(p, m, s){
h <- function(y).C("psimplex_c",
as.double(y),
as.double(m[i]),
as.double(s[i]),
as.double(1),
len=as.integer(1),
eps=as.double(1.0e-6),
pts=as.integer(5),
max=as.integer(16),
err=integer(1),
res=double(1),
## DUP=FALSE,
PACKAGE="rmutil")$res-p[i]
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0)||any(m>=1))stop("m must contain values between 0 and 1")
if(any(s<0))stop("s must be positive")
len <- max(length(p),length(m),length(s))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
interval <- c(0.000001,0.999999)
while(h(interval[1])*h(interval[2])>0)
interval <- c(interval[1]/10,1-interval[1]/10)
tmp[i] <- uniroot(h,interval)$root}
tmp}
rsimplex <- function(n=1, m, s) qsimplex(runif(n),m=m,s=s)
### two-sided power distribution
###
ptwosidedpower <- function(q, m, s=2){
if(any(q<=0)||any(q>=1))stop("q must contain values between 0 and 1")
if(any(m<=0)||any(m>=1))stop("m must contain values between 0 and 1")
if(any(s<=0))stop("s must be positive")
ifelse(q<m,m*(q/m)^s,1-(1-m)*((1-q)/(1-m))^s)}
dtwosidedpower <- function(y, m, s=2, log=FALSE){
if(any(y<=0)||any(y>=1))stop("y must contain values between 0 and 1")
if(any(m<=0)||any(m>=1))stop("m must contain values between 0 and 1")
if(any(s<=0))stop("s must be positive")
tmp <- log(s)+(s-1)*ifelse(y<m,log(y)-log(m),log(1-y)-log(1-m))
if(!log)tmp <- exp(tmp)
tmp}
qtwosidedpower <- function(p, m, s=2){
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0)||any(m>=1))stop("m must contain values between 0 and 1")
if(any(s<0))stop("s must be positive")
ifelse(p<m,m*(p/m)^(1/s),1-(1-m)*((1-p)/(1-m))^(1/s))}
rtwosidedpower <- function(n, m, s=2) qtwosidedpower(runif(n),m=m,s=s)
### continuous three-parameter distributions
###
### Box-Cox distribution
### normed to make it a real distribution with y > 0
###
pboxcox <- function(q, m, s=1, f=1){
if(any(q<=0))stop("q must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
norm <- sign(f)*pnorm(0,m,sqrt(s))
(pnorm(q^f/f,m,sqrt(s))-(f>0)*norm)/(1-(f<0)-norm)}
dboxcox <- function(y, m, s=1, f=1, log=FALSE){
if(any(y<=0))stop("y must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
norm <- sign(f)*pnorm(0,m,sqrt(s))
tmp <- (f-1)*log(y)+dnorm(y^f/f,m,sqrt(s),log=TRUE)-log(1-(f<0)-norm)
if(!log)tmp <- exp(tmp)
tmp}
qboxcox <- function(p, m, s=1, f=1){
h <- function(y){
norm <- sign(f[i])*pnorm(0,m[i],sqrt(s[i]))
(pnorm(y^f[i]/f[i],m[i],sqrt(s[i]))-(f[i]>0)*norm)/
(1-(f[i]<0)-norm)-p[i]}
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0))stop("m must be positive")
if(any(s<0))stop("s must be positive")
len <- max(length(p),length(m),length(s),length(f))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
if(length(f)!=len){
if(length(f)==1)f <- rep(f,len)
else stop("length of f incorrect")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
interval <- c(.Machine$double.xmin,20)
while(h(interval[1])*h(interval[2])>0)interval <- 2*interval
tmp[i] <- uniroot(h,interval)$root}
tmp}
rboxcox <- function(n=1, m, s=1, f=1) qboxcox(runif(n),m=m,s=s,f=f)
### Burr distribution
###
pburr <- function(q, m, s, f){
if(any(q<=0))stop("q must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
1-(1+(q/m)^s)^-f}
dburr <- function(y, m, s, f, log=FALSE){
if(any(y<=0))stop("y must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
y1 <- y/m
tmp <- log(f*s)+(s-1)*log(y1)-log(m)-(f+1)*log(1+y1^s)
if(!log)tmp <- exp(tmp)
tmp}
qburr <- function(p, m, s, f){
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
((1-p)^(-1/f)-1)^(1/s)*m}
rburr <- function(n=1, m, s, f){
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
((1-runif(n))^(-1/f)-1)^(1/s)*m}
### generalized extreme value distribution
### normed to make it a real distribution with y > 0
###
pgextval <- function(q, s, m, f){
if(any(q<=0))stop("q must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
norm <- sign(f)*exp(-m^-s)
ind <- f>0
(pweibull(exp(q^f/f),s,m)-ind+ind*norm)/(1-ind+norm)}
dgextval <- function(y, s, m, f, log=FALSE){
if(any(y<=0))stop("y must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
norm <- sign(f)*exp(-m^-s)
y1 <- exp(y^f/f)
tmp <- (f-1)*log(y)+log(y1)+dweibull(y1,s,m,log=TRUE)-log(1-(f>0)+norm)
if(!log)tmp <- exp(tmp)
tmp}
qgextval <- function(p, s, m, f){
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
norm <- sign(f)*exp(-m^-s)
ind <- f>0
(f*log(qweibull(p*(1-ind+norm)+ind-ind*norm,s,m)))^(1/f)}
rgextval <- function(n=1, s, m, f){
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
norm <- sign(f)*exp(-m^-s)
ind <- f>0
(f*log(qweibull(runif(n)*(1-ind+norm)+ind-ind*norm,s,m)))^(1/f)}
### generalized gamma distribution
###
pggamma <- function(q, s, m, f){
if(any(q<=0))stop("q must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
pgamma(q^f,s,scale=(m/s)^f)}
dggamma <- function(y, s, m, f, log=FALSE){
if(any(y<=0))stop("y must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
tmp <- log(f)+(f-1)*log(y)+dgamma(y^f,s,scale=(m/s)^f,log=TRUE)
if(!log)tmp <- exp(tmp)
tmp}
qggamma <- function(p, s, m, f) {
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
qgamma(p,s,scale=(m/s)^f)^(1/f)}
rggamma <- function(n=1, s, m, f){
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
qgamma(runif(n),s,scale=(m/s)^f)^(1/f)}
### generalized inverse Gaussian distribution
###
pginvgauss <- function(q, m, s, f){
if(any(q<=0))stop("q must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
len <- max(length(q),length(m),length(s))
if(length(q)!=len){
if(length(q)==1)q <- rep(q,len)
else stop("length of q incorrect")}
if(length(m)!=len){
if(length(m)!=1)stop("m has incorrect length")
else m <- rep(m,len)}
if(length(s)!=len){
if(length(s)!=1)stop("s has incorrect length")
else s <- rep(s,len)}
if(length(f)!=len){
if(length(f)!=1)stop("f has incorrect length")
else f <- rep(f,len)}
z <- .C("pginvgauss_c",
as.double(q),
as.double(m),
as.double(s),
as.double(f),
len=as.integer(len),
eps=as.double(1.0e-6),
pts=as.integer(5),
max=as.integer(16),
err=integer(1),
res=double(len),
## DUP=FALSE,
PACKAGE="rmutil")
if(z$err==1)warning("Unable to allocate memory for integration")
if(z$err==2)warning("Division by zero in integration")
else if(z$err==3)warning("No convergence in integration")
z$res}
dginvgauss <- function(y, m, s, f, log=FALSE){
if(any(y<=0))stop("y must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
tmp <- (f-1)*log(y)-(1/y+y/m^2)/(2*s)-f*log(m)-log(2*besselK(1/(s*m),abs(f)))
if(!log)tmp <- exp(tmp)
tmp}
qginvgauss <- function(p, m, s, f){
h <- function(y).C("pginvgauss_c",
as.double(y),
as.double(m[i]),
as.double(s[i]),
as.double(f[i]),
len=as.integer(1),
eps=as.double(1.0e-6),
pts=as.integer(5),
max=as.integer(16),
err=integer(1),
res=double(1),
## DUP=FALSE,
PACKAGE="rmutil")$res-p[i]
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
len <- max(length(p),length(m),length(s),length(f))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
if(length(f)!=len){
if(length(f)==1)f <- rep(f,len)
else stop("length of f incorrect")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
interval <- c(.Machine$double.xmin,20)
while(h(interval[1])*h(interval[2])>0)interval <- 2*interval
tmp[i] <- uniroot(h,interval)$root}
tmp}
rginvgauss <- function(n=1, m, s, f) qginvgauss(runif(n),m=m,s=s,f=f)
### generalized logistic distribution
###
pglogis <- function(q, m=0, s=1, f=1){
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
(1+exp(-sqrt(3)*(q-m)/(s*pi)))^-f}
dglogis <- function(y, m=0, s=1, f=1, log=FALSE) {
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
y1 <- exp(-sqrt(3)*(y-m)/(s*pi))
tmp <- log(3)/2+log(f*y1)-log(pi*s)-(f+1)*log(1+y1)
if(!log)tmp <- exp(tmp)
tmp}
qglogis <- function(p, m=0, s=1, f=1){
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
-log(p^(-1/f)-1)*s*pi/sqrt(3)+m}
rglogis <- function(n=1, m=0, s=1, f=1){
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
-log(runif(n)^(-1/f)-1)*s*pi/sqrt(3)+m}
### generalized Weibull distribution
###
pgweibull <- function(q, s, m, f){
if(any(q<=0))stop("q must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
(1-exp(-(q/m)^s))^f}
dgweibull <- function(y, s, m, f, log=FALSE){
if(any(y<=0))stop("y must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
y1 <- exp(-(y/m)^s)
tmp <- log(s*f)+(s-1)*log(y)+(f-1)*log(1-y1)+log(y1)-s*log(m)
if(!log)tmp <- exp(tmp)
tmp}
qgweibull <- function(p, s, m, f){
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
m*(-log(1-p^(1/f)))^(1/s)}
rgweibull <- function(n=1, s, m, f){
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
m*(-log(1-runif(n)^(1/f)))^(1/s)}
### Hjorth distribution
###
phjorth <- function(q, m, s, f){
if(any(q<=0))stop("q must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
1-(1+s*q)^(-f/s)*exp(-(q/m)^2/2)}
dhjorth <- function(y, m, s, f, log=FALSE){
if(any(y<=0))stop("y must contain positive values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
tmp <- -(f/s)*log(1+s*y)-(y/m)^2/2+log(y/m^2+f/(1+s*y))
if(!log)tmp <- exp(tmp)
tmp}
qhjorth <- function(p, m, s, f){
h <- function(y) 1-(1+s[i]*y)^(-f[i]/s[i])*exp(-(y/m[i])^2/2)-p[i]
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0))stop("m must be positive")
if(any(s<0))stop("s must be positive")
len <- max(length(p),length(m),length(s),length(f))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
if(length(f)!=len){
if(length(f)==1)f <- rep(f,len)
else stop("length of f incorrect")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
interval <- c(.Machine$double.xmin,20)
while(h(interval[1])*h(interval[2])>0)interval <- 2*interval
tmp[i] <- uniroot(h,interval)$root}
tmp}
rhjorth <- function(n=1, m, s, f) qhjorth(runif(n),m=m,s=s,f=f)
### power exponential distribution
###
ppowexp <- function(q, m=0, s=1, f=1){
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
len <- max(length(q),length(m),length(s))
if(length(q)!=len){
if(length(q)==1)q <- rep(q,len)
else stop("length of q incorrect")}
if(length(m)!=len){
if(length(m)!=1)stop("m has incorrect length")
else m <- rep(m,len)}
if(length(s)!=len){
if(length(s)!=1)stop("s has incorrect length")
else s <- rep(s,len)}
if(length(f)!=len){
if(length(f)!=1)stop("f has incorrect length")
else f <- rep(f,len)}
z <- .C("ppowexp_c",
as.double(q),
as.double(m),
as.double(s),
as.double(f),
len=as.integer(len),
eps=as.double(1.0e-6),
pts=as.integer(5),
max=as.integer(16),
err=integer(1),
res=double(len),
## DUP=FALSE,
PACKAGE="rmutil")
if(z$err==1)warning("Unable to allocate memory for integration")
if(z$err==2)warning("Division by zero in integration")
else if(z$err==3)warning("No convergence in integration")
ifelse(q-m>0,0.5+z$res,0.5-z$res)}
dpowexp <- function(y, m=0, s=1, f=1, log=FALSE){
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
s <- sqrt(s)
b <- 1+1/(2*f)
tmp <- -(abs(y-m)/s)^(2*f)/2-log(s)-lgamma(b)-b*log(2)
if(!log)tmp <- exp(tmp)
tmp}
qpowexp <- function(p, m=0, s=1, f=1){
h <- function(y) {
z <- .C("ppowexp_c",
as.double(y),
as.double(m[i]),
as.double(s[i]),
as.double(f[i]),
len=as.integer(1),
eps=as.double(1.0e-6),
pts=as.integer(5),
max=as.integer(16),
err=integer(1),
res=double(1),
## DUP=FALSE,
PACKAGE="rmutil")$res
if(y-m[i]>0) 0.5+z-p[i] else 0.5-z-p[i]}
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(s<0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
len <- max(length(p),length(m),length(s),length(f))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
if(length(f)!=len){
if(length(f)==1)f <- rep(f,len)
else stop("length of f incorrect")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
interval <- m[i]+s[i]*c(-2,2)
while(h(interval[1])*h(interval[2])>0)interval <- 2*interval
tmp[i] <- uniroot(h,interval)$root}
tmp}
rpowexp <- function(n=1, m=0, s=1, f=1) qpowexp(runif(n),m=m,s=s,f=f)
### skew Laplace distribution
###
pskewlaplace <- function(q, m=0, s=1, f=1){
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
u <- (q-m)/s
ifelse(u>0,1-exp(-f*abs(u))/(1+f^2),f^2*exp(-abs(u)/f)/(1+f^2))}
dskewlaplace <- function(y, m=0, s=1, f=1, log=FALSE){
if(any(s<=0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
tmp <- log(f)+ifelse(y>m,-f*(y-m),(y-m)/f)/s-log((1+f^2)*s)
if(!log)tmp <- exp(tmp)
tmp}
qskewlaplace <- function(p, m=0, s=1, f=1){
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(s<0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
ifelse(p<0.5,f*s*log((1+f^2)*p/f^2)+m,-s*log((1+f^2)*(1-p))/f+m)}
rskewlaplace <- function(n=1, m=0, s=1, f=1){
if(any(s<0))stop("s must be positive")
if(any(f<=0))stop("f must be positive")
q <- runif(n)
ifelse(q<0.5,f*s*log((1+f^2)*q/f^2)+m,-s*log((1+f^2)*(1-q))/f+m)}
### discrete (overdispersed) two-parameter distributions
###
### beta-binomial distribution
###
pbetabinom <- function(q, size, m, s){
##BEGIN from @hennerw in issue #5
# if (any(q > size)){
# # Updating to correctly deal with this siduation
# # stop("q must be <= size")
# if(all(q>size)) return(rep(1,len))
# Val <- q<= size
# out <- rep(1,len)
# out[Val] <- pbetabinom_c(q[Val],size[Val],m[Val],s[Val])
# return(out)
# }
# if (any(q < 0)){
# # stop("q must contain non-negative values")
# if(all(q < 0)) return(rep(0,len))
# Val <- q>=0
# out <- rep(0,len)
# out[Val] <- pbetabinom_c(q[Val],size[Val],m[Val],s[Val])
# return(out)
# }
##END from @hennerw in issue #5
if(any(q<0)){message("Negative values of q detected. `pbetabinom` returns 0 for such values.")}#; q[q<0] <- 0}
if(any(size<0))stop("size must contain non-negative values")
if(any(m<=0)||any(m>=1))stop("m must lie between 0 and 1")
if(any(s<=0))stop("s must be positive")
len <- max(length(q),length(m),length(s),length(size))
if(length(q)!=len){
if(length(q)==1)q <- rep(q,len)
else stop("length of q incorrect")}
if(length(size)!=len){
if(length(size)==1)size <- rep(size,len)
else stop("size must be the same length as q")}
if(any(q>size)){message("Elements of q that were greater than size detected. `pbetabinom` returns 1 for such values.")}#; q[q>size] <- size[q>size]}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("m and q must have the same length")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("s and q must have the same length")}
t <- s*m
u <- s*(1-m)
res <- vector("numeric",length(q))
for(i in 1:length(q)){
if(q[i] < 0){ res[i] <- 0
} else if(q[i]>size[i]){ res[i] <- 1
}else{
qq <- 0:q[i]
res[i] <- sum(exp(lbeta(qq+t[i],size[i]-qq+u[i])-
lbeta(t[i],u[i])+lchoose(size[i],qq)))
}
}
res}
dbetabinom <- function(y, size, m, s, log=FALSE){
if(any(y<0)){message("Negative values of y replaced with -1 to yield 0 value density"); y[y<0]<- -1}
if(any(size<0))stop("size must contain non-negative values")
if(any(m<=0)||any(m>=1))stop("m must lie between 0 and 1")
if(any(s<=0))stop("s must be positive")
ly <- max(length(y),length(m),length(s),length(size))
if(length(y)!=ly){
if(length(y)==1)y <- rep(y,ly)
else stop("length of y incorrect")}
if(length(size)!=ly){
if(length(size)==1)size <- rep(size,ly)
else stop("size must be the same length as y")}
if(any(y>size)){message("Elements of y exist that are greater than size");}
if(length(m)!=ly){
if(length(m)==1)m <- rep(m,ly)
else stop("m and y must have the same length")}
if(length(s)!=ly){
if(length(s)==1)s <- rep(s,ly)
else stop("s and y must have the same length")}
t <- s*m
u <- s*(1-m)
tmp <- lbeta(y+t,size-y+u)-lbeta(t,u)+lchoose(size,y)
if(!log)tmp <- exp(tmp)
tmp}
qbetabinom <- function(p, size, m, s){
h <- function(y){
t <- s[i]*m[i]
u <- s[i]*(1-m[i])
pp <- 0:y
sum(exp(lbeta(pp+t,size[i]-pp+u)-
lbeta(t,u)+lchoose(size[i],pp)))-p[i]}
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0)||any(m>=1))stop("m must lie between 0 and 1")
if(any(s<0))stop("s must be positive")
len <- max(length(p),length(m),length(s),length(size))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(size)!=len){
if(length(size)==1)size <- rep(size,len)
else stop("length of size incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
interval <- c(0,size[i])
tmp[i] <- if(h(interval[1])*h(interval[2])>0)0
else uniroot(h,interval)$root}
round(tmp)}
rbetabinom <- function(n=1, size, m, s) qbetabinom(runif(n),size=size,m=m,s=s)
### double binomial distribution
###
pdoublebinom <- function(q, size, m, s){
if(any(q<0))stop("q must contain non-negative values")
if(any(size<0))stop("n must contain non-negative values")
if(any(m<=0)||any(m>=1))stop("m must lie between 0 and 1")
if(any(s<=0))stop("s must be positive")
len <- max(length(q),length(m),length(s),length(size))
if(length(q)!=len){
if(length(q)==1)q <- rep(q,len)
else stop("length of q incorrect")}
if(length(size)!=len){
if(length(size)==1)size <- rep(size,len)
else stop("size must be the same length as q")}
if(any(q>size))stop("q must be <= size")
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("m and q must have the same length")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("s and q must have the same length")}
.C("pdb",
as.integer(q),
as.integer(size),
as.double(m),
as.double(s),
as.integer(length(q)),
res=double(length(q)),
## DUP=FALSE,
PACKAGE="rmutil")$res}
ddoublebinom <- function(y, size, m, s, log=FALSE){
if(any(y<0))stop("y must contain non-negative values")
if(any(size<0))stop("size must contain non-negative values")
if(any(m<=0)||any(m>=1))stop("m must lie between 0 and 1")
if(any(s<=0))stop("s must be positive")
ly <- max(length(y),length(m),length(s),length(size))
if(length(y)!=ly){
if(length(y)==1)y <- rep(y,ly)
else stop("length of y incorrect")}
if(length(size)!=ly){
if(length(size)==1)size <- rep(size,ly)
else stop("size must be the same length as y")}
if(any(y>size))stop("y must be <= size")
if(length(m)!=ly){
if(length(m)==1)m <- rep(m,ly)
else stop("m and y must have the same length")}
if(length(s)!=ly){
if(length(s)==1)s <- rep(s,ly)
else stop("s and y must have the same length")}
tmp <- .C("ddb",
as.integer(y),
as.integer(size),
as.double(m),
as.double(s),
as.integer(length(y)),
as.double(rep(1,length(y))),
res=double(length(y)),
## DUP=FALSE,
PACKAGE="rmutil")$res
if(!log)tmp <- exp(tmp)
tmp}
qdoublebinom <- function(p, size, m, s){
h <- function(y) .C("pdb",
as.integer(y),
as.integer(size[i]),
as.double(m[i]),
as.double(s[i]),
as.integer(1),
res=double(1),
## DUP=FALSE,
PACKAGE="rmutil")$res-p[i]
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0)||any(m>=1))stop("m must lie between 0 and 1")
if(any(s<=0))stop("s must be positive")
len <- max(length(p),length(m),length(s),length(size))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(size)!=len){
if(length(size)==1)size <- rep(size,len)
else stop("length of size incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
interval <- c(0,size[i])
tmp[i] <- if(h(interval[1])*h(interval[2])>0)0
else uniroot(h,interval)$root}
round(tmp)}
rdoublebinom <- function(n=1, size, m, s)
qdoublebinom(runif(n),size=size,m=m,s=s)
### multiplicative binomial distribution
###
pmultbinom <- function(q, size, m, s){
if(any(q<0))stop("q must contain non-negative values")
if(any(size<0))stop("size must contain non-negative values")
if(any(m<=0)||any(m>=1))stop("m must lie between 0 and 1")
if(any(s<=0))stop("s must be positive")
len <- max(length(q),length(m),length(s),length(size))
if(length(q)!=len){
if(length(q)==1)q <- rep(q,len)
else stop("length of q incorrect")}
if(length(size)!=len){
if(length(size)==1)size <- rep(size,len)
else stop("size must be the same length as q")}
if(any(q>size))stop("q must be <= size")
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("m and q must have the same length")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("s and q must have the same length")}
.C("pmb",
as.integer(q),
as.integer(size),
as.double(m),
as.double(s),
as.integer(length(q)),
res=double(length(q)),
## DUP=FALSE,
PACKAGE="rmutil")$res}
dmultbinom <- function(y, size, m, s, log=FALSE){
if(any(y<0))stop("y must contain non-negative values")
if(any(size<0))stop("size must contain non-negative values")
if(any(m<=0)||any(m>=1))stop("m must lie between 0 and 1")
if(any(s<=0))stop("s must be positive")
ly <- max(length(y),length(m),length(s),length(size))
if(length(y)!=ly){
if(length(y)==1)y <- rep(y,ly)
else stop("length of y incorrect")}
if(length(size)!=ly){
if(length(size)==1)size <- rep(size,ly)
else stop("size must be the same length as y")}
if(any(y>size))stop("y must be <= size")
if(length(m)!=ly){
if(length(m)==1)m <- rep(m,ly)
else stop("m and y must have the same length")}
if(length(s)!=ly){
if(length(s)==1)s <- rep(s,ly)
else stop("s and y must have the same length")}
tmp <- .C("dmb",
as.integer(y),
as.integer(size),
as.double(m),
as.double(s),
as.integer(length(y)),
as.double(rep(1,length(y))),
res=double(length(y)),
## DUP=FALSE,
PACKAGE="rmutil")$res
if(!log)tmp <- exp(tmp)
tmp}
qmultbinom <- function(p, size, m, s){
h <- function(y).C("pmb",
as.integer(y),
as.integer(size[i]),
as.double(m[i]),
as.double(s[i]),
as.integer(1),
res=double(1),
## DUP=FALSE,
PACKAGE="rmutil")$res-p[i]
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0)||any(m>=1))stop("m must lie between 0 and 1")
if(any(s<=0))stop("s must be positive")
len <- max(length(p),length(m),length(s),length(size))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(size)!=len){
if(length(size)==1)size <- rep(size,len)
else stop("length of size incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
interval <- c(0,size[i])
tmp[i] <- if(h(interval[1])*h(interval[2])>0)0
else uniroot(h,interval)$root}
round(tmp)}
rmultbinom <- function(n=1, size, m, s)
qmultbinom(runif(n),size=size,m=m,s=s)
### double Poisson distribution
###
pdoublepois <- function(q, m, s){
if(any(q<0))stop("q must contain non-negative values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
len <- max(length(q),length(m),length(s))
if(length(q)!=len){
if(length(q)==1)q <- rep(q,len)
else stop("length of q incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("m and q must have the same length")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("s and q must have the same length")}
.C("pdp",
as.integer(q),
as.integer(3*max(c(q,100))),
as.double(m),
as.double(s),
as.integer(length(q)),
res=double(length(q)),
## DUP=FALSE,
PACKAGE="rmutil")$res}
ddoublepois <- function(y, m, s, log=FALSE){
if(any(y<0))stop("y must contain non-negative values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
ly <- max(length(y),length(m),length(s))
if(length(y)!=ly){
if(length(y)==1)y <- rep(y,ly)
else stop("length of y incorrect")}
if(length(m)!=ly){
if(length(m)==1)m <- rep(m,ly)
else stop("m and y must have the same length")}
if(length(s)!=ly){
if(length(s)==1)s <- rep(s,ly)
else stop("s and y must have the same length")}
tmp <- .C("ddp",
as.integer(y),
as.integer(3*max(c(y,100))),
as.double(m),
as.double(s),
as.integer(length(y)),
as.double(rep(1,length(y))),
res=double(length(y)),
## DUP=FALSE,
PACKAGE="rmutil")$res
if(!log)tmp <- exp(tmp)
tmp}
qdoublepois <- function(p, m, s){
h <- function(y).C("pdp",
as.integer(y),
as.integer(3*max(c(y,100))),
as.double(m[i]),
as.double(s[i]),
as.integer(1),
res=double(1),
## DUP=FALSE,
PACKAGE="rmutil")$res-p[i]
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<0))stop("m must be positive")
if(any(s<0))stop("s must be positive")
len <- max(length(p),length(m),length(s))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
interval <- c(0,20)
if(h(interval[1])*h(interval[2])>0&&h(interval[1])>0)tmp[i] <- 0
else {
while(h(interval[1])*h(interval[2])>0)interval <- 2*interval
tmp[i] <- uniroot(h,interval)$root}}
round(tmp)}
rdoublepois <- function(n=1, m, s) qdoublepois(runif(n),m=m,s=s)
### multiplicative Poisson distribution
###
pmultpois <- function(q, m, s){
if(any(q<0))stop("q must contain non-negative values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0|s>1))stop("s must be positive and <= 1")
len <- max(length(q),length(m),length(s))
if(length(q)!=len){
if(length(q)==1)q <- rep(q,len)
else stop("length of q incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("m and q must have the same length")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("s and q must have the same length")}
.C("pmp",
as.integer(q),
as.integer(3*max(c(q,100))),
as.double(m),
as.double(s),
as.integer(length(q)),
res=double(length(q)),
## DUP=FALSE,
PACKAGE="rmutil")$res}
dmultpois <- function(y, m, s, log=FALSE){
if(any(y<0))stop("y must contain non-negative values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0|s>1))stop("s must be positive and <= 1")
ly <- max(length(y),length(m),length(s))
if(length(y)!=ly){
if(length(y)==1)y <- rep(y,ly)
else stop("length of y incorrect")}
if(length(m)!=ly){
if(length(m)==1)m <- rep(m,ly)
else stop("m and y must have the same length")}
if(length(s)!=ly){
if(length(s)==1)s <- rep(s,ly)
else stop("s and y must have the same length")}
tmp <- .C("dmp",
as.integer(y),
as.integer(3*max(c(y,100))),
as.double(m),
as.double(s),
as.integer(length(y)),
as.double(rep(1,length(y))),
res=double(length(y)),
## DUP=FALSE,
PACKAGE="rmutil")$res
if(!log)tmp <- exp(tmp)
tmp}
qmultpois <- function(p, m, s){
h <- function(y).C("pmp",
as.integer(y),
as.integer(3*max(c(y,100))),
as.double(m[i]),
as.double(s[i]),
as.integer(1),
res=double(1),
## DUP=FALSE,
PACKAGE="rmutil")$res-p[i]
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<0))stop("m must be positive")
if(any(s<=0|s>1))stop("s must be positive and <= 1")
len <- max(length(p),length(m),length(s))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
interval <- c(0,20)
if(h(interval[1])*h(interval[2])>0&&h(interval[1])>0)tmp[i] <- 0
else {
while(h(interval[1])*h(interval[2])>0)interval <- 2*interval
tmp[i] <- uniroot(h,interval)$root}}
round(tmp)}
rmultpois <- function(n=1, m, s) qmultpois(runif(n),m=m,s=s)
### power variance function Poisson distribution
###
ppvfpois <- function(q, m, s, f){
if(any(q<0))stop("q must contain non-negative values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
len <- max(length(q),length(m),length(s))
if(length(q)!=len){
if(length(q)==1)q <- rep(q,len)
else stop("length of q incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("m and q must have the same length")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("s and q must have the same length")}
if(length(f)!=len){
if(length(f)==1)f <- rep(f,len)
else stop("f and q must have the same length")}
.C("ppvfp",
as.integer(q),
as.double(m),
as.double(s),
as.double(f),
as.integer(length(q)),
res=double(length(q)),
## DUP=FALSE,
PACKAGE="rmutil")$res}
dpvfpois <- function(y, m, s, f, log=FALSE){
if(any(y<0))stop("y must contain non-negative values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
if(any(f>=1))stop("f must be < 1")
ly <- max(length(y),length(m),length(s))
if(length(y)!=ly){
if(length(y)==1)y <- rep(y,ly)
else stop("length of y incorrect")}
if(length(m)!=ly){
if(length(m)==1)m <- rep(m,ly)
else stop("m and y must have the same length")}
if(length(s)!=ly){
if(length(s)==1)s <- rep(s,ly)
else stop("s and y must have the same length")}
if(length(f)!=ly){
if(length(f)==1)f <- rep(f,ly)
else stop("f and y must have the same length")}
tmp <- log(.C("dpvfp",
as.integer(y),
as.double(m),
as.double(s),
as.double(f),
as.integer(length(y)),
as.double(rep(1,length(y))),
res=double(length(y)),
## DUP=FALSE,
PACKAGE="rmutil")$res)
if(!log)tmp <- exp(tmp)
tmp}
qpvfpois <- function(p, m, s, f){
h <- function(y).C("ppvfp",
as.integer(y),
as.double(m[i]),
as.double(s[i]),
as.double(f[i]),
as.integer(1),
res=double(1),
## DUP=FALSE,
PACKAGE="rmutil")$res-p[i]
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<0))stop("m must be positive")
if(any(s<0))stop("s must be positive")
len <- max(length(p),length(m),length(s))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
if(length(f)!=len){
if(length(f)==1)f <- rep(f,len)
else stop("f and p must have the same length")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
if(f[i]==0)tmp[i] <- qnbinom(p[i],s[i],m[i]/(m[i]+s[i]))
else {
interval <- c(0,20)
if(h(interval[1])*h(interval[2])>0&&h(interval[1])>0)
tmp[i] <- 0
else {
while(h(interval[1])*h(interval[2])>0)
interval <- 2*interval
tmp[i] <- uniroot(h,interval)$root}}}
round(tmp)}
rpvfpois <- function(n=1, m, s, f) qpvfpois(runif(n),m=m,s=s,f=f)
### gamma count distribution
###
pgammacount <- function(q, m, s){
if(any(q<0))stop("q must contain non-negative values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
1-pgamma(m*s,(q+1)*s,1)}
dgammacount <- function(y, m, s, log=FALSE){
if(any(y<0))stop("y must contain non-negative values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
tmp <- ifelse(y==0,pgamma(m*s,(y+1)*s,1,log.p=TRUE,lower.tail=FALSE),
log(pgamma(m*s,y*s+(y==0),1)-pgamma(m*s,(y+1)*s,1)))
if(!log)tmp <- exp(tmp)
tmp}
qgammacount <- function(p, m, s){
h <- function(y) 1-pgamma(m[i]*s[i],(y+1)*s[i],1)-p[i]
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0))stop("m must be positive")
if(any(s<0))stop("s must be positive")
len <- max(length(p),length(m),length(s))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
interval <- c(0,20)
if(h(interval[1])*h(interval[2])>0&&h(interval[1])>0)tmp[i] <- 0
else {
while(h(interval[1])*h(interval[2])>0)interval <- 2*interval
tmp[i] <- uniroot(h,interval)$root}}
round(tmp)}
rgammacount <- function(n=1, m, s) qgammacount(runif(n),m=m,s=s)
### Consul generalized Poisson distribution
###
pconsul <- function(q, m, s){
if(any(q<0))stop("q must contain non-negative values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
len <- max(length(q),length(m),length(s))
if(length(q)!=len){
if(length(q)==1)q <- rep(q,len)
else stop("length of q incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(any(s<0))stop("s must be positive")
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
res <- vector("numeric",length(q))
for(i in 1:length(q)){
qq <- 0:q[i]
res[i] <- sum(m[i]*exp(-(m[i]+qq*(s[i]-1))/s[i])*
(m[i]+qq*(s[i]-1))^(qq-1)/(s[i]^qq*gamma(qq+1)))}
res}
dconsul <- function(y, m, s, log=FALSE){
if(any(y<0))stop("y must contain non-negative values")
if(any(m<=0))stop("m must be positive")
if(any(s<=0))stop("s must be positive")
tmp <- log(m)-(m+y*(s-1))/s+(y-1)*log(m+y*(s-1))-y*log(s)-lgamma(y+1)
if(!log)tmp <- exp(tmp)
tmp}
qconsul <- function(p, m, s){
h <- function(y) {
pp <- 0:y
sum(m[i]*exp(-(m[i]+pp*(s[i]-1))/s[i])*(m[i]+pp*(s[i]-1))^(pp-1)/
(s[i]^pp*gamma(pp+1)))-p[i]}
if(any(p<0|p>1))stop("p must lie between 0 and 1")
if(any(m<=0))stop("m must be positive")
if(any(s<0))stop("s must be positive")
len <- max(length(p),length(m),length(s))
if(length(p)!=len){
if(length(p)==1)p <- rep(p,len)
else stop("length of p incorrect")}
if(length(m)!=len){
if(length(m)==1)m <- rep(m,len)
else stop("length of m incorrect")}
if(length(s)!=len){
if(length(s)==1)s <- rep(s,len)
else stop("length of s incorrect")}
tmp <- vector(mode="numeric",len)
for (i in 1:len){
interval <- c(0,20)
if(h(interval[1])*h(interval[2])>0&&h(interval[1])>0)tmp[i] <- 0
else {
while(h(interval[1])*h(interval[2])>0)interval <- 2*interval
tmp[i] <- uniroot(h,interval)$root}}
round(tmp)}
rconsul <- function(n=1, m, s) qconsul(runif(n),m=m,s=s)
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.