brian-j-smith / Mamba.jl

Markov chain Monte Carlo (MCMC) for Bayesian analysis in julia
Other
253 stars 52 forks source link

Model Especification - Error #136

Closed lucasfreitas1988 closed 6 years ago

lucasfreitas1988 commented 6 years ago

Friends, good afternoon!

I need help specifying a model in Mamba. My model results in the error "MethodError: no method matching Distributions.Binomial (:: Int64, :: Array {Float64,1})", although it was based on the example Seeds: Random Effect Logistic Regression, available in the Mamba documentation (http://mambajl.readthedocs.io/en/latest/examples/seeds.html). It seems that the Binomial distribution can not receive an input array, which really makes sense, but that is exactly what occurs in the example cited.

I really do not know where there is error, follow my code:

Dados = Dict{Symbol, Any}( :y => [14, 11, 16, 20, 18, 11, 30, 31, 21, 21, 6, 30, 19, 40, 21, 23, 23, 22, 21, 21, 30, 22, 22, 23, 23, 36, 38, 39, 33, 32],

  :x1 =>
  [1.93, 1.87, 2.00, 2.01, 1.98,
   1.75, 1.99, 2.03, 2.20, 1.95,
   1.90, 2.20, 2.00, 1.90, 1.91,
   1.79, 2.05, 2.10, 2.20, 1.96,
   1.98, 2.04, 2.07, 2.10, 1.97,
   1.99, 2.09, 1.89, 2.08, 1.79],

 :x2 =>
  [1, 2, 2, 3, 1,
   1, 1, 2, 3, 3,
   1, 2, 2, 1, 1,
   1, 2, 2, 2, 1,
   1, 3, 2, 2, 1,
   3, 2, 1, 2, 1],

 :x3 =>
  [1, 1, 2, 10, 23,
  12, 5, 6, 32, 14,
  11, 25, 10, 9, 1,
   1, 5, 13, 13, 1,
  15, 20, 25, 4, 1,
  4, 30, 1, 23, 16],

  :n =>
    [12, 10, 13, 11, 12,
     10, 14, 13, 11, 11,
     11, 11, 11, 11, 11,
     11, 13, 12, 11, 11,
     10, 12, 12, 13, 13,
     16, 18, 19, 13, 12]

)

Dados[:N] = length(Dados[:y])

model = Model( y = Stochastic(1, (alpha, beta1, beta2, beta3, x1, x2, x3, N, n) -> UnivariateDistribution[ begin p = invlogit(alpha + beta1 x1[i] + beta2 x2[i] + beta3 * x3[i])
Binomial(n[i], p) end for i in 1:N ], false ),

alpha = Stochastic( () -> Normal(0, 1000) ),

beta1 = Stochastic(1, s21 -> Normal(0, sqrt(s21)), false ),

beta2 = Stochastic(1, s22 -> Normal(0, sqrt(s22)), false ),

beta3 = Stochastic(1, s23 -> Normal(0, sqrt(s23)), false ),

s21 = Stochastic( () -> InverseGamma(0.001, 0.001) ),

s22 = Stochastic( () -> InverseGamma(0.001, 0.001) ),

s23 = Stochastic( () -> InverseGamma(0.001, 0.001) )

)

Could it be some initial value problem or scheme?

Thanks in advance, this is my second question in less than 10 days. Now, I'm trying to implement some models I've used in JAGS.

L

bdeonovic commented 6 years ago

In the seeds example the probability of success for an observation is defined as

p = invlogit(alpha0 + alpha1 * x1[i] + alpha2 * x2[i] +
                       alpha12 * x1[i] * x2[i] + b[i])

where alpha1, alpha2, and alpha12 are (stochastic) scalars. In your example code you have

p = invlogit(alpha + beta1 * x1[i] + beta2 * x2[i] + beta3 * x3[i])

where you defined beta1, beta2, and beta3 as (stochastic) vectors e.g.:

...
beta1 = Stochastic(1,
...

note that the 1 in the Stochastic call above means its is a vector of dimension 1. If you omit a number (like you did for alpha you get a stochastic scalar.

lucasfreitas1988 commented 6 years ago

Thank you very much!!!

You're right Benjamin. I am teaching Julia to my students by comparing them with JAGS-adjusted models, where I am more familiar. In fact, I and my students are learning Julia together and this is great. Everyone thinks the syntax is better and the codes simpler, especially compared to R and Python.

Best Regards!

bdeonovic commented 6 years ago

I'm glad to see your enthusiasm and excitement for julia; it is quite a wonderful language. Learning and teaching it at the same time sounds a bit intimidating. Good luck!

lucasfreitas1988 commented 6 years ago

It's funny, I confess rs. Starting from scratch in a new programming language is challenging, especially when the whole class is native to other languages. I'm trying to set a good example and we've solved all the doubts together, with your help, here from GitHub.

L.