R的代码从均匀随机数中获得服从卡方分布的直方图

Question

我的 text-book（用日语编写）中有一个代码可以从均匀分布生成具有 3 个自由度的卡方分布。我对此进行了改进并创建了一个代码来获得一个遵循具有 4 个自由度的卡方分布的直方图。这与 R 的分布函数非常吻合，所以我认为它可能正确工作（见下面的框 1）。

我试图进一步优化 Box1 的代码，以获得遵循具有指定自由度的卡方分布的直方图，但它没有工作，但有很多错误。（见方框 2）

我的问题：

The Box2's code to generate a chi-square distribution from a uniform distribution does not work well.
Please help me to fix the errors of the Box2's Code.

可能是“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)

Box2代码中生成y的部分被截取到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)

Answer 1

你不需要 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)
}

R的代码从均匀随机数中获得服从卡方分布的直方图

R's code to obtain a histogram following a chi-square distribution from uniform random numbers

statistics

r

chi-squared