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'
}
>>> [email protected] 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'
>>> [email protected] 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/
