probit model (gibbs sampling)

#probit model (Gibbs sampling)

library(dplyr)

n=100;m=3

X1=array(0,dim=c(n,m))

for(j in 1:m){

X1[,j]=rpois(n,sample(1:2,1))  

}

X2=array(0,dim=c(n,m))

for(j in 1:m){

X2[,j]=rpois(n,sample(20:30,1))  

}

y=c(rep(1,n),rep(0,n))

X=rbind(X1,X2)

data=cbind(X,y)

n=2*n


z=rpois(n,7)

library(mvtnorm)

ite=10000

precision=0.1

precision_vec=c()

for(i in 1:ite){


if(precision<0.98){  

beta=rmvnorm(1,solve(t(X)%*%X)%*%t(X)%*%z,solve(t(X)%*%X))

pi_z=exp(-(z-X%*%t(beta))^2/2)*ifelse((y>0 & z>=0)|(y==0 & z<0),1,0)

mat=X%*%t(beta)

z_past=c();

for(j in 1:n){

z_past=c(z_past,rnorm(1,mat[j],1)*ifelse(mat[j]>=0,1,0)*ifelse(y[j]==1,1,0)+rnorm(1,mat[j],1)*ifelse(mat[j]<0,1,0)*ifelse(y[j]==0,1,0))  

}


z=z_past

#print(beta)

precision=sum(ifelse((y==1 & z>=0)|(y==0 & z<0),1,0)
)/n

print(sum(ifelse((y==1 & z>=0)|(y==0 & z<0),1,0)
)/n)

precision_vec=c(precision_vec,precision)  

}else{

print("stop")  

}


}


precision_data=data.frame(times=1:length(precision_vec),precision=precision_vec)

plot(precision_data$times,precision_data$precision,type="p",col=2)