Iterative Calculation in R

Question

Let's say we have some variables:

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

How would you set up r to solve this? In Excel, you can make a target cell, let's call it y to zero via Goal Seek or Solver.. and this involves an iterative calculation.

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

Excel would calculate px by reducing y to zero.

Here px = 9.45945

How would r handle this?

Answer 1

First let us clarify that the question is asking to solve a system of three variables in three unknowns. That is:

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.

This system of equations is actually easy to solve by hand (eliminate sx by replacing sx in the first equation with s - px and then subtract the two remaining equations to eliminate g. Finally we are left with one equation in one unknown which is readily solved). Since the objective here is to use R we instead show several different ways to use R.

1) linear algebra Since these are linear equations we can formulate it as a matrix problem and solve it. The matrix m is formed from the coefficients of the equations above. The first column of m is the coefficients of px, the next column is the coefficients of sx and the final column is the coefficients of g. For example, the coefficients of px are 0, 1, and 1 in the three equations respectively. The right hand side vector is rhs and solve can be used to conceptually invert m and multiply it by rhs.

No packages are used, no explicit iterative calculation is involved and no manual substitution of one formula into another is required. We deal with the three equations as is and let the computer solve them.

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

The value in the above output vector are the values of px, sx and g.

2) fixed point iteration If we have an estimate of px then from the second equation we could get an estimate of sx and using that from the first equation we could get an estimate of g and from that we could get a new estimate of px using the third equation. Writing that as the function f we start off with an initial estimate of 1 and then iterate it until the values converge. No packages are used.

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 We could alternately use optimize (in the base of R) to minimize the distance between f(x) and x given bounds for px:

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

2b) uniroot or use uniroot (also base R) to find the root of 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 Below we set this up using the Ryacas package which allows direct symbolic specification of the three equations.

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

giving:

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

The above are exact values as symbolic expressions but we can evaluate them to get floating point approximations like this:

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

Answer 2

Here's a base R solution that should work as long as it is something that can be resolved to linear equation.

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

Answer 3

After some algebra, the solution in closed form of px - s/k-1 * r - g == 0 is

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

A solution using the nleqslv package:

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
# ...