R非线性最小二乘法（nls）模型拟合

Question

我试图从我的数据的G函数拟合到以下数学模式：y = A /（（1 +（B^2）*（x^2））^（（C + 1）/ 2））。该图的形状可以在这里看到：

http://www.wolframalpha.com/input/?i=y+%3D+1%2F+%28%281+%2B+%282%5E2%29*%28x%5E2%29%29%5E%28%282%2B1%29%2F2%29%29

这里有一个基本的例子就是我一直在做：

data(simdat) 

library(spatstat) 

simdat.Gest <- Gest(simdat) #Gest is a function within spatstat (explained below) 

Gvalues <- simdat.Gest$rs 

Rvalues <- simdat.Gest$r 

GvsR_dataframe <- data.frame(R = Rvalues, G = rev(Gvalues)) 

themodel <- nls(rev(Gvalues) ~ (1/(1 + (B^2)*(R^2))^((C+1)/2)), data = GvsR_dataframe, start = list(B=0.1, C=0.1), trace = FALSE) 

“葛思德”是“spatstat”内发现的功能图书馆。它是G函数或最近邻函数，它显示独立轴上粒子间的距离，以及在从属轴上寻找最近​​邻粒子的概率。因此，它从y = 0开始，并在y = 1时达到饱和点。

如果你绘制simdat.Gest，你会注意到曲线是's'形的，这意味着它从y = 0开始，到y = 1结束。因此，我认识了向量G值，这是因变量。因此，信息处于正确的方向以符合上述模型。

您可能还会注意到，我已经自动设置A = 1。这是因为G（r）总是饱和为1，所以我没有把它留在公式中。

我的问题是，我不断收到错误。对于上面的例子中，我得到这个错误：

Error in nls(rev(Gvalues) ~ (1/(1 + (B^2) * (R^2))^((C + 1)/2)), data = GvsR_dataframe, : 
    singular gradient 

我也一直在收到此错误：

Error in nls(Gvalues1 ~ (1/(1 + (B^2) * (x^2))^((C + 1)/2)), data = G_r_dataframe, : 
    step factor 0.000488281 reduced below 'minFactor' of 0.000976562 

我没有一个线索，其中的第一个错误的来源。然而，第二，我认为是因为我没有为B和C选择合适的起始值而发生的。

我希望有人能帮我弄清楚第一个错误来自哪里。另外，选择起始值以避免第二个错误的最有效方法是什么？

谢谢！

Answer 1

如上所述，您的问题很可能是起始值。有两种策略可以使用：

1. 使用蛮力来查找起始值。请参阅软件包nls2以了解执行此操作的功能。
2. 尝试对起始值进行合理的猜测。 根据您的值，可以线性化模型。

G = (1/(1 + (B^2)*(R^2))^((C+1)/2))

ln(G)=-(C+1)/2*ln(B^2*R^2+1)

If B^2*R^2 is large, this becomes approx. ln(G) = -(C+1)*(ln(B)+ln(R)), which is linear.

If B^2*R^2 is close to 1, it is approx. ln(G) = -(C+1)/2*ln(2), which is constant.

（请检查错误，这是昨天深夜，由于足球比赛。）

编辑其他信息已提供后： 数据看起来像它遵循的累积分配功能。如果它像鸭子一样嘎嘎叫，它很可能是一只鸭子。实际上?Gest表示估计CDF。

library(spatstat) 
data(simdat) 
simdat.Gest <- Gest(simdat) 
Gvalues <- simdat.Gest$rs 
Rvalues <- simdat.Gest$r 
plot(Gvalues~Rvalues) 

#let's try the normal CDF 
fit <- nls(Gvalues~pnorm(Rvalues,mean,sd),start=list(mean=0.4,sd=0.2)) 
summary(fit) 
lines(Rvalues,predict(fit)) 
#Looks not bad. There might be a better model, but not the one provided in the question.