우도 함수와 시뮬레이션을 기반으로 한 프로 비트 시뮬레이션을 만들었습니다. 이들 모두는 아래 코드로 복제 할 수 있습니다.시뮬레이션 데이터에서 플롯하기위한 신뢰 구간 추가 R
my.probit <- function(y,x) {
# use OLS to get start values
par <- lm(y~x)$coefficients
ytilde <- 2*y-1
# Run optim
res <- optim(par,probit.ll,hessian=TRUE,ytilde=ytilde,x=x)
# Return point estimates and SE based on the inverse of Hessian
names(res$par) <- c('a','b')
se=sqrt(diag(solve(res$hessian)))
names(se) <- c('a','b')
return(list(par=res$par,se=se,cov=solve(res$hessian)))
}
그리고 이것은 시뮬레이션 모델을 생성하는 기능입니다 :
probit.data <- function(N=100,a=1,b=1) {
x <- rnorm(N)
y.star <- a + b*x + rnorm(N)
y <- (y.star > 0)
return(as.data.frame(cbind(y,x,y.star)))
}
을
probit.ll <- function(par,ytilde,x) {
a <- par[1]
b <- par[2]
return( -sum(pnorm(ytilde*(a + b*x),log=TRUE)))
}
이 추정을 할 수있는 기능은 다음과 같습니다
이
는 우도 함수n 크기를 100으로 동일하게 시뮬레이션합니다.
,probit.data100 <- function(N=100,a=2,b=1) {
x <- rnorm(N)
y.star <- a + b*x + rnorm(N)
y <- (y.star > 0)
return(as.data.frame(cbind(y,x,y.star)))
}
#predicted value
se.probit.phat100 <- function(x, par, V) {
z <- par[1] + par[2] * x
# Derivative of q w.r.t. alpha and beta
J <- c(dnorm(z), dnorm(z)*par[2])
return(sqrt(t(J) %*% V %*% J) )
}
dat100 <- probit.data100()
res100 <- my.probit(dat100$y,dat100$x)
res100
이 기능은 이하의 비 - 파라 부트 스트랩 방식에 기초하여 상기 신뢰 구간을 산출한다 (샘플 기능을 통지 사용되고)
N <- dim(probit.data(N=100, a=1, b=1))[1]
npb.par <- matrix(NA,100,2)
colnames(npb.par) <- c("alpha","beta")
npb.eystar <- matrix(NA,100,N)
for (t in 1:100) {
thisdta <- probit.data(N=100, a=1, b=1)[sample(1:N,N,replace=TRUE),]
npb.par[t,] <- my.probit(thisdta$y,thisdta$x)$par
}
아래 사항이 기능은 부트 스트랩 출력 및 신뢰도를 정리
내가이 같은 그래프 (아래 하나) 플롯하지만 processres에 따라 신뢰 구간을 추가 할processres <- function(simres) {
z <- t(apply(simres,2,function(x) { c(mean(x),median(x),sd(x),quantile(x,c(0.05,0.95))) }))
rownames(z) <- colnames(simres)
colnames(z) <- c("mean","median","sd","5%","95%")
z
}
processres(npb.par)
위 기능 : 간격 내가 플롯 싶은 것이 있습니다. 이 신뢰 구간을 어떻게 플롯에 추가 할 수 있습니까?
x <- seq(-5,5,length=100)
plot(x, pnorm(1 - 0.5*x), ty='l', lwd=2, bty='n', xlab='x', ylab="Pr(y=1)")
rug(dat100$x)
다른 플롯 코드 및/또는 패키지에 대해서도 열려 있습니다. 그래프를 추가하고 신뢰 구간을 추가하면됩니다.
감사합니다. 이 이제 (즉, 평균 알파 & 베타 값을 사용하여) 예측 곡선을 나타내는, 정확하게 rnorm
이러한 수단을 통과 :
UPDATE :
,이 멋지다! –