# Re: [R] integrate and quadratic forms

From: jim holtman <jholtman_at_gmail.com>
Date: Fri 19 Jan 2007 - 00:42:38 GMT

The do give the same answer, unfortunately the examples you sent were not the same. Integral2 was missing a 'sin'. So I would assume there might be something else wrong with your functions. You might want to try breaking it down into smaller steps so you can see what you are doing. It definitely is hard to read in both cases.

> 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) {
+ sin((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)) )))
+ }
```

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

On 1/18/07, rdporto1 <rdporto1@terra.com.br> 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.022537
> > 1/2+(1/pi)*integrate(integral2,0,Inf)\$value
>  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@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

```--
Jim Holtman
Cincinnati, OH
+1 513 646 9390

What is the problem you are trying to solve?

[[alternative HTML version deleted]]

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help