这个程序的初衷是在知道总体平均值的情况下，测试sigma square 的区间，大概的意思就是，每次生成10个正态分布， mean=4, sd=4。然后重复一百次，这样这一百次，每次的平均值就符合正态分布。然后sum(Xi-mu)^2/chisquare(alfa/2,df=10)和sum(Xi-mu)^2/chisquare(1-alfa/2,df=10)就应该是他的上下区间。然后计算在100个CI中，有多少个CI包括true value.
程序如下>
ci=c(0.99,0.95,0.90) #计算alfa
alfa1=(1-ci)/2
alfa2=1-(1-ci)/2
t.var=NULL
CI=function(var,x.var){ #子程序
for (i in 1:3){
n=0
M1=x.var/qchisq(alfa1,10) #计算上限
M2=x.var/qchisq(alfa2,10) #计算下限
if (var>=M2&var<=M1){
n=n+1}
cat("1-alfa",ci,"\n")
cat("in CI",n,"\n")
cat("proporation",n/100,"\n")}}
#normal
mean=4
x.mean=NULL
var=16
x.var=0
for (p in 1:100){
T=rnorm(10,mean=4,sd=4) #生成
T.mean=mean(T)
x.mean=cbind(x.mean,T.mean)
x.var=x.var+(T.mean-mean)^2} #计算sum(xi-mu)^2
CI(var=var,x.var=x.var)
我提供一下命令作为参考。
计算sigma的CI，在知道总体平均值 mu的情况下
ci=c(0.99,0.95,0.90)
alfa1=(1-ci)/2
alfa2=1-(1-ci)/2
n=c(5,10,50,100)
for (p in 1:4){
T=rnorm(n[p],mean=4,sd=4)
t.mean=mean(T)
mean=4
t.var=0
for (q in 1:n[p]){
t.var=t.var+(T[q]-mean)^2}
cat("when n=",n[p],"\n")
for (i in 1:3){
M1=t.var/qchisq(alfa1,n[p])
M2=t.var/qchisq(alfa2,n[p])
cat("1-alfa",ci)
cat("CI is:[",M2,M1,"]","\n")}}
这个 程序可以用的。但是按照道理就改一下就可以了，但是却不行。很郁闷 。
计算mu 的CI 在知道 sigma的情况下。
ci=c(0.99,0.95,0.90)
alfa=1-(1-ci)/2
CI=function(x.mean,x.var,mean){
for (i in 1:3){
n=0
for (j in 1:100){
M=qnorm(alfa)*sqrt(x.var/10)
if (abs(x.mean[j]-mean)<=M){
n=n+1}}
cat("1-alfa",ci,"\n")
cat("in CI",n,"\n")
cat("proporation",n/p,"\n")}}
#normal
x.mean=NULL
x.var=16
for (q in 1:100){
T=rnorm(10,mean=4,sd=4)
x.mean=cbind(x.mean,mean(T))}
CI(x.mean=x.mean,x.var=x.var,mean=4)
大家帮忙看看谢谢了 。
扫码加好友，拉您进群
全部版块
我的主页
收藏
