2015-04-21 55 views
2

我想繪製一個數據幀(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") 

Plot

到目前爲止好。但我想繪製的也是擬合函數fit的第一個導數data$Y'。我感興趣的是擬合函數的梯度。

我的問題:我怎樣才能得到fit的微分函數? 我假設我可以「預測」事後繪圖的絕對值。正確?

+2

您能否提供一個可重複的例子,例如通過執行'dput(data)'並在結果中進行編輯?這將使演示解決方案變得更加容易。 (順便說一句,你可能想把'Y'顯示爲'geom_point'而不是'geom_line'!) –

+0

你可以用[這個答案](http://math.stackexchange.com/a)逐點估計導數/ 304664)。或者,由於您知道擬合的形式(32階多項式,係數爲'coef(fit)'),所以您可以編寫一個簡單的函數來手動獲取對多項式非常簡單的導數。 –

回答

2

首先,我將創建一些測試數據說,「那種」看起來像你

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") 

enter image description here

然後,因爲你有一種特殊形式的多項式,我們可以通過將每個係數乘以冪並將所有項移向下來來容易地計算導數。這裏是一個輔助函數來calcualte新的係數

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") 

enter image description here

我們可以看到拐點對應於導數爲零的地方。

+0

謝謝!這就是我一直在尋找的! – Oliver

4

您可以通過首先對數據進行相對於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() 
+0

謝謝大衛!這是可行的,並且是解決手頭問題的實用解決方案!而且,據我瞭解,MrFlick的解決方案在數學上是正確的。 – Oliver

1

據稱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

注意:所有這一切都說,它並不適用於我和我的數據。

+0

在我的情況下,它的工作原理 – Ben