Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Numerical instability in rwald #17

Open
bonStats opened this issue Feb 3, 2019 · 2 comments
Open

Numerical instability in rwald #17

bonStats opened this issue Feb 3, 2019 · 2 comments

Comments

@bonStats
Copy link

bonStats commented Feb 3, 2019

Great package, thanks for publishing. I think I have found problem in the rng_wald for very large lambda. I am getting zero and negative values from rwald, for example try rwald(n=100, lambda = 1, mu = 1e20).

Line 72 and 73 of wald-distribution.cpp should be using something like Horner's method for numerical stability. I think the following should be acceptable:

x = mu * (1 + (mu/(2.0*lambda)) * ( y - sqrt(4.0*lambda*y/mu+(y*y))));

versus the original

x = mu + (mu*mu*y)/(2.0*lambda) - mu/(2.0*lambda) * sqrt(4.0*mu*lambda*y+(mu*mu)*(y*y));

I can issue a pull request late next week if that helps. Might need to factor out lambda too.

@bonStats
Copy link
Author

bonStats commented Feb 3, 2019

Been checking this out in R... I think this is better:

 q <- 1 - sqrt( (4 / (y * mu)) + 1 )
 q[q < 0] <- 0
 x <- mu * (1 +  (mu / 2) * ( y * q ) )

@twolodzko
Copy link
Owner

twolodzko commented Feb 3, 2019

@bonStats Thanks. Sure, issue PR, I'll be happy to look at it.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants