Dear Students,
Please find below the R code I have written on the fly last Friday:
Best regards,
Jean-Michel Marin
########
x <- runif(100,-10,10)
# beta1=2 and beta2=1
px <- 1/(1+exp(-2-x))
y <- rep(0,100)
for (i in 1:100)
y[i] <- sample(c(0,1),1,prob=c(1-px[i],px[i]))
X <- cbind(rep(1,100),x)
betahat <- solve(t(X)%*%X)%*%t(X)%*%y
n <- 10000
x <- runif(n,-10,10)
px <- pnorm(1+0.5*x)
y <- rep(0,n)
for (i in 1:n)
y[i] <- sample(c(0,1),1,prob=c(1-px[i],px[i]))
plot(x,y)
curve(pnorm(1+0.5*x),-10,10,add=TRUE,lwd=2,col="blue")
logit.model <- glm(y~1+x,family="binomial")
logit.model$coefficients
f <- function(x) 1/(1+exp(-logit.model$coefficients[1]-logit.model$coefficients[2]*x))
curve(f,-10,10,add=TRUE,lwd=2,col="green")
X <- cbind(rep(1,n),x)
betahat <- solve(t(X)%*%X)%*%t(X)%*%y
abline(a=betahat[1],b=betahat[2],lwd=2,col="red")
u <- runif(1000,0,1)
sum(0.1<=u & u<=0.9)
u1 <- runif(1000,0,1)
u2 <- runif(1000,0,1)
sum(0.1<=u1 & u1<=0.9 & 0.1<=u2 & u2<=0.9)
d <- 40
n <- 10000
u <- matrix(runif(n*d,0,1),n,d)
sum(apply(0.1<=u & u<=0.9,1,sum)==d)/n
########