我已经编写了一个简单的数值示例(稍后我将使其更加复杂).你有两种物质A和B.在每一轮(即一个过程的迭代)中,它们以r的速率相互消耗. 我有一些东西是有效的(见下面的介绍),但它相当难看.首先,它显式地依赖于一个循环(也许有一个使用MAP的解决方案),然后有两个变量包含初始状态,它们被覆盖.

最后,还有100条:物质A和B不能取负值,所以一旦出现负值,这个过程就必须停止. 请看邮件末尾的复印件. 任何改进建议都是受欢迎的.

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

A <- 2000
B <- 1000 ## initial quantities of A and B

r <- 1/10 ## consumption rate

n_rounds <- 10 ## how many iterations

consumption_round <- function(ab_vec,r){ ## function determining how much
 ## A and B decreas at each round

    A1 <- ab_vec[1]-ab_vec[2]*r
    B1 <- ab_vec[2]-ab_vec[1]*r

    res <- c(A1, B1)

    return(res)
}


AB <- c(A,B) ## initial state


ini <- AB

list_res <- AB  ## other variables storing the initial state which will be overwritten ## I would like a better way to do this


for (i in seq(n_rounds)){

    print("i is,")
    print(i)
    
    output <- consumption_round(ini, r)

    ini <- output

    list_res <- rbind(list_res, output)
    
}
#> [1] "i is,"
#> [1] 1
#> [1] "i is,"
#> [1] 2
#> [1] "i is,"
#> [1] 3
#> [1] "i is,"
#> [1] 4
#> [1] "i is,"
#> [1] 5
#> [1] "i is,"
#> [1] 6
#> [1] "i is,"
#> [1] 7
#> [1] "i is,"
#> [1] 8
#> [1] "i is,"
#> [1] 9
#> [1] "i is,"
#> [1] 10

sim_res <- list_res |>
    as_tibble(.name_repair="unique") |>
    rename("A"="...1",
           "B"="...2") |>
    mutate(B=if_else(B>=0, B, 0)) |> ## B cannot be negative
    slice(1:min(which(B==0))) ## As soon as B is zero, the simulation ends.
#> New names:
#> • `` -> `...1`
#> • `` -> `...2`

sim_res
#> # A tibble: 7 × 2
#>       A      B
#>   <dbl>  <dbl>
#> 1 2000  1000  
#> 2 1900   800  
#> 3 1820   610  
#> 4 1759   428  
#> 5 1716.  252. 
#> 6 1691.   80.5
#> 7 1683.    0

创建于2024-01-18年第reprex v2.0.2

推荐答案

Iteration With Termination

如果您想要有一些终止条件,您可以try

AB <- cbind(A, B)
r <- 1 / 10
alpha <- (1 + r) * diag(2) + matrix(-r, 2, 2)
res <- AB
repeat {
    if (min(AB) == 0) break
    AB <- AB %*% alpha
    res <- rbind(res, pmax(AB, 0))
}

或者用svd进行奇异值分解

AB <- cbind(A, B)
r <- 1 / 10
alpha <- (1 + r) * diag(2) + matrix(-r, 2, 2)
X <- svd(alpha)
res <- AB
k <- 1
repeat {
    v <- AB %*% with(X, u %*% diag(d^k) %*% t(v))
    res <- rbind(res, pmax(v, 0))
    if (min(v) <= 0) {
        break
    }
    k <- k + 1
}

以至于

> res
            A       B
[1,] 2000.000 1000.00
[2,] 1900.000  800.00
[3,] 1820.000  610.00
[4,] 1759.000  428.00
[5,] 1716.200  252.10
[6,] 1690.990   80.48
[7,] 1682.942    0.00

Iteration Without Termination

您可以try 如下所示的Reduce,其中使用矩阵alpha迭代地更新值AB

AB <- cbind(A, B)
r <- 1 / 10
alpha <- (1 + r) * diag(2) + matrix(-r, 2, 2)
do.call(
    rbind,
    Reduce(`%*%`, rep(list(alpha), 10), init = AB, accumulate = TRUE)
)

这给了我们

             A         B
 [1,] 2000.000 1000.0000
 [2,] 1900.000  800.0000
 [3,] 1820.000  610.0000
 [4,] 1759.000  428.0000
 [5,] 1716.200  252.1000
 [6,] 1690.990   80.4800
 [7,] 1682.942  -88.6190
 [8,] 1691.804 -256.9132
 [9,] 1717.495 -426.0936
[10,] 1760.105 -597.8431
[11,] 1819.889 -773.8536

R相关问答推荐

使用ggcorrplot在相关性矩阵上标注supertitle和index标签

从嵌套列表中智能提取线性模型系数

ggplot geom_smooth()用于线性回归虚拟变量-没有回归线

如何根据条件计算时差(天)

如何在格子中添加双曲曲线

如何得到R中唯一的组合群?

当月份额减go 当月份额

在ggplot2中更改小提琴情节的顺序

R中边际效应包中Logistic回归的交互作用风险比

当我们有多个反斜杠和/特殊字符时使用Gsubing

函数可以跨多个列搜索多个字符串并创建二进制输出变量

在R中,我如何使用滑动窗口计算位置,然后进行过滤?

根据约束随机填充向量的元素

手动指定从相同数据创建的叠加图的 colored颜色

需要一个函数来在第一行创建一个新变量,然后用新变量替换一个不同的变量(对于多行)

R将函数参数传递给ggploy

如何使用循环从R中的聚合函数创建列,而不会在名称中给出&q;$&q;?

在直方图中显示两个变量

使用dplyr删除具有条件的行

如何根据每个子框架中分类因子的唯一计数来过滤子框架列表?