Michael O. - 1 year ago 100
R Question

Create a function from polynomial coefficients

I've found polynomial coefficients from my data:

``````R <- c(0.256,0.512,0.768,1.024,1.28,1.437,1.594,1.72,1.846,1.972,2.098,2.4029)
Ic <- c(1.78,1.71,1.57,1.44,1.25,1.02,0.87,0.68,0.54,0.38,0.26,0.17)
NN <- 3
ft <- lm(Ic ~ poly(R,NN,raw = TRUE))
pc <- coef(ft)
``````

so I can create a polynomial function:

``````f1 <- function(x) pc[1]+pc[2]*x+pc[3]*x^2+pc[4]*x^3
``````

and, for example, take a derivative:

``````g1 <- Deriv(f1)
``````

How to create a universal function so that it doesn't have to be rewritten for every new polynomial degree NN?

My original answer may not be what you really want, as it was numerical rather symbolic. Here is the symbolic solution.

``````## use `"x"` as variable name
## taking polynomial coefficient vector `pc`
## can return a string, or an expression by further parsing (mandatory for `D`)
f <- function (pc, expr = TRUE) {
stringexpr <- paste("x", seq_along(pc) - 1, sep = " ^ ")
stringexpr <- paste(stringexpr, pc, sep = " * ")
stringexpr <- paste(stringexpr, collapse = " + ")
if (expr) return(parse(text = stringexpr))
else return(stringexpr)
}

## using R base's `D` (requiring expression)
D(f(1:4 / 10, TRUE), name = "x")
# 0.2 + 2 * x * 0.3 + 3 * x^2 * 0.4

## using `Deriv::Deriv`
library(Deriv)

Deriv(f(1:4 / 10, TRUE), x = "x")
# expression(0.2 + x * (0.6 + 1.2 * x))

Deriv(f(1:4 / 10, FALSE), x = "x")
# [1] "0.2 + x * (0.6 + 1.2 * x)"
``````

If we use expression, we will be able to evaluate the result later. For example,

``````## cubic polynomial
cubic <- f(1:4 / 10, TRUE)
eval(cubic, envir = list(x = 1:4))
# [1]  1.0  4.9 14.2 31.3

## first derivative
dcubic <- D(cubic, name = "x")  ## or `dcubic <- Deriv(cubic, x = "x")`
eval(dcubic, envir = list(x = 1:4))
# [1]  2.0  6.2 12.8 21.8
``````

A polynomial function

We can use `outer` or `poly` to help us set up ordinary polynomial basis (in fact `poly` is calling `outer`, too):

``````## use `poly`, but remember `poly` does not generate intercept column
f1 <- function (x, pc) {
X <- poly(x, degree = length(pc) - 1L, raw = TRUE)
drop(X %*% pc[-1]) + pc[1]  ## add intercept
}

## use `outer`
f2 <- function (x, pc) {
X <- outer(x, seq_along(pc) - 1L, FUN = "^")
drop(X %*% pc)  ## use `drop` to flatten matrix to vector
}
``````

Some small example:

``````x <- seq(0, 1, by = 0.25)
pc <- c(-1, 0.3, -0.7, 0.9)  ## 4 coefficients, so cubic polynomial

f1(x, pc)
# [1] -1.0000000 -0.9546875 -0.9125000 -0.7890625 -0.5000000

f2(x, pc)
# [1] -1.0000000 -0.9546875 -0.9125000 -0.7890625 -0.5000000
``````

`f1` and `f2` are both vectorized functions. We can pass in a single value of `x`, or a vector value of `x` (which is done in the example).

First derivative of a polynomial function

For polynomials, derivatives are explicitly known so we can write it ourselves. It is more convenient to use `outer` than `poly` here.

``````## use `poly`
g <- function (x, pc) {
## note, degree is `length(pc) - 1` for 1st derivative
X <- outer(x, seq_len(length(pc) - 1L) - 1L, FUN = "^")
## not we update coefficients
beta <- pc[-1] * seq_len(length(pc) - 1L)
## return
drop(X %*% beta)
}
``````

Still using previous `x` and `pc`,

``````g(x, pc)
# [1] 0.30000 0.11875 0.27500 0.76875 1.60000
``````
Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download