华中科技大学 陈菲 企业统计 Project.doc

企业统计项目小组成员：指导老师：陈非个人报告姓名： 学号： 指导老师： 陈非 1.打开hitters.post1975.txt，选择数据：一个球员-赛季要求至少有100次击球和一次本垒打，即AB=100&HR=1，然后对平均本垒打（HR/AB）画图，并从中选出平均本垒打（HR/AB）最高和最低的球员-赛季。程序如下：Data-read.table(hitters.post1975.txt,header=TRUE,sep=,)data=100&Data$HR=1,dim(Data)dim(data)输出：AvgHR=data$HR/data$ABAvgHRhist(AvgHR,main=Histogram of AvgHR,xlab=AvgHR,ylab=Freq)Max.AvgHR=max(AvgHR)Min.AvgHR=min(AvgHR)Max.AvgHRMin.AvgHRdatawhich(AvgHR=Max.AvgHR),c(2,3),datawhich(AvgHR=Min.AvgHR),c(2,3),输出：结果，平均本垒打最高的是bondsba01 2001赛季。平均本垒打最低的是millafe01 1975赛季。2.给出参数为a、b的似然函数和对数似然函数，然后用contour（）函数给出似然函数的等值线图。程序如下：# graphing likelihood over range of alpha and beta x-seq(0,0.15,0.000015) density.1 - dbeta(x,2,3) density.2 - dbeta(x,4,30) density.3 - dbeta(x,7,50) density.4 - dbeta(x,3,60) mindensity - min(density.1,density.2,density.3,density.4) maxdensity - max(density.1,density.2,density.3,density.4) hist(AvgHR,freq=F) lines(x,density.1,type=l,col=black,lwd=2, ,ylim=c(mindensity,maxdensity)lines(x,density.2,type=l,col=green,lwd=2)lines(x,density.3,type=l,col=blue,lwd=2)lines(x,density.4,type=l,col=red,lwd=2)legend(0.4,25,c(Beta(1,1),Beta(1/2,1/2),Beta(1/3,1/3),Beta(0,0),col=c(black,blue,red,green),lwd=2) #loglikBeta functionn - length(AvgHR)numgrid - 100loglikBeta - function(a,b) loglik 0 & b 0) loglik - 0 for (i in 1:n) loglik - loglik + a*log(AvgHRi)+b*log(1-AvgHRi) loglik - loglik -log(beta(a, b) loglikalpharange - ppoints(numgrid)*0.2+1.9 # a between 1.9 and 2.1betarange - ppoints(numgrid)*7+67 # b between 67 and 75full - matrix(NA,nrow=numgrid,ncol=numgrid)for (i in 1:numgrid) for (j in 1:numgrid) fulli,j - loglikBeta(alpharangei,betarangej) full - exp(full - max(full)full - full/sum(full)contour(alpharange,betarange,full,xlab=a,ylab=b,drawlabels=F,main=Contour Plot)运行的结果：上图为参数为a、b的似然函数的等值线图，a的范围是1.92.1，b的范围是6775。3.用格点搜寻方法找出使似然函数最大的a和b。# calculating probabilities for grid sampler:alphamarginal - rep(NA,numgrid)for (i in 1:numgrid) alphamarginali - sum(fulli,)betaconditional - matrix(NA,nrow=numgrid,ncol=numgrid)for (i in 1:numgrid) for (j in 1:numgrid) betaconditionali,j - fulli,j/sum(fulli,) # plotting marginal distribution of alphapar(mfrow=c(1,1)plot(alpharange,alphamarginal,type=l,main=marginal dist. of alpha)# plotting marginal distribution of beta given alphaalpharange20alpharange50alpharange80par(mfrow=c(3,1) #将三幅图按行排列显示在一个窗口中plot(betarange,betaconditional25,type=l,main=dist. of beta for alpha = 24.9)plot(betarange,betaconditional50,type=l,main=dist. of beta for alpha = 29.9)plot(betarange,betaconditional75,type=l,main=dist. of beta for alpha = 34.9)# sampling grid values:alpha.samp - rep(NA,10000)beta.samp - rep(NA,10000)for (m in 1:10000) a - sample(1:100,size=1,replace=F,prob=alphamarginal) b - sample(1:100,size=1,replace=F,prob=betaconditionala,) alpha.sampm - alpharangea beta.sampm - betarangebpar(mfrow=c(1,1)contour(alpharange,betarange,full,xlab=alpha,ylab=beta,drawlabels=F,col=2)points(alpha.samp,beta.samp)# calculating posterior means/intervals for alpha and betapar(mfrow=c(2,1)hist(alpha.samp,main=Alpha Samples)hist(beta.samp,main=Beta Samples)a-mean(alpha.samp)b-mean(beta.samp) 输出：var0-var(AvgHR) #AvgHR的真实方差var1-a*b/(a+b)2)*(a+b+1) #在参数条件下方差var0var1输出：alpha.sampsort - sort(alpha.samp)beta.sampsort = 0)/10000par(mfrow=c(1,1)plot(t,AvgHR)4. 用New-Raphson方法求使似然函数最大的a和b。beta.2D.NR - function(olda,oldb,AvgHR) n - length(AvgHR) Jaa - n*trigamma(olda+oldb)-trigamma(olda) #对a求二阶导 Jab - n*trigamma(olda+oldb) #先对a 求导，再对b 求导 Jbb - n*trigamma(olda+oldb)-trigamma(oldb) #对b求二阶导 J - rbind(c(Jaa,Jab),c(Jab,Jbb) Jinv - solve(J) #求J的逆矩阵 ga - sum(log(AvgHR)-n*digamma(olda)+n*digamma(olda+oldb) #a的一阶导 gb - sum(log(1-AvgHR)-n*digamma(oldb)+n*digamma(olda+oldb) gvec - t(t(c(ga,gb) oldvec - t(t(c(olda,oldb) newvec - oldvec - Jinv%*%gvec newvecold - c(5,10) # starting valueitersmat - NULLitersmat - rbind(itersmat,old)new - c(-1,-1)diff 0.00000001) new - beta.2D.NR(old1,old2,AvgHR) itersmat - rbind(itersmat,c(new1,new2) diff - sum(abs(new-old) old - new itersmatnumiters - length(itersmat,1)alpha.hat - itersmatnumiters,1beta.hat - itersmatnumiters,2alpha.hatbeta.hat输出：abline(v=alpha.hat,col=red)abline(h=beta.hat,col=red)5. 用正态分布代替精英球员，截掉一半的正态分布代替非精英球员# function for calculating log-likelihood:loglik.mix - function(alpha,mu0,x,sigma0,sigma1) phi0 - dnorm(x,mu0,1) #Elite player phi1 - dnorm(x,0,1) #Not elite player loglik - sum(log(alpha*phi0 + (1-alpha)*phi1) loglik# Expectation function:Estep - function(alpha,mu0,sigsq0,sigsq1) ind - rep(NA,n) for (i in 1:n) prob1 - 2*(1-alpha)*dnorm(AvgHRi,mean=0,sd=sqrt(sigsq1) prob0 - alpha*dnorm(AvgHRi,mean=mu0,sd=sqrt(sigsq0) indi - prob0/(prob0+prob1) ind# Maximization functionMstep - function(ind) alpha - sum(ind)/n mu0 - sum(ind*AvgHR)/sum(ind) mu1 - sum(1-ind)*AvgHR)/sum(1-ind) sigsq0 - sum(ind*(AvgHR-mu0)2)/sum(ind) sigsq1 - sum(1-ind)*(AvgHR)2)/sum(1-ind) c(alpha,mu0,sigsq0,sigsq1)# Starting values for EM algorithm:curalpha - 0.5curmu0 - 1cursigsq0 - 2cursigsq1 - 3itermat - c(curalpha,curmu0,cursigsq0,cursigsq1)# Running EM algorithmdiff - 1numiters 0.000001 | numiters = 100) numiters - numiters+1 curind - Estep(curalpha,curmu0,cursigsq0,cursigsq1) curparam - Mstep(curind) curalpha - curparam1 curmu0 - curparam2 cursigsq0 - curparam3 cursigsq1 - curparam4 itermat - rbind(itermat,curparam) diff - max(abs(itermatnumiters,-itermatnumiters-1,) print (numiters)par(mfrow=c(2,3)for (i in 1:4) plot(1:length(itermat,1),itermat,i)lastiter - length(itermat,1)itermatlastiter,6. # Different Starting values: curalpha - 0.5curmu0 - 1cursigsq0 - 2cursigsq1 - 2itermat - c(curalpha,curmu0,cursigsq0,cursigsq1)diff - 1numiters 0.000001 | numiters = 100) numiters - numiters+1 curind - Estep(curalpha,curmu0,cursigsq0,cursigsq1) curparam - Mstep(curind) curalpha - curparam1 curmu0 - curparam2 cursigsq0 - curparam3 cursigsq1 - curparam4 itermat - rbind(itermat,curparam) diff - max(abs(itermatnumiters,-itermatnumiters-1,) print (numiters)par(mfrow=c(2,3)for (i in 1:4) plot(1:length(itermat,1),itermat,i)lastiter - length(itermat,1)itermatlastiter,#graphing histx-seq(0,0.20,0.000020) par(mfrow=c(1,1) hist(AvgHR,ylim=range(0,30),freq=F) dens-curalpha*dnorm(x,curmu0,sqrt(cursigsq0)+ (1-curalpha)*2*dnorm(x,0,sqrt(cursigsq1) lines(x,dens,col = red)7. abreu-datadata$name=abreubo01,AvgHR01-abreu$HR/abreu$AB p0- cumprod(pnorm(AvgHR01,curmu0,sqrt(cursigsq0) p0- p08 p1- cumprod(pnorm(AvgHR01,0,sqrt(cursigsq1) p1- p18 p - (curalpha*p0)/(curalpha*p0+(1-curalpha)*p1) p8. # function for calculating log-likelihood:loglik.mix - function(alpha,mu0,x,sigma0,a1,b1) phi0 - dnorm(x,mu0,sigma0) #Elite player phi1 - dgamma(x,a1,b1) #Not elite player loglik - sum(log(alpha*phi0 + (1-alpha)*phi1) loglik# Expectation function:Estep - function(alpha,mu0,sigma0,a1,b1) ind - rep(NA,n) for (i in 1:n) prob1 - 2*(1-alpha)*dgamma(AvgHRi,a1,b1) prob0 - alpha*dnorm(AvgHRi,mean=mu0,sd=sqrt(sigma0) indi - prob0/(prob0+prob1) ind# Maximization functionMstep - function(ind) alpha - sum(ind)/n mu0 - sum(ind*AvgHR)/sum(ind) a1 - sum(1-ind)*AvgHR)/sum(1-ind) sigsq0 - sum(ind*(AvgHR-mu0)2)/sum(ind) b1 - sum(1-ind)*(AvgHR)2)/sum(1-ind) c(alpha,mu0,sigsq0,a1,b1)# Starting values for EM algorithm:curalpha - 0.5curmu0 - 1cura1 - 2cursigsq0 - 2curb1 - 3itermat - c(curalpha,curmu0,cura1,cursigsq0,curb1)# Running EM algorithmdiff - 1numiters 0.000001 | numiters = 100) numiters - numiters+1 curind - Estep(curalpha,curmu0,cura1,cursigsq0,curb1) curparam - Mstep(curind) curalpha - curparam1 curmu0 - curparam2 cura1 - curparam3 cursigsq0 - curparam4 curb1 - curparam5 itermat - rbind(itermat,curparam) diff - max(abs(itermatnumiters,-itermatnumiters-1,) print (numiters)par(mfrow=c(2,3)for (i in 1:5) plot(1:length(itermat,1),itermat,i)lastiter - length(itermat,1)itermatlastiter,#change starting values# Starting values for EM algorithm:curalpha - 0.5curmu0 - 4cura1 - 4cursigsq0 - 1curb1 - 2itermat - c(curalpha,curmu0,cura1,cursigsq0,curb1)# Running EM algorithmdiff - 1num