r - 任意阶多项式的函数(首选符号方法)

Question

我从我的数据中找到了多项式系数:

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)

所以我可以创建一个多项式函数:

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

例如，取一个导数:

g1 <- Deriv(f1)

如何创建一个通用函数，以便不必为每个新的多项式次数重写它 NN ?

Answer 1

我最初的答案可能不是你真正想要的，因为它是数字而不是象征性的。这是象征性的解决方案。

## 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)
  }

## an example cubic polynomial with coefficients 0.1, 0.2, 0.3, 0.4
cubic <- f(pc = 1:4 / 10, TRUE)

## using R base's `D` (requiring expression)
dcubic <- D(cubic, name = "x")
# 0.2 + 2 * x * 0.3 + 3 * x^2 * 0.4

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

dcubic <- Deriv(cubic, x = "x", nderiv = 1L)
# expression(0.2 + x * (0.6 + 1.2 * x))

Deriv(f(1:4 / 10, FALSE), x = "x", nderiv = 1L)  ## use string, get string
# [1] "0.2 + x * (0.6 + 1.2 * x)"

当然， Deriv使高阶导数更容易得到。我们可以简单地设置 nderiv .对于 D但是，我们必须使用递归(参见 ?D 的示例)。

Deriv(cubic, x = "x", nderiv = 2L)
# expression(0.6 + 2.4 * x)

Deriv(cubic, x = "x", nderiv = 3L)
# expression(2.4)

Deriv(cubic, x = "x", nderiv = 4L)
# expression(0)

如果我们使用表达式，我们将能够稍后评估结果。例如，

eval(cubic, envir = list(x = 1:4))  ## cubic polynomial
# [1]  1.0  4.9 14.2 31.3

eval(dcubic, envir = list(x = 1:4))  ## its first derivative
# [1]  2.0  6.2 12.8 21.8

上面意味着我们可以为一个函数包装一个表达式。使用函数有几个优点，一个是我们可以使用 curve 来绘制它。或 plot.function .

fun <- function(x, expr) eval.parent(expr, n = 0L)

注意，成功 fun需要 expr以符号表示 x .如 expr是根据 y 定义的例如，我们需要定义 fun与 function (y, expr) .现在让我们使用 curve绘制 cubic和 dcubic , 范围 0 < x < 5 :

curve(fun(x, cubic), from = 0, to = 5)  ## colour "black"
curve(fun(x, dcubic), add = TRUE, col = 2)  ## colour "red"

最方便的方式当然是定义单个函数 FUN而不是做 f + fun组合。这样，我们也不必担心 f 使用的变量名的一致性。和 fun .

FUN <- function (x, pc, nderiv = 0L) {
  ## check missing arguments
  if (missing(x) || missing(pc)) stop ("arguments missing with no default!")
  ## expression of polynomial
  stringexpr <- paste("x", seq_along(pc) - 1, sep = " ^ ")
  stringexpr <- paste(stringexpr, pc, sep = " * ")
  stringexpr <- paste(stringexpr, collapse = " + ")
  expr <- parse(text = stringexpr)
  ## taking derivatives
  dexpr <- Deriv::Deriv(expr, x = "x", nderiv = nderiv)
  ## evaluation
  val <- eval.parent(dexpr, n = 0L)
  ## note, if we take to many derivatives so that `dexpr` becomes constant
  ## `val` is free of `x` so it will only be of length 1
  ## we need to repeat this constant to match `length(x)`
  if (length(val) == 1L) val <- rep.int(val, length(x))
  ## now we return
  val
  }

假设我们要计算系数为 pc <- c(0.1, 0.2, 0.3, 0.4) 的三次多项式及其衍生产品 x <- seq(0, 1, 0.2) ，我们可以简单地做:

FUN(x, pc)
# [1] 0.1000 0.1552 0.2536 0.4144 0.6568 1.0000

FUN(x, pc, nderiv = 1L)
# [1] 0.200 0.368 0.632 0.992 1.448 2.000

FUN(x, pc, nderiv = 2L)
# [1] 0.60 1.08 1.56 2.04 2.52 3.00

FUN(x, pc, nderiv = 3L)
# [1] 2.4 2.4 2.4 2.4 2.4 2.4

FUN(x, pc, nderiv = 4L)
# [1] 0 0 0 0 0 0

现在绘图也很容易:

curve(FUN(x, pc), from = 0, to = 5)
curve(FUN(x, pc, 1), from = 0, to = 5, add = TRUE, col = 2)
curve(FUN(x, pc, 2), from = 0, to = 5, add = TRUE, col = 3)
curve(FUN(x, pc, 3), from = 0, to = 5, add = TRUE, col = 4)