x = seq(-1, 1, 0.01);  y <- x;
a = -0.5; b=6;
th = exp(a+b*x)/(1+exp(a+b*x));
k = 1; # k = 5;
for (j in 1:length(x)){
  y[j] = rbinom(1, k, th[j]);
}
plot(x, y)

z <- y/k;
mod <- glm(z ~ x, family = binomial)
summary(mod)

zfit <- mod$fitted.values
yfit <- k * zfit
lines(x, yfit)
lines(x, k*th, col='red')