Как вычислить интервалы прогнозирования для окружности, подходящей для R

Я хочу вычислить интервал предсказания радиуса из круга, подходящего по формуле > r² = (xh)²+(yk)². 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. (Я мог бы вычислить это вручную, но у меня возникли проблемы с созданием моей матрицы дизайна).`

Итак, вот что я попробовал дальше:

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

Вопрос 1а: Как 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 ссылки на Википедию для уточнения используемой терминологии: (см. ниже)

доверительный интервал

интервал прогнозирования

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

Edit4: добавлен рабочий пример для линейной регрессии.


r regression predict nls
person Toby    schedule 06.08.2013    source источник

Ответы (3)


arrow_upward
2
arrow_downward

Мне трудно понять, что вы хотите сделать. Позвольте мне проиллюстрировать, как выглядят данные, и кое-что о «прогнозе».

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

введите здесь описание изображения

person IRTFM    schedule 06.08.2013

arrow_upward
1
arrow_downward

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

person ben    schedule 07.08.2013
comment
Ok. я думал, что прогнозирование работает и для нелинейных. Как небрежно с моей стороны. Я просматривал симуляции MCMC, выбирая значения из моего файла vcov. Я еще не пробовал кодировку. Будет сообщение как можно скорее. - person Toby; 07.08.2013
comment
Чтобы было ясно, не линейность и нелинейность функции определяет, существуют ли доверительные интервалы; это то, описывает ли функция определенное распределение вероятностей. - person ben; 07.08.2013

arrow_upward
0
arrow_downward

Вот решение для поиска 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]

Запустите эту команду в REPL, чтобы наблюдать за локальными минутами:

o=optim(f=tFun,par=runif(3),method="CG");o$par
person wespiserA    schedule 06.08.2013
comment
Найти h,k и r не проблема. Это уже было частью результата, указанного в коде плаката. - person IRTFM; 07.08.2013