pymc - Drawing from a 2-D prior that is only available as samples in pymc2 -

i'm trying play around bayesian updating, , have situation in using posterior previous runs prior. 2d prior on alpha , beta, have traces, alphatrace , betatrace. stack them , use code adopted https://gist.github.com/jcrudy/5911624 make kde based stochastic.

#from https://gist.github.com/jcrudy/5911624 def kernelsmoothing(name, dataset, bw_method=none,  observed=false, value=none):     '''create pymc node distribution comes kernel smoothing density estimate.'''     density = gaussian_kde(dataset, bw_method)      def logp(value):         #print "val", value         d = density(value)         if d == 0.0:             return float('-inf')         return np.log(d)      def random():         result = none         sample=density.resample(1)         #print sample, sample.shape         result = sample[0][0],sample[1][0]         return result      if value == none:         value = random()      dtype = type(value)      result = pymc.stochastic(logp = logp,                              doc = 'a kernel smoothing density node.',                              name = name,                              parents = {},                              random = random,                              trace = true,                              value = none,                              dtype = dtype,                              observed = observed,                              cache_depth = 2,                              plot = true,                              verbose = 0)     return result 

note critical thing here obtain 2-values joint prior: why need 2-d prior , not 2 1-d priors.

the model so:

ctrace=np.vstack((alphatrace, betatrace)) cnew=kernelsmoothing("cnew", ctrace) @pymc.deterministic def alphanew(cnew=cnew, name='alphanew'):     return cnew[0] @pymc.deterministic def betanew(cnew=cnew, name='betanew'):     return cnew[1] newtheta=pymc.beta("newtheta", alphanew, betanew) newexp = pymc.binomial('newexp', n=[14], p=[newtheta], value=[4], observed=true) model3=pymc.model([cnew, alphanew, betanew, newtheta, newexp]) mcmc3=pymc.mcmc(model3) mcmc3.sample(20000,5000,5) 

in case wondering, 71st experiment in hierarchical rat tumor example in chapter 5 in gelman's bda. "prior" using posterior on alpha , beta after 70 experiments.

but, when sample, things blow error:

valueerror: maximum competence reported stochastic cnew <= 0... may need write custom step method class. 

its not cnew care updating stochastic, rather alphanew , betanew. how ought structuring code make error go away?

edit: initial model gave me posteriors wish use prior:

tumordata="""0 20  0 20  0 20  0 20  0 20  0 20  0 20  0 19  0 19  0 19  0 19  0 18  0 18  0 17  1 20  1 20  1 20  1 20  1 19  1 19  1 18  1 18  3 27  2 25  2 24  2 23  2 20  2 20  2 20  2 20  2 20  2 20  1 10  5 49  2 19  5 46  2 17  7 49  7 47  3 20  3 20  2 13  9 48  10 50  4 20  4 20  4 20  4 20  4 20  4 20  4 20  10 48  4 19  4 19  4 19  5 22  11 46  12 49  5 20  5 20  6 23  5 19  6 22  6 20  6 20  6 20  16 52  15 46  15 47  9 24  """ tumortuples=[e.strip().split() e in tumordata.split("\n")] tumory=np.array([np.int(e[0].strip()) e in tumortuples if len(e) > 0]) tumorn=np.array([np.int(e[1].strip()) e in tumortuples if len(e) > 0]) n = tumorn.shape[0]  mu = pymc.uniform("mu",0.00001,1., value=0.13) nu = pymc.uniform("nu",0.00001,1., value=0.01)  @pymc.deterministic def alpha(mu=mu, nu=nu, name='alpha'):     return mu/(nu*nu)  @pymc.deterministic def beta(mu=mu, nu=nu, name='beta'):     return (1.-mu)/(nu*nu)   thetas=pymc.container([pymc.beta("theta_%i" % i, alpha, beta) in range(n)]) deaths = pymc.binomial('deaths', n=tumorn, p=thetas, value=tumory, size=n, observed=true) 

i use joint-posterior model on alpha, beta input "new model" @ top. begs question if ought including theta1..theta70 in model @ top update along alpha , beta new data binomial n=14, y=4. cant little model prior 2d sample array working :-(

i found question since ran similar proble. according documentation of pymc.stepmethod.competence, problem none of built-in samplers handle dtype associated stochastic variable.

i not sure needs done resolve that. maybe 1 of sampler methods can extended handle special types?

hopefully more pymc mojo can shine light on needs done..

def competence(s):     """     function used sampler determine step method class     should used handle stochastic variables.      return value should competence     score 0 3, assigned follows:      0:  can't handle variable.     1:  can handle variable, i'm generalist ,         shouldn't top choice (metropolis         , friends fall category).     2:  i'm designed type of situation,         more specialized.     3:  made situation, let me handle variable.      in order eligible inclusion in registry, sampling     method's init method must work single argument,     stochastic object.      if want exclude particular step method     consideration handling variable, this:      competence functions must called 'competence' , decorated     '@staticmethod' decorator. example:          @staticmethod         def competence(s):             if isinstance(s, mystochasticsubclass):                 return 2             else:                 return 0      :seealso: pick_best_methods, assign_method     """