如何计算R中的二重积分
How to calculate the double integration in R
这是我的 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))
}
现在我要计算误差项。
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.
我知道 r 中有函数可以做 integrate
但我不确定如何在这里实现它。
我正在尝试做下面提到的功能数据分析问题:
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)
以下是我的做法(针对一个版本的 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
这是我的 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))
}
现在我要计算误差项。
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.
我知道 r 中有函数可以做 integrate
但我不确定如何在这里实现它。
我正在尝试做下面提到的功能数据分析问题:
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)
以下是我的做法(针对一个版本的 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