1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
# NUMERICAL INTEGRATION
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="statmod")
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="statmod")
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)
}
|
