Python MCMC.isample方法代码示例

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import isample [as 别名]
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)

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import isample [as 别名]
def test_interactive():
    S = MCMC(disaster_model)
    S.isample(
        200,
        100,
        2,
        out=open(
            'testresults/interactive.log',
            'w'),
        progress_bar=0)

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import isample [as 别名]
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)

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import isample [as 别名]
    # Print comparison.
    signed_errors = numpy.zeros([len(typed_molecules)], numpy.float64)
    for (i, molecule) in enumerate(typed_molecules):
        # Get metadata.
        name = OEGetSDData(molecule, 'name').strip()
        try:
            dg_exp           = float(OEGetSDData(molecule, 'dG(exp)')) * units.kilocalories_per_mole
            signed_errors[i] = energies[molecule] / units.kilocalories_per_mole - dg_exp / units.kilocalories_per_mole
        except Exception as exception:
            # We couldn't find an experimental dG in the SDF file---ditch this molecule.
            print "Couldn't find dG(exp) for molecule '%s'; discarding it." % name
            typed_molecules.remove(molecule)
            continue

        # Form output.
        outstring = "%48s %8.3f %8.3f" % (name, dg_exp / units.kilocalories_per_mole, energies[molecule] / units.kilocalories_per_mole)

        print outstring

    print "Initial RMS error %8.3f kcal/mol" % (signed_errors.std())

    # Create MCMC model.
    model = create_model(typed_molecules, parameters)

    # Sample models.
    from pymc import MCMC
    sampler = MCMC(model, db='txt', name=mcmcDbName)
    sampler.isample(iter=mcmcIterations, burn=0, save_interval=1, verbose=True)
    #sampler.sample(iter=mcmcIterations, burn=0, save_interval=1, verbose=True, progress_bar=True)
    sampler.db.close()

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import isample [as 别名]
# PyMC implementation of Panel 6.4 from Royle & Dorazio (2008) pp. 217
# MA MacNeil - 04.03.14

import Mbht
import sys
import os
import pdb
from pymc import MCMC, BinaryMetropolis, Metropolis, AdaptiveMetropolis
from pymc import Matplot as mp
import pdb


M = MCMC(Mbht)
#M = MCMC(models,db='sqlite',dbname='xx_dbase')
xex = 6
M.isample(10**xex, 10**xex-10**(xex-1), thin=10**(xex-4), verbose=2)
#M.isample(100000, 80000, thin=10, verbose=2)

try:
    os.mkdir('Outputs')
except OSError:
    pass
os.chdir('Outputs')
M.write_csv("zz_results.csv")
mp.plot(M)

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import isample [as 别名]
    
    # Print comparison.
    for molecule in typed_molecules:
        # Get metadata.
        name = OEGetSDData(molecule, 'name').strip()
        try:
            dg_exp = float(OEGetSDData(molecule, 'dG(exp)')) * units.kilocalories_per_mole
        except Exception as exception:
            # We couldn't find an experimental dG in the SDF file---ditch this molecule.
            print "Couldn't find dG(exp) for molecule '%s'; discarding it." % name
            typed_molecules.remove(molecule)
            continue

        # Form output.
        outstring = "%48s %8.3f %8.3f" % (name, dg_exp / units.kilocalories_per_mole, energies[molecule] / units.kilocalories_per_mole)
        
        print outstring

    # Create MCMC model.
    model = create_model(typed_molecules, parameters)

    # Sample models.
    from pymc import MCMC
    sampler = MCMC(model, db='txt')
    sampler.isample(iter=10, burn=0, save_interval=1, verbose=True)
    
    
    
    
    

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import isample [as 别名]
def test_interactive():
    S = MCMC(DisasterModel)
    S.isample(200,100,2,  out=open('testresults/interactive.log', 'w'))

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import isample [as 别名]
def r(s=s, e=e, l=l):
    """ Concatenate Poisson means """
    out = np.empty(len(disasters_array))
    out[:s] = e
    out[s:] = l
    return out



from pymc.examples import DisasterModel
from pymc import MCMC
M = MCMC(DisasterModel)



M.isample(iter=10000, burn=1000, thin=10)



M.trace('s')[:]
#array([41, 40, 40, ..., 43, 44, 44])



from pylab import hist, show
hist(M.trace('l')[:])
#(array([   8,   52,  565, 1624, 2563, 2105, 1292,  488,  258,   45]),
#array([ 0.52721865,  0.60788251,  0.68854637,  0.76921023,  0.84987409,
#   0.93053795,  1.01120181,  1.09186567,  1.17252953,  1.25319339]),
#<a list of 10 Patch objects>)
show()

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import isample [as 别名]

#.........这里部分代码省略.........
        #@pm.deterministic
        #def C(eval_fun = gp.matern.euclidean, diff_degree=diff_degree, amp=amp, scale=scale):
        #    return gp.NearlyFullRankCovariance(eval_fun, diff_degree=diff_degree, amp=amp, scale=scale)
                    
        self.goal_rate = np.empty(N,dtype=object)
        self.def_rate = np.empty(N,dtype=object)
        self.goal_var = np.empty(N,dtype=object)
        self.def_var = np.empty(N,dtype=object)
        self.match_rate = np.empty(len(league.games)*2,dtype=object)
        self.outcome_future = np.empty(len(league.games),dtype=object)
        self.match_goals_future = np.empty(len(league.future_games)*2,dtype=object)
        self.home_adv = Uniform(name = 'home_adv',lower=0.,upper=2.0)
        self.league = league
        
        fmesh = np.arange(0.,league.n_days+2.)            

        for t in league.teams.values():
            # Prior parameters of C
            diff_degree_g = pm.Uniform('diff_degree_g_%i'%t.team_id, 1., 3)
            amp_g = pm.Uniform('amp_g_%i'%t.team_id, .01, 2.)
            scale_g = pm.Uniform('scale_g_%i'%t.team_id, 1., 10.)
            diff_degree_d = pm.Uniform('diff_degree_d_%i'%t.team_id, 1., 3)
            amp_d = pm.Uniform('amp_d_%i'%t.team_id, .01, 2.)
            scale_d = pm.Uniform('scale_d_%i'%t.team_id, 1., 10.)
            
            @pm.deterministic(name='C_d%i'%t.team_id)
            def C_d(eval_fun = gp.matern.euclidean, diff_degree=diff_degree_d, amp=amp_d, scale=scale_d):
                return gp.NearlyFullRankCovariance(eval_fun, diff_degree=diff_degree, amp=amp, scale=scale)
            
            @pm.deterministic(name='C_g%i'%t.team_id)
            def C_g(eval_fun = gp.matern.euclidean, diff_degree=diff_degree_g, amp=amp_g, scale=scale_g):
                return gp.NearlyFullRankCovariance(eval_fun, diff_degree=diff_degree, amp=amp, scale=scale)
            
            
            self.goal_rate[t.team_id] = Exponential('goal_rate_%i'%t.team_id,beta=1)
            self.def_rate[t.team_id] = Exponential('def_rate_%i'%t.team_id,beta=1)
            
            @pm.deterministic(name='M_d%i'%t.team_id)
            def M_d(eval_fun = linfun, c=self.def_rate[t.team_id]):
                return gp.Mean(eval_fun, c=c)
            @pm.deterministic(name='M_g%i'%t.team_id)
            def M_g(eval_fun = linfun, c=self.goal_rate[t.team_id]):
                return gp.Mean(eval_fun, c=c)
            
            self.def_var[t.team_id] = gp.GPSubmodel('smd_%i'%t.team_id,M_d,C_d,fmesh)
            self.goal_var[t.team_id] = gp.GPSubmodel('smg_%i'%t.team_id,M_g,C_g,fmesh)


        for game in range(len(league.games)):
            gd = int(game/(league.n_teams/2))
            assert(gd<league.n_days)
            self.match_rate[2*game] = Poisson('match_rate_%i'%(2*game),
                    mu=Deterministic(eval=clip_rate,
                                     parents={'val':
                                        self.goal_var[league.games[game].hometeam.team_id].f_eval[gd] - 
                                        self.def_var[league.games[game].awayteam.team_id].f_eval[gd] + self.home_adv},
                                     doc='clipped goal rate',name='clipped_h_%i'%game),
                    value=league.games[game].homescore, observed=True)
            self.match_rate[2*game+1] = Poisson('match_rate_%i'%(2*game+1),
                    mu=Deterministic(eval=clip_rate,
                                     parents={'val':
                                        self.goal_var[league.games[game].awayteam.team_id].f_eval[gd] - 
                                        self.def_var[league.games[game].hometeam.team_id].f_eval[gd]},
                                     doc='clipped goal rate',name='clipped_a_%i'%game),
                    value=league.games[game].awayscore, observed=True)


        for game in range(len(league.future_games)):
            gd = league.n_days
            self.match_goals_future[2*game] = Poisson('match_goals_future_%i_home'%game,
                    mu=Deterministic(eval=clip_rate,
                                     parents={'val':
                                        self.goal_var[league.future_games[game][0].team_id].f_eval[gd] - 
                                        self.def_var[league.future_games[game][1].team_id].f_eval[gd] + self.home_adv},
                                     doc='clipped goal rate',name='clipped_fut_h_%i'%game))

            self.match_goals_future[2*game+1] = Poisson('match_goals_future_%i_away'%game,
                    mu=Deterministic(eval=clip_rate,
                                     parents={'val':
                                        self.goal_var[league.future_games[game][1].team_id].f_eval[gd] - 
                                        self.def_var[league.future_games[game][0].team_id].f_eval[gd]},
                                     doc='clipped goal rate',name='clipped_fut_a_%i'%game))

            self.outcome_future[game] = Deterministic(eval=outcome_eval,parents={
                'home':self.match_goals_future[2*game],
                'away':self.match_goals_future[2*game+1]},name='match_outcome_future_%i'%game,
                dtype=int,doc='The outcome of the match'
                )
            
    def run_mc(self,nsample = 30000,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=30,verbose=verbose)
        else:
            self.M.sample(iter=nsample, burn=1000, thin=30,verbose=verbose)
        if doplot:
            from pymc.Matplot import plot
            plot(self.M)

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import isample [as 别名]
class LeagueBasicModel(object):
    """MCMC model of a football league."""
    
    #TODO: optimal Kelly Bettor
    #TODO: refine model
    #TODO: identify columns for autotesting
    
    def __init__(self, fname, playedto=None):
        super(LeagueBasicModel, self).__init__()
        league = League(fname,playedto)

        N = len(league.teams)
        def outcome_eval(home=None,away=None):
            if home > away:
                return 1
            if home < away:
                return -1
            if home == away:
                return 0
            
        def clip_rate(val):
            if val>0.2:return val
            else: return 0.2
        def linfun(x, c):
            return 0.*x+ c
                    
        self.goal_rate = np.empty(N,dtype=object)
        self.match_rate = np.empty(len(league.games)*2,dtype=object)
        self.outcome_future = np.empty(len(league.games),dtype=object)
        self.match_goals_future = np.empty(len(league.future_games)*2,dtype=object)
        self.home_adv = Uniform(name = 'home_adv',lower=0.,upper=2.0)
        self.league = league
        

        for t in league.teams.values():
            
            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=Deterministic(eval=clip_rate,
                                     parents={'val':
                                        self.goal_rate[league.games[game].hometeam.team_id] + self.home_adv},
                                     doc='clipped goal rate',name='clipped_h_%i'%game),
                    value=league.games[game].homescore, observed=True)
            self.match_rate[2*game+1] = Poisson('match_rate_%i'%(2*game+1),
                    mu=Deterministic(eval=clip_rate,
                                     parents={'val':
                                        self.goal_rate[league.games[game].awayteam.team_id]},
                                     doc='clipped goal rate',name='clipped_a_%i'%game),
                    value=league.games[game].awayscore, observed=True)


        for game in range(len(league.future_games)):
            self.match_goals_future[2*game] = Poisson('match_goals_future_%i_home'%game,
                    mu=Deterministic(eval=clip_rate,
                                     parents={'val':
                                        self.goal_rate[league.future_games[game][0].team_id] + self.home_adv},
                                     doc='clipped goal rate',name='clipped_fut_h_%i'%game))

            self.match_goals_future[2*game+1] = Poisson('match_goals_future_%i_away'%game,
                    mu=Deterministic(eval=clip_rate,
                                     parents={'val':
                                        self.goal_rate[league.future_games[game][1].team_id]},
                                     doc='clipped goal rate',name='clipped_fut_a_%i'%game))

            self.outcome_future[game] = Deterministic(eval=outcome_eval,parents={
                'home':self.match_goals_future[2*game],
                'away':self.match_goals_future[2*game+1]},name='match_outcome_future_%i'%game,
                dtype=int,doc='The outcome of the match'
                )
            
    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)

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import isample [as 别名]
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)

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import isample [as 别名]
def bayesian_model():
    model = MCMC(disastermodel)
    model.isample(iter=10000, burn=1000, thin=10)
    print model.trace('l')[:]
    plot(model) 
    return model