mychallenge <- function (Nsim = 10^4) { library(mcsm) data(challenger) x <- challenger[, 2] y <- challenger[, 1] n <- length(x) a <- b <- rep(0, Nsim) # m <- mean(x) # x <- x - m res <- glm(y ~ x, binomial()) print(res$coef) a[1] <- res$coef[1] b[1] <- res$coef[2] for (t in 2:Nsim) { uni <- runif(n) * exp(y * (a[t - 1] + b[t - 1] * x))/(1 + exp(a[t - 1] + b[t - 1] * x)) z <- log(uni/(1-uni)) mina <- max(( z-b[t-1]*x)[y==1]) maxa <- min((-z-b[t-1]*x)[y==0]) a[t] <- runif(1, mina, maxa) minb <- -Inf maxb <- Inf if (sum((y==1)&(x>0))>0) minb <- max(minb, max((( z-a[t])/x)[(y==1) & (x>0)])) if (sum((y==0)&(x<0))>0) minb <- max(minb, max(((-z-a[t])/x)[(y==0) & (x<0)])) if (sum((y==1)&(x<0))>0) maxb <- min(maxb, min((( z-a[t])/x)[(y==1) & (x<0)])) if (sum((y==0)&(x>0))>0) maxb <- min(maxb, min(((-z-a[t])/x)[(y==0) & (x>0)])) b[t] <- runif(1, minb, maxb) } # data.frame(a=a-b*m, b=b) data.frame(a=a,b=b) }