def test_fit(self):
        p = self._build_parent()
        s = MyStochastic(self.STOCHASTIC_NAME, p)

        mcmc = MCMC({p, s})

        mcmc.sample(100, burn=10, thin=2)

 def test_nd(self):
     M = MCMC([self.NDstoch()], db=self.name, dbname=os.path.join(testdir, 'ND.'+self.name), dbmode='w')
     M.sample(10, progress_bar=0)
     a = M.trace('nd')[:]
     assert_equal(a.shape, (10,2,2))
     db = getattr(pymc.database, self.name).load(os.path.join(testdir, 'ND.'+self.name))
     assert_equal(db.trace('nd')[:], a)

    def test_simple(self):

        # Priors
        mu = Normal('mu', mu=0, tau=0.0001)
        s = Uniform('s', lower=0, upper=100, value=10)
        tau = s ** -2

        # Likelihood with missing data
        x = Normal('x', mu=mu, tau=tau, value=m, observed=True)

        # Instantiate sampler
        M = MCMC([mu, s, tau, x])

        # Run sampler
        M.sample(10000, 5000, progress_bar=0)

        # Check length of value
        assert_equal(len(x.value), 100)
        # Check size of trace
        tr = M.trace('x')()
        assert_equal(shape(tr), (5000, 2))

        sd2 = [-2 < i < 2 for i in ravel(tr)]

        # Check for standard normal output
        assert_almost_equal(sum(sd2) / 10000., 0.95, decimal=1)

def test_interactive():
    if 'sqlite' not in dir(pymc.database):
        raise nose.SkipTest
    M=MCMC(disaster_model,db='sqlite',
           dbname=os.path.join(testdir, 'interactiveDisaster.sqlite'),
           dbmode='w')
    M.isample(10, out=open('testresults/interactivesqlite.log', 'w'), progress_bar=0)

    def test_fit_with_sibling(self):
        p = self._build_parent()
        s = MyStochastic(self.STOCHASTIC_NAME, p)
        sib = MyStochastic(self.SIBLING_NAME, p)

        mcmc = MCMC({p, s, sib})

        mcmc.sample(100, burn=10, thin=2)

    def test_pymc_model(self):
        """ Tests sampler """

        sampler = MCMC(model_omm.pymc_parameters)
        self.assert_(isinstance(model_omm, TorsionFitModelOMM))
        self.assert_(isinstance(sampler, pymc.MCMC))

        sampler.sample(iter=1)

def mcmc(prob, nsample=100, modulename = 'model' ):
    try:
        mystr = "from " + modulename + " import model"
        exec(mystr)
    except:
        print 'cannot import', modulename
    M = MCMC( model(prob) )
    M.sample(nsample)
    return M

def test_interactive():
    S = MCMC(disaster_model)
    S.isample(
        200,
        100,
        2,
        out=open(
            'testresults/interactive.log',
            'w'),
        progress_bar=0)

 def test_zcompression(self):
     db = pymc.database.hdf5.Database(dbname=os.path.join(testdir, 'DisasterModelCompressed.hdf5'),
                                      dbmode='w',
                                      dbcomplevel=5)
     S = MCMC(DisasterModel, db=db)
     S.sample(45,10,1)
     assert_array_equal(S.e.trace().shape, (35,))
     S.db.close()
     db.close()
     del S

 def test_zcompression(self):
     with warnings.catch_warnings():
         warnings.simplefilter('ignore')
         db = pymc.database.hdf5.Database(dbname=os.path.join(testdir, 'disaster_modelCompressed.hdf5'),
                                          dbmode='w',
                                          dbcomplevel=5)                                 
         S = MCMC(disaster_model, db=db)
         S.sample(45,10,1, progress_bar=0)
         assert_array_equal(S.trace('early_mean')[:].shape, (35,))
         S.db.close()
         db.close()
         del S

    def MCMC( self, nruns=10000, burn=1000, init_error_std=1., max_error_std=100., verbose=1 ):
        ''' Perform Markov Chain Monte Carlo sampling using pymc package

            :param nruns: Number of MCMC iterations (samples)
            :type nruns: int
            :param burn: Number of initial samples to burn (discard)
            :type burn: int
            :param verbose: verbosity of output
            :type verbose: int
            :param init_error_std: Initial standard deviation of residuals
            :type init_error_std: fl64
            :param max_error_std: Maximum standard deviation of residuals that will be considered
            :type max_error_std: fl64
            :returns: pymc MCMC object
        '''
        if max_error_std < init_error_std:
            print "Error: max_error_std must be greater than or equal to init_error_std"
            return
        try:
            from pymc import Uniform, deterministic, Normal, MCMC, Matplot
        except ImportError as exc:
            sys.stderr.write("Warning: failed to import pymc module. ({})\n".format(exc))
            sys.stderr.write("If pymc is not installed, try installing:\n")
            sys.stderr.write("e.g. try using easy_install: easy_install pymc\n")
        def __mcmc_model( self, init_error_std=1., max_error_std=100. ):
            #priors
            variables = []
            sig = Uniform('error_std', 0.0, max_error_std, value=init_error_std)
            variables.append( sig )
            for nm,mn,mx in zip(self.parnames,self.parmins,self.parmaxs):
                evalstr = "Uniform( '" + str(nm) + "', " +  str(mn) + ", " +  str(mx) + ")"
                variables.append( eval(evalstr) )
            #model
            @deterministic()
            def residuals( pars = variables, p=self ):
                values = []
                for i in range(1,len(pars)):
                    values.append(float(pars[i]))
                pardict = dict(zip(p.parnames,values))
                p.forward(pardict=pardict, reuse_dirs=True)
                return numpy.array(p.residuals)*numpy.array(p.obsweights)
            #likelihood
            y = Normal('y', mu=residuals, tau=1.0/sig**2, observed=True, value=numpy.zeros(len(self.obs)))
            variables.append(y)
            return variables

        M = MCMC( __mcmc_model(self, init_error_std=init_error_std, max_error_std=max_error_std) )
        M.sample(iter=nruns,burn=burn,verbose=verbose)
        return M

def analizeMwm():
	masked_values = np.ma.masked_equal(x, value=None)
	print("m v: ", masked_values)

	print("dmwm da: ", dmwm.disasters_array)

	Mwm = MCMC(dmwm)
	Mwm.sample(iter=10000, burn=1000, thin=10)

	print("Mwm t: ", Mwm.trace('switchpoint')[:])

	hist(Mwm.trace('late_mean')[:])
	# show()

	plot(Mwm)

    def test_zcompression(self):

        original_filters = warnings.filters[:]
        warnings.simplefilter("ignore")
        try:
            db = pymc.database.hdf5.Database(dbname=os.path.join(testdir, 'disaster_modelCompressed.hdf5'),
                                             dbmode='w',
                                             dbcomplevel=5)
            S = MCMC(disaster_model, db=db)
            S.sample(45,10,1, progress_bar=0)
            assert_array_equal(S.trace('early_mean')[:].shape, (35,))
            S.db.close()
            db.close()
            del S
        finally:
            warnings.filters = original_filters

def estimate_failures(samples, #samples from noisy labelers
                      n_samples=10000, #number of samples to run MCMC for
                      burn=None, #burn-in. Defaults to n_samples/2
                      thin=10, #thinning rate. Sample every k samples from markov chain 
                      alpha_p=1, beta_p=1, #beta parameters for true positive rate
                      alpha_e=1, beta_e=10 #beta parameters for noise rates
                      ):

  if burn is None:
    burn = n_samples / 2

  S,N = samples.shape
  p = Beta('p', alpha=alpha_p, beta=beta_p) #prior on true label
  l = Bernoulli('l', p=p, size=S)
  e_pos = Beta('e_pos', alpha_e, beta_e, size=N) # error rate if label = 1
  e_neg = Beta('e_neg', alpha_e, beta_e, size=N) # error rate if label = 0

  @deterministic(plot=False)
  def noise_rate(l=l, e_pos=e_pos, e_neg=e_neg):
    #probability that a noisy labeler puts a label 1
    return np.outer(l, 1-e_pos) + np.outer(1-l, e_neg)

  noisy_label = Bernoulli('noisy_label', p=noise_rate, size=samples.shape, value=samples, observed=True)
  variables = [l, e_pos, e_neg, p, noisy_label, noise_rate]
  model = MCMC(variables, verbose=3)
  model.sample(iter=n_samples, burn=burn, thin=thin)
  model.write_csv('out.csv', ['p', 'e_pos', 'e_neg'])
  p = np.median(model.trace('p')[:])
  e_pos = np.median(model.trace('e_pos')[:],0)
  e_neg = np.median(model.trace('e_neg')[:],0)
  return p, e_pos, e_neg

def bimodal_gauss(data,pm):
    '''run MCMC to get regression on bimodal normal distribution'''
    m1 = np.mean(data[pm])/2.
    m2 = np.mean(data[pm])*2.
    dm = m2 - m1
    size = len(data[pm])

    ### set up model
    p = Uniform( "p", 0.2 , 0.8) #this is the fraction that come from mean1 vs mean2
    # p = distributions.truncated_normal_like('p', mu=0.5, tau=0.001, a=0., b=1.)
    # p = Normal( 'p', mu=(1.*sum(comp0==1))/size, tau=1./0.1**2 ) # attention: wings!, tau = 1/sig^2
    # p = Normal( 'p', mu=0.5, tau=1./0.1**2 ) # attention: wings!, tau = 1/sig^2
    
    ber = Bernoulli( "ber", p = p, size = size) # produces 1 with proportion p
    precision = Gamma('precision', alpha=0.01, beta=0.01)
    
    dmu = Normal( 'dmu', dm, tau=1./0.05**2 ) # [PS] give difference between means, finite
    # dmu = Lognormal( 'dmu', 0.3, tau=1./0.1**2)
    
    mean1 = Normal( "mean1", mu = m1,          tau = 1./0.1**2 ) # better to use Normals versus Uniforms,
                                                                 # if not truncated
    mean2 = Normal( "mean2", mu = mean1 + dmu, tau = 1./0.1**2 ) # tau is 1/sig^2
    
    @deterministic
    def mean( ber = ber, mean1 = mean1, mean2 = mean2):
        return ber*mean1 + (1-ber)*mean2

    
    obs = Normal( "obs", mean, precision, value = data[pm], observed = True)
    model = Model( {"p":p, "precision": precision, "mean1": mean1, "mean2":mean2, "obs":obs} )
    
    from pymc import MCMC, Matplot



    M = MCMC(locals(), db='pickle', dbname='metals.pickle')
    iter = 3000; burn = 2000; thin = 10
    M.sample(iter=iter, burn=burn, thin=thin)
    M.db.commit()

    mu1 = np.mean(M.trace('mean1')[:])
    mu2 = np.mean(M.trace('mean2')[:])
    p   = np.mean(M.trace('p')[:])
    return p, mu1, 0.1, mu2, 0.1, M

    def bayesian_regression(self, Methodology):
        
        fit_dict                = OrderedDict()
        
        fit_dict['methodology'] = r'Inference $\chi^{2}$ model'
        
        #Initial guess for the fitting:
        Np_lsf                  = polyfit(self.x_array, self.y_array, 1)
        m_0, n_0                = Np_lsf[0], Np_lsf[1]
                
        MCMC_dict               = self.lr_ChiSq(self.x_array, self.y_array, m_0, n_0)
        
        myMCMC                  = MCMC(MCMC_dict)
        
        myMCMC.sample(iter=10000, burn=1000)

        fit_dict['m'], fit_dict['n'], fit_dict['m_error'], fit_dict['n_error'] = myMCMC.stats()['m']['mean'], myMCMC.stats()['n']['mean'], myMCMC.stats()['m']['standard deviation'], myMCMC.stats()['n']['standard deviation']
        
        return fit_dict

class LeagueModel(object):
    """MCMC model of a football league."""
    def __init__(self, fname):
        super(LeagueModel, self).__init__()
        league = fuba.League(fname)

        N = len(league.teams)
        #dummy future games
        future_games = [[league.teams["Werder Bremen"],league.teams["Dortmund"]]]

        self.goal_rate = np.empty(N,dtype=object)
        self.match_rate = np.empty(len(league.games)*2,dtype=object)
        self.match_goals_future = np.empty(len(future_games)*2,dtype=object)
        self.home_adv = Normal(name = 'home_adv',mu=0,tau=10.)

        for t in league.teams.values():
            print t.name,t.team_id
            self.goal_rate[t.team_id] = Exponential('goal_rate_%i'%t.team_id,beta=1)

        for game in range(len(league.games)):
            self.match_rate[2*game] = Poisson('match_rate_%i'%(2*game),
                    mu=self.goal_rate[league.games[game].hometeam.team_id] + self.home_adv,
                    value=league.games[game].homescore, observed=True)
            self.match_rate[2*game+1] = Poisson('match_rate_%i'%(2*game+1),
                    mu=self.goal_rate[league.games[game].hometeam.team_id],
                    value=league.games[game].homescore, observed=True)

        for game in range(len(future_games)):
            self.match_goals_future[2*game] = Poisson('match_goals_future_%i'%(2*game),
                    mu=self.goal_rate[future_games[game][0].team_id] + self.home_adv)
            self.match_goals_future[2*game+1] = Poisson('match_goals_future_%i'%(2*game+1),
                    mu=self.goal_rate[future_games[game][1].team_id])

    def run_mc(self,nsample = 10000,interactive=False):
        """run the model using mcmc"""
        from pymc.Matplot import plot
        from pymc import MCMC
        self.M = MCMC(self)
        if interactive:
            self.M.isample(iter=nsample, burn=1000, thin=10)
        else:
            self.M.sample(iter=nsample, burn=1000, thin=10)
        plot(self.M)

def fit_std_curve_by_pymc(i_vals, i_sds, dpx_concs):
    import pymc
    from pymc import Uniform, stochastic, deterministic, MCMC
    from pymc import Matplot
    # Define prior distributions for both Ka and Kd
    ka = Uniform('ka', lower=0, upper=1000)
    kd = Uniform('kd', lower=0, upper=1000)

    @stochastic(plot=True, observed=True)
    def quenching_model(ka=ka, kd=kd, value=i_vals):
        pred_i = quenching_func(ka, kd, dpx_concs)
        # The first concentration in dpx_concs should always be zero
        # (that is, the first point in the titration should be the
        # unquenched fluorescence), so we assert that here:
        assert dpx_concs[0] == 0
        # The reason this is necessary is that in the likelihood calculation
        # we skip the error for the first point, since (when the std. err
        # is calculated by well) the error is 0 (the I / I_0 ratio is
        # always 1 for each well, the the variance/SD across the wells is 0).
        # If we don't skip this first point, we get nan for the likelihood.
        # In addition, the model always predicts 1 for the I / I_0 ratio
        # when the DPX concentration is 0, so it contributes nothing to
        # the overall fit.
        return -np.sum((value[1:] - pred_i[1:])**2 / (2 * i_sds[1:]**2))

    pymc_model = pymc.Model([ka, kd, quenching_model])
    mcmc = MCMC(pymc_model)
    mcmc.sample(iter=155000, burn=5000, thin=150)
    Matplot.plot(mcmc)

    plt.figure()
    num_to_plot = 1000
    ka_vals = mcmc.trace('ka')[:]
    kd_vals = mcmc.trace('kd')[:]
    if num_to_plot > len(ka_vals):
        num_to_plot = len(ka_vals)
    for i in range(num_to_plot):
        plt.plot(dpx_concs, quenching_func(ka_vals[i], kd_vals[i], dpx_concs),
                 alpha=0.01, color='r')
    plt.errorbar(dpx_concs, i_vals, yerr=i_sds, linestyle='', marker='o',
            color='k', linewidth=2)

    return (ka_vals, kd_vals)

 def run_mc(self,nsample = 10000,interactive=False):
     """run the model using mcmc"""
     from pymc.Matplot import plot
     from pymc import MCMC
     self.M = MCMC(self)
     if interactive:
         self.M.isample(iter=nsample, burn=1000, thin=10)
     else:
         self.M.sample(iter=nsample, burn=1000, thin=10)
     plot(self.M)

 def run_mc(self,nsample = 10000,interactive=False,doplot=False,verbose=0):
     """run the model using mcmc"""
     from pymc import MCMC
     self.M = MCMC(self)
     if interactive:
         self.M.isample(iter=nsample, burn=1000, thin=10,verbose=verbose)
     else:
         self.M.sample(iter=nsample, burn=1000, thin=10,verbose=verbose)
     if doplot:
         from pymc.Matplot import plot
         plot(self.M)