/
mcmc.jl
75 lines (61 loc) · 1.89 KB
/
mcmc.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
using DataFrames
using GLM
# function initialize_mcmc(rollcall,colnames)
# meanrollcall=copy(rollcall)
# for c in colnames:
# meanrollcall[c] = meanrollcall[c]-mean(meanrollcall[c])
# end
# iterator=eachrow(rollcall)
# for i in 1:length(iterator)
# sum=0
# for j in 1:length(iterator[i])
# sum=sum+iterator[i][j]
# end
# meanrow=sum/length(iterator[i])
# for k in colnames
# meanrow[i,k] = df[i,k] - meanrow
# end
# end
# end
function initialize_x0(rollcall)
# subtract out column and row means, as instructed
meanrollcall=copy(rollcall)
for j in 1:size(meanrollcall,2)
meanrollcall[:,j] = meanrollcall[:,j] - mean(meanrollcall[:,j])
end
for i in 1:size(meanrollcall,1)
meanrollcall[i,:] = meanrollcall[i,:] - mean(meanrollcall[i,:])
end
trans = transpose(meanrollcall)
# get nxn correlation matrix so that we can find eigenvalues for initial x0
corr = meanrollcall*trans
eigs = eig(corr)
max_eig_loc = findmax(eigs[1])[2]
x0 = eigs[2][:,max_eig_loc]
x0
end
function initialize_beta(rollcall,init_x)
#implement probit regression model
# for each bill (j = 1:m)
# perform regression using GLM on the pertaining col/ of the roll call matrix
# and all the x0s
# this will give one beta (x coeff) and one alpha (intercept) (= B_j and a_j)
beta0 = Array(Float64, size(rollcall,2))
alpha0 = Array(Float64, size(rollcall,2))
for j in 1:size(rollcall,2)
data = DataFrame(X=init_x,Y=rollcall[:,j])
result = glm(Y~X,data,Binomial(),ProbitLink())
beta0[j] = coef(result)[2]
alpha0[j] = coef(result)[1]
end
#return initial beta and alpha
(beta0, alpha0)
end
# function for iterative MCMC process
function update_vars(rollcall, x, beta, alpha)
# first sample y* from Normal distribution N(x(t-1)beta(t-1) - alpha(t-1)),1)
# truncated based on observed 1 or 0 below or above 0
# for abstentions we sample y* from untruncated mu(t-1)
# next sample beta
(x, beta, alpha)
end