Inspired by Chuck's elegant solution to the Poisson tail area problem, here is a simple solution to the normal tail area expectation that does not use integrate().
Gu.k <- function(k) {1/sqrt(2*pi) * exp(-k*k/2) - k * pnorm(k, lower=FALSE)} > k <- 1:10 > Gu.k(k) [1] 8.331547e-02 8.490703e-03 3.821543e-04 7.145258e-06 5.346166e-08 1.563570e-10 1.760326e-13 7.550262e-17 1.224779e-20 7.474560e-25 > > sapply(k, function(k)integrate(function(x) (x-k)*dnorm(x),lower=k,upper=Inf)$val) [1] 8.331547e-02 8.490702e-03 3.821543e-04 7.145259e-06 5.346163e-08 1.563570e-10 1.760326e-13 7.550262e-17 1.224779e-20 7.474560e-25 > Ravi. ---------------------------------------------------------------------------- ------- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html ---------------------------------------------------------------------------- -------- -----Original Message----- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Ravi Varadhan Sent: Thursday, June 14, 2007 2:08 PM To: 'Charles C. Berry' Cc: 'kavindra malik'; r-help@stat.math.ethz.ch Subject: Re: [R] Normal and Poisson tail area expectations in R Perfect, Chuck. I got a closed-form solution after some algebraic labor, but your solution is simple and elegant. Ravi. ---------------------------------------------------------------------------- ------- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html ---------------------------------------------------------------------------- -------- -----Original Message----- From: Charles C. Berry [mailto:[EMAIL PROTECTED] Sent: Thursday, June 14, 2007 1:36 PM To: Ravi Varadhan Cc: 'kavindra malik'; r-help@stat.math.ethz.ch Subject: RE: [R] Normal and Poisson tail area expectations in R Ravi, This looks simple to me. km_G <- function(lambda,k) lambda*ppois(k-1,lambda,lower=FALSE) - k*ppois(k,lambda,lower=FALSE) Am I confused here? Chuck On Wed, 13 Jun 2007, Ravi Varadhan wrote: > > More interesting is the Poisson convolution. I don't know if there is an > analytic solution to this. I looked at Jolley's "Summation of Series" and > Abramowitz and Stegun, but no help there. It seems that discrete FFT > technique should work. Does anyone know the answer? > > Ravi. > ---------------------------------------------------------------------------- > ------- > > Ravi Varadhan, Ph.D. > > Assistant Professor, The Center on Aging and Health > > Division of Geriatric Medicine and Gerontology > > Johns Hopkins University > > Ph: (410) 502-2619 > > Fax: (410) 614-9625 > > Email: [EMAIL PROTECTED] > > Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html > > > > ---------------------------------------------------------------------------- > -------- > > -----Original Message----- > From: [EMAIL PROTECTED] > [mailto:[EMAIL PROTECTED] On Behalf Of kavindra malik > Sent: Wednesday, June 13, 2007 5:45 PM > To: Charles C. Berry > Cc: r-help@stat.math.ethz.ch > Subject: Re: [R] Normal and Poisson tail area expectations in R > > Thank you very much. This solves the problem I was trying to solve. I am new > to R and am learning. A great lesson in the power of R... > > "Charles C. Berry" <[EMAIL PROTECTED]> wrote: On Wed, 13 Jun 2007, > kavindra malik wrote: > >> I am interested in R functions for the following integrals / sums > (expressed best I can in text) - >> >> Normal: G_u(k) = Integration_{Lower limit=k}^{Upper limit=infinity} [(u > -k) f(u) d(u)], where where u is N(0,1), and f(u) is the density function. >> >> Poisson: G(lambda,k) = Sum_{Lower limit=k}^{Upper limit=infinity} [(x-k) > p(x, lambda)] where P(x,lambda) is the Poisson prob function with parameter > lambda. >> >> The Normal expression is very commonly used in inventory management to >> determine safety stocks (and its tabular values can be found in some >> texts) - and I am also looking for Poisson and/or Gamma as that'd fit >> the situation better. >> >> I am wondering if there are standard functions in R that might allow me to > get these values, instead of needing to do the numerical integration, etc. > myself. > > Not that I know of, but it is not difficult to do the integration: > >> k <- 1.1 # for example >> integrate(function(x) (x-k)*dnorm(x),lower=k,upper=Inf) > 0.06861951 with absolute error < 5.5e-07 >> > > see > > ?integrate > ?qnorm > ?qpois > ?qgamma > >> Thank you very much. >> >> >> >> --------------------------------- >> Sucker-punch spam with award-winning protection. >> >> [[alternative HTML version deleted]] >> >> ______________________________________________ >> R-help@stat.math.ethz.ch mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. >> > > Charles C. Berry (858) 534-2098 > Dept of Family/Preventive > Medicine > E mailto:[EMAIL PROTECTED] UC San Diego > http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 > > > > > > --------------------------------- > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > Charles C. Berry (858) 534-2098 Dept of Family/Preventive Medicine E mailto:[EMAIL PROTECTED] UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.