2013-08-06 5 views
2

수식> r² = (x-h) ² + (y-k) ² 인 원 맞춤에서 반경의 예측 간격을 계산하고 싶습니다. r- 원의 반경, 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() 그러므로 내가 nlsLM를 사용하여 원의 적합을 계산하기 위해 노력하고, nls.lm 작동하지 않습니다. (나는 손으로 계산하지만, 내 Designmatrix을 만드는 문제를 가질 수) .`

을 그래서 내가 뭘하려 다음 :

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

질문 A : 결과

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

: 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 : 사용 된 용어에 대한 명확한 설명을 위해 2 개의 위키 백과 링크 추가 : (c.f. 아래)

confidence interval

prediction interval

EDIT3 : 질문 (들)의 일부 다른 표현

Edit4 : 선형 회귀를위한 작업 예를

답변

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

enter image description here

0

찾을 수있는 솔루션입니다 h, k, r은 기본 R의 최적화 함수를 사용합니다. 근본적으로 최적화 할 데이터가 포함 된 클로저 인 비용 함수를 만듭니다. 나는 RSS 가치에, 그렇지 않으면 우리는 -Inf에 갈 것입니다. 로컬 optima 문제가 있으므로 몇 번 실행해야합니다 ...REPL에서이 명령을

# 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 
+1

h, k 및 r을 찾는 것이 문제가 아니 었습니다. 그것은 이미 포스터의 코드에서 "out"이라는 결과의 일부였습니다. –

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 샘플링 또는 유사한 방법을 사용하여 사후 신뢰 구간을 얻으십시오.

+0

확인. 나는 비선형 적으로도 예측이 가능하다고 생각했다. 얼마나 조잡한가. 저는 내 vcov에서 값을 선택하는 MCMC 시뮬레이션을 살펴 보았습니다. 나는 코딩에 아직 도전하지 않았다. 빨리 게시 할 것입니다. – Toby

+0

명확성 - 신뢰 구간이 존재하는지 여부를 결정하는 것은 함수의 선형성과 비선형 성이 아닙니다. 함수가 정의 된 확률 분포를 설명하는지 여부입니다. – ben