Michael O. - 1 year ago 100

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?

Recommended for you: Get network issues from **WhatsUp Gold**. **Not end users.**

Answer Source

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