Revision 936fa1723eb80b0f9d4b12a1e50e5c542b549b35 authored by Mark Clements on 03 March 2021, 16:10:02 UTC, committed by cran-robot on 03 March 2021, 16:10:02 UTC
1 parent 3344f17
gaussquad.R
## 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 Guassian 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)
}
Computing file changes ...