Mamba - logistic regression (bernoulli response).

New to Julia
1

Please can someone post a really simple logistic regression implemented in Mamba, i.e. something like this but that actually works. I have been unable to figure it out and have not found another reference that might help. Thanks.

using Distributed, Mamba, LinearAlgebra
@everywhere using Mamba
using DataFrames

dat = Dict{Symbol, Any}(
    :y => [1 0 1 1 0 0 0 0 1 0 1 0 1 1 0 0 0 0 0 1 1 1 1 1 0 1 1 1 0 1 1 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 1],
    :trt => [0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1],
    :remote => [1 1 0 1 1 1 0 1 1 1 1 1 1 1 1 1 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 1 0 1 0 0 1 1 0 1 0 1 0 0 0 1 1 0 0]
  )
dat[:N] = length(dat[:y])
dat[:xmat] = [ones(dat[:N]) transpose(dat[:trt]) transpose(dat[:remote])]

inits = [
    Dict{Symbol, Any}(
      :y => dat[:y],
      :b => rand(Normal(0, 1), size(dat[:xmat], 2))
    )
    for i in 1:3 # 3 chains
    ]

scheme1 = [NUTS(:b)]

model = Model(

    y = Stochastic(1,
    (p, N) -> UnivariateDistribution[(
      Bernoulli(p[i])) for i in 1:N], 
false 
    ),

    p = Logical(1,
      (xmat, b) -> invlogit.(xmat * b),
      false
    ),

    b = Stochastic(1,
      () -> MvNormal(3, sqrt(100))
    )
)

setsamplers!(model, scheme1)
sim1 = mcmc(model, dat, inits,
          1000,
          burnin=250,
          thin=2,
          chains=3)

Aside - the above gives:

MethodError: no method matching Array{Float64,1}(::Array{Int64,2})

You might have used a 2d row vector where a 1d column vector was required.
Note the difference between 1d column vector [1,2,3] and 2d row vector [1 2 3].
You can convert to a column vector with the vec() function.
Closest candidates are:
  Array{Float64,1}(::AbstractArray{S,N}) where {T, N, S} at array.jl:497
  Array{Float64,1}() where T at boot.jl:413
  Array{Float64,1}(!Matched::UndefInitializer, !Matched::Int64) where T at boot.jl:394
  ...
convert(::Type{Array{Float64,1}}, ::Array{Int64,2}) at array.jl:489
setinits!(::ArrayStochastic{1}, ::Model, ::Array{Int64,2}) at dependent.jl:168
setinits!(::Model, ::Dict{Symbol,Any}) at initialization.jl:11
setinits!(::Model, ::Array{Dict{Symbol,Any},1}) at initialization.jl:24
#mcmc#23(::Int64, ::Int64, ::Int64, ::Bool, ::Function, ::Model, ::Dict{Symbol,Any}, ::Array{Dict{Symbol,Any},1}, ::Int64) at mcmc.jl:30
(::getfield(Mamba, Symbol("#kw##mcmc")))(::NamedTuple{(:burnin, :thin, :chains),Tuple{Int64,Int64,Int64}}, ::typeof(mcmc), ::Model, ::Dict{Symbol,Any}, ::Array{Dict{Symbol,Any},1}, ::Int64) at none:0
top-level scope at none:0
include_string(::Module, ::String, ::String) at loading.jl:1005
include_string(::Module, ::String, ::String, ::Int64) at eval.jl:30
(::getfield(Atom, Symbol("##114#119")){String,Int64,String})() at eval.jl:94
withpath(::getfield(Atom, Symbol("##114#119")){String,Int64,String}, ::String) at utils.jl:30
withpath at eval.jl:46 [inlined]
#113 at eval.jl:93 [inlined]
with_logstate(::getfield(Atom, Symbol("##113#118")){String,Int64,String}, ::Base.CoreLogging.LogState) at logging.jl:397
with_logger at logging.jl:493 [inlined]
#112 at eval.jl:92 [inlined]
hideprompt(::getfield(Atom, Symbol("##112#117")){String,Int64,String}) at repl.jl:85
macro expansion at eval.jl:91 [inlined]
(::getfield(Atom, Symbol("##111#116")))(::Dict{String,Any}) at eval.jl:86
handlemsg(::Dict{String,Any}, ::Dict{String,Any}) at comm.jl:164
(::getfield(Atom, Symbol("##19#21")){Array{Any,1}})() at task.jl:259
2

Answering own question - but wondering why this does not work when I wrap it in a function… :roll_eyes:

using Distributed, Mamba, LinearAlgebra
# @everywhere using Mamba
using DataFrames

dat = Dict{Symbol, Any}(
    :y => [1.,0,1,1,0,0,0,0,1,0,1,0,1,1,0,0,0,0,0,1,1,1,1,1,0,1,1,1,0,1,1,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,1],
    :trt => [0.,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1],
    :remote => [1,1,0,1,1,1,0,1,1,1,1,1,1,1,1,1,0,0,0,1,0,0,1,0,0,0,0,0,0,1,0,0,1,0,1,0,0,1,1,0,1,0,1,0,0,0,1,1,0,0]
  )

dat[:N] = length(dat[:y])
# dat[:xmat] = [ones(dat[:N]) transpose(dat[:trt]) transpose(dat[:remote])]

inits = [
    Dict{Symbol, Any}(
      :y => dat[:y],
      :b1 => rand(Normal(0, 1)),
      :b2 => rand(Normal(0, 1))

    )
    for i in 1:3
    ]

scheme1 = [NUTS(:b1),
            NUTS(:b2)]

model = Model(

    y = Stochastic(1,
    (b1, b2, trt, remote, N) ->
      UnivariateDistribution[(
        p = invlogit(b1 + b2 * trt[i]);
        Bernoulli(p)) for i in 1:N
      ],
    false
    ),

    b1 = Stochastic(
      () -> Normal(0, 1000)
    ),
    b2 = Stochastic(
      () -> Normal(0, 1000)
    )
)

setsamplers!(model, scheme1)
sim1 = mcmc(model, dat, inits,
          1000,
          burnin=250,
          thin=2,
          chains=3)

sim1


describe(sim1)
3

OK, so the bizarre thing is when I wrap this in a function as in:

using Distributed, Mamba, LinearAlgebra
# @everywhere using Mamba
using DataFrames

function mcmctester()
  dat = Dict{Symbol, Any}(
      :y => [1.,0,1,1,0,0,0,0,1,0,1,0,1,1,0,0,0,0,0,1,1,1,1,1,0,1,1,1,0,1,1,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,1],
      :trt => [0.,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1],
      :remote => [1,1,0,1,1,1,0,1,1,1,1,1,1,1,1,1,0,0,0,1,0,0,1,0,0,0,0,0,0,1,0,0,1,0,1,0,0,1,1,0,1,0,1,0,0,0,1,1,0,0]
    )

  dat[:N] = length(dat[:y])
  # dat[:xmat] = [ones(dat[:N]) transpose(dat[:trt]) transpose(dat[:remote])]

  inits = [
      Dict{Symbol, Any}(
        :y => dat[:y],
        :b1 => rand(Normal(0, 1)),
        :b2 => rand(Normal(0, 1))

      )
      for i in 1:3
      ]

  scheme1 = [NUTS(:b1),
              NUTS(:b2)]

  model = Model(

      y = Stochastic(1,
      (b1, b2, trt, remote, N) ->
        UnivariateDistribution[(
          p = invlogit(b1 + b2 * trt[i]);
          Bernoulli(p)) for i in 1:N
        ],
      false
      ),

      b1 = Stochastic(
        () -> Normal(0, 1000)
      ),
      b2 = Stochastic(
        () -> Normal(0, 1000)
      )
  )

  setsamplers!(model, scheme1)
  sim1 = mcmc(model, dat, inits,
            1000,
            burnin=250,
            thin=2,
            chains=3)
  return sim1
end


mysim = mcmctester()

I get an error saying

MethodError: no method matching (::getfield(Main, Symbol("##105#106")))(::Model)
The applicable method may be too new: running in world age 25213, while current world is 25218.
Closest candidates are:
  #105(::Model) at none:0 (method too new to be called from this world context.)
setinits!(::ScalarStochastic, ::Model, ::Float64) at dependent.jl:163
setinits!(::Model, ::Dict{Symbol,Any}) at initialization.jl:11
setinits!(::Model, ::Array{Dict{Symbol,Any},1}) at initialization.jl:24
#mcmc#23(::Int64, ::Int64, ::Int64, ::Bool, ::Function, ::Model, ::Dict{Symbol,Any}, ::Array{Dict{Symbol,Any},1}, ::Int64) at mcmc.jl:30
#mcmc at none:0 [inlined]
mcmctester() at julia_tmp3.jl:49
top-level scope at none:0
include_string(::Module, ::String, ::String) at loading.jl:1005
include_string(::Module, ::String, ::String, ::Int64) at eval.jl:30
(::getfield(Atom, Symbol("##114#119")){String,Int64,String})() at eval.jl:94
withpath(::getfield(Atom, Symbol("##114#119")){String,Int64,String}, ::String) at utils.jl:30
withpath at eval.jl:46 [inlined]
#113 at eval.jl:93 [inlined]
with_logstate(::getfield(Atom, Symbol("##113#118")){String,Int64,String}, ::Base.CoreLogging.LogState) at logging.jl:397
with_logger at logging.jl:493 [inlined]
#112 at eval.jl:92 [inlined]
hideprompt(::getfield(Atom, Symbol("##112#117")){String,Int64,String}) at repl.jl:85
macro expansion at eval.jl:91 [inlined]
macro expansion at dynamic.jl:24 [inlined]
(::getfield(Atom, Symbol("##111#116")))(::Dict{String,Any}) at eval.jl:86
handlemsg(::Dict{String,Any}, ::Dict{String,Any}) at comm.jl:164
(::getfield(Atom, Symbol("##19#21")){Array{Any,1}})() at task.jl:259

If anyone can shed some light that would be greatly appreciated!

UPDATE

OK, so this exists as a reported issue in GitHub: Calling mcmc() from within a function · Issue #149 · brian-j-smith/Mamba.jl · GitHub.

4

This does not answer your question directly, but I added a logistic regression example for DynamicHMC.jl:

https://github.com/tpapp/DynamicHMCExamples.jl/blob/master/src/example_logistic_regression.jl

1 Like
5

Thank you, that’s great. I will have a look tomorrow.