Re: [R] integrate and quadratic forms

From: Thomas Lumley <tlumley_at_u.washington.edu>
Date: Fri 19 Jan 2007 - 17:48:33 GMT

As the documentation for integrate() says, the function must be vectorized

   f: an R function taking a numeric first argument and returning a

           numeric vector of the same length.

so you can't use sum(). You need matrix operations or an explicit loop to add up the terms.

         -thomas

On Thu, 18 Jan 2007, rdporto1 wrote:

> 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.
>

Thomas Lumley			Assoc. Professor, Biostatistics
tlumley@u.washington.edu	University of Washington, Seattle

______________________________________________
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 Sat Jan 20 04:53:22 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.