[R] integrate and quadratic forms

From: rdporto1 <rdporto1_at_terra.com.br>
Date: Fri 19 Jan 2007 - 00:17:00 GMT


Hi all.

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

The parameters are:

lambda=c(.6,.3,.1)
h=c(2,2,2)
sigma=c(0,0,0)
q=3

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)/ (u*(prod((1+lambda^2*u^2)^(h/4))*
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.

Rogerio.



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. Received on Fri Jan 19 11:21:02 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Fri 19 Jan 2007 - 18:30:27 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.