How to calculate the double integration in R

Question

This is my r code to calculate beta values for each case which is pretty simple

data =data.frame(
  "t" = seq(0, 1, 0.001)
)


B3t <- function(t){
  t**3 - 1.6*t**2 +0.76*t+1
}

B2t <- function(t){

  ifelse(t >= 0 & t < 0.342,
         ((t-0.5)^2-0.025),

         ifelse( data$t >=  0.342 & data$t <= 0.658, 
                 0,

                 ifelse(t >  0.658 & t <= 1, 
                        (-(t-0.5)^2+0.025),
                        0
                 )))

}

B1t <- function(t){
  0
}


X1t <- function(t){
  a0 = rnorm(1)
  a1 = rnorm(1)
  a2 = rnorm(1)
  a3 = rnorm(1)

  return(a0 + a1*t + a2*(t^2) + a3*(t^3))
}

X2t <- function(t){
  a0 = rnorm(1)
  a1 = rnorm(1)
  a2 = rnorm(1)
  a3 = rnorm(1)
  a4 = rnorm(1)

  return(a0 + a1 * sin(2*pi*t) + a2 * cos(2*pi*t) + a3 * sin(4*pi*t) + a4 * cos(4*pi*t))
}

Now I want to calculate the error term.

I have one issue: Can anyone help me with this question?

How do I solve the double integration in order to calculate the error term.

I know there are functions in r to do integrate but I am not sure how do I implement it here.

I am trying to do functional data analysis problem mentioned below:

What I don't know is how to find the variance in order to find the error term which follows normal distribution N(0, variance)

You gan generate the standard normal distribution using the function `rnorm()`. Also, the `integral` function works with functions and limits. You don't need to calculate the values of the function, just define it: `B3t <- function(t) (t**3 - 1.6*t**2 +0.76*t+1 )` . Then the single integral of `B3t` from 0 to 1 would be `integrate(B3t,0,1)$value` — Rohit, May 22 '19 at 06:15
@jenesaisquoi I created the functions as directed. No idea what to do next. — Stupid_Intern, May 22 '19 at 06:57
@Rohit any idea what is `s` in the double integration function? Is it same as `t`? — Stupid_Intern, May 22 '19 at 06:59
For double integration, you can check out other questions on the topic on the site like : https://stackoverflow.com/questions/8913603/calculating-double-integrals-in-r-quickly — Rohit, May 22 '19 at 07:30

score 7 · Accepted Answer · answered May 28 '19 at 21:05

Here is how I would do it (for one version of beta and X). Note that a double integral is just two single integrals nested into one another. A parameter n defines the number of random samples I am using for estimating the expectation in the integral.

beta <- function(t){
    return(t*t*t-1.6*t*t+0.76*t+1)
}

myX <- function(a,t){
    pt <- c(1,t,t*t,t*t*t)
    return(sum(a*pt))
}

## computes the expectation by averaging over n samples
myE <- function(n,s,t){
    samp <- sapply(seq(n),function(x){
          a <- rnorm(4)
          myX(a,s)*myX(a,t)})
    return(mean(samp,na.rm=T))
}

## funtion inside the first integral
myIntegrand1 <- function(s,t,n){
   return(beta(s)*myE(n,s,t))
}

## function inside the second integral
myIntegrand2 <- function(t,n){
    v <- integrate(myIntegrand1,0,1,t=t,n=n)
    return(beta(t)*v$value)
}

## computes sigma
mySig <- function(n){
  v <- integrate(myIntegrand2,0,1,n=n)
  return( 0.25*v$value)
}
## tests various values of n (number of samples drawn to compute the expectation)
sapply(seq(3),function(x)
             c("100"=mySig(100),"1000"=mySig(1000),"10000"=mySig(10000)))
## output shows you the level of precision you may expect:
##       [,1]     [,2]     [,3]
##  100  48.61876 47.85445 58.2094
## 1000  52.95681 50.61860 50.61702
## 10000 54.88292 53.02073 54.48635

How to calculate the double integration in R

1 Answers1