고려대학교 세종캠퍼스 수시 합격률 예측(몬테칼로 추정) 논문

소개글

고려대학교 세종캠퍼스 수시 합격률 예측(몬테칼로 추정) 논문에 대한 보고서 자료입니다.

본문내용

',theta,'|x)')))
plot(theta,dbeta(theta,a+x,b+n-x),type='l',xlim = c(0.1,0.27),xlab=expression(theta),ylab=expression(paste(pi,'(',theta,'|x)')))
abline(v=c(out[1],out[2]))
abline(h=k,lty=2,col=2)
title('95% HPD region for theta')
mtext(round(k,3),side=2,at=k,col=2)
abline(v=c(qu[1],qu[2]),lty=2,col=3)
#legend(0.5,9,c('HPD Interval','Classical Interval'),lty=1:2,col=2:3)
legend(0.2,9,c('HPD Interval','Classical Interval'),lty=1:2,col=2:3)
#두 신뢰구간의 길이 비교
(hpd.len=out[2]-out[1]) #HPD Interval Length
(classic.len=qu[2]-qu[1]) #Classical Interval Length
#사전사후 비교#
#posterior with beta prior
y=seq(0,1,length=100)
par(mfrow=c(1,1))
plot(y,dbeta(y,1/2,1/2),xlab='x',ylab='f(x)',ylim=c(0,15),type='l',lty=2,main='Beta(1/2,1/2)')
#plot(y,dbeta(y,1/2,1/2),xlab='x',ylab='f(x)',ylim=c(0,15),xlim=c(0.15,0.25),type='l',lty=2,main='Beta(1/2,1/2)')
lines(y,dbeta(y,x+1/2,n-x+1/2),lty=1)
legend(0.6,13,c('prior','posterior'), lty=2:1)
abline(v=(a+x-1)/(a+b+n-2),col=2, lty=2)
abline(v=x/n,col=3, lty=3)
#legend(0.2,10,c('0.1804(Bayes)','0.1823(MLE)'),lty=2:3, col=2:3)
legend(0.6,10,c('0.1804(Bayes)','0.1823(MLE)'),lty=2:3, col=2:3)
#예측분포 11~14연도의자료바탕으로 15년도 예측
#predictive Distribution
par(mfrow=c(1,1))
m=56; z=0:56
pred.z=(gamma(m+1)/(gamma(z+1)*gamma(m-z+1)))*beta(a+x+z,b+n-x+m-z)/beta(a+x,b+n-x)
plot(z,pred.z,xlab='z',ylab='probability',type='h')
title('Predictive Distribution: a=1/2, b=1/2, n=170, x=31, m=56')
#예측 평균
(m.z<-sum(z*pred.z))
#예측 분산
sum((z-m.z)^2*pred.z)
#Monte Carlo method
a=1/2; b=1/2
n=170; x=31; m=56
N=10000; theta=1:N; z=1:N
for (i in 1:N) {
theta[i]=rbeta(1,a+x,b+n-x)
z[i]=rbinom(1,m,theta[i])
}
plot(table(z)/N,xlab='z',ylab='predictive density',type='h',main='')
title('Monte Carlo method: a=1/2, b=1/2, n=170, x=31, m=56')
mean(z); var(z)
#예측분포 15연도의자료바탕으로 16년도 예측(정시정원 6명)
par(mfrow=c(1,1))
n=56;x=4;
m=30; z=0:30
pred.z=(gamma(m+1)/(gamma(z+1)*gamma(m-z+1)))*beta(a+x+z,b+n-x+m-z)/beta(a+x,b+n-x)
plot(z,pred.z,xlab='z',ylab='probability',type='h')
title('Predictive Distribution: a=1/2, b=1/2, n=56, x=4, m=30')
#예측 평균
(m.z<-sum(z*pred.z))
#예측 분산
sum((z-m.z)^2*pred.z)
#Monte Carlo method
N=10000; theta=1:N; z=1:N
for (i in 1:N) {
theta[i]=rbeta(1,a+x,b+n-x)
z[i]=rbinom(1,m,theta[i])
}
plot(table(z)/N,xlab='z',ylab='predictive density',type='h',main='')
title('Monte Carlo method: a=1/2, b=1/2, n=56, x=4, m=30')
mean(z); var(z)
#예측분포 11~15연도의자료바탕으로 16년도 예측(정시정원 6명)
par(mfrow=c(1,1))
n=226;x=35;
m=30; z=0:30
pred.z=(gamma(m+1)/(gamma(z+1)*gamma(m-z+1)))*beta(a+x+z,b+n-x+m-z)/beta(a+x,b+n-x)
plot(z,pred.z,xlab='z',ylab='probability',type='h')
title('Predictive Distribution: a=1/2, b=1/2, n=226, x=35, m=30')
#예측 평균
(m.z<-sum(z*pred.z))
#예측 분산
sum((z-m.z)^2*pred.z)
#Monte Carlo method
a=1/2; b=1/2
n=226; x=35; m=30
N=10000; theta=1:N; z=1:N
for (i in 1:N) {
theta[i]=rbeta(1,a+x,b+n-x)
z[i]=rbinom(1,m,theta[i])
}
plot(table(z)/N,xlab='z',ylab='predictive density',type='h',main='')
title('Monte Carlo method: a=1/2, b=1/2, n=226, x=35, m=30')
mean(z); var(z)

키워드

합격률,   분석,   통계,   몬테칼로