r - R中的迭代计算

Question

假设我们有一些变量:

s = 95
k = 100
r = .05
g = sx * .1
sx = s - px
px = need to solve for this

你会如何设置 r 来解决这个问题？在 Excel 中，您可以创建一个目标单元格，让我们通过 Goal Seek 或 Solver 将其称为 y 为零。这涉及到迭代计算。

y = px - s/k-1 * r - g

Excel 会通过将 y 减为零来计算 px。

这里px = 9.45945

你会如何处理这个？

Answer 1

首先让我们澄清一下，问题是要求求解一个包含三个未知数的三个变量的系统。即:

Solve this system of 3 equations in the 3 unknowns px, sx and g:

       0.1 * sx      - g  = 0 
 px +        sx           = s 
 px                  - g  = s/k-1 * r

where s, k and r are known constants given in the question and we have 
rearranged the equations to put the constant terms on the right hand
side and have vertically aligned the variables on the left hand side.

这个方程组实际上很容易手工求解(通过将第一个方程中的 sx 替换为 s - px 来消除 sx，然后减去剩余的两个方程以消除 g。最后我们留下一个未知数中的一个方程很容易解决)。由于此处的目标是使用 R，因此我们展示了几种不同的使用 R 的方法。

1) 线性代数 由于这些是线性方程，我们可以将其表示为矩阵问题并求解。矩阵 m 由上述等式的系数形成。 m的第一列是px的系数，下一列是sx的系数，最后一列是的系数>g。例如，px的系数在三个方程中分别为0、1、1。右侧向量是 rhs，solve 可用于在概念上反转 m 并将其乘以 rhs。

不使用包，不涉及显式迭代计算，也不需要手动将一个公式替换为另一个公式。我们按原样处理这三个方程，让计算机求解。

s <- 95; k <- 100; r <- .05
m <- matrix(c(0, 1, 1,   0.1, 1, 0,   -1, 0, -1), 3)
rhs <- c(0, s, s/k-1 * r)

solve(m, rhs)
## [1]  9.454545 85.545455  8.554545

上述输出向量中的值为px、sx和g的值。

2) 定点迭代 如果我们有一个 px 的估计，那么从第二个方程我们可以得到 sx 的估计值并使用它第一个方程，我们可以得到 g 的估计值，从中我们可以使用第三个方程得到 px 的新估计值。将其写为函数 f，我们从初始估计值 1 开始，然后迭代它直到值收敛。没有使用包。

f <- function(px) {
  sx <- s - px
  g <- .1 * sx
  s/k-1 * r + g 
}

px <- 1
for(i in 1:50) {
  prev <- px
  px <- f(px)
  cat(i, px, "\n")
  if (abs(px - prev) < 1e-5) break
}

## 1 10.3 
## 2 9.37 
## 3 9.463 
## 4 9.4537 
## 5 9.45463 
## 6 9.454537 
## 7 9.454546 

2a) optimize 我们可以交替使用 optimize(在 R 的基础上)来最小化 f(x) 和 x px 的给定边界:

optimize(function(px) (f(px) - px)^2, c(0, 100))
## $minimum
## [1] 9.454545
##
## $objective
## [1] 3.155444e-30

2b) uniroot 或使用 uniroot(也是基于 R)求 f(px) - px 的根:

uniroot(function(px) f(px) - px, c(0, 100))
## $root
## [1] 9.454545
## 
## $f.root
## [1] 0
##
## $iter
## [1] 1
##
## $init.it
## [1] NA
##
## $estim.prec
## [1] 90.54545

3) Ryacas 下面我们使用 Ryacas 包进行设置，它允许直接对三个方程进行符号说明。

# devtools::install_github("mikldk/ryacas0")

library(Ryacas0)

s <- 95
k <- 100
r <- .05

sx <- Sym("sx")
px <- Sym("px")
g <- Sym("g")

res <- Solve(List(g == sx * 0.1, sx == s - px, px - s/k-1 * r - g == 0), 
   List(px, sx, g))

给予:

res
## Yacas matrix:
##      [,1]            [,2]           [,3]        
## [1,] px == 1.05/0.11 sx == 9.4/0.11 g == 9.4/1.1

以上是作为符号表达式的精确值，但我们可以计算它们以获得浮点近似值，如下所示:

Eval(res)
## "( px == 9.54545454545454 )" "( sx == 85.4545454545455 )" "( g == 8.54545454545454 )"