1

我有我想使用 R 拟合以下等式的数据：

Z(u,w)=z0*F(w)*[1-exp((-b*u)/F(w))]

其中 z0 和 b 是常数，F(w), w=0,...,9 是一个递减阶跃函数，它取决于 w 且 F(0)=1 和 u=1,...,50。

Z(u,w) 是以 50x10 矩阵的形式观察到的一组数据（u=50,...,1 在行的一侧，w=0,...,9 沿列）。例如，由于我没有很好地解释，Z(42,3) 将是第 9 行向下和第 4 列中的元素。

使用 F(0)=1 我能够仅使用第一列（即 w=0）和代码来获得 b 和 z0 的估计值：

n0=nls(zuw~z0*(1-exp(-b*u)),start=list(z0=283,b=0.03),options(digits=10))

然后，我通过遍历每一列并使用我找到的 b 和 z0 的值，找到了 w=1,...,9 的 F(w)。

但是，我想找到一种方法来一次估计所有 12 个参数（b、z0 和 F(w) 的 10 个值），因为 b 和 z0 应该适合所有数据，而不仅仅是第一列。

有谁知道这样做的任何方法？所有帮助将不胜感激！

谢谢詹姆斯

4

1 回答 1

0

这可能是nls(...)函数的公式接口对您不利的情况。作为替代方案，您可以nls.lm(...)在minpack.lm包中使用以编程方式定义的函数执行非线性回归。为了证明这一点，首先我们创建了一个人工数据集，该数据集按照您的功能形式设计，并添加了随机误差（误差 ~ N[0,1]）。

u  <- 1:50
w  <- 0:9
z0 <- 100
b  <- 0.02
F  <- 10/(10+w^2)
# matrix containing data, in OP's format: rows are u, cols are w
m  <- do.call(cbind,lapply(w,function(w)
                       z0*F[w+1]*(1-exp(-b*u/F[w+1]))+rnorm(length(u),0,1)))

所以现在我们有一个矩阵m，它相当于你的数据集。该矩阵采用所谓的“宽”格式 - 不同值的响应w位于不同列中。我们需要“长”格式：所有响应都在一个列中，并有一个单独的列标识u和w。我们melt(...)在reshape2包中使用。

# prepend values of u
df.wide <- data.frame(u=u, m)
library(reshape2)
# reshape to long format: col1 = u, col2=w, col3=z
df   <- melt(df.wide,id="u",variable.name="w", value.name="z")
df$w <- as.numeric(substr(df$w,2,4))-1

现在我们有一个包含、和df列的数据框。该函数采用（至少）4 个参数：是拟合参数的初始估计向量，是在每一步计算残差的函数，是因变量 ( )，是包含自变量的向量或矩阵变量（）。uwznls.lm(...)parfnobservedzxxu, v

接下来我们定义一个函数 ，f(par, xx)其中par是一个11 元素向量。前两个元素包含z0和的估计值b。接下来的 9 个包含 的估计值F(w), w=1:9。这是因为您声明F(0)已知为 1.xx是一个包含两列的矩阵：分别为u和的值w。f(par,xx)然后计算给定参数估计值和z的所有值的响应估计值。uw

library(minpack.lm)
# model function
f <- function(pars, xx) {
  z0 <- pars[1]
  b  <- pars[2]
  F  <- c(1,pars[3:11])
  u  <- xx[,1]
  w  <- xx[,2]
  z  <- z0*F[w+1]*(1-exp(-b*u/F[w+1]))
  return(z)
}
# residual function
resids <- function(p, observed, xx) {observed - f(p,xx)}

接下来，我们使用 执行回归nls.lm(...)，它使用高度鲁棒的拟合算法 (Levenberg-Marquardt)。因此，我们可以将par参数（包含 的初始估计z0, b, and F）设置为 all 1's，这与创建数据集时使用的值（“实际”值）相距甚远。nls.lm(...)返回一个包含多个组件的列表（请参阅文档）。该par组件包含拟合参数的最终估计值。

# initial parameter estimates; all 1's
par.start <- c(z0=1, b=1, rep(1,9))   
# fit using Levenberg-Marquardt algorithm
nls.out <- nls.lm(par=par.start, 
                  fn = resids, observed = df$z, xx = df[,c("u","w")], 
                  control=nls.lm.control(maxiter=10000, ftol=1e-6, maxfev=1e6))
par.final <- nls.out$par
results <- rbind(predicted=c(par.final[1:2],1,par.final[3:11]),actual=c(z0,b,F))
print(results,digits=5)
#               z0        b                                                                          
# predicted 102.71 0.019337 1 0.90456 0.70788 0.51893 0.37804 0.27789 0.21204 0.16199 0.13131 0.10657
# actual    100.00 0.020000 1 0.90909 0.71429 0.52632 0.38462 0.28571 0.21739 0.16949 0.13514 0.10989

所以回归在恢复“实际”参数值方面做得很好。最后，我们绘制结果ggplot只是为了确保这一切都是正确的。我不能过分强调绘制最终结果的重要性。

df$pred <- f(par.final,df[,c("u","w")])
library(ggplot2)
ggplot(df,aes(x=u, color=factor(w)))+
  geom_point(aes(y=z))+ geom_line(aes(y=pred))

于 2014-03-22T16:53:45.033 回答