0
set.seed(1)

n=100 
s=1.5 
xbar=18 

ll=2.9307
ul=2.9768


fun = function(w) { 
  d=w^((n/2)-1)
  e=(t-xbar-(((n-1)*s^2)/(2*w)))^2
  g=n/((n-1)*s^2)
  h=(e*g+1)*w
  
  d*exp(-h/2)
}

Int_w = Vectorize(function(t) { integrate(fun, 0, Inf)$value })
integrate(Int_w , ll, ul)

#> Error in t - xbar : non-numeric argument to binary operator

I am trying to perform a double integral in R. fun is my function here. w and t are my variables. Integrations are done in the order: dw, dt. Why am I getting the error at the end (last line of the posted code above)?

Phil
  • 7,287
  • 3
  • 36
  • 66
Dihan
  • 311
  • 1
  • 7
  • `t` is not defined in the function nor outside it. The inner integral's function argument is `w` not `t`. Are you looking for [this](https://stackoverflow.com/questions/43189512/double-integral-in-r)? – Rui Barradas Sep 03 '22 at 17:05
  • Int_w is the function of t, isn't it? – Dihan Sep 03 '22 at 17:13
  • No, it is not, it's a function of `w`. When I changed the `t` in Vectorize to `w` and the same in the function there's no error, I got `6.438021e-61 with absolute error < 7.1e-75`. – Rui Barradas Sep 03 '22 at 17:26
  • May I see your code? – Dihan Sep 03 '22 at 17:36

1 Answers1

0

Here is a solution with package cubature, that has functions for multidimensional integration over hypercubes. Note that the integrand is a function of one variable, in this case a 2-dim vector.

library(cubature)

fun <- function(x) {
  t <- x[1]   # outer integral
  w <- x[2]   # inner integral
  d <- w^((n/2)-1)
  e <- (t - xbar - (((n - 1)*s^2)/(2*w)))^2
  g <- n/((n - 1)*s^2)
  h <- (e*g + 1)*w
  d*exp(-h/2)
}

n <- 100 
s <- 1.5 
xbar <- 18 

ll <- 2.9307
ul <- 2.9768

hcubature(fun, c(ll, 0), c(ul, Inf))
#> $integral
#> [1] 0
#> 
#> $error
#> [1] 0
#> 
#> $functionEvaluations
#> [1] 85
#> 
#> $returnCode
#> [1] 0

pcubature(fun, c(ll, 0), c(ul, Inf))
#> $integral
#> [1] NaN
#> 
#> $error
#> [1] 0
#> 
#> $functionEvaluations
#> [1] 9
#> 
#> $returnCode
#> [1] 0

Created on 2022-09-03 by the reprex package (v2.0.1)

Rui Barradas
  • 70,273
  • 8
  • 34
  • 66
  • No. The function *fun* has two variables *w* and *t*. How can you write the function in terms of *w* only? You just can not replace *t* with *w*. – Dihan Sep 03 '22 at 21:35
  • @Dihan OK, but then, you only have one value of each `ll` and `ul`. They are the limits of what integral, the inner or the outer one? And the other integral? – Rui Barradas Sep 03 '22 at 21:52
  • outer integral limits are **ll** and **ul**. 0 and *Inf* are the limits for inner integrals. – Dihan Sep 03 '22 at 22:02
  • @Dihan Done, see now. – Rui Barradas Sep 03 '22 at 22:23