Advertisement
Guest User

pymc3 slow snippet

a guest
Jun 9th, 2017
91
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.     with pm.Model() as model_context:
  2.         print "Setting hyper priors..."
  3.         H1_mu = pm.Normal('h1_mu', numpy.array([100, 1]), numpy.array([1.0/1000, 1.0/1000]), shape=2)
  4.         H1_precision = pm.Gamma('h1_precision', numpy.array([0.1, 0.1]), numpy.array([0.1, 0.1]), shape=2)
  5.  
  6.         H1_alpha = pm.Gamma('h1_a', numpy.array([0.1, 0.1]), numpy.array([0.1, 0.1]), shape=2)
  7.         H1_beta = pm.Gamma('h1_b', numpy.array([0.1, 0.1]), numpy.array([0.1, 0.1]), shape=2)
  8.         w1 = pm.Beta('w1', 15, 85)
  9.  
  10.         print "Setting priors and likelihood..."
  11.         Z = numpy.empty(nGenes, dtype=object)
  12.         MU = numpy.empty(nGenes, dtype=object)
  13.         ALPHA = numpy.empty(nGenes, dtype=object)
  14.         for i,gene in enumerate(fullG):
  15.             Z[i] = pm.Categorical('Z_%d' % i, tt.stack([1.0-w1,w1]))
  16.             MU[i] = pm.Normal('mu_%d' % i, H1_mu[Z[i]], H1_precision[Z[i]])
  17.             ALPHA[i] = pm.Gamma('alpha_%d' % i, alpha=H1_alpha[Z[i]], beta=H1_beta[Z[i]])
  18.             COUNTS[i] = pm.NegativeBinomial('counts_%d' % i, mu=MU[i], alpha=ALPHA[i], observed=gene.reads.flatten())
  19.  
  20.     ### ABOVE IS SLOW
  21.  
  22.     ### BELOW IS EVEN SLOWER ###
  23.     print "Performing sampling..."
  24.     with model_context:
  25.         #start = pm.find_MAP()
  26.         #step1 = pm.Metropolis(vars=[COUNTS, MU, ALPHA, H1_mu, H1_precision, H1_alpha, H1_beta, w1])
  27.         #step2 = pm.ElemwiseCategoricalStep(Z, np.arange(2))
  28.  
  29.         db = Text('%s_negbinom_pymc3.trace' % label)
  30.         trace = pm.sample(SAMPLES, njobs=2, random_seed=123, trace=db)
Advertisement
RAW Paste Data Copied
Advertisement