galen2.jl

using Omega

# Add uncertainty to this if you want
region_size = 30_000 

# The number of groups
ngroups = ciid(ω -> Int(floor(uniform(ω, 1, region_size))))

# Produced random Boolean correlation matrix
function correlation_model(rng, nprobes = 10)
  # Sample ids without replacement
  neuron_ids = Int[]

  while length(neuron_ids) < nprobes
    id = Int(floor(uniform(rng, 1, region_size)))
    if id ∉ neuron_ids
      push!(neuron_ids, id)
    end
  end

  # Split the ids [1, 2, 3] into groups: [[1,2], ]
  region_id(i) = Int(floor(nprobes*i/region_size))
  issamegroup(i, j) = region_id(i) == region_id(j)

  correlated = [issamegroup(i,j) for i in neuron_ids, j in neuron_ids]

  # Return summary
  count(correlated)
end

# THis is just omega stuff to turn function into random variable
correlation_model_ = ciid(correlation_model)

rand(correlation_model_)

data =  [  true  false  false  false  false  false  false  false  false  false
           false   true  false  false  false  false   true  false  false  false
           false  false   true  false  false  false  false  false  false  false
           false  false  false   true  false  false  false  false   true  false
           false  false  false  false   true  false  false  false  false  false
           false  false  false  false  false   true  false  false  false  false
           false   true  false  false  false  false   true  false  false  false
           false  false  false  false  false  false  false   true  false  false
           false  false  false   true  false  false  false  false   true  false
           false  false  false  false  false  false  false  false  false   true]

# Draw condition samples using replica exchange MCMC 
samples = rand(ngroups, correlation_model_ ==ₛ count(data), 1000; alg = Replica)

using UnicodePlots

# Plot a histogram, take only last 400 samples for burn in
histogram(samples[500:end])