绘制函数和导函数答案

【问题标题】：Plotting a function and a derivative function绘制函数和导函数
【发布时间】：2015-04-21 17:35:14
【问题描述】：

我想绘制一个数据框 (X,Y) data 以及一个拟合函数和拟合函数的导数。

fit <- lm(data$Y ~ poly(data$X,32,raw=TRUE))
data$fitted_values <- predict(fit, data.frame(x=data$X))

据我所知，这给了我一个 32 次多项式函数fit，我用它来计算函数值并将它们存储在data$fitted 中。绘制这些系列就像ggplot2 的魅力。

ggplot(data, aes(x=X)) + 
    geom_line(aes(y = Y), colour="red") + 
    geom_line(aes(y = predict), colour="blue")

到目前为止一切顺利。但我也想绘制拟合函数fit 的一阶导数data$Y'。我感兴趣的是拟合函数的梯度。

我的问题：如何获得fit的导数函数？我假设我可以“预测”之后绘制的绝对值。对吗？

【问题讨论】：

您能否提供一个可重现的示例，例如通过执行dput(data) 并在结果中进行编辑？这将使演示解决方案变得更加容易。（顺便说一句，您可能希望将Y 显示为geom_point 而不是geom_line！）
您可以使用this answer 逐点估算导数。或者，由于您知道拟合的形式（系数为 coef(fit) 的 32 次多项式），您可以编写一个简单的函数来手动求导，这对于多项式来说非常简单。

标签： r plot derivative

【解决方案1】：

首先，我将创建一些“有点”像你的测试数据

set.seed(15)
rr<-density(faithful$eruptions)
dd<-data.frame(x=rr$x)
dd$y=rr$y+ runif(8,0,.05)

fit <- lm(y ~ poly(x,32,raw=TRUE), dd)
dd$fitted <- fitted(fit)

ggplot(dd, aes(x=x)) + 
    geom_line(aes(y = y), colour="red") + 
    geom_line(aes(y = fitted), colour="blue")

然后，因为您有一个特殊形式的多项式，我们可以通过将每个系数乘以幂并将所有项向下移动来轻松计算导数。这是一个计算新系数的辅助函数

deriv_coef<-function(x) {
    x <- coef(x)
    stopifnot(names(x)[1]=="(Intercept)")
    y <- x[-1]
    stopifnot(all(grepl("^poly", names(y))))
    px <- as.numeric(gsub("poly\\(.*\\)","",names(y)))
    rr <- setNames(c(y * px, 0), names(x))
    rr[is.na(rr)] <- 0
    rr
}

我们可以使用...

dd$slope <- model.matrix(fit) %*% matrix(deriv_coef(fit), ncol=1)

现在我可以绘图了

ggplot(dd, aes(x=x)) + 
    geom_line(aes(y = y), colour="red") + 
    geom_line(aes(y = fitted), colour="blue") + 
    geom_line(aes(y = slope), colour="green")

我们可以看到拐点对应于导数为零的地方。

【讨论】：

谢谢！这就是我要找的东西！

【解决方案2】：

您可以通过首先根据X 对数据进行排序，然后找出每对连续值之间的差异来近似导数。

data <- d[order(d$X), ]
data$derivative = c(diff(d$fitted_values) / diff(d$X), NA)

（请注意我是如何在末尾添加NA 的，因为考虑差异会使它稍微短一些）。之后你可以绘制这个：

ggplot(data, aes(X, derivative)) + geom_line()

【讨论】：

谢谢大卫！这是可行的，并且是解决手头问题的务实解决方案！而且，据我了解，MrFlick 的解决方案是数学上正确的方法。

【解决方案3】：

据称 quantchem 包可以通过导数功能做到这一点。

说明

计算给定 x 多项式的导数。

用法

导数(obj, x)

参数

obj：'lm' 类的对象，适合 y ~ x + I(x^2) + I(x^3) + ... 方式。

x：x 值的向量

示例

x = 1:10 y = 抖动(x+x^2)

拟合 = lm(y~x+I(x^2))

导数（拟合，1:10）

Source

注意：所有这些都说了，它对我和我的数据都不起作用。

【讨论】：