search for: abramowitz

Hi anyone coded up Stirling numbers in R? [I need unsigned Stirling numbers of the first kind] cheers Robin -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743

...using the pgamma() function. Is Prof. Ripley's message (the first above) saying that upperIncompletegamma(x,\rate) = pgamma(x,\rate) * gamma(\rate) ? I have not yet been able to see why this would be true. I keep tripping myself up on the differing notation of R, mathworld, wikipaedia and Abramowitz and Stegun's "Handbook of mathematical functions". Would someone please confirm that this is mathematically correct? 3) Adding an igamma() function or documentation for how to generate it. -------------------------------------------------------------------- Since I am at least th...

I am using pnorm() with the log.p=T argument to get approximations to ln \Phi(x) and qnorm with the log.p=T argument to get estimates of \Phi^{-1}(exp(x)). What approximations are used in these two functions (I noticed in the source pnorm.c it doesn't look like Abramowitz and Stegen) and where can I find the citation? Thanks, Richard Morey

...el; } - if (df > 4e5 || del*del > 2*M_LN2*(-DBL_MIN_EXP)) { + if (df > 4e5 || del*del > 2*M_LN2*(-(DBL_MIN_EXP))) { /*-- 2nd part: if del > 37.62, then p=0 below FIXME: test should depend on `df', `tt' AND `del' ! */ /* Approx. from Abramowitz & Stegun 26.7.10 (p.949) */ -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not...

..., /*lower_tail*/0, log_p) : pbeta (n / (n + x * x), n / 2.0, 0.5, /*lower_tail*/1, log_p); thus avoiding the loss of accuracy in the pbeta routine when 1-1./nx is calculated. It also makes the if (n > 4e5) { /*-- Fixme(?): test should depend on `n' AND `x' ! */ ??? /* Approx. from? Abramowitz & Stegun 26.7.8 (p.949) */ ??? val = 1./(4.*n); ??? return pnorm(x*(1. - val)/sqrt(1. + x*x*2.*val), 0.0, 1.0, ???????? lower_tail, log_p); } code unneccessary. Ian Smith Personally, I think the code should also guard against the possible overflow of the x * x expressions. _____________...

...es of x by allowing for exponential scaling, but there is no facility for dealing with large \nu. What would work for me would be an lbesselK function in the manner of lgamma which returned the log of K_\nu(x) for large \nu. Does anybody have any leads on this? Note that I have trawled through Abramowitz and Stegun and found 9.7.8 which doesn't work for me because of the complication in the definition of the x argument. I have also seen a result of Ismail (1977) reported by Barndorff-Nielsen and Blaesild which has the other problem, the treatment of the x argument is too simple. To do the...

Full_Name: Morten Welinder Version: 2 OS: Solaris/space/gcc2.95.2 Submission from: (NULL) (65.213.85.217) I changed src/nmath/standalone/test.c to read: --------------------------------------------------------------------------------- #define MATHLIB_STANDALONE 1 #include <Rmath.h> #include <stdio.h> int main() { double x; for (x = 99990; x <= 100009; x++) printf

...ions include elliptic integrals, Airy functions, hypergeometric functions, and so on. Most functions optionally return error estimates. This library is a little odd in that documentation is limited to a pointer to the GSL reference manual, and occasional re-productions of tables and figures from Abramowitz and Stegun. It passes R CMD check under Fedora Linux, and MacOSX 10.3 if the GSL (version 1.5) is installed in the default locations. A windows version is available at http://www.stats.ox.ac.uk/pub/RWin/2.0.0/ (also, a big Thank You to Kurt Hornik for advice on configure.ac) -- Robin Ha...

Hi List I'm thinking of putting together some R functionality (maybe a package) of elliptic functions. An elliptic function is one that is doubly periodic on the complex plane and I think they are cool. One objective of mine is to reproduce Abramowitz and Stegun's beautiful diagram on p643 (figure 18.5). The gsl package does sn, cn, dn et seq, but before I put a huge amount of effort in, has anyone else coded up any related material? I'd be interested to hear if anyone has anything on Weierstrass's P function (in particular, co...

Hi the following gave me a shock: > atan(2) [1] 1.107149 > atan(2+0i) [1] -0.4636476+0i > or, perhaps more of a gotcha: > atan(1.0001+0i) [1] -0.7853482+0i > atan(0.9999+0i) [1] 0.7853482+0i > evidently atan()'s branch cuts aren't where I thought they were. Where do I look for documentation on this? -- Robin Hankin Uncertainty Analyst National

...ions include elliptic integrals, Airy functions, hypergeometric functions, and so on. Most functions optionally return error estimates. This library is a little odd in that documentation is limited to a pointer to the GSL reference manual, and occasional re-productions of tables and figures from Abramowitz and Stegun. It passes R CMD check under Fedora Linux, and MacOSX 10.3 if the GSL (version 1.5) is installed in the default locations. A windows version is available at http://www.stats.ox.ac.uk/pub/RWin/2.0.0/ (also, a big Thank You to Kurt Hornik for advice on configure.ac) -- Robin Ha...

...-------------------------------------- CC| |tony3@garlandconsulting.us ------- Comment #12 from tony3@garlandconsulting.us 2008-10-02 12:46 CST ------- The --bwlimit work-around did not fix the problem for me. Instead, I found this page: http://marc-abramowitz.com/archives/2007/10/14/solving-rsync-hangs-with-cygwin/ which described running rsync on the client against a TCP daemon service (rather than over ssh) which solved the problem for me. After setting up an rsyncd.conf file and installing the rsyncd daemon as a service on the remote machine: cygr...

...-------------------------- --- src/nmath/pt.c.~1~ Wed Feb 3 12:21:46 1999 +++ src/nmath/pt.c Wed Mar 10 18:15:42 1999 @@ -41,6 +41,11 @@ if(!finite(n)) return pnorm(x, 0.0, 1.0); #endif + if (n > 4e5) { /*-- Fixme(?): test should depend on `n' AND `x' ! */ + /* Approx. from Abramowitz & Stegun 26.7.8 (p.949) */ + val = 1./(4.*n); + return pnorm(x*(1. - val)/sqrt(1. + x*x*2.*val), 0.0, 1.0); + } val = 0.5 * pbeta(n / (n + x * x), n / 2.0, 0.5); return (x > 0.0) ? 1 - val : val; } -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-...

...expansion is incorrect. However, since this approximation just gives an initial solution, the end result of the function is still correct. To fix the approximation, in the snippet below the line gamma = sigma; should be replaced by gamma = 1.0/sigma; /* the skewness */ The reference is Abramowitz and Stegun 'Handbook of Mathmatical Functions' pages 935 and 928 Mikael double qpois(double p, double lambda, int lower_tail, int log_p) { double mu, sigma, gamma, z, y; #ifdef IEEE_754 if (ISNAN(p) || ISNAN(lambda)) return p + lambda; #endif if(!R_FINITE(lambda)) ML_ERR_r...

Hi all. I may be wrong, (and often am), but in trying to determine how to calculate the erf function, the documentation for 'pnorm' states: ## if you want the so-called 'error function' erf <- function(x) 2 * pnorm(x * sqrt(2)) - 1 ## and the so-called 'complementary error function' erfc <- function(x) 2 * pnorm(x * sqrt(2), lower=FALSE) Should, instead, it read:

Hello all, I'm hoping to find a way to evaluate the following sort of integral in R. \int_a^b \int_{g(y)}^Inf f(x,y) dx dy. The integral has no closed form and so must be evaluated numerically. The "adapt" package provides for multidimensional integration but does not appear to allow the limits of integration to be a function. I need to evaluate a number of integrals of this

Dear R users, i try to integrate lgamma from 0 to Inf. But here i get the message "roundoff error is detected in the extrapolation table", if i use 1.0e120 instead of Inf the computation works, but this is against the suggestion of integrates help information to use Inf explicitly. Using stirlings approximation doesnt bring the solution too. ## Stirlings approximation lgammaApprox

...for (i=0; i<=3; i++) pn[i] /= IBETA_SMALL; } while (fabs(next-prev) > DOUBLE_EPS*prev); factor = a*log(x) + (b-1)*log(1-x); factor -= LGAMMA(a+1) + LGAMMA(b) - LGAMMA(a+b); I = exp(factor) * next; return(flipped ? 1-I : I); } /* * Incomplete gamma function. * Reference: Abramowitz and Stegun. * Assumptions: x >= 0; df > 0. */ double igamma(x, df) double x, df; { double factor, term, gintegral, pn[6], rn, ak, bk; double increment, df1; int i, count, k; if (x <= 0.0) return(0.0); if (df < 1.0) { increment = exp(df*log(x) - x - LGAMMA(df + 1.0)); df...

Apologies for the off-topic question; as usual I'm trying to draw upon the unparalleled knowledge and sagacity of the r-help list. Please reply off-list if you can help me out. A collaborator of mine found a formula we need, on sheets which he had photocopied out of a book, some years ago. He cannot remember which book (he's getting to be as senile and forgetful as I am, poor bloke!).

Hi, R guys: I'm using L-BFGS-B method of optim for minimization problem. My function called besselI function which need non-negative parameter and the besselI will overflow if the parameter is too large. So I set the constraint box which is reasonable for my problem. But the point outside the box was test, and I got error. My program and the error follows. This program depends on CircStats