※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。

# n = 19 (nonsmoker)
# n = 49 (former smoker)
# n = 32 (smoker)
smoke.num <- c(rep(1, 19), rep(2, 49), rep(3, 32))
label.smoke <- c("Nonsmoker","Former smoker","Smoker")
smoke <- factor(smoke.num, levels=1:3, labels=label.smoke)

mu <- mean(log(smoke.num))
sg <- sqrt(var(log(smoke.num)))

# Empirical distribution
Pi <- cumsum(as.numeric(table(smoke)) / length(smoke))
Pi <- Pi[-length(Pi)]

# Continuous critical value
CrV <- exp(mu + sg * qnorm(Pi))

# Simulation
set.seed(27011)
nsim <- 1000
smoke.sim.continuous <- exp(rnorm(nsim, mean=mu, sd=sg))
smoke.sim <- cut(smoke.sim.continuous, breaks=c(0,CrV,max(smoke.sim.continuous)), labels=label.smoke)

table(smoke.sim)

# Plot
hist(smoke.sim.continuous, freq=F, xlab='Smoke', main='')
lines(density(smoke.sim.continuous))
abline(v=CrV, col=4)