From: Thomas Lumley <tlumley_at_u.washington.edu>

Date: Fri 19 Jan 2007 - 17:48:33 GMT

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

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