Michael O. - 3 months ago 32

R Question

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
```

**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
```