using DelimitedFiles, Distributions, LinearAlgebra
affair = readdlm("affair.csv", ',')
y = affair[:,1];
X = affair[:,2:end];
n, k = size(X);
#find the MLE and the information matrix
S = ones(k,1); #score
betat = (X'*X)\(X'*y); #initial guess
e = 10^(-5);  #tolerance level
while sum(abs.(S)) > e  #stopping criterion
    Xbetat = X*betat;
    phi = pdf.(Normal(0,1),Xbetat);
    Phi = cdf.(Normal(0,1),Xbetat);
    global S = sum(repeat(y.*phi./Phi-(1 .- y).*phi./(1 .-Phi),outer = [1,k]).*X,dims=1)';
    d = phi.^2 ./ (Phi.*(1 .-Phi));
    global IM = X'*diagm(vec(d))*X; # information matrix
    global betat = betat + IM\S;
end

burnin = 500;
nloop = 5000+burnin;
V = IM\I;  #scale matrix for the proposal
B = cholesky(Hermitian(V)).L;
nu = 5; #df for the proposal
b0 = zeros(k,1); #prior mean
V0 = 10*I; #prior covariance
# log-posterior density
logf(b)= (y'*log.(cdf.(Normal(0,1),X*b)) + (1 .-y)'*log.(cdf.(Normal(0,1),-X*b))- 0.5*(b-b0)'*(V0\(b-b0)))[1];
# log-proposal density
logprop(b)  = (-0.5*(k+nu)*log.(1 .+ (b-betat)'*(V\(b-betat))/nu))[1];
store_beta = zeros(nloop,k);
beta = betat;
for i = 1:nloop
# candidate draw from the t proposal
    global beta
    betac = betat + B*randn(k,1)*sqrt(nu/rand(Gamma(nu/2,2)));
    rho = logf(betac)-logf(beta) + logprop(beta)-logprop(betac);
    if exp(rho) > rand()
        beta = betac;
    end
    store_beta[i,:] = beta';
end
store_beta = store_beta[burnin+1:end,:]; # discard the burn-in
println(mean(store_beta,dims=1))
println(std(store_beta,dims=1))