R: 統計檢定動態 Demo

寫完了信賴區間的動態 Demo, 乾脆又加寫了統計假設剪定的動態 Demo

 

?View Code LANGUAGE
###################################################################
# Test
###################################################################
#require(BHH2)
 
n = 30
nRepeat=500
mu0 = 60
mu1 = 63
sd = 5
alpha = 0.05
# delay time between each point
delaytime = 0.01
###############################
alpha2 = alpha/2
z = qnorm(1-alpha2)
 
xbarsd = sd/sqrt(n)
 
limL0 = mu0 - 3*xbarsd
limH0 = mu0 + 3*xbarsd
 
limL1 = mu1 - 3*xbarsd
limH1 = mu1 + 3*xbarsd
 
xbars = replicate(nRepeat,mean(rnorm(n,mu0,sd)))
xbar.p = z*xbarsd + mu0
 
ymax = max(dnorm(mu0,mu0,xbarsd),dnorm(mu1,mu1,xbarsd))
 
plot(mids,dnorm(0)*runif(length(mids)),
	xlim=c(limL0,limH1),ylim=c(0,ymax),
    xlab = expression(bar(X)),ylab = "counts",
    type="n",main=eval(substitute(expression(paste("模擬產生的",
		bar(X)," 檢定值數目:",nR)),list(nR=nRepeat)))
	)
 
fH0 = function(x) dnorm(x,mu0,xbarsd)
curve(fH0,limL0,limH0,col=4,lwd=2,add=T)
 
fH1 = function(x) dnorm(x,mu1,xbarsd)
curve(fH1,limL1,limH1,col=2,lwd=2,add=T)
abline(v=xbar.p,col=3)
abline(v=mu0,col=4,lty=2)
abline(v=mu1,col=2,lty=2)
text(xbar.p-0.5,ymax*2/3,paste("臨界值",round(xbar.p,2)),col=3)
text(limL0,ymax/2,expression(H[0]),col=4,cex=2)
text(limH1,ymax/2,expression(H[1]),col=2,cex=2)
 
legend("topleft",legend=c(
				eval(substitute(expression(H[0]:mu == mu0),list(mu0=mu0))) , 
			    eval(substitute(expression(H[1]:mu == mu1),list(mu1=mu1)))	),
				lty=c(1,1),col=c(4,2),lwd=c(2,2))
 
#----------------------------------------------
# H0
#----------------------------------------------
y = hist(xbars,breaks=50,plot=F)
mids = y$mids
xcounts = y$counts
maxcount = max(xcounts)
 
xs = rep(mids,xcounts)
s2 = function(x){
  if (x == 0)
    return(NULL)
  else 
    return(seq(x))
}
ys = unlist(sapply(xcounts,s2))		
 
n2 = length(xs)
pidx = sample(n2)
for (i in pidx)
{
  Sys.sleep(delaytime)
  points(xs[i],dnorm(mu0,mu0,xbarsd)*ys[i]/maxcount,col=4)
}
 
cat("超過臨界值 ",round(xbar.p,2)," 的檢定值數量: ",sum(xbars >= xbar.p),"個檢定值\n")
cat("超過臨界值 ",round(xbar.p,2)," 的檢定值比例: ",round(sum(xbars >= xbar.p)/nRepeat,2),"\n")
 
#--------------------------------------------
# H1
#--------------------------------------------
 
xbars2 = replicate(nRepeat,mean(rnorm(n,mu1,sd)))
 
y = hist(xbars2,breaks=50,plot=F)
mids = y$mids
xcounts = y$counts
maxcount = max(xcounts)
 
xs = rep(mids,xcounts)
s2 = function(x){
  if (x == 0)
    return(NULL)
  else 
    return(seq(x))
}
 
ys = unlist(sapply(xcounts,s2))			
 
n2 = length(xs)
pidx = sample(n2)
for (i in pidx)
{
  Sys.sleep(delaytime)
  points(xs[i],dnorm(mu1,mu1,xbarsd)*ys[i]/maxcount,col=2)
}
 
cat("超過臨界值 ",round(xbar.p,2)," 的檢定值數量: ",sum(xbars2 <= xbar.p),"個檢定值\n")
cat("超過臨界值 ",round(xbar.p,2)," 的檢定值比例: ",round(sum(xbars2 <= xbar.p)/nRepeat,2),"\n")