Как вычислить интервалы предсказания для окружности, подходящей в 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) результаты теперь отличаются вне и OUT1 ...

Edit2: Добавлено 2 ссылки для википедии для выяснения урожая по терминологии: (c.f. ниже)

confidence interval

prediction interval

Edit3: Некоторые перефразировка вопрос (ы)

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. Существует локальная проблема с оптимизацией, поэтому вам нужно запустить это несколько раз ...

# 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] 

Выполнить эту команду в РЕПЛ соблюдать местные мин:

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

Answer 3

Я думаю, что этот вопрос не отвечает в его нынешнем виде. Любая функция 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 или аналогичный, чтобы получить ваши задние доверительные интервалы.