s-news
[Top] [All Lists]

Another question about double integration

To: <s-news@lists.biostat.wustl.edu>
Subject: Another question about double integration
From: "Dr. Takashi Kikuchi" <takashi.kikuchi@st-hughs.oxford.ac.uk>
Date: Sat, 7 Jul 2001 11:55:58 +0100
References: <Pine.SOL.3.96.1010705120432.997K-100000@foster> <001501c10570$54753cd0$761401a3@stats.ox.ac.uk> <001001c105a3$ef903ca0$a6e501a3@ox.ac.uk>
Dear All
Please tell me the meaning of an error message written in the bottom and the
solution. Thanks!
Takashi

 benefit2 <- function(n, B = 1, a = 1, g = 1, w = 5, u = 5)
{
 B <<- B
 a <<- a
 g <<- g
 w <<- w
 u <<- u
 K2 <- (gamma((g + n + 1)/2)/(gamma(1/2) * gamma((g + n)/2)) *
(a^(g/2)))/(gamma((n - 1)/2)) * sqrt(n/
  (1 + n * w))
 fun <- function(x)
 unlist(lapply(x, function(x)
 {
  integrate(f = function(y, x)
  {
   K2 * ((u + n * w * y)/(1 + n * w) * (x^((n - 3)/2)) * (a + x + (n * ((y -
u)^2))/(1 +
    n * w))^( - (g + n)/2))
  }
  , lower = (-16 * B * n * w + 8 * B * g * n * w + 8 * B * (n^2) * w - 2 * n
* u * w - 8 * g *
   n * u * w - 8 * (n^2) * u * w - 16 * B * (n^2) * (w^2) + 8 * B * g *
(n^2) * (w^2) +
   8 * B * (n^3) * (w^2) + 6 * sqrt(n * w * (1 + n * w)) * sqrt(-8 * (B^2) +
4 * (B^2) *
   g + 4 * (B^2) * n + 16 * B * u - 8 * B * g * u - 8 * B * n * u - 8 *
(u^2) + 4 * g * (
   u^2) + 4 * n * (u^2) - 9 * a * w - 8 * (B^2) * n * w + 4 * (B^2) * g * n
* w + 4 * (
   B^2) * (n^2) * w - 9 * x * w + 16 * B * n * u * w - 8 * B * g * n * u *
w - 8 * B * (
   n^2) * u * w - 8 * n * (u^2) * w + 4 * g * n * (u^2) * w + 4 * (n^2) *
(u^2) * w - 8 *
   a * n * (w^2) + 4 * a * g * n * (w^2) + 4 * a * (n^2) * (w^2) - 8 * n * x
* (w^2) +
   4 * g * n * x * (w^2) + 4 * (n^2) * x * (w^2)))/(2 * (-9 * n * w - 8 *
(n^2) * (w^2) +
   4 * g * (n^2) * (w^2) + 4 * (n^3) * (w^2))), upper = (1/(2 * (-49 * n *
w - 8 * (n^2
   ) * (w^2) + 4 * g * (n^2) * (w^2) + 4 * (n^3) * (w^2))) * (-16 * B * n *
w + 8 * B *
   g * n * w + 8 * B * (n^2) * w - 82 * n * u * w - 8 * g * n * u * w - 8 *
(n^2) * u *
   w - 16 * B * (n^2) * (w^2) + 8 * B * g * (n^2) * (w^2) + 8 * B * (n^3) *
(w^2) + 14 *
   sqrt(n * w * (1 + n * w)) * sqrt(-8 * (B^2) + 4 * (B^2) * g + 4 * (B^2) *
n + 16 * B *
   u - 8 * B * g * u - 8 * B * n * u - 8 * (u^2) + 4 * g * (u^2) + 4 * n *
(u^2) - 49 *
   a * w - 8 * (B^2) * n * w + 4 * (B^2) * g * n * w + 4 * (B^2) * (n^2) *
w - 49 * x *
   w + 16 * B * n * u * w - 8 * B * g * n * u * w - 8 * B * (n^2) * u * w -
8 * n * (u^
   2) * w + 4 * g * n * (u^2) * w + 4 * (n^2) * (u^2) * w - 8 * a * n *
(w^2) + 4 * a *
   g * n * (w^2) + 4 * a * (n^2) * (w^2) - 8 * n * x * (w^2) + 4 * g * n * x
* (w^2) +
   4 * (n^2) * x * (w^2)))), x = x)$integral
 }
 ))
 integrate(f = fun, lower = 0, upper = Inf)$integral
}



> benefit2(10)
Error in integrate: Missing value where logical needed: if(abs(lower) != Inf
&& abs(upper) != Inf) c(integrate.f(f, lower, upper, subdivisions, rel.tol,
abs.tol,
  keep.xy, aux, ...), list(call = match.call())) else c(integrate.i(f,
lower, upper,
 . . .



<Prev in Thread] Current Thread [Next in Thread>