如何计算 R 中的圆拟合的预测区间

Question

我希望用公式 > r² = (x-h)²+(y-k)² 计算半径的预测区间。 r-圆的半径，x，y，是高斯坐标，h，k，标记拟合圆的中心。

# data
x <- c(1,2.2,1,2.5,1.5,0.5,1.7)
y <- c(1,1,3,2.5,4,1.7,0.8)
# using nls.lm from minpack.lm (minimising the sum of squared residuals)
library(minpack.lm)

residFun <- function(par,x,y) {
  res <- sqrt((x-par$h)^2+(y-par$k)^2)-par$r
  return(res)
}
parStart <- list("h" = 1.5, "k" = 2.5, "r" = 1.7)
out <- nls.lm(par = parStart, x = x, y = y, lower =NULL, upper = NULL, residFun)

问题是，predict() 不适用于 nls.lm，因此我尝试使用 nlsLM 计算圆拟合。 （我可以手动计算，但在创建我的 Designmatrix 时遇到了麻烦）。`

这就是我接下来尝试的：

dat = list("x" = x,"y" = y)
out1 <- nlsLM(y ~ sqrt(-(x-h)^2+r^2)+k, start = parStart )

导致：

Error in stats:::nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

问题 1a：nlsLM() 如何处理圆形拟合？ （优点是通用的predict() 可用。 问题 1b：如何获得我的圆拟合的预测区间？

线性回归的例子（这就是我想要的圆形回归）

attach(faithful)     
eruption.lm = lm(eruptions ~ waiting) 
newdata = data.frame(waiting=seq(45,90, length = 272)) 
# confidence interval
conf <- predict(eruption.lm, newdata, interval="confidence") 
# prediction interval
pred <- predict(eruption.lm, newdata, interval="predict")
# plot of the data [1], the regression line [1], confidence interval [2], and prediction interval [3]
plot(eruptions ~ waiting)
lines(conf[,1] ~ newdata$waiting, col = "black") # [1]
lines(conf[,2] ~ newdata$waiting, col = "red") # [2]
lines(conf[,3] ~ newdata$waiting, col = "red") # [2]
lines(pred[,2] ~ newdata$waiting, col = "blue") # [3]
lines(pred[,3] ~ newdata$waiting, col = "blue") # [3]

亲切的问候

编辑摘要：

Edit1：重新排列 nlsLM 中的公式，但参数 (h,k,r) 结果现在在 out 和 out1 中有所不同...

Edit2：添加了 2 个 wikipedia 链接，用于澄清所用术语的 puprose：（参见下文）

confidence interval

prediction interval

Edit3：问题的一些改写

Edit4：添加了线性回归的工作示例

Answer 1

我认为这个问题以目前的形式无法回答。任何基于线性模型的predict() 函数都要求预测变量是输入设计矩阵的线性函数。 r^2 = (x-x0)^2 + (y-y0)^2 不是设计矩阵的线性函数（类似于 [x0 x y0 y]，所以我认为您无法找到可以给您置信区间的线性模型拟合。如果有人比我聪明的人有办法做到这一点，不过，我很想听听。

解决这类问题的一般方法是创建一个分层非线性模型，其中您的超参数将是x0 和y0（您的 h 和 k）在您的搜索空间，然后 r^2 将分布 ~N((x-x0)^2+(y-y0)^2, \sigma)。然后，您将使用 MCMC 抽样或类似方法来获得您的后验置信区间。

Answer 2

我很难弄清楚你想做什么。让我来说明一下数据是什么样的以及关于“预测”的一些内容。

plot(x,y, xlim=range(x)*c(0, 1.5), ylim=range(y)*c(0, 1.5))
lines(out$par$h+c(-1,-1,1,1,-1)*out$par$r, # extremes of x-coord
      out$par$k+c(-1,1,1,-1 ,-1)*out$par$r, # extremes of y-coord
      col="red")

那么我们所说的“预测区间”是什么？ （我确实意识到你在考虑一个圆圈，如果你只是想在这个背景上画一个圆圈，那也很容易。）

lines(out$par$h+cos(seq(-pi,pi, by=0.1))*out$par$r, #center + r*cos(theta)
      out$par$k+sin(seq(-pi,pi, by=0.1))*out$par$r, #center + r*sin(theta)
      col="red")

Answer 3

这是一个使用基数 R 的 optim 函数找到 h,k,r 的解决方案。您实际上创建了一个成本函数，它是一个包含您希望优化的数据的闭包。我必须 RSS 值，否则我们会去 -Inf。有一个局部最优问题，所以你需要运行几次...

# data
x <- c(1,2.2,1,2.5,1.5,0.5,1.7)
y <- c(1,1,3,2.5,4,1.7,0.8)

residFunArg <- function(xVector,yVector){

  function(theta,xVec=xVector,yVec=yVector){
  #print(xVec);print(h);print(r);print(k)
    sum(sqrt((xVec-theta[1])^2+(yVec-theta[2])^2)-theta[3])^2
  }
}

rFun = residFunArg(x,y);

o = optim(f=rFun,par=c(0,0,0))


h = o$par[1]
k = o$par[2]
r = o$par[3]

在 REPL 中运行此命令以观察本地分钟数：

o=optim(f=tFun,par=runif(3),method="CG");o$par