2017-05-30 123 views
1

我是編程和使用R軟件的新手,所以我非常感謝您對我正在嘗試解決的當前問題的反饋。用不同的累積分佈擬合實驗數據點使用R

因此,我必須擬合一些函數(兩/三參數函數)的累積分佈。這似乎是非常直接的任務,但我現在一直在嗡嗡聲一段時間。

讓我告訴你什麼是我的變量:

x=c(0.01,0.011482,0.013183,0.015136,0.017378,0.019953,0.022909,0.026303,0.0302,0.034674,0.039811,0.045709,0.052481,0.060256,0.069183,0.079433,0.091201,0.104713,0.120226,0.138038,0.158489,0.18197,0.20893,0.239883,0.275423,0.316228,0.363078,0.416869,0.47863,0.549541,0.630957,0.724436,0.831764,0.954993,1.096478,1.258925,1.44544,1.659587,1.905461,2.187762,2.511886,2.884031,3.311311,3.801894,4.365158,5.011872,5.754399,6.606934,7.585776,8.709636,10,11.481536,13.182567,15.135612,17.378008,19.952623,22.908677,26.30268,30.199517,34.673685,39.810717,45.708819,52.480746,60.255959,69.183097,79.432823,91.201084,104.712855,120.226443,138.038426,158.489319,181.970086,208.929613,239.883292,275.42287,316.227766,363.078055,416.869383,478.630092,549.540874,630.957344,724.43596,831.763771,954.992586,1096.478196) 
    y=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00044816,0.00127554,0.00221488,0.00324858,0.00438312,0.00559138,0.00686054,0.00817179,0.00950625,0.01085188,0.0122145,0.01362578,0.01514366,0.01684314,0.01880564,0.02109756,0.0237676,0.02683182,0.03030649,0.0342276,0.03874555,0.04418374,0.05119304,0.06076553,0.07437854,0.09380666,0.12115065,0.15836926,0.20712933,0.26822017,0.34131335,0.42465413,0.51503564,0.60810697,0.69886817,0.78237651,0.85461023,0.91287236,0.95616228,0.98569093,0.99869001,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999) 

這是我設置x軸的對劇情: Raw data

經過一番研究,我試圖與乙狀結腸功能,如在其中一篇帖子中找到的(由於我的聲望不夠高,我無法添加鏈接)。這是代碼:

# sigmoid function definition 
sigmoid = function(params, x) { 
    params[1]/(1 + exp(-params[2] * (x - params[3]))) 
} 

# fitting code using nonlinear least square 
fitmodel <- nls(y~a/(1 + exp(-b * (x-c))), start=list(a=1,b=.5,c=25)) 

# get the coefficients using the coef function 
params=coef(fitmodel) 

# asigning to y2 sigmoid function 
y2 <- sigmoid(params,x) 

# plotting y2 function 
plot(y2,type="l") 

# plotting data points 
points(y) 

這使我得到了一些很好的擬合結果(我不知道如何量化這個)。但是,當我觀察Sigmuid擬合函數的圖時,我不明白爲什麼S形現在發生在x值從40到7的範圍內(查看S形應該是從x值10至200)。

Raw data

因爲我無法解釋這種現象,我想嘗試的擬合方程韋伯,但到目前爲止,我不能讓代碼運行。

綜上所述:

  1. 你有什麼想法,爲什麼是乙狀結腸給我說,奇怪的配件?
  2. 你知道這個擬合方法有更好的兩個或三個參數方程嗎?
  3. 我怎麼能確定合適的善良?像r^2那樣的東西?
+2

它繪製數組的索引,因爲你不提供的x值。嘗試'plot(x,y2,type =「l」)'和'points(x,y)'。 – Lyngbakr

+0

@Lyngbakr謝謝。這解決了我的第一個問題。爲了更好地看到這個S曲率,我輸入了plot(x,y,type =「l」,log =「x」)''。但這只是證實了這個合適的看起來不太好。 – numb

+0

我最初的想法是嘗試'a + b * tanh(x/c)',但是這也給了蹩腳的結果... – Lyngbakr

回答

0

通過不同的功能和不同的數據集,我發現了最好的解決方案,給出了我發佈的所有問題的答案。

的代碼,因爲它遵循有關規定的數據集:

df <- data.frame(x=c(0.01,0.011482,0.013183,0.015136,0.017378,0.019953,0.022909,0.026303,0.0302,0.034674,0.039811,0.045709,0.052481,0.060256,0.069183,0.079433,0.091201,0.104713,0.120226,0.138038,0.158489,0.18197,0.20893,0.239883,0.275423,0.316228,0.363078,0.416869,0.47863,0.549541,0.630957,0.724436,0.831764,0.954993,1.096478,1.258925,1.44544,1.659587,1.905461,2.187762,2.511886,2.884031,3.311311,3.801894,4.365158,5.011872,5.754399,6.606934,7.585776,8.709636,10,11.481536,13.182567,15.135612,17.378008,19.952623,22.908677,26.30268,30.199517,34.673685,39.810717,45.708819,52.480746,60.255959,69.183097,79.432823,91.201084,104.712855,120.226443,138.038426,158.489319,181.970086,208.929613,239.883292,275.42287,316.227766,363.078055,416.869383,478.630092,549.540874,630.957344,724.43596,831.763771,954.992586,1096.478196), 
     y=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00044816,0.00127554,0.00221488,0.00324858,0.00438312,0.00559138,0.00686054,0.00817179,0.00950625,0.01085188,0.0122145,0.01362578,0.01514366,0.01684314,0.01880564,0.02109756,0.0237676,0.02683182,0.03030649,0.0342276,0.03874555,0.04418374,0.05119304,0.06076553,0.07437854,0.09380666,0.12115065,0.15836926,0.20712933,0.26822017,0.34131335,0.42465413,0.51503564,0.60810697,0.69886817,0.78237651,0.85461023,0.91287236,0.95616228,0.98569093,0.99869001,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999)) 


library(drc) 
fm <- drm(y ~ x, data = df, fct = G.3()) #The Gompertz model G.3() 
plot(fm) 

#Gompertz Coefficients and residual standard error 

summary(fm) 

情節after fitting

1
# Data 
df <- data.frame(x=c(0.01,0.011482,0.013183,0.015136,0.017378,0.019953,0.022909,0.026303,0.0302,0.034674,0.039811,0.045709,0.052481,0.060256,0.069183,0.079433,0.091201,0.104713,0.120226,0.138038,0.158489,0.18197,0.20893,0.239883,0.275423,0.316228,0.363078,0.416869,0.47863,0.549541,0.630957,0.724436,0.831764,0.954993,1.096478,1.258925,1.44544,1.659587,1.905461,2.187762,2.511886,2.884031,3.311311,3.801894,4.365158,5.011872,5.754399,6.606934,7.585776,8.709636,10,11.481536,13.182567,15.135612,17.378008,19.952623,22.908677,26.30268,30.199517,34.673685,39.810717,45.708819,52.480746,60.255959,69.183097,79.432823,91.201084,104.712855,120.226443,138.038426,158.489319,181.970086,208.929613,239.883292,275.42287,316.227766,363.078055,416.869383,478.630092,549.540874,630.957344,724.43596,831.763771,954.992586,1096.478196), 
      y=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00044816,0.00127554,0.00221488,0.00324858,0.00438312,0.00559138,0.00686054,0.00817179,0.00950625,0.01085188,0.0122145,0.01362578,0.01514366,0.01684314,0.01880564,0.02109756,0.0237676,0.02683182,0.03030649,0.0342276,0.03874555,0.04418374,0.05119304,0.06076553,0.07437854,0.09380666,0.12115065,0.15836926,0.20712933,0.26822017,0.34131335,0.42465413,0.51503564,0.60810697,0.69886817,0.78237651,0.85461023,0.91287236,0.95616228,0.98569093,0.99869001,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999)) 

# sigmoid function definition 
sigmoid = function(x, a, b, c) { 
    a * exp(-b * exp(-c * x)) 
} 

# fitting code using nonlinear least square 
fitmodel <- nls(y ~ sigmoid(x, a, b, c), start=list(a=1,b=.5,c=-2), data = df) 

# plotting y2 function 
plot(df$x, predict(fitmodel),type="l", log = "x") 

# plotting data points 
points(df) 

enter image description here

我使用的功能是Gompertz functionthis blog post解釋了爲什麼R 2不應與非線性擬合,使用,並提供了一種替代。

+0

謝謝。我沒有意識到Gompertz函數的存在。它甚至適用於我的其他數據集。關於R^2,我還在某處推薦使用它來進行非線性迴歸,並且您用迴歸的標準誤差發送的這個鏈接似乎是最好的方法。 – numb

+0

順便說一句,還有[R包](https://cran.r-project.org/web/packages/sigmoid/index.html)專門用於sigmoid函數。 – Lyngbakr

+0

太棒了!這個軟件包有一些我在[growthmodels]找不到的模型(https://cran.r-project.org/web/packages/growthmodels/growthmodels.pdf)。你真的幫了我很多! – numb