你是要把这个过程在R里实现吗

嗯嗯，是的，我自己写了一个代码，但是结果很差

set.seed(1234)
N<-1:30###重复的次数
n<-sum(N)
n1<-2000
b<-rep(1,3)###原始系数
x<-matrix(rnorm(length(N)*(length(b)-1)),ncol = length(b)-1)
X<-cbind(rep(1,length(N)),x)
colnames(X)<-c('x0','x1','x2')
p<-exp(X%*%b)/(1+exp(X%*%b))
y<-matrix(NA,length(N))
for(i in 1:length(N)){
  y[i,]<-sum(rbinom(N[i],1,p[i]))
}##生成的y
cf<-function(X,y){
  u<-list(NULL)
  length(u)<-length(N)
  for(i in 1:length(N)){
    for(j in 1:n1)
      u[[i]][[j]]<-sort(runif(N[i]))
  }
  k<-list(NULL)
  length(k)<-length(N)
  for(i in 1:length(N)){
    for(l in 1:n1){
      if(N[i]==1){
        k[[i]]<-mean(u[[i]])
      }else{
        h<-u[[i]][[1]]
        for(j in 1:n1) {
          if(j<=n1){
            h<-u[[i]][[j]]+h
          }}
        h<-h-u[[i]][[1]]
        k[[i]]<-h/n1
      }
    }
  }
  p1<-matrix(NA,length(N))
  for(i in 1:length(N)){
    if(y[i]==0){
      p1[i]<-0.5*k[[i]]
    }else if(y[i]==N[i]){
      p1[i]<-0.5*(k[[i]][as.numeric(y[i])]+1)
    }else{
      p1[i]<-0.5*(k[[i]][as.numeric(y[i])]+k[[i]][as.numeric(y[i])+1])
    }
  }
  library('MASS')##p1就是p*
  B1<-ginv(t(X)%*%X)%*%t(X)%*%log(p1/(1-p1))
  return(B1)
}
B1<-cf(X,y)
B1###模拟得出的系数
大神们看看我的程序哪里有错吗？

set.seed(1234)
N<-1:30###重复的次数
n<-sum(N)
n1<-2000
b<-rep(1,3)###原始系数
x<-matrix(rnorm(length(N)*(length(b)-1)),ncol = length(b)-1)
X<-cbind(rep(1,length(N)),x)
colnames(X)<-c('x0','x1','x2')
p<-exp(X%*%b)/(1+exp(X%*%b))
y<-matrix(NA,length(N))
for(i in 1:length(N)){
  y[i,]<-sum(rbinom(N,1,p))
}##生成的y
cf<-function(X,y){
  u<-list(NULL)
  length(u)<-length(N)
  for(i in 1:length(N)){
    for(j in 1:n1)
      u[][[j]]<-sort(runif(N))
  }
  k<-list(NULL)
  length(k)<-length(N)
  for(i in 1:length(N)){
    for(l in 1:n1){
      if(N==1){
        k[]<-mean(u[])
      }else{
        h<-u[][[1]]
        for(j in 1:n1) {
          if(j<=n1){
            h<-u[][[j]]+h
          }}
        h<-h-u[][[1]]
        k[]<-h/n1
      }
    }
  }
  p1<-matrix(NA,length(N))
  for(i in 1:length(N)){
    if(y==0){
      p1<-0.5*k[]
    }else if(y==N){
      p1<-0.5*(k[][as.numeric(y)]+1)
    }else{
      p1<-0.5*(k[][as.numeric(y)]+k[][as.numeric(y)+1])
    }
  }
  library('MASS')##p1就是p*
  B1<-ginv(t(X)%*%X)%*%t(X)%*%log(p1/(1-p1))
  return(B1)
}
B1<-cf(X,y)
B1###模拟得出的系数
想问一问我的程序哪里有错误，为什么模拟的结果比glm的差很多