添加置信區間從R中的仿真數據繪圖

Question

我已經創建了基於似然函數和仿真的probit仿真，所有這些都可以用下面的代碼複製。

這是似然函數：

probit.ll <- function(par,ytilde,x) { 
    a <- par[1] 
    b <- par[2] 
return( -sum(pnorm(ytilde*(a + b*x),log=TRUE))) 
} 

這是做估值的功能：

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

這模擬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 <- 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) 

我想繪製這樣的圖表（下一個），但添加一個基於processres置信區間上方運行。如何將這些置信區間添加到情節中？

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) 

我也開放給不同的劇情代碼和/或包。我只是想要一個基於這個模擬的圖表增加置信區間。

謝謝！

Answer 1

這裏的基礎上添加的模擬結果的陰影CI的方式：

UPDATE：這個現在（即，使用平均阿爾法&的β值）繪製的預期曲線，並正確地傳遞這些手段rnorm。

x <- seq(-5,5,length=100) 
plot(x, pnorm(1 - 0.5*x), ty='n', lwd=2, bty='n', xlab='x', ylab="Pr(y=1)", 
    xaxs = 'i', ylim=c(0, 1)) 

params <- processres(npb.par) 
sims <- 100000 
sim.mat <- matrix(NA, ncol=length(x), nrow=sims) 
for (i in 1:sims) { 
    alpha <- rnorm(1, params[1, 1], params[1, 3]) 
    beta <- rnorm(1, params[2, 1], params[2, 3]) 
    sim.mat[i, ] <- pnorm(alpha - beta*x) 
} 

CI <- apply(sim.mat, 2, function(x) quantile(x, c(0.05, 0.95))) 
polygon(c(x, rev(x)), c(CI[1, ], rev(CI[2, ])), col='gray', border=NA) 
lines(x, pnorm(params[1, 1] - params[2, 1]*x), lwd=2) 
rug(dat100$x) 
box()