import matplotlib
matplotlib.use('TKAgg', warn=False, force=True)
import pymc3 as pm
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy as sp
from scipy import stats
import seaborn as sbn
#%matplotlib inline
import sys
print 'version: {}'.format(repr(sys.version_info))
from subprocess import Popen, PIPE, STDOUT
p = Popen('pip freeze', shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True)
output = p.stdout.read()
print output


gene_data = {
    'gene1': {
        'class1': [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 10, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
        'class2': [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 1, 1, 2, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0]
    },
    'gene2': {
        'class1': [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 2, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        'class2': [3, 1, 0, 1, 0, 0, 0, 1, 1, 2, 1, 0, 0, 0, 2, 0, 0, 2, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 2, 0, 0, 0, 0, 3, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 2, 1, 1, 0, 1, 1, 2, 1, 1, 1, 0, 1, 2, 0, 0, 1, 2, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 2, 0, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 2, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 2, 0, 2, 0, 0, 0, 1]
    },
    'gene3': {
        'class1': [13, 3, 9, 12, 24, 6, 22, 14, 5, 14, 18, 18, 7, 9, 13, 20, 9, 24, 12, 14, 4, 9, 5, 13, 5, 7, 16, 4, 13, 15, 2, 15, 8, 4, 4, 5, 8, 12, 2, 18, 11, 7, 7, 12, 25, 9, 8, 12, 2, 2, 14, 1, 2, 7, 9, 21, 7, 7, 8, 6, 5, 19, 36, 9, 6, 12, 8, 7, 3, 7, 7, 7, 6, 9, 7, 0, 12, 22, 13, 7, 5, 17, 3, 9, 8, 9, 3, 6, 20, 8, 15, 7, 8, 9, 13, 5, 4, 13, 6, 4, 10, 9, 6, 18, 13, 16, 16, 8, 10, 18, 8, 11, 5, 13, 22, 13, 5, 4, 3, 7, 6, 6, 8, 27, 9, 20, 5, 9, 12, 8, 19, 11, 8, 8, 23, 4, 4, 6, 9, 8, 4, 14, 6, 14, 2, 12, 3, 10, 11, 6, 5, 17, 20, 5, 4, 5, 1, 14, 16, 8, 13, 10, 12, 9, 3, 28, 24, 8, 5, 5, 2, 6, 19, 4, 15, 3, 4, 2, 9, 8, 4, 13, 11, 30, 10, 2, 8, 8, 11, 17, 17, 14, 6, 6, 3, 24, 7, 22, 9, 33, 6, 9, 12, 22, 10, 13, 12, 31, 4, 3, 11, 1, 2, 9, 7, 0, 2, 10, 7, 4, 9, 5, 9, 7, 18, 15, 5, 14, 8, 2, 13, 4, 7, 3, 14, 9, 6, 7, 6, 15, 4, 14, 4, 2, 8, 21, 10, 9, 12, 7, 3, 12, 3, 5, 4, 13, 4, 14, 48, 6, 7, 7, 9, 4, 8, 8, 8, 15, 5, 7, 8, 10, 9, 3, 6, 10, 6, 5, 8, 1, 2, 5, 7, 8, 14, 2, 11, 5, 7, 7, 7, 10, 17, 7, 7, 8, 12, 8, 9, 2, 11, 14, 7, 8, 9, 10, 15, 5, 6, 4, 16, 9, 19, 14, 12, 16, 6, 10, 7, 8, 11, 7, 4, 9, 3, 12, 5, 2, 32, 17, 6, 6, 12, 7, 15, 4, 7, 10, 4, 5, 8, 6, 6, 3, 10, 11, 14, 22, 11, 9, 9, 11, 18, 10, 3, 9, 11, 6, 9, 6, 5, 14, 6, 9, 8, 8, 13, 12, 19, 13, 8, 7, 15, 4, 17, 15, 10, 3, 5, 9, 14, 4, 7, 8, 11, 16, 6, 7, 4, 9, 12, 7, 5, 8, 17, 7, 1, 1, 11, 16],
        'class2': [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    },
    'gene4': {
        'class1': [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        'class2': [1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 2, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 2, 0, 0, 2, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0]
    }
}

def beta_params(data):
    params = sp.stats.beta.fit(data, loc=0, scale=1)
    return params[0], params[1]

# grab the per-gene beta parameters
beta_pars = dict()
for gene in gene_data.keys():
    c1_vec = np.array(gene_data[gene]['class1'])/5000.
    c2_vec = np.array(gene_data[gene]['class2'])/5000.
    beta_pars[gene] = beta_params(np.concatenate((c1_vec, c2_vec)))
    

laplace_FC = 0.25
nbinom_n = 25.
nbinom_mu = 5000.
nbinom_p = 1.-nbinom_mu/(nbinom_mu + nbinom_n)
nbinom_alpha = nbinom_n

with pm.Model() as model:
    # shared priors -- read depth per class
    c1_reads_prior = pm.distributions.discrete.NegativeBinomial('c1_num_reads', mu=nbinom_mu, alpha=nbinom_alpha)
    c2_reads_prior = pm.distributions.discrete.NegativeBinomial('c2_num_reads', mu=nbinom_mu, alpha=nbinom_alpha)
    # per-gene priors -- logFC, alpha, beta
    gene_priors, gene_pnames, gene_obs, gene_onames = [], [], [], []
    for gene in gene_data.keys():
        lfcn, c2_an, c2_bn = '{}_log_FC'.format(gene), '{}_c2_alpha'.format(gene), '{}_c2_beta'.format(gene)
        lfc_prior = pm.distributions.continuous.Laplace(lfcn, mu=0, b=laplace_FC)
        c2_a_prior = pm.distributions.continuous.Gamma(c2_an, alpha=(beta_pars[gene][0] ** 2)/1e-2, beta=beta_pars[gene][0]/1e-2)
        c2_b_prior = pm.distributions.continuous.Gamma(c2_bn, alpha=(beta_pars[gene][1] ** 2)/2, beta=beta_pars[gene][1]/2)
        # per-gene intermediates -- rates
        c2_pn, c1_pn = '{}_c2_prob'.format(gene), '{}_c1_prob'.format(gene)
        c2_proba = pm.distributions.continuous.Beta(c2_pn, c2_a_prior, c2_b_prior)
        c1_proba = pm.Deterministic(c1_pn, c2_proba * np.exp(lfc_prior))
        gene_priors.append((lfc_prior, c2_a_prior, c2_b_prior, c2_proba, c1_proba))
        gene_pnames.append((lfcn, c2_an, c2_bn, c2_pn, c1_pn))
        # per-gene observations
        c1_cn, c2_cn = '{}_c2_counts'.format(gene), '{}_c1_counts'.format(gene)
        class2_counts_obs = pm.distributions.discrete.Binomial(c2_cn, n=c2_reads_prior, p=c2_proba, observed=gene_data[gene]['class2'])
        class1_counts_obs = pm.distributions.discrete.Binomial(c1_cn, n=c1_reads_prior, p=c1_proba, observed=gene_data[gene]['class1'])
        gene_obs.append((class1_counts_obs, class2_counts_obs))
        gene_onames.append((c1_cn, c2_cn))
    for RV in model.basic_RVs:
        print(RV.name, RV.logp(model.test_point), model.test_point[RV.name] if RV.name in model.test_point else 'NA')
        
    trace = pm.sample(8000, progressbar=False, tune=5000, step=pm.SMC(), seed=314159, chains=200)
   

# STDOUT
"""
version: sys.version_info(major=2, minor=7, micro=14, releaselevel='final', serial=0)
You are using pip version 9.0.3, however version 18.1 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.
absl-py==0.1.10
alabaster==0.7.10
anaconda-client==1.6.5
anaconda-navigator==1.6.9
anaconda-project==0.8.0
appnope==0.1.0
appscript==1.0.1
asn1crypto==0.22.0
astroid==1.5.3
astropy==2.0.2
Babel==2.5.0
backports-abc==0.5
backports.functools-lru-cache==1.5
backports.shutil-get-terminal-size==1.0.0
backports.ssl-match-hostname==3.5.0.1
backports.weakref==1.0rc1
beautifulsoup4==4.6.0
bitarray==0.8.1
blaze==0.11.3
bleach==1.5.0
bokeh==0.12.10
boto==2.48.0
bottle==0.12.13
Bottleneck==1.2.1
BrainGenesetAnnotation==0.0.0
cdecimal==2.3
certifi==2018.1.18
cffi==1.10.0
chardet==3.0.4
click==6.7
cloudpickle==0.5.2
clyent==1.2.2
colorama==0.3.9
conda==4.3.34
configparser==3.5.0
contextlib2==0.5.5
CProfileV==1.0.7
cryptography==2.0.3
ctc-model==0.5
cycler==0.10.0
Cython==0.26.1
cytoolz==0.8.2
dask==0.16.0
datashape==0.5.4
decorator==4.1.2
distributed==1.20.0
docutils==0.14
editdistance==0.4
entrypoints==0.2.3
enum34==1.1.6
et-xmlfile==1.0.1
fastcache==1.0.2
Flask==0.12.2
Flask-Cors==3.0.3
funcsigs==1.0.2
functools32==3.2.3.post2
futures==3.1.1
gevent==1.2.2
gmpy2==2.0.8
google-images-download==1.1.1
greenlet==0.4.12
grin==1.2.1
h5py==2.7.1
heapdict==1.0.0
html5lib==0.9999999
idna==2.6
imageio==2.2.0
imagesize==0.7.1
ipaddress==1.0.18
ipykernel==4.6.1
ipython==5.4.1
ipython-genutils==0.2.0
ipywidgets==7.0.5
isort==4.2.15
itsdangerous==0.24
jdcal==1.3
jedi==0.10.2
Jinja2==2.9.6
joblib==0.12.5
jsonschema==2.6.0
jupyter-client==5.1.0
jupyter-console==5.2.0
jupyter-core==4.4.0
Keras==2.1.4
kiwisolver==1.0.1
lazy-object-proxy==1.3.1
llvmlite==0.20.0
locket==0.2.0
loompy==2.0.8
lxml==4.1.0
Mako==1.0.7
Markdown==2.6.11
MarkupSafe==1.0
matplotlib==2.2.3
mccabe==0.6.1
mistune==0.8.1
mock==2.0.0
mpmath==0.19
msgpack-python==0.4.8
multipledispatch==0.4.9
natsort==5.2.0
navigator-updater==0.1.0
nbconvert==5.3.1
nbformat==4.4.0
networkx==1.11
nltk==3.2.5
nose==1.3.7
notebook==5.2.1
numba==0.35.0+0.g24edca2.dirty
numexpr==2.6.2
numpy==1.15.2
numpydoc==0.7.0
odo==0.5.1
olefile==0.44
openpyxl==2.4.8
packaging==16.8
pandas==0.23.4
pandocfilters==1.4.2
partd==0.3.8
path.py==10.3.1
pathlib2==2.3.0
patsy==0.4.1
pbr==3.1.1
pep8==1.7.0
pexpect==4.2.1
pickleshare==0.7.4
Pillow==4.2.1
ply==3.10
prompt-toolkit==1.0.15
protobuf==3.4.1
psutil==5.4.0
ptyprocess==0.5.2
py==1.5.2
PyAudio==0.2.7
pycodestyle==2.3.1
pycosat==0.6.3
pycparser==2.18
pycrypto==2.6.1
pycurl==7.43.0
pyflakes==1.6.0
Pygments==2.2.0
pygpu==0.6.9
pylint==1.7.4
pymc3==3.5
pynipt==2.0.1
pyodbc==4.0.21
pyOpenSSL==17.2.0
pyparsing==2.2.1
pysam==0.14.1
PySocks==1.6.7
pytest==3.2.5
python-dateutil==2.7.3
pytz==2018.5
PyWavelets==0.5.2
PyYAML==3.12
pyzmq==16.0.3
QtAwesome==0.4.4
qtconsole==4.3.1
QtPy==1.3.1
requests==2.18.4
rope==0.10.7
ruamel-yaml==0.11.14
scandir==1.6
scikit-image==0.13.1
scikit-learn==0.19.1
scipy==1.1.0
seaborn==0.9.0
selenium==3.11.0
simplegeneric==0.8.1
singledispatch==3.4.0.3
six==1.11.0
snowballstemmer==1.2.1
sortedcollections==0.5.3
sortedcontainers==1.5.7
Sphinx==1.6.3
sphinxcontrib-websupport==1.0.1
spyder==3.2.4
SQLAlchemy==1.1.13
statsmodels==0.8.0
subprocess32==3.5.2
sympy==1.1.1
tables==3.4.2
tblib==1.3.2
tensorflow==1.5.0
tensorflow-tensorboard==1.5.1
terminado==0.6
testpath==0.3.1
Theano==1.0.3
toolz==0.8.2
tornado==4.5.2
tqdm==4.28.1
traitlets==4.3.2
typing==3.6.2
unicodecsv==0.14.1
urllib3==1.22
wcwidth==0.1.7
webencodings==0.5.1
Werkzeug==0.12.2
widgetsnbextension==3.0.8
wrapt==1.10.11
xlrd==1.1.0
XlsxWriter==1.0.2
xlwings==0.11.4
xlwt==1.2.0
zict==0.1.3

('c1_num_reads', array(-7.832604), array(5000))
('c2_num_reads', array(-7.832604), array(5000))
('gene1_log_FC', array(0.69314718), array(0.))
('gene1_c2_alpha_log__', array(-0.73011356), array(-2.06277612))
('gene1_c2_beta_log__', array(4.79035285), array(6.05586589))
('gene1_c2_prob_logodds__', array(-2.3913697), array(-8.11864201))
('gene2_log_FC', array(0.69314718), array(0.))
('gene2_c2_alpha_log__', array(-0.60908696), array(-1.95173937))
('gene2_c2_beta_log__', array(4.32940381), array(5.59491823))
('gene2_c2_prob_logodds__', array(-2.30479137), array(-7.54665761))
('gene3_log_FC', array(0.69314718), array(0.))
('gene3_c2_alpha_log__', array(1.25839589), array(-0.12418241))
('gene3_c2_beta_log__', array(2.57729241), array(3.84288108))
('gene3_c2_prob_logodds__', array(-1.08162627), array(-3.96706349))
('gene4_log_FC', array(0.69314718), array(0.))
('gene4_c2_alpha_log__', array(0.6264108), array(-0.75347523))
('gene4_c2_beta_log__', array(3.76310741), array(5.02862668))
('gene4_c2_prob_logodds__', array(-1.45897885), array(-5.7821019))
('gene1_c1_counts', array(-584.80881175), 'NA')
('gene1_c2_counts', array(-594.13276173), 'NA')
('gene2_c1_counts', array(-940.90331163), 'NA')
('gene2_c2_counts', array(-983.62540712), 'NA')
('gene3_c1_counts', array(-37490.37652495), 'NA')
('gene3_c2_counts', array(-26386.14570245), 'NA')
('gene4_c1_counts', array(-5982.7069653), 'NA')
('gene4_c2_counts', array(-6116.68814086), 'NA')
/usr/local/anaconda2/lib/python2.7/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters
/usr/local/anaconda2/lib/python2.7/site-packages/scipy/optimize/minpack.py:163: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last ten iterations.
  warnings.warn(msg, RuntimeWarning)
Adding model likelihood to RVs!
Init new trace!
Sample initial stage: ...
Beta: 0.000000 Stage: 0
Initializing chain traces ...
Sampling ...
Beta: 0.000028 Stage: 1
Initializing chain traces ...
Sampling ...
"""
# STDERR
"""
Traceback (most recent call last):
  File "foo.oob.py", line 84, in <module>
    trace = pm.sample(8000, progressbar=False, tune=5000, step=pm.SMC(), seed=314159, chains=200)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/sampling.py", line 340, in sample
    **kwargs)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/step_methods/smc.py", line 518, in sample_smc
    _iter_parallel_chains(**sample_args)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/step_methods/smc.py", line 694, in _iter_parallel_chains
    for _ in p:
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/backends/smc_text.py", line 58, in paripool
    yield function(work_item)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/step_methods/smc.py", line 655, in _work_chain
    return _sample(*work)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/step_methods/smc.py", line 603, in _sample
    for strace in sampling:
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/step_methods/smc.py", line 636, in _iter_sample
    point, out_list = step.step(point)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/backends/smc_text.py", line 103, in step
    apoint, alist = self.astep(self.bij.map(point))
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/step_methods/smc.py", line 212, in astep
    varlogp = self.check_bnd(q)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/theano/compile/function_module.py", line 917, in __call__
    storage_map=getattr(self.fn, 'storage_map', None))
  File "/usr/local/anaconda2/lib/python2.7/site-packages/theano/gof/link.py", line 325, in raise_with_op
    reraise(exc_type, exc_value, exc_trace)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/theano/compile/function_module.py", line 903, in __call__
    self.fn() if output_subset is None else\
IndexError: index out of bounds
Apply node that caused the error: Subtensor{int64}(inarray, ScalarFromTensor.0)
Toposort index: 26
Inputs types: [TensorType(float64, vector), Scalar(int64)]
Inputs shapes: [(2,), ()]
Inputs strides: [(8,), ()]
Inputs values: [array([5806, 3882]), 3]
Outputs clients: [[Elemwise{Composite{Switch(Cast{int8}(GE(Cast{int64}(i0), i1)), ((i2 + scalar_gammaln((i3 + Cast{int64}(i0))) + (i4 * Cast{int64}(i0))) - scalar_gammaln((i5 + Cast{int64}(i0)))), i6)}}(Subtensor{int64}.0, TensorConstant{0}, TensorConstant{-187.36735574088576}, TensorConstant{25.0}, TensorConstant{-0.0049875..7418644428}, TensorConstant{1}, TensorConstant{-inf})]]

Backtrace when the node is created(use Theano flag traceback.limit=N to make it longer):
  File "foo.oob.py", line 84, in <module>
    trace = pm.sample(8000, progressbar=False, tune=5000, step=pm.SMC(), seed=314159, chains=200)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/step_methods/smc.py", line 167, in __init__
    self.check_bnd = logp_forw([model.varlogpt], vars, shared)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/step_methods/smc.py", line 728, in logp_forw
    out_list, inarray0 = join_nonshared_inputs(out_vars, vars, shared)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/theanof.py", line 241, in join_nonshared_inputs
    for var, slc, shp, dtyp in ordering.vmap}
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/theanof.py", line 241, in <dictcomp>
    for var, slc, shp, dtyp in ordering.vmap}
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/theanof.py", line 254, in reshape_t
    return x[0]
  File "foo.oob.py", line 84, in <module>
    trace = pm.sample(8000, progressbar=False, tune=5000, step=pm.SMC(), seed=314159, chains=200)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/step_methods/smc.py", line 167, in __init__
    self.check_bnd = logp_forw([model.varlogpt], vars, shared)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/step_methods/smc.py", line 728, in logp_forw
    out_list, inarray0 = join_nonshared_inputs(out_vars, vars, shared)
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/theanof.py", line 241, in join_nonshared_inputs
    for var, slc, shp, dtyp in ordering.vmap}
  File "/usr/local/anaconda2/lib/python2.7/site-packages/pymc3/theanof.py", line 241, in <dictcomp>
    for var, slc, shp, dtyp in ordering.vmap}

HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.
"""