Michael O. Michael O. - 1 month ago 20
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?

Answer

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

Original Answer (for numerical function)

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