if (FALSE) {
##  Polynomial minimax approximation of data points
##  (see the Remez algorithm)
n <- 10; m <- 101           # polynomial of degree 10; no. of data points
xi <- seq(-1, 1, length = m)
yi <- 1 / (1 + (5*xi)^2)    # Runge's function
pval <- function(p, x)      # Horner scheme
    outer(x, (length(p) - 1):0, "^") %*% p
pfit <- function(x, y, n)   # polynomial fitting of degree n
    qr.solve(outer(x, seq(n, 0), "^"), y)
fn1 <- function(p)           # objective function
    max(abs(pval(p, xi) - yi))
pf <- pfit(xi, yi, 10)      # start with a least-squares fitting
sol1 <- pureCMAES(pf, fn1, rep(-200, 11), rep(200, 11))
zapsmall(sol1$xmin)
# [1]  -50.24826    0.00000  135.85352    0.00000 -134.20107    0.00000
# [7]   59.19315    0.00000  -11.55888    0.00000    0.93453
print(sol1$fmin, digits = 10)
# [1] 0.06546780411
##  Polynomial fitting in the L1 norm
##  (or use LP or IRLS approaches)
fn2 <- function(p)
    sum(abs(pval(p, xi) - yi))
sol2 <- pureCMAES(pf, fn2, rep(-100, 11), rep(100, 11))
zapsmall(sol2$xmin)
# [1] -21.93238   0.00000  62.91083   0.00000 -67.84847   0.00000 
# [7]  34.14398   0.00000  -8.11899   0.00000   0.84533
print(sol2$fmin, digits = 10)
# [1] 3.061810639
}
Run the code above in your browser using DataLab