Python pymc.MCMC类代码示例

    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