Thomas Lumley tlumley at u.washington.edu
Fri Jan 19 18:48:33 CET 2007

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

-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.
>
>
> Rogerio.
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help