2014-02-12 22 views
0

我想自定義公式,以適應數據和我收到警告味精並沒有得到執行它..這是我的R-代碼..錯誤在執行定製適合與NLS功能

x <- c(3, 33, 146, 227, 342, 351, 353, 444, 556, 571, 709, 759, 836, 
860, 968, 1056, 1726, 1846, 1872, 1986, 2311, 2366, 2608, 2676, 
3098, 3278, 3288, 4434, 5034, 5049, 5085, 5089, 5089, 5097, 5324, 
5389, 5565, 5623, 6080, 6380, 6477, 6740, 7192, 7447, 7644, 7837, 
7843, 7922, 8738, 10089, 10237, 10258, 10491, 10625, 10982, 11175, 
11411, 11442, 11811, 12559, 12559, 12791, 13121, 13486, 14708, 
15251, 15261, 15277, 15806, 16185, 16229, 16358, 17168, 17458, 
17758, 18287, 18568, 18728, 19556, 20567, 21012, 21308, 23063, 
24127, 25910, 26770, 27753, 28460, 28493, 29361, 30085, 32408, 
35338, 36799, 37642, 37654, 37915, 39715, 40580, 42015, 42045, 
42188, 42296, 42296, 45406, 46653, 47596, 48296, 49171, 49416, 
50145, 52042, 52489, 52875, 53321, 53443, 54433, 55381, 56463, 
56485, 56560, 57042, 62551, 62651, 62661, 63732, 64103, 64893, 
71043, 74364, 75409, 76057, 81542, 82702, 84566, 88682) 

y <- c(1:136) 
df <- data.frame(x,y) 

fit <- nls(y ~ a*(1-exp(-x/b))^c, data=df, start = list(a=100,b=1000,c=0.5), 
    algorithm="port",lower=list(a=100,b=100,c=0.5),upper=list(a=200,b=10000,c=2)) 

警告信息:

1: In min(x) : no non-missing arguments to min; returning Inf 
2: In max(x) : no non-missing arguments to max; returning -Inf 

如何將這個數據代入這個公式的自定義,以及如何找到R平方值..

在此先感謝..

+1

您的代碼適用於我,沒有任何警告。 – thelatemail

+0

我可以改變算法,以「線性」.. – Chammu

+0

如果你相信自己,你可以做任何事情。 – thelatemail

回答

1

首先,要直接回答您的問題,對於非線性擬合,R 不是一個有意義的概念。對於線性模型,R 是由該模型解釋的數據集中總變異性的分數。此計算僅在SST = SSR + SSE時有效,對所有線性模型均如此,但對於非線性模型不一定適用。請參閱this question以獲取更完整的說明和一些其他參考。

所以,雖然可以用於線性模型訪問R 作爲

rsq <- summary(lm(...))$r.squared 

用於非線性模型的summary(...)方法不返回的R 值。

其次,您確實需要養成根據您的數據繪製擬合曲線的習慣。

plot(x,y) 
lines(x,predict(fit)) 

只是有沒有辦法來解釋這是一個 「好」 的選擇。如果我們重新運行不會對A,B的約束模型,以及c我們得到了一個更好的契合:

fit.2 <- nls(y ~ a*(1-exp(-x/b))^c, data=df, start = list(a=100,b=1000,c=0.5), 
      algorithm="port") 

par(mfrow=c(1,2)) 
plot(x,y, main="Constrained Model",cex=0.5) 
lines(x,predict(fit), col="red") 
plot(x,y, main="UNconstrained Model",cex=0.5) 
lines(x,predict(fit.2), col="red") 

顯然,這種模式是「更好」,但是,這並不意味着這很好」。除其他外,我們需要查看殘差。對於一個合適的模型,殘差不應該依賴於x。因此,讓我們來看看:

plot(y,residuals(fit),main="Residuals: 1st Model") 
plot(y,residuals(fit.2),main="Residuals: 2ndd Model") 

在第二個模型的殘差小得多,不遵循一個趨勢,雖然有似乎是一些底層結構。這是需要考慮的事情 - 這意味着可能由於實際效應或可能由於數據收集方法而導致一些低振幅振盪。另外,迴歸建模的基本原理(開始假設)是殘差正態分佈且方差不變(例如,方差不依賴於x)。我們可以用一個Q-Q圖來檢查這個圖,它繪製了N(0,1)分位數的(標準化)殘差的分位數。如果殘差是正態分佈的,這應該是一條直線。

se <- summary(fit)$sigma 
qqnorm(residuals(fit)/se, main="Q-Q Plot, 1st Model") 
qqline(residuals(fit)/se,probs=c(0.25,0.75)) 
se.2 <- summary(fit.2)$sigma 
qqnorm(residuals(fit.2)/se.2, main="Q-Q Plot, 2nd Model") 
qqline(residuals(fit.2)/se.2,probs=c(0.25,0.75)) 

正如你所看到的,從第2模型的殘差非常接近正態分佈。海事組織這一點,比其他任何事情都更顯示第二種模式是一種「好」模式。

最後,您的數據在我看來就像是具有2個峯值的分佈的cdf。一個非常簡單的方法來類似數據模型爲:

Y〜一個 ×(1 - EXP(-x/B ))+一個 ×(1 - EXP( -x/b ))

當我嘗試這種模式,它比你的模型的約束版本更好一點。

+0

最小平方我得到..如何mle(最大似然估計)方法可以應用於此** y〜a *(1-exp(-x/b))^ c **方程。基本上我的項目是比較通過對Musa數據集的這個方程的最小二乘法和mle方法輸出兩個估計參數(a,b,c)。 – Chammu

+0

該方法告訴天氣哪種方法估計得好。最小二乘法或最大似然估計 – Chammu