Merge 456b0dd into c014311

pymc3/distributions/discrete.py

@@ -389,7 +389,7 @@ def _random(self, lower, upper, size=None):
         # as array-like.
         samples = stats.uniform.rvs(lower, upper - lower - np.finfo(float).eps,
                                     size=size)
-        return np.floor(samples).astype('int32')
+        return np.floor(samples).astype('int64')
 
     def random(self, point=None, size=None, repeat=None):
         lower, upper = draw_values([self.lower, self.upper], point=point)

pymc3/examples/gelman_bioassay.py

@@ -1,5 +1,5 @@
 import pymc3 as pm
-from numpy import ones, array
+from numpy import ones, array, random
 
 # Samples for each dose level
 n = 5 * ones(4, dtype=int)
@@ -18,14 +18,13 @@
     # Data likelihood
     deaths = pm.Binomial('deaths', n=n, p=theta, observed=[0, 1, 3, 5])
 
-    step = pm.NUTS()
-
-
 def run(n=1000):
     if n == "short":
         n = 50
     with model:
-        trace = pm.sample(n, step)
+        random.seed(42)
+        trace = pm.sample(n, init='random')
+        pm.summary(trace, varnames=['alpha', 'beta'])
 
 if __name__ == '__main__':
     run()

pymc3/sampling.py

@@ -108,6 +108,7 @@ def sample(draws, step=None, init='advi', n_init=200000, start=None,
         * NUTS : Run NUTS to estimate starting points and covariance matrix. If
         njobs > 1 it will sample starting points from the estimated posterior,
         otherwise it will use the estimated posterior mean.
+        * random : Draw random values from priors
         * None : Do not initialize.
     n_init : int
         Number of iterations of initializer
@@ -158,6 +159,10 @@ def sample(draws, step=None, init='advi', n_init=200000, start=None,
             start = start_
     else:
         step = assign_step_methods(model, step)
+        if init=='random':
+            start = [random_init(model) for _ in range(njobs)]
+            if njobs == 1:
+                start = start[0]
 
     if njobs is None:
         import multiprocessing as mp
@@ -402,7 +407,19 @@ def sample_ppc(trace, samples=None, model=None, vars=None, size=None, random_see
                                                          size=size))
 
     return {k: np.asarray(v) for k, v in ppc.items()}
-
+
+def random_init(model):
+
+    def trans_sample(v): 
+        return v.distribution.transform_used.backward(v.distribution.dist.random()).eval()
+
+    var_dict = {v.name:v.distribution.random() for v in model.vars 
+                                if not v.name.endswith('_')}
+    trans_var_dict = {v.name: trans_sample(v) for v in model.vars 
+                                if v.name.endswith('_')}
+    var_dict.update(trans_var_dict)
+
+    return var_dict
 
 def init_nuts(init='ADVI', njobs=1, n_init=500000, model=None,
               random_seed=-1, **kwargs):
@@ -421,6 +438,7 @@ def init_nuts(init='ADVI', njobs=1, n_init=500000, model=None,
         * ADVI_MAP: Initialize ADVI with MAP and use MAP as starting point.
         * MAP : Use the MAP as starting point.
         * NUTS : Run NUTS and estimate posterior mean and covariance matrix.
+        * random : Draw random values from prior.
     njobs : int
         Number of parallel jobs to start.
     n_init : int
@@ -472,6 +490,11 @@ def init_nuts(init='ADVI', njobs=1, n_init=500000, model=None,
         start = np.random.choice(init_trace, njobs)
         if njobs == 1:
             start = start[0]
+    elif init == 'random':    
+        start = [random_init(model) for _ in range(njobs)]
+        if njobs == 1:
+            start = start[0]
+        cov = None
     else:
         raise NotImplemented('Initializer {} is not supported.'.format(init))