probit-Stan.R

# Simulando dados de probit

set.seed(1928506)
x1 <- rnorm(5000,0,1)
P <- pnorm(1.5 + 1.2*x1) #- 2.3*x2 + 1.9*x3 + rnorm(5000))
y <- rbinom(5000, 1,P)


# Dá uma olhada. Aqui dá o resultado certinho.
glm(y ~ x1 , family=binomial(link="probit")) 




# Agora vamos para o Bayes maldito!
data <- list(N=nrow(X), x1=x1, y=y)
stanstr <-
'
data {
int<lower=1> N;
vector[N] x1;
int<lower=0,upper=1> y[N];
}

transformed data {


}
parameters {
real beta;
real alpha;

}
model {

y[N] ~ bernoulli(Phi(alpha + beta * x1[N]));
}
generated quantities { // Generated quantities block. Not used presently.
}
'

fit <- stan(model_code = stanstr, data=data, iter=12000, warmup=2000, thin=10, chains=3)