R help - Jan 2007 - integrate and quadratic forms

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.

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] 1.022537
> 1/2+(1/pi)*integrate(integral2,0,Inf)$value
[1] 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] 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.
>


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

What is the problem you are trying to solve?

	[[alternative HTML version deleted]]

Thanks Jim. But a narrower problem still remains.
Please, look at this new code:

lambda=c(.6,.3)
integral = function(u) {
  theta = (atan(.6*u) + atan(.3*u))/2 - .1*u/2
  rho = (1+.6^2*u^2)^(1/4) * (1+.3^2*u^2)^(1/4)
  integrand = sin(theta)/(u*rho)
}
> integrate(integral,0,Inf)$value
[1] 1.222688

If I replace
  theta = sum(atan(lambda*u))/2 - .1*u/2
I get
> integrate(integral,0,Inf)$value
[1] 1.517134

It seems there is a problem between the sum() and integrate()
functions. 

Is there a way to solve this problem?

Thanks in advance.

Rogerio.

---------- Part of the Original Message -----------
> They 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.
>
> Jim Holtman
> Cincinnati, OH
> +1 513 646 9390
> 
> What is the problem you are trying to solve?

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 at 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 at u.washington.edu	University of Washington, Seattle