Python MCMC.use_step_method方法代码示例

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import use_step_method [as 别名]
def analizeM():
	M = MCMC(dm)
	print("M: ", M)

	M.sample(iter=10000, burn=1000, thin=10)
	print("M t: ", M.trace('switchpoint')[:])

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

	plot(M)
	# show()

	print("M smd dm sp: ", M.step_method_dict[dm.switchpoint])
	print("M smd dm em: ", M.step_method_dict[dm.early_mean])
	print("M smd dm lm: ", M.step_method_dict[dm.late_mean])

	M.use_step_method(Metropolis, dm.late_mean, proposal_sd=2.)

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

#.........这里部分代码省略.........
  mcmcParams=args['mcmcString']
  surveyParams=args['surveyString']
  priorParams=args['priorsString']

  maxIter=int(mcmcParams[0])
  burnIter=int(mcmcParams[1])
  thinFactor=int(mcmcParams[2])

  if surveyParams[5] == 'Inf':
    magLim = np.Inf
  else:
    magLim = float(surveyParams[5])
  S=U.UniformDistributionSingleLuminosity(int(surveyParams[0]), float(surveyParams[1]),
      float(surveyParams[2]), float(surveyParams[3]), float(surveyParams[4]),
      surveyLimit=magLim)
  #S.setRandomNumberSeed(53949896)
  S.generateObservations()
  lumCalModel=L.UniformSpaceDensityGaussianLFBook(S,float(surveyParams[1]), float(surveyParams[2]),
      float(priorParams[0]), float(priorParams[1]), float(priorParams[2]), float(priorParams[3]))

  class SurveyData(IsDescription):
    """
    Class that holds the data model for the data from the simulated parallax survey. Intended for use
    with the HDF5 files through the pytables package.
    """
    trueParallaxes = Float64Col(S.numberOfStarsInSurvey)
    absoluteMagnitudes = Float64Col(S.numberOfStarsInSurvey)
    apparentMagnitudes = Float64Col(S.numberOfStarsInSurvey)
    parallaxErrors = Float64Col(S.numberOfStarsInSurvey)
    magnitudeErrors = Float64Col(S.numberOfStarsInSurvey)
    observedParallaxes = Float64Col(S.numberOfStarsInSurvey)
    observedMagnitudes = Float64Col(S.numberOfStarsInSurvey)

  baseName="LumCalSimSurvey-{0}".format(S.numberOfStars)+"-{0}".format(S.minParallax)
  baseName=baseName+"-{0}".format(S.maxParallax)+"-{0}".format(S.meanAbsoluteMagnitude)
  baseName=baseName+"-{0}".format(S.varianceAbsoluteMagnitude)

  h5file = openFile(baseName+".h5", mode = "w", title = "Simulated Survey")
  group = h5file.createGroup("/", 'survey', 'Survey parameters, data, and MCMC parameters')
  parameterTable = h5file.createTable(group, 'parameters', SurveyParameters, "Survey parameters")
  dataTable = h5file.createTable(group, 'data', SurveyData, "Survey data")
  mcmcTable = h5file.createTable(group, 'mcmc', McmcParameters, "MCMC parameters")

  surveyParams = parameterTable.row
  surveyParams['kind']=S.__class__.__name__
  surveyParams['numberOfStars']=S.numberOfStars
  surveyParams['minParallax']=S.minParallax
  surveyParams['maxParallax']=S.maxParallax
  surveyParams['meanAbsoluteMagnitude']=S.meanAbsoluteMagnitude
  surveyParams['varianceAbsoluteMagnitude']=S.varianceAbsoluteMagnitude
  surveyParams['parallaxErrorNormalizationMagnitude']=S.parallaxErrorNormalizationMagnitude
  surveyParams['parallaxErrorSlope']=S.parallaxErrorSlope
  surveyParams['parallaxErrorCalibrationFloor']=S.parallaxErrorCalibrationFloor
  surveyParams['magnitudeErrorNormalizationMagnitude']=S.magnitudeErrorNormalizationMagnitude
  surveyParams['magnitudeErrorSlope']=S.magnitudeErrorSlope
  surveyParams['magnitudeErrorCalibrationFloor']=S.magnitudeErrorCalibrationFloor
  surveyParams['apparentMagnitudeLimit']=S.apparentMagnitudeLimit
  surveyParams['numberOfStarsInSurvey']=S.numberOfStarsInSurvey
  surveyParams.append()
  parameterTable.flush()

  surveyData = dataTable.row
  surveyData['trueParallaxes']=S.trueParallaxes
  surveyData['absoluteMagnitudes']=S.absoluteMagnitudes
  surveyData['apparentMagnitudes']=S.apparentMagnitudes
  surveyData['parallaxErrors']=S.parallaxErrors
  surveyData['magnitudeErrors']=S.magnitudeErrors
  surveyData['observedParallaxes']=S.observedParallaxes
  surveyData['observedMagnitudes']=S.observedMagnitudes
  surveyData.append()
  dataTable.flush()

  mcmcParameters = mcmcTable.row
  mcmcParameters['iterations']=maxIter
  mcmcParameters['burnIn']=burnIter
  mcmcParameters['thin']=thinFactor
  mcmcParameters['minMeanAbsoluteMagnitude']=float(priorParams[0])
  mcmcParameters['maxMeanAbsoluteMagnitude']=float(priorParams[1])
  mcmcParameters['priorTau']="OneOverX"
  mcmcParameters['tauLow']=float(priorParams[2])
  mcmcParameters['tauHigh']=float(priorParams[3])
  mcmcParameters.append()
  dataTable.flush()

  h5file.close()

  # Run MCMC and store in HDF5 database
  baseName="LumCalResults-{0}".format(S.numberOfStars)+"-{0}".format(S.minParallax)
  baseName=baseName+"-{0}".format(S.maxParallax)+"-{0}".format(S.meanAbsoluteMagnitude)
  baseName=baseName+"-{0}".format(S.varianceAbsoluteMagnitude)

  M=MCMC(lumCalModel.pyMCModel, db='hdf5', dbname=baseName+".h5", dbmode='w', dbcomplevel=9,
      dbcomplib='bzip2')
  M.use_step_method(Metropolis, M.priorParallaxes)
  M.use_step_method(Metropolis, M.priorAbsoluteMagnitudes)
  start=now()
  M.sample(iter=maxIter, burn=burnIter, thin=thinFactor)
  finish=now()
  print "Elapsed time in seconds: %f" % (finish-start)
  M.db.close()

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import use_step_method [as 别名]
# M.sample(iter=400000, burn=50000, thin=10,verbose=0)
# np.save('mc_data/nm_rho_s.npy',M.trace('rho_s')[:])
# np.save('mc_data/nm_alpha.npy',M.trace('alpha')[:])
# np.save('mc_data/nm_beta.npy',M.trace('beta')[:])

# np.save('mc_data/nm_ia.npy',M.trace('interaction_angle')[:])
# np.save('mc_data/nm_il.npy',M.trace('interaction_length')[:])
# np.save('mc_data/nm_ig.npy',M.trace('ignore_length')[:])

# network model with alignment
M = MCMC(networkModelAlignMC)
M.use_step_method(
    pymc.AdaptiveMetropolis,
    [
        networkModelAlignMC.interaction_angle,
        networkModelAlignMC.interaction_length,
        networkModelAlignMC.align_weight,
        networkModelAlignMC.ignore_length,
    ],
    delay=1000,
)
M.sample(iter=400000, burn=50000, thin=10, verbose=0)
np.save("mc_data/nma_rho_s.npy", M.trace("rho_s")[:])
np.save("mc_data/nma_alpha.npy", M.trace("alpha")[:])
np.save("mc_data/nma_beta.npy", M.trace("beta")[:])

np.save("mc_data/nma_ia.npy", M.trace("interaction_angle")[:])
np.save("mc_data/nma_il.npy", M.trace("interaction_length")[:])
np.save("mc_data/nma_aw.npy", M.trace("align_weight")[:])
np.save("mc_data/nma_ig.npy", M.trace("ignore_length")[:])

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import use_step_method [as 别名]
@deterministic(plot=False)
def r(s=s, e=e, l=l):
    """Allocate appropriate mean to time series"""
    out = np.empty(len(disasters_array))
    # Early mean prior to switchpoint
    out[:s] = e
    # Late mean following switchpoint
    out[s:] = l
    return out

# Where the value of x is None, the value is taken as missing.
D = Impute('D', Poisson, x, mu=r)



M.step_method_dict[DisasterModel.s]
#[<pymc.StepMethods.DiscreteMetropolis object at 0x3e8cb50>]

M.step_method_dict[DisasterModel.e]
#[<pymc.StepMethods.Metropolis object at 0x3e8cbb0>]

M.step_method_dict[DisasterModel.l]
#[<pymc.StepMethods.Metropolis object at 0x3e8ccb0>]



from pymc import Metropolis
M.use_step_method(Metropolis, DisasterModel.l, proposal_sd=2.)

# 需要导入模块: from pymc import MCMC [as 别名]
# 或者: from pymc.MCMC import use_step_method [as 别名]
#            
#            if np.log(random.random()) > logp_p - logp:
#                
#                self.reject()
                
            
            
            
    
    return locals()
    
if __name__ == '__main__':
    
    model = make_bma()
    M = MCMC(model)
    M.use_step_method(model['ModelMetropolis'], model['regression_model'])
    M.sample(iter=5000, burn=1000, thin=1)
    
    model_chain = M.trace("regression_model")[:]
    
    from collections import Counter
    
    counts = Counter(model_chain).items()
    counts.sort(reverse=True, key=lambda x: x[1])
    
    for f in counts[:10]:
        columns = unpack(f[0])
        print('Visits:', f[1])
        print(np.array([1. if i in columns else 0 for i in range(0,M.rank)]))
        print(M.coefficients.flatten())
        X = sm.add_constant(M.X[:, columns], prepend=True)