2016-01-14 26 views
3

我在2D圖上有點。我想找到適合這個模型的最好的第三個多項式並且得到它的一階導數。但我無法獲得D函數的工作。下面是簡單的例子:使用`D`來獲得我的迴歸多項式的導數時出錯

a <- 0:10 
b <- c(2, 4, 5, 8, 9, 12, 15, 16, 18, 19, 20) 
plot(a, b) 
m1 <- lm(b ~ a + I(a^2) + I(a^3)) 
s <- coef(m1) 

## try to get 1st derivative of the regression polynomial 
D(expression(s[1] + s[2] * a + (a^2) * s[3] + (a^3) * s[4]), "a") 

錯誤d(式(S [1] + S [2] * A +(A^2)* S [3] +(A^3)* S [4]):

功能「[」是不是在衍生品表

我想避免由差分計算數值導感謝您的幫助

+0

工作正常,thx,但是我怎樣才能繪製派生? – Bobesh

回答

2

你看到錯誤消息「! 函數'['不在派生表中「是因爲D只能識別符號操作的某一組函數。您可以在?D找到他們:

The internal code knows about the arithmetic operators ‘+’, ‘-’, 
‘*’, ‘/’ and ‘^’, and the single-variable functions ‘exp’, ‘log’, 
‘sin’, ‘cos’, ‘tan’, ‘sinh’, ‘cosh’, ‘sqrt’, ‘pnorm’, ‘dnorm’, 
‘asin’, ‘acos’, ‘atan’, ‘gamma’, ‘lgamma’, ‘digamma’ and 
‘trigamma’, as well as ‘psigamma’ for one or two arguments (but 
derivative only with respect to the first). (Note that only the 
standard normal distribution is considered.) 

雖然"["實際上是R A功能(讀?Extract?"[")。

爲了證明類似的行爲,可以考慮:

s <- function (x) x 

D(expression(s(x) + x^2), name = "x") 
# Error in D(expression(s(x) + x^2), name = "x") : 
# Function 's' is not in the derivatives table 

在這裏,即使s已被定義爲一個簡單的功能,D可以做到與它無關。

您的問題已通過我最近的回答Function for derivatives of polynomials of arbitrary order (symbolic method preferred)解決。在我的三個答案中提供了三種方法,其中沒有一種方法基於數​​值導數。我個人更喜歡the one using outer(LaTeX數學公式的唯一答案),因爲對於多項式,一切都是確切的。

要使用該解決方案,使用功能g那裏,並要評估導數值指定參數x(比如0:10),並pc你多項式迴歸係數s。默認情況下,nderiv = 0L因此多項式本身返回,就好像調用了predict.lm(m1, newdata = list(a = 0:10))一樣。但是一旦指定了nderiv,就可以得到迴歸曲線的精確導數。

a <- 0:10 
b <- c(2, 4, 5, 8, 9, 12, 15, 16, 18, 19, 20) 
plot(a, b) 
m1 <- lm(b ~ a + I(a^2) + I(a^3)) 
s <- coef(m1) 
#(Intercept)   a  I(a^2)  I(a^3) 
# 2.16083916 1.17055167 0.26223776 -0.02020202 

## first derivative at your data points 
g(0:10, s, nderiv = 1) 
# [1] 1.1705517 1.6344211 1.9770785 2.1985237 2.2987568 2.2777778 2.1355866 
# [8] 1.8721834 1.4875680 0.9817405 0.3547009 

其他備註:建議您使用poly(a, degree = 3, raw = TRUE)而非I()。他們在這裏也是這樣做的,但是poly更簡潔,並且如果你想要交互就更容易了,就像How to write interactions in regressions in R?