在R中使用NLS拟合渐近线

pdkcd3nj  于 2023-09-27  发布在  其他
关注(0)|答案(1)|浏览(121)

我正试图用nls来拟合一个模型

data <- structure(list(x = c(0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 
                     0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 
                     0.19, 0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 
                     0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 
                     0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 
                     0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61, 0.62, 
                     0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, 0.72, 0.73, 
                     0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83, 0.84, 
                     0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 
                     0.96, 0.97, 0.98, 0.99, 1), 
               y = c(0, 0.485920598224734, 0.579084967320261, 
                                                       0.629311870308531, 0.66413220309958, 0.688711933099295, 0.708245568969946, 
                                                       0.724788081171333, 0.738425093472615, 0.749807346519394, 0.760647315694837, 
                                                       0.769951765276708, 0.778616890715529, 0.786471444472986, 0.793875046379542, 
                                                       0.800645032394326, 0.806792819019893, 0.812866398378857, 0.818266404087108, 
                                                       0.823723492308131, 0.828706795673146, 0.833244855438536, 0.837942746239689, 
                                                       0.842309558466764, 0.846539372663185, 0.850717812598111, 0.854388218169364, 
                                                       0.85805291548934, 0.862094357393613, 0.865770471216143, 0.86910408996204, 
                                                       0.872403459200274, 0.875571538659132, 0.878585495333505, 0.881485286982333, 
                                                       0.884584867425864, 0.887478950823415, 0.890293118703085, 0.893267117618518, 
                                                       0.895978536975198, 0.89874703884465, 0.901190170391301, 0.903810257727545, 
                                                       0.906407512058681, 0.908776436338728, 0.911179610126438, 0.913645574678197, 
                                                       0.916151496988897, 0.918457630504895, 0.920558266974912, 0.922852983988355, 
                                                       0.924925079201986, 0.926900134143905, 0.928966521106259, 0.930821702771356, 
                                                       0.932751091703057, 0.934743271398807, 0.936695493335617, 0.938704797785199, 
                                                       0.940520021691355, 0.942449410623056, 0.944150469503667, 0.945914319148329, 
                                                       0.947643919285327, 0.949436310186374, 0.951114536061877, 0.952821303193767, 
                                                       0.954670776607586, 0.956343294231812, 0.957907355081771, 0.959505665439393, 
                                                       0.961024060279134, 0.962559579872706, 0.964140765476496, 0.965573536547079, 
                                                       0.966994891115107, 0.968530410708679, 0.970071638553529, 0.971550075634329, 
                                                       0.972954305448526, 0.974472700288267, 0.975876930102463, 0.977189827896224, 
                                                       0.978691097982133, 0.980049661786112, 0.981516682364358, 0.982886662670891, 
                                                       0.984171019208266, 0.985409709735423, 0.986694066272797, 0.988006964066558, 
                                                       0.989222821588606, 0.990358763592773, 0.99155749636099, 0.992716271370266, 
                                                       0.993983503153809, 0.995324942203956, 0.996375260438965, 0.997591117961013, 
                                                       0.998909724006051, 1)), row.names = c(NA, -101L), class = c("tbl_df", 
                                                                                                                   "tbl", "data.frame"))

The data looks like this
我试图拟合一个渐近模型,由于数据的性质,我们知道它应该收敛于一个点,它收敛的值是有用的信息。
我尝试使用ax^B+c拟合模型,但似乎不起作用。

> model =  nls(y~a*x^(b),
+              data = data,
+              start = list(a=-10, b=-0.2))
Error in numericDeriv(form[[3L]], names(ind), env, central = nDcentral) : 
  Missing value or an infinity produced when evaluating the model

从我收集到的信息来看,这些模型与数据不符。这是真的吗?我在不同的时间段从另一个来源生成数据,其他时间段生成的数据我可以很好地建模,数据遵循类似的形状,例如,

data <- structure(list(x = c(0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 
                     0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 
                     0.19, 0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 
                     0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 
                     0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 
                     0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61, 0.62, 
                     0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, 0.72, 0.73, 
                     0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83, 0.84, 
                     0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 
                     0.96, 0.97, 0.98, 0.99, 1), y = c(0.387696586431811, 0.50631394232361, 
                                                       0.573625136575837, 0.620130450617488, 0.653935039565606, 0.680528424328709, 
                                                       0.702201768036288, 0.720147005264378, 0.73551633943648, 0.748852762970566, 
                                                       0.760487368804423, 0.771175048836208, 0.780975399794722, 0.78935205112075, 
                                                       0.796907591961063, 0.804330695626262, 0.811118100850909, 0.81746184153892, 
                                                       0.823328808396517, 0.828639539118631, 0.833526470880376, 0.838598814687283, 
                                                       0.843386418567692, 0.848240241035659, 0.852398768334271, 0.856808926265603, 
                                                       0.860755554084031, 0.864616097738635, 0.868152170314207, 0.871648511737245, 
                                                       0.875105122007748, 0.878171042611661, 0.881203853921796, 0.884395589842069, 
                                                       0.887593947621097, 0.890540674767407, 0.893235771280999, 0.896036817534682, 
                                                       0.898765023342052, 0.901360791974307, 0.90402277919412, 0.906492732510015, 
                                                       0.909022282554713, 0.911326689401715, 0.913558255802404, 0.915922259378207, 
                                                       0.917922060722445, 0.92012713968811, 0.922193159619905, 0.924504188325663, 
                                                       0.926484124093633, 0.92848392543787, 0.930477104923352, 0.932159057047313, 
                                                       0.934238320696619, 0.936125550442009, 0.938045889481177, 0.939840413203986, 
                                                       0.941654802503063, 0.943429460649604, 0.945190875078635, 0.947071482965268, 
                                                       0.948746813230474, 0.950534715094527, 0.952263020229778, 0.95381915703738, 
                                                       0.955467999867563, 0.957037380392676, 0.958514054895209, 0.960143032149124, 
                                                       0.961705790815482, 0.963202330894282, 0.964698870973082, 0.966195411051882, 
                                                       0.967658841836904, 0.969162003774459, 0.970506241101877, 0.972049134191968, 
                                                       0.973373505943118, 0.974658146541734, 0.976075224315465, 0.977406217925372, 
                                                       0.978863026851637, 0.980081448862696, 0.981425686190114, 0.982551402178591, 
                                                       0.983955236234811, 0.985067708505778, 0.986425189550707, 0.987656855279277, 
                                                       0.98886865543158, 0.990080455583882, 0.991530642651392, 0.992802039532497, 
                                                       0.994033705261067, 0.995066715226964, 0.996291759096778, 0.997523424825349, 
                                                       0.998702115683872, 1)), row.names = c(NA, -100L), class = c("tbl_df", 
                                                                                                                   "tbl", "data.frame"))

model2 =  nls(y~a*x^(b)+c,
+              data = data,
+              start = list(a=-10, b=-0.2, c=1.5))

这个很好用。
我需要一个能同时装下这两个的模型。我是否使用了错误的公式,或者是否有比NLS更好的替代方案?谢谢

uwopmtnx

uwopmtnx1#

如果使用更好的起始值,该模型似乎可以很好地处理第一个数据集。对于这个模型,你可以通过做一个对数-对数回归来得到合适的初始值:

m1 <- lm(log(y) ~ log(x), data1, subset = x > 0)
model =  nls(y~a*x^(b),
             data = data1,
             start = list(a=exp(coef(m1)[1]), b=coef(m1)[2]))

“更好的起始值”通常是这类问题的答案。它还可以帮助为nls提供函数梯度,您可以按如下方式执行(在这种情况下,我们必须再次删除x=0数据点,因为导数是NaN...)

power_law <- deriv(~a*x^b, c("a", "b"),
                   function.arg = c("a", "b", "x"))

model =  nls(y~power_law(a, b, x),
             data = subset(data1, x>0),
             start = list(a=exp(coef(m1)[1]), b=coef(m1)[2]))

正如我确信我在其他地方评论过的那样,您也可以通过

model <- glm(y~log(x), family = gaussian(link="log"), data = data1)

(尽管截距项将是log(a)而不是a

相关问题