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

我想用公式>r²=（x-h）2 +（y-k）2来计算半径的预测间隔。圆的半径，x，y，是高斯坐标，h，k，标记拟合圆的中心。如何计算拟合R的圆的预测间隔

# 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）的结果现在是在出并OUT1不同...

编辑2：增加了两个维基百科链接，用于澄清puprose上使用的术语：（c.f.下文）

confidence interval

prediction interval

EDIT3：问题（S）的一些改述

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

来源

2013-08-06 Toby

我有一个很难搞清楚你想做什么。让我来说明数据的外观和“预测”的内容。

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")

enter image description here

来源

2013-08-06 22:03:42

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

# 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

来源

2013-08-06 22:07:16 wespiserA

找到h，k和r不是问题。这已经是海报代码中名为“out”的结果的一部分。 –

我觉得这个问题是不是在目前的形式交代。任何基于线性模型的函数都需要预测变量为输入设计矩阵的线性函数。 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采样或类似的方法得到后验置信区间。

来源

2013-08-07 14:19:16 ben

好的。我认为预测也适用于非线性。我多么sl。。我一直在寻找MCMC模拟从我的vcov中选择值。我还没有挑战编码。将尽快发布。 – Toby

要清楚的是 - 函数的线性与非线性决定了置信区间是否存在;这是函数是否描述了一个确定的概率分布。 – ben

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

回答

相关问题