R：改进实施简单模型

Question

我编写了一个简单的数值示例（稍后我将使其变得更复杂）。你有两种物质 A 和 B。在每一轮（即过程的迭代），它们以速率 r 相互消耗。 我有一些有效的东西（参见下面的代表），但它相当丑陋。首先，它明确依赖于循环（也许有一个带有映射的解决方案），然后有两个包含被覆盖的初始状态的变量。最后，没有停止条件：物质 A 和 B 不能取负值，因此一旦发生这种情况，过程就必须停止。 请查看帖子末尾的 reprex。 欢迎任何改进建议。

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}

Answer 1

终止迭代

如果你想要一些终止条件，你可以尝试

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

这样

> 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

---

不断迭代

您可以尝试如下所示的

Reduce

，其中矩阵

alpha

用于迭代更新

A

和

B

的值

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