尝试在循环中递归地定义、计算和积分r函数

cl25kdpy 于 2023-02-14 发布在其他

关注(0)|答案(2)|浏览(128)

This is the intergral I'm trying to evaluate
我要做的是把g_k（x，s）放到一个列表中，这样我就可以在循环中递归地使用它们。

n = 500
k = 3
g_1 <- 1
rho <- sqrt( ((k-1)*(n-k))/(k*(n-k+1)) )
tau2 <- (sqrt( n / (k*(n-k+1))))

g_lst <- c()

#Coded discretely for k=1,2,3,4

g1 <- function(x,s){ return( 1 ) }
g_lst <- c(g_lst,g1)

g2 <- function(x,s){ return( integrate(function(y){g_lst[[1]](y,s)*(dnorm(y,rho*x,tau2)+dnorm(y,-rho*x,tau2))},0,s)[[1]] ) }
g_lst <- c(g_lst,g2)

g3 <- function(x,s){ return( integrate(function(y){g_lst[[2]](y,s)*(dnorm(y,rho*x,tau2)+dnorm(y,-rho*x,tau2))},0,s)[[1]] ) }
g_lst <- c(g_lst,g3)

g4 <- function(x,s){ return( integrate(function(y){g_lst[[3]](y,s)*(dnorm(y,rho*x,tau2)+dnorm(y,-rho*x,tau2))},0,s)[[1]] ) }
g_lst <- c(g_lst,g4)

#Trying to generalise to k=1,2,...,k
g_lst2 <- c()
g1 <- function(x,s){ return( 1 ) }
g_lst2 <- c(g_lst2,g1)

for (i in 1:3){
  
  i <- force(i)
  gk <- function(x,s){ return( integrate(function(y){g_lst[[i]](y,s)*(dnorm(y,rho*x,tau2)+dnorm(y,-rho*x,tau2))},0,s)[[1]] ) }
  force(gk)
  g_lst2 <- c(g_lst2,gk)
}

下面是我从列表中评估函数得到的相应值。g_lst给了我正确的值，而g_lst2（对于所有i〉=1）g_lsti给了我g_lst2[4]的值。从我在堆栈交换中找到的线程中，我觉得我需要使用force（）函数，但我使用它的方式没有帮助。

g_lst2[[2]](2,2)
#[1] 1.811419
g_lst2[[3]]](2,2)
#[1] 1.811419
g_lst2[[4]](2,2)
#[1] 1.811419

g_lst[[2]](2,2)
#[1] 0.7380149
g_lst[[3]](2,2)
#[1] 1.156224
g_lst[[4]](2,2)
#[1] 1.811419

来源：https://stackoverflow.com/questions/75413304/trying-to-define-evaluate-and-integrate-r-functions-recursively-in-loop

2条答案

按热度按时间

tvz2xvvm1#

在for循环中使用force()实际上没有什么帮助，因为for不会创建一个新的作用域，您需要创建一个生成函数来构建列表。

#Trying to generalise to k=1,2,...,k
g_lst2 <- c()
g1 <- function(x,s){ return( 1 ) }
g_lst2 <- c(g_lst2,g1)

for (i in 1:3){
  make_gk <- function(i) {
    force(i)
    function(x,s){ return( integrate(function(y){g_lst2[[i]](y,s)*(dnorm(y,rho*x,tau2)+dnorm(y,-rho*x,tau2))},0,s)[[1]] ) }
  }
  gk <- make_gk(i)
  g_lst2 <- c(g_lst2,gk)
}

注意，我们在函数内部调用force()，因为这会创建一个新的作用域。

g_lst2[[2]](2,2)
# [1] 0.7380149
g_lst2[[3]](2,2)
# [1] 1.156224
g_lst2[[4]](2,2)
# [1] 1.811419

或者，您可以使用lapply，它将为您处理强制

g_lst2 <- c()
g1 <- function(x,s){ return( 1 ) }
g_lst2 <- c(g_lst2,g1)

g_lst2 <- c(g_lst2, lapply(1:3, function(i) 
  function(x,s){ return( integrate(function(y){g_lst2[[i]](y,s)*(dnorm(y,rho*x,tau2)+dnorm(y,-rho*x,tau2))},0,s)[[1]] ) }
))

赞(0）回复(0）举报 2023-02-14

4bbkushb2#

你可以很容易地使用递归来泛化。下面是两种方法

n <- 500
k <- 3
rho <- sqrt( ((k-1)*(n-k))/(k*(n-k+1)) )
tau2 <- (sqrt( n / (k*(n-k+1))))


g <- function(x,s,k){
  if(k == 1) 1
  else 
    integrate(function(y)g(y,s, k-1)*(
      dnorm(y, x*rho, tau2) + dnorm(y, -x*rho, tau2)),0,s)[[1]]
}

g(2,2,1)
#> [1] 1
g(2,2,2)
#> [1] 0.7380149
g(2,2,3)
#> [1] 1.156224
g(2,2,4)
#> [1] 1.811419

第二种方法是返回一个函数，然后在所需的点计算该函数

g1 <- function(k){
  if(k == 1) function(...) 1
  else  function(x,s)
    integrate(function(y)g1(k-1)(y,s)*(
      dnorm(y, x*rho, tau2) + dnorm(y, -x*rho, tau2)),0,s)[[1]]
}
g1(1)(2,2)
#> [1] 1
g1(2)(2,2)
#> [1] 0.7380149
g1(3)(2,2)
#> [1] 1.156224
g1(4)(2,2)
#> [1] 1.811419

创建于2023年2月10日，使用reprex v2.0.2
当然，这个方法是缓慢的，因为它不保存中间结果。即如果你计算g(2,2,10)，然后如果你需要g(2,2,8)，你不应该计算它，而是从已经计算过的g(2,2,10)表中读取它。但上面提供的方法确实再次计算该值。我们可以跳过这个过程使用memoization：

赞(0）回复(0）举报 2023-02-14

我来回答

尝试在循环中递归地定义、计算和积分r函数

2条答案

相关问题

热门标签

最新问答