brian-j-smith / Mamba.jl

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

Problem with symmetric matrix #128

Open EnzoHG opened 6 years ago

EnzoHG commented 6 years ago

This is my very first time in julia, and i can't seem to follow the traceback and crack the problem. i have visited number of forums and i haven't found a similar issue. i checked out the birats example as a guide for this model2 distribution2

and i think i have all the syntax right( mu = R should be mu_cero). this is my code: `

 using Mamba
 using Distributions
 ##Data
 data=Dict{Symbol,Any}(
 :z => [1.0 1.0 0.0 0.0 1.0 1.0 0.0 1.0 0.0 0.0
  0.0 0.0 1.0 1.0 1.0 1.0 0.0 1.0 1.0 0.0
  0.0 1.0 1.0 0.0 0.0 1.0 1.0 1.0 0.0 1.0
  0.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0
  1.0 1.0 0.0 1.0 0.0 1.0 1.0 1.0 1.0 0.0],

 :y => [0.0 1.0 1.0 1.0 1.0 0.0 1.0 0.0 0.0 0.0
    1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 1.0
    1.0 0.0 0.0 1.0 1.0 1.0 0.0 1.0 1.0 0.0
    1.0 0.0 0.0 1.0 0.0 0.0 1.0 1.0 0.0 1.0
    1.0 1.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 1.0],

  :x => [1.0 1.0 0.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0
    1.0 0.0 1.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0
    1.0 0.0 1.0 0.0 0.0 1.0 0.0 1.0 1.0 1.0
    1.0 1.0 0.0 1.0 0.0 1.0 1.0 1.0 0.0 0.0
    1.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 0.0 1.0]
 )
data[:sig] = eye(5)
data[:J] = 10
data[:I] = 5

##Model
model=Model(

z = Stochastic(2,
 (b, x, y,J,I) ->
    UnivariateDistribution[
    begin
        p=(1+exp(transpose(b[j])*(y[i, j]-x[i, j])))^(-1)
        Bernoulli(p)
    end
    for i in 1:I, j in 1:J
    ],
 false
),

b=Stochastic(2,
 (mu, sig, J) ->
    MultivariateDistribution[
        MvNormal(mu,sig)
        for j in 1:J
    ]
),

mu=Stochastic(1,
 (mu0, sig0)->
    MultivariateNormal(mu0, sig0)
),

mu0=Stochastic(1,
    ()->Uniform(0,1),
 false
),

sig0=Stochastic(2,
    ()->Uniform(0,1),
 false
)

)

##Initial Values
inits=[
   Dict(:z => data[:z], :b =>zeros(5,10), :mu => rand(5), :mu0 => rand(5), :sig0 => rand(5,5)),
   Dict(:z => data[:z], :b =>rand(5,10), :mu => ones(5), :mu0 => ones(5), :sig0 => rand(5,5))
]

#Sample Scheme
 scheme = [AMWG(:mu, 1.0), AMWG(:b, 1.0),  Slice(:sig0, 10.0)]
 setsamplers!(model,scheme)

 ##MCMC simulation
 sim = mcmc(model, data, inits, 10000, burnin=2500, thin=2, chains=2)
 describe(sim)`

where the error and the traceback are:

LoadError: ArgumentError: matrix is not symmetric/Hermitian. This error can be avoided by calling cholfact(Hermitian(A)) which will ignore either the upper or lower triangle of the matrix. while loading C:\Users\Enzo Hernandez\Downloads\tarea2_mamba.jl, in expression starting on line 80 in at base\ in #mcmc#29 at Mamba\src\model\mcmc.jl:29 in setinits! at Mamba\src\model\initialization.jl:24 in setinits! at Mamba\src\model\initialization.jl:11 in setinits! at Mamba\src\model\dependent.jl:169 in at base\ in Distributions.MvNormal at Mamba\src\distributions\constructors.jl:29 in cholfact at base\linalg\cholesky.jl:349

I really want to learn to use mamba, so any piece of advice would be very much appreciate it.

bdeonovic commented 6 years ago

@EnzoHG sorry for getting back to you so late.

Problem is with sig0. Currently you are setting it so that all the values are uniform(0,1) which produces a nonsymmetric covariance matrix. This does not match your model specification where sig0 = A'A

EnzoHG commented 6 years ago

@bdeonovic no problem, at this point i'm doing for my self. i come to this solution ``

 A=Logical(2,
     sig0 -> transpose(sig0)*sig0
  ),

sig0=Stochastic(2,
 (I,J) ->UnivariateDistribution[
        Uniform(0,1)
        for i in 1:I, j in 1:J
    ],
 false
 )    ´´

but this doesn't seem to work, it's no longer the symmetric issue but this:

LoadError: DimensionMismatch("incompatible distribution for stochastic node") while loading C:\Users\Enzo Hernandez\Downloads\tarea2_mamba.jl, in expression starting on line 89 in at base\ in #mcmc#29 at Mamba\src\model\mcmc.jl:29 in setinits! at Mamba\src\model\initialization.jl:24 in setinits! at Mamba\src\model\initialization.jl:11 in setinits! at Mamba\src\model\dependent.jl:171

i followed the example in the birats model with the normal distribution matrix.

i really appreciate that you take the time to help me out 👍 thanks again

bdeonovic commented 6 years ago

Your model is still wrong. You cannot have sig0 as a covariance matrix because a matrix filled with random values is not guaranteed to be positive definite. Also your new code still doesn't follow your model that you have a picture of in original post.

your model is something like this:

 A=Stochastic(1,
     ()-> Uniform(0,1)
  ),

sig0=Logical(2,
 (A) -> A'A,
 false
 )