Sampler for discrete distributions

[ghost]

The samplers currently available work for continuous distributions, but I think there is not a sampler for discrete distributions.

I want to test a Bayesian variable selection algorithm that requires a Bernoulli distribution as a prior at the bottom of the hierarchical model.

As far as I know, Mamba allows you to create your own sampler, although it seems to be focused on conditional distributions for Gibbs sampling. I'm not sure if this can be used to create arbitrary samplers.

[brian-j-smith]

You are correct. The package does not currently have a discrete sampler. However, all the infrastructure needed to implement one is there. It should be pretty straightforward and easier than the continuous samplers. For discrete parameters with finite support, one could write a Sampler that would compute the probability mass function of the conditional distribution and sample from it directly at each step of the MCMC algorithm.

There would be some implementation details to work out and development, testing, and documentation would take time, but it is doable. It's something that I could probably get done in November.

Note that the performance of such an MCMC algorithm is going to depend greatly on the number of variables being selected.

[ghost]

Thank you. Great to know it's possible and (relatively) straightforward.

Should we consider this a feature request?

[goedman]

Brian, are you thinking MH and/or a Gibbs sampler?

[brian-j-smith]

Here's what I'm thinking...

The general algorithm implemented in the package does block sampling of parameters from their full conditional distributions. For a block that is all discrete (finite support) parameters, one can determine all possible values that the parameters could take on. Mamba provides a function that could then be used to compute the unnormalized pdf at each value. Divide the unnormalized pdfs by their sum and you have the exact full conditional probabilities. That allows for a draw to be generated directly (Gibbs sampling) from the full conditional for the block.

A performance consideration is that the pdfs will need to be computed at each value at each iteration of the overall MCMC algorithm. The number of evaluations could be very large if there are a large number of discrete parameters in a block.

That's the first thought that popped into my head anyway. Would love to hear any other ideas that you might have.

[goedman]

Thanks Brian,

I’ll study some of the samplers first to get a better understanding of Mamba’s structure.

My question was triggered by https://www.stat.tamu.edu/~fliang/papers/2005/equation.pdf .

But I’m not even sure if this is applicable at all in this case. My level of understanding is basically William Bolstad’s 2 (great) books, but clearly not on par with you!

Rob J. Goedman goedman@icloud.com

On Oct 22, 2014, at 6:38 PM, Brian J Smith notifications@github.com wrote:

Here's what I'm thinking...

The general algorithm implemented in the package does block sampling of parameters from their full conditional distributions. For a block that is all discrete (finite support) parameters, one can determine all possible values that the parameters could take on. Mamba provides a function that could then be used to compute the unnormalized pdf at each value. Divide the unnormalized pdfs by their sum and you have the exact full conditional probabilities. That allows for a draw to be generated directly (Gibbs sampling) from the full conditional for the block.

A performance consideration is that the pdfs will need to be computed at each value at each iteration of the overall MCMC algorithm. The number of evaluations could be very large if there are a large number of discrete parameters in a block.

That's the first thought that popped into my head anyway. Would love to hear any other ideas that you might have.

— Reply to this email directly or view it on GitHub.

[binarybana]

@goedman, This is definitely off-topic here, but I've actually implemented another of Dr. Liang's MCMC techniques in Julia here (SAMC).

I evaluated adding this method to Mamba, but it's sufficiently different from most MCMC techniques that (at the time) it made more sense to implement it by itself (along with MH, and AMWG for comparison (the later being inspired by Mamba's implementation!)).

[brian-j-smith]

@binarybana, nice job on the SAMC julia package.

[goedman]

@binarybana Hi Jason, thanks for your update/contribution.

I tried to get your README example to work but ran into a couple of issues. As I am working on OS X,

I can't have 2 files named samc.jl and SAMC.jl. I renamed samc.jl to same_record.jl and updated SAMC.jl accordingly (the 'include' statement).

I stored the provided README example in a little script 'readme.jl' and tried to run it:

''' julia> include("/Users/rob/Projects/Julia/Rob/SAMC/Examples/readme.jl") ERROR: MySampler has no method matching MySampler(::Float64, ::Float64) in include at /Applications/Julia-0.3.4.app/Contents/Resources/julia/lib/julia/sys.dylib in include_from_node1 at /Applications/Julia-0.3.4.app/Contents/Resources/julia/lib/julia/sys.dylib while loading /Users/rob/Projects/Julia/Rob/SAMC/Examples/readme.jl, in expression starting on line 34 '''

Shouldn't the data matrix being added as an argument?

I tried giving it a dummy matrix, like rand(5,5), that makes it fail 2 lines further down:

''' julia> mh MHRecord(MySampler(0.0,0.0,5x5 Array{Float64,2}: 0.567589 0.858875 0.604416 0.0232634 0.502679 0.0762462 0.573267 0.515966 0.145528 0.907228 0.479726 0.910614 0.0764755 0.99086 0.792443 0.0435808 0.166115 0.0707901 0.933435 0.966226 0.168827 0.872621 0.556068 0.993351 0.611206),MySampler(0.0,0.0,5x5 Array{Float64,2}: 0.567589 0.858875 0.604416 0.0232634 0.502679 0.0762462 0.573267 0.515966 0.145528 0.907228 0.479726 0.910614 0.0764755 0.99086 0.792443 0.0435808 0.166115 0.0707901 0.933435 0.966226 0.168827 0.872621 0.556068 0.993351 0.611206),Inf,{},0,0,1,0,1,Dict{Int64,Int64}(),Dict{Int64,Int64}())

julia> include("/Users/rob/Projects/Julia/Rob/SAMC/Examples/readme.jl") ERROR: unrecognized keyword argument "burn" in include at /Applications/Julia-0.3.4.app/Contents/Resources/julia/lib/julia/sys.dylib in include_from_node1 at /Applications/Julia-0.3.4.app/Contents/Resources/julia/lib/julia/sys.dylib while loading /Users/rob/Projects/Julia/Rob/SAMC/Examples/readme.jl, in expression starting on line 36 '''

I'm probably on a completely wrong track here :-)

[binarybana]

@goedman , thanks for the bug report! All those are mistakes on my part. I've since cleaned up the API, and made the README and tests the same to make sure they run properly (although proper doctest/literate programming support would make my copy/pasting unnecessary but I digress). So please try again and file bug reports/feature requests on the SAMC tracker. You can now wear your "First official user of the SAMC.jl" badge with pride!

And sorry for the repeated hijacking of this thread!

[goedman]

Brian,

Looks like I totally missed this, but DGS() looks like what you described above?

The Eyes example works fine, but I have been trying to warp it into below (Jags) example from Wagenmaker's BCM book:

model{

Observed Returns

for (i in 1:m){ k[i] ~ dbin(theta,n) }

Priors on Rate Theta and Number n

theta ~ dbeta(1,1) n ~ dcat(p[]) for (i in 1:nmax){ p[i] <- 1/nmax } }

I believe that should be possible with the DGS() sampler for :p, but I can't get it to work. Any chance, if this is possible in Mamba, you can give me a little push in the right direction?

[brian-j-smith]

@goedman: Yes, DGS is the discrete sampler describe in my earlier post (10/22). A quick tip for implementing the Jags example... the stochastic nodes are k (the response), theta, and n. Use DGS to sample n and one of the continuous samplers for theta.

[goedman]

@brian-j-smith Thanks a lot Brian!

[brian-j-smith]

Long-Term Update:

In addition to the Discrete Gibbs Sampler, several binary samplers have been implemented recently in Mamba by @bdeonovic, including:

• Binary Hamiltonian Monte Carlo (BHMC)
• Binary Metropolised Gibbs (BMG)
• Binary Modified Metropolised Gibbs (BMMG)

[ghost]

@brian-j-smith Great news! Thank you for the update.