# Case 010: constructed data; not empirical observations. Base R only.
set.seed(20261007)
x <- rep(0:1, each=100)
e <- rep(seq(-3,3,length.out=100),2)
y_full <- 10 + 4*x + e
observed <- rep(FALSE,200)
for(g in 0:1) {
  k <- if(g==0) 80L else 40L
  ix <- which(x==g)
  observed[ix[c(seq_len(k/2),seq.int(101-k/2,100))]] <- TRUE
}
y <- y_full; y[!observed] <- NA_real_
means <- sapply(0:1,function(g) mean(y[x==g],na.rm=TRUE))
counts <- sapply(0:1,function(g) sum(observed & x==g))
variances <- sapply(0:1,function(g) var(y[x==g],na.rm=TRUE))
# Saturated two-group normal observed-data likelihood, fixed target weights 1/2.
fiml_mean <- mean(means)
fiml_se <- sqrt(sum(.25*variances*(counts-1)/counts^2))
# Proper normal-model MI with independent reference priors p(mu,sigma^2) ~ 1/sigma^2.
M <- 2000L; Q <- U <- numeric(M)
for(i in seq_len(M)) {
  completed <- y
  for(g in 0:1) {
    j <- g+1L
    variance_draw <- (counts[j]-1)*variances[j]/rchisq(1,counts[j]-1)
    mean_draw <- rnorm(1,means[j],sqrt(variance_draw/counts[j]))
    missing <- x==g & !observed
    completed[missing] <- rnorm(sum(missing),mean_draw,sqrt(variance_draw))
  }
  Q[i] <- mean(completed)
  U[i] <- sum(sapply(0:1,function(g) .25*var(completed[x==g])/100))
}
W <- mean(U); B <- var(Q); T <- W+(1+1/M)*B
mi_se <- sqrt(T)
# Large-complete-sample Rubin df; no small-sample adjustment here.
df <- (M-1)*(1+W/((1+1/M)*B))^2
interval <- mean(Q)+c(-1,1)*qt(.975,df)*mi_se
print(counts)
print(c(deletion_mean=mean(y,na.rm=TRUE),fiml_mean=fiml_mean,fiml_SE=fiml_se,
        MI_mean=mean(Q),MI_SE=mi_se,within=W,between=B,total=T))
print(interval)
stopifnot(identical(as.integer(counts),c(80L,40L)),
          abs(mean(y,na.rm=TRUE)-34/3)<1e-10,abs(fiml_mean-12)<1e-10,
          T>W,abs(mean(Q)-12)<.1)
print(R.version.string)