JuliaStats/MCMC.jl — pull request #19 by dmbates (2013-09-18 17:32:30+00:00)
Title: Fix logit link and add const to global objects

The inverse of the logit link should have a negative sign in the
exponent. This is present in other examples but not in these two.

It might be a good idea to contrast the MCMC results, at least in the
documentation, with those from fitting a GLM. Right now it is kind of
clumsy to fit a GLM with a given model matrix but that could be
changed.

With the data from `example/binomial.jl` and the definition of the coefficients
```jl
julia> m = model(ex, vars=zeros(nbeta), gradient=true)
MCMCLikelihoodModel(##ll#319,# function,nothing,nothing,##ll#340,nothing,nothing,[:vars=>PDims(1,(10,))],10,[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0])

julia> res = run(m * RWM(0.05), steps=1000:10000)
10 parameters, 9001 samples (per parameter), 4.1 sec.

julia> res = run(m * RWM(0.05), steps=1000:10000)
10 parameters, 9001 samples (per parameter), 3.3 sec.

julia> head(res.samples)
6x10 DataFrame:
vars.1 vars.2 vars.3 vars.4 vars.5 vars.6 vars.7 vars.8 vars.9 vars.10
[1,] -0.21506 -0.443009 -0.0923568 -0.765506 -1.50518 1.45962 -0.166601 -1.21317 -0.510049 -0.470377
[2,] -0.180667 -0.460005 -0.0881983 -0.730407 -1.50652 1.42352 -0.257553 -1.22122 -0.427078 -0.489092
[3,] -0.172506 -0.516425 -0.0515133 -0.704989 -1.46919 1.46088 -0.209047 -1.23566 -0.464969 -0.438487
[4,] -0.172506 -0.516425 -0.0515133 -0.704989 -1.46919 1.46088 -0.209047 -1.23566 -0.464969 -0.438487
[5,] -0.172506 -0.516425 -0.0515133 -0.704989 -1.46919 1.46088 -0.209047 -1.23566 -0.464969 -0.438487
[6,] -0.172506 -0.516425 -0.0515133 -0.704989 -1.46919 1.46088 -0.209047 -1.23566 -0.464969 -0.438487

julia> [colwise(mean, res.samples) beta0 ] # mean sample vs original coefs
10x2 Array{Any,2}:
[-0.21482366575183443] -0.0717078
[-0.4725536125628145] -0.250912
[-0.020901613339271534] -0.100819
[-0.5815990199432688] -0.512176
[-1.405229813825649] -1.19319
[1.3544778459403048] 1.22552
[-0.2580605077816288] -0.246421
[-1.0614956381864762] -1.00204
[-0.423581899678947] -0.59563
[-0.5593022590885337] -0.681701

julia> res = run(m * HMC(2, 0.1), steps=1000:10000)
10 parameters, 9001 samples (per parameter), 10.0 sec.

julia> [colwise(mean, res.samples) beta0 ] # mean sample vs original coefs
10x2 Array{Any,2}:
[-0.2046847981894676] -0.0717078
[-0.462927351516865] -0.250912
[-0.024959659831684532] -0.100819
[-0.5805002542280536] -0.512176
[-1.413034972454002] -1.19319
[1.358553312357889] 1.22552
[-0.2584049999023584] -0.246421
[-1.058018576520552] -1.00204
[-0.42451432264431094] -0.59563
[-0.5546427132866298] -0.681701

julia> res = run(m * NUTS(), steps=1000:10000)
10 parameters, 9001 samples (per parameter), 17.4 sec.

julia> [colwise(mean, res.samples) beta0 ] # mean sample vs original coefs
10x2 Array{Any,2}:
[-0.20488317740629458] -0.0717078
[-0.465418897468125] -0.250912
[-0.024554388672484415] -0.100819
[-0.5822459997386241] -0.512176
[-1.4147503896216957] -1.19319
[1.3579861510776314] 1.22552
[-0.2573139709672702] -0.246421
[-1.0575141534744146] -1.00204
[-0.42458208644443585] -0.59563
[-0.5552305354105672] -0.681701

julia> fr = DataFrame(y = Y,x1=X[:,2],x2=X[:,3],x3=X[:,4],x4=X[:,5],x5=X[:,6],x6=X[:,7],x7=X[:,8],x8=X[:,9],x9=X[:,10]);

julia> using GLM

julia> glm(:(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9), fr, Binomial())

Formula: y ~ :(+(x1,x2,x3,x4,x5,x6,x7,x8,x9))

Coefficients:

10x4 DataFrame:
Estimate Std.Error z value Pr(>|z|)
[1,] -0.206247 0.0878745 -2.34707 0.0189218
[2,] -0.467542 0.0924892 -5.0551 4.30171e-7
[3,] -0.0246491 0.0844057 -0.292031 0.770263
[4,] -0.58779 0.0903566 -6.50522 7.75793e-11
[5,] -1.4265 0.117335 -12.1575 5.23478e-34
[6,] 1.37016 0.112457 12.1839 3.78862e-34
[7,] -0.259542 0.0864448 -3.0024 0.00267858
[8,] -1.06695 0.105179 -10.1441 3.51962e-24
[9,] -0.430233 0.0900796 -4.77614 1.78689e-6
[10,] -0.56062 0.0927221 -6.04624 1.48267e-9
```