如何计算R中的双重积分

Posted 2023-03-24

【中文标题】如何计算R中的双重积分

【英文标题】：How to calculate the double integration in R

【发布时间】：2019-10-08 13:02:04

【问题描述】：

这是我的 r 代码，用于计算每种情况的 beta 值，非常简单

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))

现在我要计算误差项。

我有一个问题： 谁能帮我解答这个问题？

如何解决双重积分以计算误差项。

我知道 r 中有函数可以做integrate，但我不确定如何在这里实现它。

我正在尝试做下面提到的功能数据分析问题：

我不知道如何找到方差以找到服从正态分布 N(0, variance) 的误差项

【问题讨论】：

您可以使用函数rnorm() 生成标准正态分布。此外，integral 函数适用于函数和限制。您不需要计算函数的值，只需定义它：B3t <- function(t) (t**3 - 1.6*t**2 +0.76*t+1 )。那么B3t从0到1的单积分就是integrate(B3t,0,1)$value

@jenesaisquoi 我按照指示创建了函数。不知道下一步该做什么。

@Rohit 知道双重积分函数中的s 是什么吗？和t一样吗？

对于双重集成，您可以在网站上查看有关该主题的其他问题，例如：***.com/questions/8913603/…

【参考方案1】：

这就是我的做法（对于 beta 和 X 的一个版本）。请注意，双积分只是两个相互嵌套的单积分。参数n 定义了我用于估计积分期望的随机样本数。

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