【问题标题】:do a nonlinear least square fit in r在 r 中进行非线性最小二乘拟合
【发布时间】:2015-08-05 10:13:43
【问题描述】:

我有两个向量:

y <- c(0.044924, 0.00564, 0.003848, 0.002385, 0.001448, 0.001138, 
0.001025, 0.000983, 0.00079, 0.000765, 0.000721, 0.00061, 0.000606, 
0.000699, 0.000883, 0.001069, 0.001226, 0.001433, 0.00162, 0.001685, 
0.001604, 0.001674, 0.001706, 0.001683, 0.001505, 0.001497, 0.001416, 
0.001449, 0.001494, 0.001544, 0.00142, 0.001458, 0.001544, 0.001279, 
0.00159, 0.001756, 0.001749, 0.001909, 0.001885, 0.002063, 0.002265, 
0.002137, 0.002391, 0.002619, 0.002733, 0.002957, 0.003244, 0.003407, 
0.003563, 0.003889, 0.004312, 0.004459, 0.004946, 0.005248, 0.005302, 
0.00574, 0.006141, 0.006977, 0.007386, 0.007843, 0.008473, 0.008949, 
0.010164, 0.010625, 0.011279, 0.01191, 0.012762, 0.014539, 0.01477)

x <- 0:68

我正在尝试使用非线性最小二乘函数来拟合数据,但我不断收到错误:

nlsModel(formula, mf, start, wts) 中的错误:初始参数估计处的奇异梯度矩阵

我的代码是:

a=0.00012
b=0.08436
k=0.21108
fit = nls(y ~ (a*b*exp(b*x)*k)/((k*b)+(a*(exp(b*x)-1))), start=list(a=a,b=b,k=k))

我输入的参数是我知道接近预期值的参数。有谁知道我在这里做错了什么?

我尝试了参数abk 的各种初始值,但总是出现某种错误。

【问题讨论】:

  • 当我将您的起始曲线与点进行比较时,例如使用plot(x, y); lines((a*b*exp(b*x)*k)/((k*b)+(a*(exp(b*x)-1))), col = "red"),它看起来并不那么相似——特别是没有表示接近 0 的凹凸。为什么你确信数据可以被这条曲线拟合,并且参数值很接近?
  • 我和@DavidRobinson 在一起。看起来您使用的方程无法拟合上翘的抛物线。我使用了具有一系列起始值的方程,并且 x 小于 10 和 x 大于 60 的唯一曲线看起来很奇怪,其中有锐角。
  • 非常感谢@DavidRobinson 和 Jean V. Adams。我同意你们俩。在这种情况下,y=死亡率,x=年龄,模型使用的是混合分布(脆弱模型)。 (a*b*exp(b*x)*k)/((k*b)+(a*(exp(b*x)-1))) 它的伽马分布与 Gompertz 是基线危险。我怎样才能找到合适的呢?需要你的帮助。非常感谢

标签: r nls gamma-distribution


【解决方案1】:

请改用optim()。您必须创建一个以 a、b 和 k 作为输入(作为向量收集)并返回平方误差作为结果的函数:

func <-function(pars) {
  a <- pars["a"]
  b <- pars["b"]
  k <- pars["k"]
  fitted <- (a*b*exp(b*x)*k)/((k*b)+(a*(exp(b*x)-1)))
  sum((y-fitted)^2)  
  } 

然后我们使用初始值运行 optim():

result <- optim(c(a=0.00012, b=0.08436, k=0.21108), func)

测试结果拟合:

plot(x, y)
a <- result$par["a"]
b <- result$par["b"]
k <- result$par["k"]
lines((a*b*exp(b*x)*k)/((k*b)+(a*(exp(b*x)-1))), col = "blue")

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 2012-06-29
    • 1970-01-01
    • 2017-02-24
    • 1970-01-01
    • 2012-08-29
    • 1970-01-01
    • 2018-12-27
    • 1970-01-01
    相关资源
    最近更新 更多