|
Thanks for Wolfgang Schwarz, Tom Jagger, Tim Hesterberg and
Bill Dunlap from Insightful for answers.
Tom provided pdf, cdf, quantile and random functions in
his message, and Bill pointed out that they are available
on the web:
Tim says he has put in an enhancement request to Insightful.
For my application
(fitting a convolution of the Wald and Exponential) I can get NAs in
the cdf (pinvgauss) function as the second term in the summand trades off
exp(big1)*pnorm(-big2).
Wolfgang Schwarz (in press, The Ex-Wald distribution as a
descriptive model of response times, Behaviour Research Methods, Instruments and
Computers) suggested using and approximation to pnorm expressed in terms of an
exp function due to Derenzo (1997, Mathematics of Computation, 31, 214-225),
which allows the product to be simplified,
avoiding the NAs. The approximation isn't terribly accurate unless big is very
big, and limited experiments with numerical integration of the pdf suggest
using it is only needed when exp(big1) is Inf or pnorm(-big2) is zero (again my
*limited* numerical experiments suggest the approximation is accurate to around
6 decimal places in these cases).
The following function implements the approximation. As I
am interested in the Weiner process application the parameters of this
function are mean drift rate = m and barrier = a (and drift variance is set to 1
without loss of generality). In terms of the web script parameters, the Inverse
Gaussian mean, mu = a/m, "precision", lambda = a^2, and w = q, is
a quantile:
# Wald cumulative density with protection against numerical error
pwald
<- function(w,m,a){
sqrtw <- sqrt(w); k1 <-
(m*w-a)/sqrtw; k2 <- (m*w+a)/sqrtw
p1 <- exp(2*a*m);
p2 <- pnorm(-k2)
bad <- (p1==Inf) | (p2==0); p <-
p1*p2
p[bad] <- (exp(-(k1[bad]^2)/2 -
0.94/(k2[bad]^2))/(k2[bad]*((2*pi)^.5)))
p +
pnorm(k1)
}
As an example of the problem with m = 19 and a = 19:
pinvgauss(1,19/19,19^2)
[1] NA
pwald(1,19,19)
[1] 0.5104916
For the quantile function, Tim suggested using uniroot() (slow) or approx()
(a bit faster) to solve the cdf and so get the quantile
function. The web scripts contain a much faster
algorithm contributed by Dr Paul Bagshaw, 23 Dec 98. I reproduce it below
in my parameterization (and using pdf dwald <-
function(w,m,a){ a*exp(-(a-m*w)^2/(2*w))/sqrt(2*pi*w^3)}):
qwald <- function(p,m,a)
{
thi <- a*m; rthi <-
sqrt(thi)
U <- qnorm (p)
r1 <-
1 + U/rthi + U^2/(2*thi) + U^3/(8*thi*rthi)
x <-
r1
for (i in 1:10)
{
cum <-
pwald(x,rthi,rthi)
dx <-
(cum-p)/dwald(x,rthi,rthi)
dx <-
ifelse (is.finite(dx), dx, ifelse (p >
cum,-1,1))
dx[dx < -1] <-
-1
if (all(dx == 0.))
break
x <- x -
dx
}
q <- x*a/m;
q[p==0] <- 0; q[p==1] <- Inf
q
}
Note that the second last line fixes the problem of NAs for p = 0 and p = 1
in qinvgauss. I would appreciate a reference for this function and maybe some
information on its accuracy, particularly why the limit of 10 in the loop and
why the df==0 test rather than some epsilon. For long p might it not be more
efficient to only loop on elements where dx > epsilon? It seems pretty fast
anyway.
Finally, for completeness here is my version of the web script rinvgauss
(which as Bill pointed out is way more efficient that qwald(runif) as qwald is
not explicit):
rwald <- function(n,m,a)
{
if(length(n)>1) n <-
length(n)
if(length(m)>1 && length(m)!=n) m <-
rep(m,length=n)
if(length(a)>1 && length(a)!=n)
lambda <- rep(a,length=n)
y2 <- rchisq(n,1); y2onm <-
y2/m
u <- runif(n)
r1 <- (2*a + y2onm -
sqrt(y2onm*(4*a+y2onm)))/(2*m)
r2 <- (a/m)^2/r1
ifelse(u < a/(a+m*r1), r1, r2)
}
Wolfgang sent me a Matlab version of this function which appears very
similar (mu = m in my notation)
FUNCTION generate_ig(mu,sigma,a)
LOCAL r,x,kb,kk
kk = a / mu
kb = a * mu / sigma ^ 2
x = (FN generate_nv(0,1))^2
r = RND
x = 1 + (x -
(x * x + 4 * kb * x) ^ 0.5) / (2 * kb)
IF r > 1 / (1 + x) THEN
x = 1 / x
ENDIF
RETURN kk * x
ENDFUNC
Wolfgang gave the citation: J. Dagpunar, (1988). pp. 79-80. '
Principles of Random Variate Generation.' Clarendon Press, Oxford. The S scripts
on the web cite: Chhikara and Folks, The Inverse Gaussian Distribution, Marcel
Dekker, 1989, page 53. Dont suppose anyone knows the original source (our
library doesnt have either book!)?
Thanks again to everyone for their help!
Associate Professor ANDREW HEATHCOTE
School of Behavioural
Science Telephone
Building W, University
Drive Area:
1-61-2
University of
Newcastle
Fax: 49216980
NSW, 2308,
AUSTRALIA Office: 49215952
|