r - R 的代码,用于从均匀随机数中获取遵循卡方分布的直方图
问题描述
我的教科书(用日语编写)中有一个代码,可以从均匀分布中生成具有 3 个自由度的卡方分布。我对此进行了改进并创建了一个代码来获得一个直方图,该直方图遵循具有 4 个自由度的卡方分布。这与 R 的分布函数非常吻合,所以我认为它可能工作正常(参见下面的 Box1)。
我试图进一步细化 Box1 的代码,以获得具有指定自由度的卡方分布的直方图,但它并没有出现很多错误。(见方框2)
我的问题:
Box2 从均匀分布生成卡方分布的代码效果不佳。
请帮助我修复 Box2 代码的错误。
可能“ y<-ifelse(x<0.2,1,ifelse(x<0.4,2,ifelse(x<0.6,3,ifelse(x<0.8,4,5)))) 的概括不起作用在方框 2 中。
Box1:获取遵循卡方分布的直方图的代码,具有 4 个自由度(可能工作正常)
ite <- 10000
sc <- numeric(ite) #★1
A<- c(20,20,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)
Box 1 的代码是基于以下事实设计的;如果500=sum(A)*5个均匀随机数被分成5个大小相同的房间,每个房间进入的数字的期望值为100。这里,第1个房间,第2个房间,...,第5个房间是由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) 的逻辑 (y)
A<- c(20,20,20,20,20)
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))))
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),ylim=c(0,0.35),breaks="Scott",freq=F)
生成 y 的 Box2 代码部分被剪切到 Box2' 中。如果您查看 Box2' 的 table(y),您会发现太多的 y[i] 为零。我希望 Box 2' 中 table(y) 的输出与 Box 1' 中 table(y) 的输出大致相同。
Box2'在 Box2 的代码上逐步生成均匀随机数 (x) 的逻辑 (y)
A<- c(20,20,20,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)
解决方案
您不需要 ifelse 来打破随机均匀分布,您可以使用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,breaks=seq(0,1,length.out=numb+2),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,2.5))
for(i in seq(2,12,2)){
chiq_dist_n(i,10000)
}
推荐阅读
- sql - 我需要多次查询单个表并在单个结果集中显示结果
- reactjs - ImageBackground 没有包装 stack.navigator 并且 stack.screens 不可见。(@react-navigation/stack": "^5.10.0 ) React-Native
- html - 尝试在侧边栏旁边对齐图像
- blazor - 如何使用 Blazor 检索 Web 组件自定义事件结果
- javascript - Redux:useSelector 在 useDispatch 之后没有更新状态
- python-3.x - 如何使用具有多个返回值的python赋值运算符
- android - 使用 ViewPager 的 ImageSlider 未显示在 destroyitem 中删除的图像
- python - 使用 Python 连续读取和绘制 CSV 文件
- python-3.x - 为什么 pip3 在 Ubuntu 20.04 上安装/更新软件包后要创建一个 kdewallet?
- gis - exif中使用什么CRS?