[R] integrate and quadratic forms

rdporto1 rdporto1 at terra.com.br
Fri Jan 19 01:17:00 CET 2007

Hi all.

I'm trying to numerically invert the characteristic function
of a quadratic form following Imhof's (1961, Biometrika 48)

The parameters are:


I've implemented Imhof's procedure two ways that, for me,
should give the same answer:

#more legible
integral1 = function(u) {
  o=(1/2)*sum(h*atan(lambda*u)+sigma^2*lambda*u/(1+lambda^2*u^2)) - q*u/2
  rho=prod((1+lambda^2*u^2)^(h/4))*exp( (1/2)*sum((sigma*lambda*u)^2/(1+lambda^2*u^2)) )
  integrand = sin(o)/(u*rho) 

#same as above
integral2= function(u) {
((1/2)*sum(h*atan(lambda*u)+sigma^2*lambda*u/(1+lambda^2*u^2)) - q*u/2)/
exp( (1/2)*sum((sigma*lambda*u)^2/(1+lambda^2*u^2)) )))

The following should be near 0.18. However, nor the answers are near this
value neither they agree each other!

> 1/2+(1/pi)*integrate(integral1,0,Inf)$value
[1] 1.022537
> 1/2+(1/pi)*integrate(integral2,0,Inf)$value
[1] 1.442720

What's happening? Is this a bug or OS specific? Shouldn't they give the 
same answer? Why do I get results so different from 0.18? In time:
the procedure works fine for q=.2.

I'm running R 2.4.1 in a PC with Windows XP 32bits. Other ways (in R) to
find the distribution of general quadratic forms are welcome.

Thanks in advance.


More information about the R-help mailing list