如何計算擬合R的圓的預測間隔

Question

我想用公式>r²=（x-h）2 +（y-k）2來計算半徑的預測間隔。圓的半徑，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）的結果現在是在出並OUT1不同...

編輯2：增加了兩個維基百科鏈接，用於澄清puprose上使用的術語：（c.f.下文）

confidence interval

prediction interval

EDIT3：問題（S）的一些改述

Edit4：增加了線性迴歸的工作示例

Answer 1

我有一個很難搞清楚你想做什麼。讓我來說明數據的外觀和「預測」的內容。

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 2

這裏找到一個解決方案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 

Answer 3

我覺得這個問題是不是在目前的形式交代。任何基於線性模型的函​​數都需要預測變量爲輸入設計矩陣的線性函數。 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採樣或類似的方法得到後驗置信區間。