0

I need to evaluate an integral in the following form:

\int_a^b f(x) \int_0^x g(t)(x-t)dtdx

Can you please suggest a way? I assume that this integral can't be done in the standard approach suggested in the following answer:

Standard approach

Update: Functions are added in the following image. f(x) basically represents a pdf of a uniform distribution but the g(t) is a bit more complicated. a and b can be any positive real numbers.

enter image description here

Rel_Ai
  • 581
  • 2
  • 11
  • Why do you assume it can't be done like in the answers to the question you link to? Also, can you give an example of `g()` and of the integration limits a, b? – Rui Barradas Jul 17 '22 at 06:52
  • @Rui Barradas I have added the functions f(t) and g(t). The reason I think it is different from the answer I refer to is because I believe it is not valid to transform the original function to \int_a^b \int_0^x f(x)g(t)(x-t)dtdx. May be I am wrong? – Rel_Ai Jul 17 '22 at 07:07

1 Answers1

1

The domain of integration is a simplex (triangle) with vertices (a,a), (a,b) and (b,b). Use the SimplicialCubature package:

library(SimplicialCubature)

alpha <- 3
beta <- 4
g <- function(t){
  ((beta/t)^(1/2) + (beta/t)^(3/2)) * exp(-(t/beta + beta/t - 2)/(2*alpha^2)) / 
    (2*alpha*beta*sqrt(2*pi))
}
a <- 1
b <- 2
h <- function(tx){
  t <- tx[1]
  x <- tx[2]
  g(t) * (x-t)
}

S <- cbind(c(a, a), c(a ,b), c(b, b))
adaptIntegrateSimplex(h, S)
# $integral
# [1] 0.01962547
# 
# $estAbsError
# [1] 3.523222e-08

Another way, less efficient and less reliable, is:

InnerFunc <- function(t, x) { g(t) * (x - t) }
InnerIntegral <- Vectorize(function(x) { integrate(InnerFunc, a, x, x = x)$value})
integrate(InnerIntegral, a, b)
# 0.01962547 with absolute error < 2.2e-16
Stéphane Laurent
  • 75,186
  • 15
  • 119
  • 225
  • I wonder where is the part with f(x) here. Is it missing from your solution? I am struggling to understand it. – Rel_Ai Jul 17 '22 at 18:09
  • 1
    Ah yes, just divide by **b-a**. – Stéphane Laurent Jul 17 '22 at 22:50
  • One more thing, in integrate(InnerFunc, a, x, x = x)$value, shouldn't it be 0, x instead of a,x? My understanding is, it is the inner interval – Rel_Ai Jul 18 '22 at 04:54
  • Ah yes you're right, sorry. I'll do the corrections later. – Stéphane Laurent Jul 18 '22 at 05:33
  • There is a funny thing happening. As I understand, in your second approach, you can replace x with any other letter you wish. But when I used 'l' (to be consistent with my own notations) the result comes negative. I cleared all the variables before executing but the same. I tried the letter 'v' or 'z', and it works. This is rather spooky to me and I can't seem to find any logical reason for that. – Rel_Ai Jul 18 '22 at 05:39
  • This is correct actually, because f(x)=0 for x < a. – Stéphane Laurent Jul 18 '22 at 05:40
  • If you have some time can you please reflect on this problem? https://stackoverflow.com/questions/73135559/double-integration-with-a-differentiation-inside-in-r – Rel_Ai Jul 27 '22 at 09:38