如何解决R的代码从均匀随机数中得出遵循卡方分布的直方图
我在自己的text-book(用日语写)中有一个代码,可以从均匀分布生成具有3个自由度的卡方分布。我对此进行了改进,并创建了一个代码来获得直方图,该直方图遵循具有4个自由度的卡方分布。这与R的分布函数非常吻合,所以我认为它可能正确运行(请参见下面的Box1)。
我试图进一步细化Box1的代码,以获得具有指定自由度的卡方分布后的直方图,但是它并没有很多错误。 (请参见Box2)
我的问题:
Box2的从均匀分布生成卡方分布的代码不能很好地工作。
请帮助我修复Box2密码中的错误。
大概是“ y
Box1:用于获取具有4个自由度的卡方分布的直方图的代码(可能正确运行)
ite <- 10000
sc <- numeric(ite) #★1
A<- c(20,20,20) #★2
for(i in 1:ite){
s<- runif(sum(A)*5) #★3
y<-ifelse(s<0.2,1,ifelse(s<0.4,2,ifelse(s<0.6,3,ifelse(s<0.8,4,5)))) #★4
z1 <- table(y)
z2 <- A*5
z3 <- (z1-z2)^2 /z2
sc[i] <- sum(z3)
}
hist(sc,ylim=c(0,0.35),breaks="Scott",freq=F)
curve(dchisq(x,4),add=T)
方框1的代码是基于以下事实而设计的; 如果将500 = sum(A)* 5个统一随机数分为五个大小相同的房间,则进入每个房间的数字的期望值为100。 在此,第一室,第二室,...和第五室由0≤x<0.2、0.2≤x <0.4,...以及0.8≤x≤1定义。我们可以从以下Box'1的table(y)的输出中看到这一点。当然,Box 1'的sum(table(y))始终为500。
Box1'在Box1的代码上逐步生成统一的随机数(x)的逻辑
A<- c(20,20)
s<- runif(sum(A)*5) #★3
y<-ifelse(s<0.2,5))))
table(y)
sum(table(y))
Box2:用于获得遵循自由度n(具有许多误差)的卡方分布的直方图的代码
chiq_dist_n<-function(numb,itr){
A<-numeric(numb) #★2
aa<-numeric(numb) #★4-1
for(i in 1:numb){
A[i]=20
} #★2
ntot=sum(A)
for(i in 1:numb){
if (i ==1){aa[i]= A[i]/ntot
}else{
aa[i]=aa[i-1]+(A[i]/ntot)
}
} #★4-2
sc<-numeric(itr) #★1
y<-numeric(ntot*numb) #★4-3
for(i in 1:itr){
x<-runif(ntot*numb)
for(k in 1:ntot*numb){
for(j in 1:numb){
if (x[k]<aa[numb-j+1]) {
y[k]<-j
} else {}
}
}#★3
z1<-table(y)
z2<-A*ntot
z3<-(z1-z2)^2/z2
sum(z3)
sc[i]<-sum(z3)
}
return(sc)
}
hist(chiq_dist(10,1000),freq=F)
Box2代码中生成y的部分被切成Box2'。 如果查看Box2'的表(y),您会看到太多y [i]为零。 我希望方框2'中表(y)的输出与方框1'中表(y)的输出大致相同。
Box2'在Box2的代码上逐步生成统一的随机数(x)的逻辑
A<- c(20,20)
ntot=sum(A)
numb=length(A)
aa<-numeric(numb)
for(i in 1:numb){
if (i ==1){aa[i]= A[i]/ntot
}else{
aa[i]=aa[i-1]+(A[i]/ntot)
}
} #★4-2
y<-numeric(ntot*numb)
x<-runif(ntot*numb)
for(k in 1:ntot*numb){
for(j in 1:numb){
if (x[k]<aa[numb-j+1]) {
y[k]<-j
} else {}
}
}#★3
table(y)
解决方法
您不需要费力就能打破随机均匀分布,只需使用cut()
并指定中断次数即可,例如:
set.seed(111)
v = runif(10)
[1] 0.59298128 0.72648112 0.37042200 0.51492383 0.37766322 0.41833733
[7] 0.01065785 0.53229524 0.43216062 0.09368152
cut(v,breaks=seq(0,1,length.out=numb+2),labels=1:5)
[1] 3 4 2 3 2 3 1 3 3 1
我不确定A或它的作用,但是为了模拟卡方,我假设您对标签1:(df + 1)进行了随机采样,其中df是自由度。如果将采样次数固定为500,那么我们知道每个中断的预期值为500 /(df + 1)。
因此无需更改太多代码。
chiq_dist_n<-function(numb,ite){
sc <- numeric(ite)
for(i in 1:ite){
x<- runif(500) #★3
y<- cut(x,labels=1:(numb+1))
z1 <- table(y)
z2 <- length(x)/(numb+1)
z3 <- (z1-z2)^2 /z2
sc[i] <- sum(z3)
}
hist(sc,ylim=c(0,0.35),breaks="Scott",freq=F,main=paste0("df=",numb))
curve(dchisq(x,numb),add=T)
}
我们尝试从4到9:
par(mfrow=c(3,2))
par(mar=c(2.5,2.5,2.5))
for(i in seq(2,12,2)){
chiq_dist_n(i,10000)
}
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。