R

擬合平滑樣條方程及其解析導數

  • April 22, 2015

我需要將樣條函數擬合到數據集。我試過bs,nssmooth.spline。在我的例子中,通過 獲得的曲線比和smooth.spline更好地遵循數據趨勢。bs``ns

但是,我不知道如何

  1. 獲得擬合smooth.spline函數的方程?

(即如何解釋smooth.spline 輸出的方程係數) 2. 計算的擬合函數? 3. 計算擬合函數的解析一階導數,並獲得其方程(帶有實際係數)?

數據集的一個例子可以計算為

library(stats)
x <- 1:11
y <- c(0.2,0.40, 0.6, 0.75, 0.88, 0.99, 1.1, 1.15, 1.16, 1.16, 1.16 )
plot(x, y)
spline <- smooth.spline(x,y)
lines(spline)

這適合使用截斷功率基礎的自然樣條曲線(線性尾部受限)。在此示例中,未使用默認節點(基於預測變量的分位數);相反,我們指定 4 個結。獲得擬合優度檢驗的唯一方法是假設一個比這更豐富的模型,看看它是否改進了下面擬合的模型。但是anova()通過將非線性項匯集到復合(“塊”)測試(F = 175.38)中來測試線性關係的擬合優度。

require(rms)
x <- 1:11
y <- c(0.2,0.40, 0.6, 0.75, 0.88, 0.99, 1.1, 1.15, 1.16, 1.16, 1.16 )
dd <- datadist(x); options(datadist='dd')

f <- ols(y ~ rcs(x, c(3, 5, 7, 9)))
f

Linear Regression Model

ols(formula = y ~ rcs(x, c(3, 5, 7, 9)))

               Model Likelihood     Discrimination    
                  Ratio Test           Indexes        
Obs       11    LR chi2     66.08    R2       0.998    
sigma 0.0201    d.f.            3    R2 adj   0.996    
d.f.       7    Pr(> chi2) 0.0000    g        0.383    

Residuals

     Min        1Q    Median        3Q       Max 
-0.027360 -0.011739  0.001227  0.009892  0.031166 

         Coef    S.E.   t     Pr(>|t|)
Intercept  0.0465 0.0224  2.08 0.0762  
x          0.1741 0.0072 24.18 <0.0001 
x' -0.1004 0.0311 -3.23 0.0144 
x'' 0.0542 0.0913 0.59 0.5715 

anova(f)

Analysis of Variance Response: y 

Factor d.f. Partial SS MS F P 
x 3 1.152321844 0.3841072814 946.15 <.0001
Nonlinear 2 0.142398208 0.0711991039 175.38 <.0001
REGRESSION 3 1.152321844 0.3841072814 946.15 <.0001
ERROR 7 0.002841792 0.0004059703 

ggplot(Predict(f)) + geom_point(aes(x=x, y=y), data=data.frame(x,y))
Function(f) ## if have latex installed can also use latex(f)

function(x = 6) {0.046475489+0.17411942* x-0.002790266*pmax(x- 3,0)^3+0.0015048699*pmax(x-5,0)^3+0.0053610582*pmax(x-7,0)^3-0.0040756621*pmax(x-9,0)^3 }

函數以最簡單的形式重新表達受限三次樣條。一階導數是:

function(x) 0.174 - 3 * 0.00279 * pmax(x - 3, 0) ^ 2 + 3 * 0.0015 * pmax(x - 5, 0) ^ 2 + ...

引用自:https://stats.stackexchange.com/questions/147733

comments powered by Disqus