Topic: Utilities

Topic Description:

Basic R utilities used in many packages and functions, such as date and string manipulations, statistical distributions, geometry of lines and distances. The R package date is required for the two data functions.



File: utilities/distributions.r

View File Source Download File No help file available

Function: dgammaPlusdexp

Function Description: dgammaPlusdexp

The PDF of the sum of gamma and negative exponential distribution. The shape and rate of the gamma are a and r; mean and sd are the mean and sd of the gamma. Lambda is the rate of the exponential. This is only an appoximation based on the observation that the resulting distribution is very close to a gamma. So I simply work out a new gamma whose mean is the sum of the means of the initial gamma and exponential, and likewise for the new variance. As long as gamma's rate > the exponential lambda, the distribution can be specified (using pgamma) as in dgammaMinusdexp. But if rate < lambda, this fails. The gamma approximation fails if the sd is sufficiently higher than the mean, and the mean is low. Then the gamma is absurdly skewed, and the shape of the sum is dominated by the exponential at low z, nothing like a true gamma. It appears to work reasonably well as long as the mean >= sd, or even if mean>=0.5*sd.

Function Arguments:

ArgumentDefault Value
z
mean
sd
lambda
draws10000
div.01
xrangec(0,25)
xmax4
graphitFALSE

Function Source:

dgammaPlusdexp=function(z,mean,sd,lambda,draws=10000,div=.01,xrange=c(0,25),xmax=4,graphit=FALSE)
{
r=mean/(sd^2)
a=mean*r

pdf.z=numeric()
pdf.z[z<=0]=0
part2=numeric()

if(graphit)
{
gamma=rgamma(draws,shape=a,rate=r)
expon=rexp(draws,rate=lambda)

net=gamma+expon

highvalue=max(max(z),max(gamma),max(expon),max(net))
x=seq(0,highvalue+div,by=div)
dev.set(2)
hist(gamma,breaks=x,xlim=xrange)
dev.set(3)
hist(expon,breaks=x,xlim=xrange)
dev.set(4)
hist(net,breaks=x,xlim=xrange)
}

newmean=mean+1/lambda
newvar=sd^2+1/(lambda^2)
newshape=newmean^2/newvar
newrate=newmean/newvar
pdf.z[z>0]=dgamma(z[z>0],shape=newshape,rate=newrate)
if(graphit) lines(z,draws*div*pdf.z)

# browser()
return(pdf.z)
}