Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: Integrals


From   "Jamie Griffin" <Jamie.Griffin@lshtm.ac.uk>
To   <statalist@hsphsun2.harvard.edu>
Subject   Re: st: Integrals
Date   Fri, 30 Sep 2005 11:53:19 +0100

The integrals can also be calculated in a single step:

I(r, b)=A+sign(b)^(r-1)*2^(r/2-1)/sqrt(pi)*gamma((r+1)/2)*P((r+1)/2,
b^2/2)
where P is the incomplete gamma function and A is:
A = -2^(r/2-1)*((r-1)/2)!/sqrt(pi) for r odd
A = (r-1)!/((r/2-1)!*2^(r/2)) for r even

Only sign(`b')^(`r'-1)*gammap((`r'+1)/2, (`b')^2/2) needs to be
recalculated for different values of b.

The code below checks these.

Jamie Griffin.


clear
set obs 2000
tempname b I
scalar  `b'=6*uniform()-3
di `b'

gen double x=`b'-(_n-1)/100
gen double f=normden(x)

gen double g=.
forval r=3(2)15{
	qui replace g=x^(`r')*f
	integ g x 
	scalar
`I'=-2^(`r'/2-1)/sqrt(_pi)*exp(lnfact((`r'-1)/2))+sign(`b')^(`r'-1)*2^(`r'/2-1)/sqrt(_pi)*exp(lngamma((`r'+1)/2))*gammap((`r'+1)/2,
(`b')^2/2)
	di `I'
}

forval r=2(2)16{
	qui replace g=x^(`r')*f
	integ g x 
	scalar
`I'=exp(lnfact(`r'-1))/(exp(lnfact(`r'/2-1))*2^(`r'/2))+sign(`b')^(`r'-1)*2^(`r'/2-1)/sqrt(_pi)*exp(lngamma((`r'+1)/2))*gammap((`r'+1)/2,
(`b')^2/2)
	di `I'
}




>>> Jamie.Griffin@lshtm.ac.uk 09/29/05 3:29 pm >>>
You can evaluate the integrals by integrating by parts or directly.
Let I(r, b)=integral from minus infinity to b of x^r*f(x)dx where f is
the standard normal density and r is a positive integer.
Putting x^r*f(x)=x^(r-1)*x*f(x) and integrating by parts gives 
I(r, b)=-b^(r-1)*exp(-b^2/2)/sqrt(2*pi)+(r-1)*I(r-2, b)
I(0,b)=normal(b) 
and 
I(1,b)= - normalden(b) 
so any I() for any higher integer can be found recursively.

The integral of f(x)^2 can be found by substituting w=sqrt(2)*x
J(b)=1/(2*sqrt(pi))*normal(sqrt(2)*b)

The code below checks these results numerically (in Stata 8).

Jamie Griffin

clear
set obs 2000
local b=4*uniform()-2
di `b'

gen x=`b'-(_n-1)/100
gen f=normden(x)

gen f2=f^2
integ f2 x
di 1/(2*sqrt(_pi))*norm(sqrt(2)*`b')

integ f x
local I0=norm(`b')
di `I0'

gen g1=x*f
integ g1 x
local I1=-normden(`b')
di `I1'

gen g2=x^2*f
integ g2 x
local I2=-(`b')*exp(-(`b')^2/2)/sqrt(2*_pi)+`I0'
di `I2'

gen g3=x^3*f
integ g3 x
local I3=-(`b')^2*exp(-(`b')^2/2)/sqrt(2*_pi)+2*`I1'
di `I3'

gen g4=x^4*f
integ g4 x
local I4=-(`b')^3*exp(-(`b')^2/2)/sqrt(2*_pi)+3*`I2'
di `I4'

gen g5=x^5*f
integ g5 x
local I5=-(`b')^4*exp(-(`b')^2/2)/sqrt(2*_pi)+4*`I3'
di `I5'

gen g6=x^6*f
integ g6 x
local I6=-(`b')^5*exp(-(`b')^2/2)/sqrt(2*_pi)+5*`I4'
di `I6'


>>> ben.jann@soz.gess.ethz.ch 09/29/05 2:28 pm >>>
Hi, I am interested in evaluating integrals of

 x^r * f(x)

and

 f(x)^2

from minus infinity to b, where b < infinity.
r is a positive integer and f(u) is the standard 
normal density.

Is it correct that such integrals can only be 
solved via numerical integration and that they
cannot be expressed in terms of normal() and 
normalden()?

Thanks,
ben

*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html 
*   http://www.stata.com/support/statalist/faq 
*   http://www.ats.ucla.edu/stat/stata/ 
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html 
*   http://www.stata.com/support/statalist/faq 
*   http://www.ats.ucla.edu/stat/stata/
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index