R/distns.R

In nclbayes: Functions, Demos and Examples of Bayesian Analysis and Basic MCMC

#' The Inverse Chi distribution
#'
#' Density, distribution function, quantile function and random generation for the inverse
#' chi distribution with parameters a and b.
#' @param x vector of quantities.
#' @param p vector of probabilities.
#' @param n number of observations.
#' @param a first parameter. Must be strictly positive.
#' @param b second parameter. Must be strictly positive.
#' @note If X~Inv-Chi(a,b) then it has density f(x)=2b^ax^(-2*a-1)e^(-b/x^2)/Gamma(a). Also 1/X^2~Gamma(a,b).
#' @keywords character
#' @importFrom stats dgamma
#' @examples dinvchi(1, 2, 2)
#' @export dinvchi qinvchi pinvchi rinvchi
dinvchi=function(x,a,b){2*dgamma(1./x^2,a,b)/x^3}

#' @rdname dinvchi
#' @examples pinvchi(1, 2, 2)
pinvchi=function(x,a,b){1-pgamma(1./x^2,a,b)}

#' @rdname dinvchi
#' @examples qinvchi(0.95, 2, 2)
qinvchi=function(p,a,b){1/sqrt(qgamma(1-p,a,b))}

#' @rdname dinvchi
#' @examples rinvchi(1, 2, 2)
rinvchi=function(n,a,b){1/sqrt(rgamma(n,a,b))}

#' The Generalised t distribution
#'
#' Density, distribution function, quantile function and random generation for the generalised
#' t distribution with a degrees of freedom, mean b and scale c.
#' @param x vector of quantities.
#' @param p vector of probabilities.
#' @param n number of observations.
#' @param a degrees of freedom. Must be strictly positive.
#' @param b mean.
#' @param c scale. Must be strictly positive.
#' @note If X~t_a(b,c) then it has density f(x)=(1_(x-b)^2/(ac))^(-(a+1)/2)/(sqrt(ac)*Beta(a,b)). Also (X-b)/sqrt(c)~t_a.
#' @keywords character
#' @examples dgt(1, 10, 2, 2)
#' @export dgt pgt qgt rgt
dgt=function(x,a,b,c){dt((x-b)/sqrt(c),a)/sqrt(c)}

#' @rdname dgt
#' @importFrom stats pt
#' @examples pgt(1, 10, 2, 2)
pgt=function(x,a,b,c){pt((x-b)/sqrt(c),a)}

#' @rdname dgt
#' @importFrom stats qt
#' @examples qgt(0.95, 10, 2, 2)
qgt=function(p,a,b,c){b+sqrt(c)*qt(p,a)}

#' @rdname dgt
#' @examples rgt(1, 10, 2, 2)
#' @importFrom stats rt
rgt=function(n,a,b,c){b+sqrt(c)*rt(n,a)}

#' The Normal-Gamma distribution
#'
#' Density and random generation for the normal-gamma distribution with parameters b, c, g and h.
#' Also contours and confidence regions.
#' @param mu vector of quantities.
#' @param tau vector of quantities.
#' @param b mean.
#' @param c must be strictly positive.
#' @param g must be strictly positive.
#' @param h must be strictly positive.
#' @param n sample size.
#' @note If (mu,tau)^T~NGa(b,c,g,h) then mu|tau~N(b,1/(c*tau)) and tau~Ga(g,h). Also mu~t_(2g)(b,h/(gc)).
#' @keywords character
#' @examples dnormgamma(1, 2, 1, 2, 5, 3)
#' @export dnormgamma rnormgamma
dnormgamma=function(mu,tau,b,c,g,h){dnorm(mu,b,1/sqrt(c*tau))*dgamma(tau,g,h)}

#' @rdname dnormgamma
#' @examples rnormgamma(1, 1, 2, 5, 3)
rnormgamma=function(n,b,c,g,h){
  output=matrix(ncol=2,nrow=n);colnames(output)=c("mu","tau");
  output[,2]=rgamma(n,g,h);output[,1]=rnorm(n,b,1/sqrt(c*output[,2]));
  output
}

#' @rdname dnormgamma
#' @param p probability.
#' @param ... arguments to be passed to the plot function when plotting contours.
#' @importFrom graphics contour
#' @importFrom stats rgamma rnorm
#' @examples
#' mu=seq(2,4,len=1000)
#' tau=seq(0,30,len=1000)
#' NGacontour(mu,tau,3,1,4,0.35,0.95)
#' NGacontour(mu,tau,2.5,30,20,2,0.95,add=TRUE,lty=3)
#' @export NGacontour
NGacontour <- function (mu, tau, b, c, g, h, p=NULL,...) {
  pdf=function(mu,tau){dnormgamma(mu,tau,b,c,g,h)}
  z=outer(mu,tau,pdf)
  if (length(p)==0) {
    contour(mu,tau,z,xlab=expression(mu),ylab=expression(tau),...)}
  else {
    NGaclevels=function(p,b,c,g,h){n=100000
                                   tau=rgamma(n,g,h)
                                   mu=rnorm(n,b,1/sqrt(c*tau))
                                   y=dnormgamma(mu,tau,b,c,g,h)
                                   y=sort(y)
                                   signif(y[floor(n*(1-p))],2)}
    pdflevels=NGaclevels(p,b,c,g,h)
    contour(mu,tau,z,levels=pdflevels,xlab=expression(mu),ylab=expression(tau),...)}
}

#' The Normal-InvChi distribution
#'
#' Density and random generation for the normal-inv-chi distribution with parameters b, c, g and h.
#' @param mu vector of quantities.
#' @param sigma vector of quantities.
#' @param b mean.
#' @param c must be strictly positive.
#' @param g must be strictly positive.
#' @param h must be strictly positive.
#' @param n number of observations.
#' @note If (mu,sigma)^T~NInvChi(b,c,g,h) then mu|sigma~N(b,sigma^2/c) and sigma~Inv-Chi(g,h). Also mu~t_(2g)(b,h/(gc)).
#' @keywords character
#' @examples dnorminvchi(1, 1/sqrt(2), 1, 2, 5, 3)
#' @export dnorminvchi rnorminvchi
dnorminvchi=function(mu,sigma,b,c,g,h){
  dnorm(mu,b,sigma/sqrt(c))*dinvchi(sigma,g,h)
}

#' @rdname dnorminvchi
#' @examples rnorminvchi(1, 1, 2, 5, 3)
rnorminvchi=function(n,b,c,g,h){
  output=matrix(ncol=2,nrow=n);colnames(output)=c("mu","sigma");
  output[,2]=rinvchi(n,g,h);output[,1]=rnorm(n,b,output[,2]/sqrt(c));
  output
}

#' The Bivariate Normal distribution
#'
#' Density of the bivariate normal distribution with mean vector b and covariance matrix c.
#' @param x vector of quantities.
#' @param y vector of quantities.
#' @param b mean vector (2x1).
#' @param c covariance matrix. Must be positive definite 2x2 matrix.
#' @keywords character
#' @importFrom stats dnorm
#' @examples
#' dbvnorm(1, 2, c(1,2), diag(1,2))
#' @export
dbvnorm=function(x,y,b=c(0,0),c=diag(1,2)){
    xs=(x-b[1])/sqrt(c[1,1])
    ys=(y-b[2])/sqrt(c[2,2])
    rho=c[1,2]/sqrt(c[1,1]*c[2,2])
    dnorm(ys,rho*xs,sqrt(1-rho*rho))*dnorm(xs)/sqrt(c[1,1]*c[2,2])
}

#' The Bivariate t distribution
#'
#' Density (and its contours) of the bivariate t distribution with a degrees of freedom,
#' mean vector b and covariance matrix c.
#' @param x vector of quantities.
#' @param y vector of quantities.
#' @param a degrees of freedom. Must be strictly positive.
#' @param b mean vector (2x1).
#' @param c covariance matrix. Must be a positive definite 2x2 matrix.
#' @keywords character
#' @importFrom stats dt
#' @examples
#' dbvt(1, 2, 10, c(1,2), diag(1,2))
#' @export
dbvt=function(x,y,a,b,c){
  xs=(x-b[1])/sqrt(c[1,1])
  ys=(y-b[2])/sqrt(c[2,2])
  rho=c[1,2]/sqrt(c[1,1]*c[2,2])
  dgt(ys,a+1,rho*xs,(1-rho*rho)*(a+xs^2)/(a+1))*dt(xs,a)/sqrt(c[1,1]*c[2,2])
}

#' @rdname dbvt
#' @param p probability.
#' @param ... arguments to be passed to the plot function when plotting contours.
#' @importFrom stats qf
#' @examples
#' beta1=seq(-3,3,len=100)
#' beta2=seq(-3,3,len=100)
#' tcontour(beta1,beta2,10,c(0,0),matrix(c(1,0.8,0.8,1),ncol=2))
#' @export
tcontour=function (x,y,a,b,c,p=NULL, ...) {
  pdf=function(x,y){dbvt(x,y,a,b,c)}
  z=outer(x,y,pdf)
  if (length(p)==0) {
    contour(x,y,z,...)}
  else {
    pdflevels=(gamma(a/2+1)/(sqrt(det(c))*a*pi*gamma(a/2)))*(1+2*qf(p,2,a)/a)^(-a/2-1)
    contour(x,y,z,levels=signif(pdflevels,2),...)
  }
}