Python pymc.Normal方法代码示例

# 需要导入模块: import pymc [as 别名]
# 或者: from pymc import Normal [as 别名]
def LogNormalWrapper(name, mean, stddev):
    """
    Create a PyMC Normal stochastic, automatically converting parameters to be appropriate for a `LogNormal` distribution.
    Note that the resulting distribution is Normal, not LogNormal.
    This is appropriate if we want to describe the log of a volume or a fluorescence reading.

    Parameters
    ----------
    mean : float
       Mean of exp(X), where X is lognormal variate.
    stddev : float
       Standard deviation of exp(X), where X is lognormal variate.

    Returns
    -------
    stochastic : pymc.Normal
       Normal stochastic for random variable X

    """
    # Compute parameters of lognormal distribution
    mu = np.log(mean**2 / np.sqrt(stddev**2 + mean**2))
    tau = np.sqrt(np.log(1.0 + (stddev/mean)**2))**(-2)
    stochastic = pymc.Normal(name, mu=mu, tau=tau)
    return stochastic 

# 需要导入模块: import pymc [as 别名]
# 或者: from pymc import Normal [as 别名]
def NormalWrapper(name, mean, stddev, size=1, observed=False, value=None):
    """
    Create a PyMC Normal stochastic, automatically converting parameters to be appropriate for a `Normal` distribution.

    Notes
    -----
    Everything is coerced into a 1D array.

    Parameters
    ----------
    mean : float or array-like (but not pymc variable)
       Mean of exp(X), where X is lognormal variate.
    stddev : float or array-like (but not pymc variable)
       Standard deviation of exp(X), where X is lognormal variate.
    size : list of int, optional, default=1
       Size vector
    observed : bool, optional, default=False
       If True, the stochastic is fixed to this observed value.
    value : float, optional, default=None
       If observed=True, the observed value of the real quantity

    Returns
    -------
    stochastic : pymc.Normal
       Normal stochastic for random variable X.
       Name is prefixed with log_prefix
       Deterministic encoding the exponentiated lognormal (real quantity)

    """
    @pymc.deterministic(name=name + '_mu')
    def mu(mean=mean, stddev=stddev):
        mu = mean
        return mu
    @pymc.deterministic(name=name + '_tau')
    def tau(mean=mean, stddev=stddev):
        tau = stddev**(-2)
        return tau
    stochastic = pymc.Normal(name, mu=mu, tau=tau, size=size, observed=observed, value=value)

    return stochastic

#=============================================================================================
# PyMC models
#============================================================================================= 

# 需要导入模块: import pymc [as 别名]
# 或者: from pymc import Normal [as 别名]
def run_mcmc(pymc_model, nthin=20, nburn=0, niter=20000, map=True, db='ram', dbname=None):
    """
    Sample the model with pymc. Initial values of the parameters can be chosen with a maximum a posteriori estimate.

    Parameters
    ----------
    pymc_model : pymc model
       The pymc model to sample.
    nthin: int
        The number of MCMC steps that constitute 1 iteration.
    nburn: int
        The number of MCMC iterations during the burn-in.
    niter: int
        The number of production iterations.
    map: bool
        Whether to initialize the parameters before MCMC with the maximum a posteriori estimate.
    db : str
        How to store model, default = 'ram' means not storing it. To store model use storage = 'pickle'. If not,
        supply the name of the database backend that will store the values of the stochastics and deterministics sampled
        during the MCMC loop.
   dbname : str
        name for storage object, default = None. To store model use e.g. dbname = 'my_mcmc.pickle'

    Returns
    -------
    mcmc : pymc.MCMC
       The MCMC samples.

    """
    # Find MAP:
    if map == True:
        pymc.MAP(pymc_model).fit()

    # Sample the model with pymc
    mcmc = pymc.MCMC(pymc_model, db=db, dbname=dbname, name='Sampler', verbose=True)

    step_methods = 'AdaptiveMetropolis'
    mcmc.use_step_method(pymc.Metropolis, getattr(pymc_model, 'DeltaG'), proposal_sd=0.1, proposal_distribution='Normal')
    mcmc.use_step_method(pymc.Metropolis, getattr(pymc_model, 'log_F_PL'), proposal_sd=0.1, proposal_distribution='Normal')
    mcmc.use_step_method(pymc.Metropolis, getattr(pymc_model, 'log_F_P'), proposal_sd=0.1, proposal_distribution='Normal')
    mcmc.use_step_method(pymc.Metropolis, getattr(pymc_model, 'log_F_L'), proposal_sd=0.1, proposal_distribution='Normal')
    mcmc.use_step_method(pymc.Metropolis, getattr(pymc_model, 'log_F_plate'), proposal_sd=0.1, proposal_distribution='Normal')
    mcmc.use_step_method(pymc.Metropolis, getattr(pymc_model, 'log_F_buffer'), proposal_sd=0.1, proposal_distribution='Normal')
    mcmc.use_step_method(pymc.Metropolis, getattr(pymc_model, 'log_F_buffer_control'), proposal_sd=0.1, proposal_distribution='Normal')
    if hasattr(pymc_model, 'epsilon_ex'):
        mcmc.use_step_method(pymc.Metropolis, getattr(pymc_model, 'epsilon_ex'), proposal_sd=10000.0, proposal_distribution='Normal')
    if hasattr(pymc_model, 'epsilon_em'):
        mcmc.use_step_method(pymc.Metropolis, getattr(pymc_model, 'epsilon_em'), proposal_sd=10000.0, proposal_distribution='Normal')
    # Uncomment below to use log_F_PL highly correlated with DeltaG to improve sampling somewhat.
    #mcmc.use_step_method(pymc.AdaptiveMetropolis, [pymc_model.log_F_PL, pymc_model.DeltaG], scales={ pymc_model.log_F_PL : 0.1, pymc_model.DeltaG : 0.1  })

    mcmc.sample(iter=(nburn+niter), burn=nburn, thin=nthin, progress_bar=False, tune_throughout=True)

    return mcmc 

# 需要导入模块: import pymc [as 别名]
# 或者: from pymc import Normal [as 别名]
def _create_competitive_binding_model(self, receptor_name, DeltaG_prior):
        """
        Create the binding free energy priors, binding reaction models, and list of all species whose concentrations will be tracked.

        Populates the following fields:
        * parameter_names['DeltaGs'] : parameters associated with DeltaG priors
        * reactions : reactions for GeneralBindingModel with DeltaG parameters substituted for free energies
          formatted like [ ('DeltaG1', {'RL': -1, 'R' : +1, 'L' : +1}), ('DeltaG2', {'RP' : -1, 'R' : +1, 'P' : +1}) ]
        * conservation_equations : conservation relationships for GeneralBindingModel with species names substituted with log concentrations
          formatted like [ ('receptor', {'receptor' : +1, 'receptor:ligand' : +1}), ('ligand' : {'receptor:ligand' : +1, 'ligand' : +1}) ]
        * complex_names : names of all complexes
        * all_species : names of all species (ligands, receptor, complexes)
        """
        self.parameter_names['binding affinities'] = list()
        self.complex_names = list()
        self.reactions = list() # reactions for GeneralBindingModel with pymc parameters substituted for free energies
        self.conservation_equations = list() # list of conservation relationships for GeneralBindingModel with pymc parameters substituted for log concentrations
        receptor_conservation_equation = { receptor_name : +1 } # Begin to populate receptor conservation equation
        for ligand_name in self.ligand_names:
            # Create complex name
            complex_name = receptor_name + ':' + ligand_name
            self.complex_names.append(complex_name)
            # Create the DeltaG prior
            name = 'DeltaG (%s + %s -> %s)' % (receptor_name, ligand_name, complex_name) # form the name of the pymc variable
            if DeltaG_prior == 'uniform':
                DeltaG = pymc.Uniform(name, lower=DG_min, upper=DG_max) # binding free energy (kT), uniform over huge range
            elif DeltaG_prior == 'chembl':
                DeltaG = pymc.Normal(name, mu=0, tau=1./(12.5**2)) # binding free energy (kT), using a Gaussian prior inspured by ChEMBL
            else:
                raise Exception("DeltaG_prior = '%s' unknown. Must be one of 'uniform' or 'chembl'." % DeltaG_prior)
            self.model[name] = DeltaG
            self.parameter_names['binding affinities'].append(name)
            # Create the reaction for GeneralBindingModel
            self.reactions.append( (DeltaG, {complex_name : -1, receptor_name : +1, ligand_name : +1}) )
            # Create the conservation equation for GeneralBindingModel
            self.conservation_equations.append( (ligand_name, {ligand_name : +1, complex_name : +1}) )
            # Keep track of receptor conservation.
            receptor_conservation_equation[complex_name] = +1
        # Create the receptor conservation equation for GeneralBindingModel
        self.conservation_equations.append( [receptor_name, receptor_conservation_equation] )

        # Create a list of all species that may be present in assay plate
        self.all_species = [self.receptor_name] + self.ligand_names + self.complex_names