Hierarchical QFA with pyMC

pyMC_Hierarchy.py

def hierarchy_inf(data,par,iter=250000,burn=1000,thin=100):
    priors={}
    x0=mc.Uniform('x0',par.x0_min,par.x0_max)
    tau=mc.Uniform('tau',par.tau_min,par.tau_max)

    priors["x0"]=x0
    priors["tau"]=tau
  
    r=mc.Uniform('r',par.r_min,par.r_max)
    r_delta=mc.Uniform('r_delta',0,(par.r_max-par.r_min)/2)
    
    K=mc.Uniform('K',par.K_min,par.K_max)
    K_delta=mc.Uniform('K_delta',0,(par.K_max-par.K_min)/2)
    
    priors["r"]=r
    priors["r_delta"]=r_delta
    priors["K"]=K
    priors["K_delta"]=K_delta

    grps=data.groupby("Gene")
    for grp in grps:
        grplab,reps=grp

        r_gen=mc.Uniform("r_{}".format(grplab),r-r_delta,r+r_delta)
        r_gen_delta=mc.Uniform("r_delta_{}".format(grplab),0,r_delta)
        K_gen=mc.Uniform("K_{}".format(grplab),K-K_delta,K+K_delta)
        K_gen_delta=mc.Uniform("K_delta_{}".format(grplab),0,K_delta)

        priors["r_{}".format(grplab)]=r_gen
        priors["r_delta_{}".format(grplab)]=r_gen_delta
        priors["K_{}".format(grplab)]=K_gen
        priors["K_delta_{}".format(grplab)]=K_gen_delta
        
        reps=reps.groupby("ID")
        for rep in reps:
            replab,repdf=rep
            r_rep=mc.Uniform("r_{0}_{1}".format(grplab,replab),r_gen-r_gen_delta,r_gen+r_gen_delta)
            K_rep=mc.Uniform("K_{0}_{1}".format(grplab,replab),K_gen-K_gen_delta,K_gen+K_gen_delta)
            @mc.deterministic(plot=False)
            def logisticobs(x0=x0,r=r_rep,K=K_rep):
                return(logistic(x0,r,K,repdf.ExptTime))
            obs=mc.Normal('obs_{0}_{1}'.format(grplab,replab),mu=logisticobs,tau=tau,value=repdf.Intensity,observed=True)
            priors["r_{0}_{1}".format(grplab,replab)]=r_rep
            priors["K_{0}_{1}".format(grplab,replab)]=K_rep
            priors['obs_{0}_{1}'.format(grplab,replab)]=obs
    M=mc.MCMC(priors)
    M.sample(iter=iter, burn=burn, thin=thin,progress_bar=False)
    return(M)