Итерационный расчет в R

Question

Допустим, у нас есть некоторые переменные:

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

Как бы вы настроили r для решения этой проблемы?В Excel вы можете создать целевую ячейку, давайте назовем ее y равной нулю с помощью Goal Seek или Solver ... и это потребует итерационных вычислений.

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

Excel вычислит px, уменьшив y до нуля.

Здесь px = 9.45945

Как бы это обработать?

Answer 1

После некоторой алгебры решение в закрытой форме px - s/k-1 * r - g == 0 имеет вид

px <- (10*k*r + 10*s + k*s)/(11*k)
# [1] 9.545455

Решение с использованием пакета nleqslv:

library(nleqslv)

fn <- function(px) {
  s <- 95
  k <- 100
  r <- .05
  sx <- s - px
  g <- sx * .1
  y <- px - s/k-1 * r - g
  return(y)
}

nleqslv(1, fn)

# $x
# [1] 9.545455
# ...

Answer 2

Для начала поясним, что вопрос состоит в том, чтобы решить систему из трех переменных с тремя неизвестными.То есть:

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, заменив sx в первом уравнении на s - px, а затем вычтите два оставшихся уравнения, чтобы исключить 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 )"

Answer 3

Вот базовое решение R, которое должно работать до тех пор, пока оно может быть разрешено в линейное уравнение.

s = 95
k = 100
r = .05

g = function(sx) {
    sx * 0.1
}

sx = function(px) {
    s - px
}

fn = function(px) {
    px - s/k - 1 * r - g(sx(px))
}

foo = function(fn) {
    x1 = 1
    x2 = 2
    y1 = fn(x1)
    y2 = fn(x2)
    m = (y2 - y1)/(x2 - x1)  # Line slope
    solve(-m, -m * x1 + y1)
}

foo(fn)
#> [1] 9.545455