Re: [R] Normal and Poisson tail area expectations in R

From: Ravi Varadhan <rvaradhan_at_jhmi.edu>
Date: Thu, 14 Jun 2007 15:04:55 -0400

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: rvaradhan_at_jhmi.edu

Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html  



-----Original Message-----
From: r-help-bounces_at_stat.math.ethz.ch
[mailto:r-help-bounces_at_stat.math.ethz.ch] On Behalf Of Ravi Varadhan Sent: Thursday, June 14, 2007 2:08 PM
To: 'Charles C. Berry'
Cc: 'kavindra malik'; r-help_at_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: rvaradhan_at_jhmi.edu

Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html  



-----Original Message-----
From: Charles C. Berry [mailto:cberry_at_tajo.ucsd.edu] Sent: Thursday, June 14, 2007 1:36 PM
To: Ravi Varadhan
Cc: 'kavindra malik'; r-help_at_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: rvaradhan_at_jhmi.edu
>
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
>


> --------
>
> -----Original Message-----
> From: r-help-bounces_at_stat.math.ethz.ch
> [mailto:r-help-bounces_at_stat.math.ethz.ch] On Behalf Of kavindra malik
> Sent: Wednesday, June 13, 2007 5:45 PM
> To: Charles C. Berry
> Cc: r-help_at_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" <cberry_at_tajo.ucsd.edu> 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_at_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:cberry_at_tajo.ucsd.edu UC San Diego
> http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901
>
>
>
>
>
> ---------------------------------
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help_at_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:cberry_at_tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901

R-help_at_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_at_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. Received on Thu 14 Jun 2007 - 19:18:11 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 14 Jun 2007 - 19:31:58 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.