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 <[EMAIL PROTECTED]> 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.
