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

______________________________________________
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 Fri Jan 19 11:47:51 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 - 02:30:28 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.