## NUMERICAL INTEGRATION ## Borrowed from statmod -- thanks to Gordon Smyth for excellent software gauss.quad <- function(n,kind="legendre",alpha=0,beta=0) # Calculate nodes and weights for Gaussian quadrature. # Adapted from Netlib routine gaussq.f # Gordon Smyth, Walter and Eliza Hall Institute # Suggestion from Stephane Laurent 6 Aug 2012 # Created 4 Sept 2002. Last modified 28 Aug 2016. { n <- as.integer(n) if(n<0L) stop("need non-negative number of nodes") if(n==0L) return(list(nodes=numeric(0L), weights=numeric(0L))) kind <- match.arg(kind,c("legendre","chebyshev1","chebyshev2","hermite","jacobi","laguerre")) i <- 1L:n i1 <- i[-n] switch(kind, legendre={ lnmuzero <- log(2) a <- rep_len(0,n) b <- i1/sqrt(4*i1^2-1) }, chebyshev1={ lnmuzero <- log(pi) a <- rep_len(0,n) b <- rep_len(0.5,n-1L) b[1] <- sqrt(0.5) }, chebyshev2={ lnmuzero <- log(pi/2) a <- rep_len(0,n) b <- rep_len(0.5,n-1L) }, hermite={ lnmuzero <- log(pi)/2 a <- rep_len(0,n) b <- sqrt(i1/2) }, jacobi={ ab <- alpha+beta # muzero <- 2^(ab+1) * gamma(alpha+1) * gamma(beta+1) / gamma(ab+2) lnmuzero <- (ab+1)*log(2) + lgamma(alpha+1) + lgamma(beta+1) - lgamma(ab+2) a <- i a[1] <- (beta-alpha)/(ab+2) i2 <- i[-1] abi <- ab+2*i2 a[i2] <- (beta^2-alpha^2)/(abi-2)/abi b <- i1 b[1] <- sqrt(4*(alpha+1)*(beta+1)/(ab+2)^2/(ab+3)) i2 <- i1[-1] abi <- ab+2*i2 b[i2] <- sqrt(4*i2*(i2+alpha)*(i2+beta)*(i2+ab)/(abi^2-1)/abi^2) }, laguerre={ a <- 2*i-1+alpha b <- sqrt(i1*(i1+alpha)) lnmuzero <- lgamma(alpha+1) }) b <- c(b,0) z <- rep_len(0,n) z[1] <- 1 ierr <- 0L out <- .Fortran("gausq2",n,as.double(a),as.double(b),as.double(z),ierr,PACKAGE="rstpm2") x <- out[[2]] w <- out[[4]] w <- exp(lnmuzero + 2*log(abs(w))) list(nodes=x,weights=w) } gauss.quad.prob <- function(n,dist="uniform",l=0,u=1,mu=0,sigma=1,alpha=1,beta=1) # Calculate nodes and weights for Gaussian quadrature using probability densities. # Adapted from Netlib routine gaussq.f # Gordon Smyth, Walter and Eliza Hall Institute # Corrections for n=1 and n=2 by Spencer Graves, 28 Dec 2005 # Created 4 Sept 2002. Last modified 28 Aug 2016. { n <- as.integer(n) if(n<0L) stop("need non-negative number of nodes") if(n==0L) return(list(nodes=numeric(0L), weights=numeric(0L))) dist <- match.arg(dist,c("uniform","beta1","beta2","normal","beta","gamma")) if(n==1L){ switch(dist, uniform={x <- (l+u)/2}, beta1=,beta2=,beta={x <- alpha/(alpha+beta)}, normal={x <- mu}, gamma={x <- alpha*beta} ) return(list(nodes=x, weights=1)) } if(dist=="beta" && alpha==0.5 && beta==0.5) dist <- "beta1" if(dist=="beta" && alpha==1.5 && beta==1.5) dist <- "beta2" i <- 1L:n i1 <- 1L:(n-1L) switch(dist, uniform={ a <- rep_len(0,n) b <- i1/sqrt(4*i1^2-1) }, beta1={ a <- rep_len(0,n) b <- rep_len(0.5,n-1L) b[1] <- sqrt(0.5) }, beta2={ a <- rep_len(0,n) b <- rep_len(0.5,n-1L) }, normal={ a <- rep_len(0,n) b <- sqrt(i1/2) }, beta={ ab <- alpha+beta a <- i a[1] <- (alpha-beta)/ab i2 <- 2:n abi <- ab-2+2*i2 a[i2] <- ((alpha-1)^2-(beta-1)^2)/(abi-2)/abi b <- i1 b[1] <- sqrt(4*alpha*beta/ab^2/(ab+1)) i2 <- i1[-1] # 2:(n-1) abi <- ab-2+2*i2 b[i2] <- sqrt(4*i2*(i2+alpha-1)*(i2+beta-1)*(i2+ab-2)/(abi^2-1)/abi^2) }, gamma={ a <- 2*i+alpha-2 b <- sqrt(i1*(i1+alpha-1)) }) b <- c(b,0) z <- rep_len(0,n) z[1] <- 1 ierr <- 0L out <- .Fortran("gausq2",n,as.double(a),as.double(b),as.double(z),ierr,PACKAGE="rstpm2") x <- out[[2]] w <- out[[4]]^2 switch(dist, uniform = x <- l+(u-l)*(x+1)/2, beta1=,beta2=,beta = x <- (x+1)/2, normal = x <- mu + sqrt(2)*sigma*x, gamma = x <- beta*x) list(nodes=x,weights=w) }