Source code for pymc.distributions

"""
pymc.distributions

A collection of common probability distributions. The objects associated
with a distribution called 'dist' are:

  dist_like : function
    The log-likelihood function corresponding to dist. PyMC's convention
    is to sum the log-likelihoods of multiple input values, so all
    log-likelihood functions return a single float.
  rdist : function
    The random variate generator corresponding to dist. These take a
    'size' argument indicating how many variates should be generated.
  dist_expval : function
    Computes the expected value of a dist-distributed variable.
  Dist : Stochastic subclass
    Instances have dist_like as their log-probability function
    and rdist as their random function.
"""

#-------------------------------------------------------------------
# Decorate fortran functions from pymc.flib to ease argument passing
#-------------------------------------------------------------------
# TODO: Add exponweib_expval
# TODO: categorical, mvhypergeometric
# TODO: __all__

__docformat__ = 'reStructuredText'

from . import flib, utils
import numpy as np
import scipy.stats as stats
gaussian_kde = stats.gaussian_kde
from .Node import ZeroProbability
from .PyMCObjects import Stochastic, Deterministic
from .CommonDeterministics import Lambda
from numpy import pi, inf
import itertools
import pdb
from . import utils
import warnings

from . import six
from .six import print_
xrange = six.moves.xrange


def poiscdf(a, x):
    x = np.atleast_1d(x)
    a = np.resize(a, x.shape)
    values = np.array([flib.gammq(b, y) for b, y in zip(a.ravel(), x.ravel())])
    return values.reshape(x.shape)

# Import utility functions
import inspect
import types
from copy import copy
random_number = np.random.random
inverse = np.linalg.pinv

sc_continuous_distributions = ['beta', 'cauchy', 'chi2',
                               'degenerate', 'exponential', 'exponweib',
                               'gamma', 'half_cauchy', 'half_normal',
                               'inverse_gamma', 'laplace', 'logistic',
                               'lognormal', 'noncentral_t', 'normal',
                               'pareto', 't', 'truncated_pareto', 'uniform',
                               'weibull', 'skew_normal', 'truncated_normal',
                               'von_mises']
sc_bool_distributions = ['bernoulli']
sc_discrete_distributions = ['binomial', 'betabin', 'geometric', 'poisson',
                             'negative_binomial', 'categorical', 'hypergeometric',
                             'discrete_uniform', 'truncated_poisson']

sc_nonnegative_distributions = ['bernoulli', 'beta', 'betabin', 'binomial', 'chi2', 'exponential',
                                'exponweib', 'gamma', 'half_cauchy', 'half_normal',
                                'hypergeometric', 'inverse_gamma', 'lognormal',
                                'weibull']

mv_continuous_distributions = ['dirichlet', 'mv_normal',
                               'mv_normal_cov', 'mv_normal_chol', 'wishart',
                               'wishart_cov']

mv_discrete_distributions = ['multivariate_hypergeometric', 'multinomial']

mv_nonnegative_distributions = ['dirichlet', 'wishart',
                                'wishart_cov', 'multivariate_hypergeometric',
                                'multinomial']

availabledistributions = (sc_continuous_distributions +
                          sc_bool_distributions +
                          sc_discrete_distributions +
                          mv_continuous_distributions +
                          mv_discrete_distributions)

# Changes lower case, underscore-separated names into "Class style"
# capitalized names For example, 'negative_binomial' becomes
# 'NegativeBinomial'
capitalize = lambda name: ''.join([s.capitalize() for s in name.split('_')])


# ==============================================================================
# User-accessible function to convert a logp and random function to a
# Stochastic subclass.
# ==============================================================================

# TODO Document this function
def bind_size(randfun, shape):
    def newfun(*args, **kwargs):
        try:
            return np.reshape(randfun(size=shape, *args, **kwargs), shape)
        except ValueError:
            # Account for non-array return values
            return randfun(size=shape, *args, **kwargs)
    newfun.scalar_version = randfun
    return newfun


def new_dist_class(*new_class_args):
    """
    Returns a new class from a distribution.

    :Parameters:
      dtype : numpy dtype
        The dtype values of instances of this class.
      name : string
        Name of the new class.
      parent_names : list of strings
        The labels of the parents of this class.
      parents_default : list
        The default values of parents.
      docstr : string
        The docstring of this class.
      logp : function
        The log-probability function for this class.
      random : function
        The random function for this class.
      mv : boolean
        A flag indicating whether this class represents array-valued
        variables.

      .. note::
        stochastic_from_dist provides a higher-level version.
        stochastic_from_data is suited for non-parametric distributions.

      :SeeAlso:
        stochastic_from_dist, stochastic_from_data
    """

    (dtype,
     name,
     parent_names,
     parents_default,
     docstr,
     logp,
     random,
     mv,
     logp_partial_gradients) = new_class_args

    class new_class(Stochastic):
        __doc__ = docstr

        def __init__(self, *args, **kwds):
            (dtype,
             name,
             parent_names,
             parents_default,
             docstr,
             logp,
             random,
             mv,
             logp_partial_gradients) = new_class_args
            parents = parents_default

            # Figure out what argument names are needed.
            arg_keys = [
                'name',
                'parents',
                'value',
                'observed',
                'size',
                'trace',
                'rseed',
                'doc',
                'debug',
                'plot',
                'verbose']
            arg_vals = [
                None, parents, None, False, None, True, True, None, False, None, -1]
            if 'isdata' in kwds:
                warnings.warn(
                    '"isdata" is deprecated, please use "observed" instead.')
                kwds['observed'] = kwds['isdata']
                pass

            # No size argument allowed for multivariate distributions.
            if mv:
                arg_keys.pop(4)
                arg_vals.pop(4)

            arg_dict_out = dict(zip(arg_keys, arg_vals))
            args_needed = ['name'] + parent_names + arg_keys[2:]

            # Sort positional arguments
            for i in xrange(len(args)):
                try:
                    k = args_needed.pop(0)
                    if k in parent_names:
                        parents[k] = args[i]
                    else:
                        arg_dict_out[k] = args[i]
                except:
                    raise ValueError(
                        'Too many positional arguments provided. Arguments for class ' +
                        self.__class__.__name__ +
                        ' are: ' +
                        str(args_needed))

            # Sort keyword arguments
            for k in args_needed:
                if k in parent_names:
                    try:
                        parents[k] = kwds.pop(k)
                    except:
                        if k in parents_default:
                            parents[k] = parents_default[k]
                        else:
                            raise ValueError('No value given for parent ' + k)
                elif k in arg_dict_out.keys():
                    try:
                        arg_dict_out[k] = kwds.pop(k)
                    except:
                        pass

            # Remaining unrecognized arguments raise an error.
            if len(kwds) > 0:
                raise TypeError('Keywords ' +
                                str(kwds.keys()) +
                                ' not recognized. Arguments recognized are ' +
                                str(args_needed))

        # Determine size desired for scalar variables.
        # Notes
        # -----
        # Case | init_val     | parents       | size | value.shape | bind size
        # ------------------------------------------------------------------
        # 1.1  | None         | scalars       | None | 1           | 1
        # 1.2  | None         | scalars       | n    | n           | n
        # 1.3  | None         | n             | None | n           | 1
        # 1.4  | None         | n             | n(m) | n (Error)   | 1 (-)
        # 2.1  | scalar       | scalars       | None | 1           | 1
        # 2.2  | scalar       | scalars       | n    | n           | n
        # 2.3  | scalar       | n             | None | n           | 1
        # 2.4  | scalar       | n             | n(m) | n (Error)   | 1 (-)
        # 3.1  | n            | scalars       | None | n           | n
        # 3.2  | n            | scalars       | n(m) | n (Error)   | n (-)
        # 3.3  | n            | n             | None | n           | 1
        # 3.4  | n            | n             | n(m) | n (Error)   | 1 (-)

            if not mv:

                shape = arg_dict_out.pop('size')
                shape = None if shape is None else tuple(np.atleast_1d(shape))

                init_val = arg_dict_out['value']
                init_val_shape = None if init_val is None else np.shape(
                    init_val)

                if len(parents) > 0:
                    pv = [np.shape(utils.value(v)) for v in parents.values()]
                    biggest_parent = np.argmax(
                        [(np.prod(v) if v else 0) for v in pv])
                    parents_shape = pv[biggest_parent]

                    # Scalar parents can support any shape.
                    if np.prod(parents_shape) < 1:
                        parents_shape = None

                else:
                    parents_shape = None

                def shape_error():
                    raise ValueError(
                        'Shapes are incompatible: value %s, largest parent %s, shape argument %s' %
                        (shape, init_val_shape, parents_shape))

                if init_val_shape is not None and shape is not None and init_val_shape != shape:
                    shape_error()

                given_shape = init_val_shape or shape
                bindshape = given_shape or parents_shape

                # Check consistency of bindshape and parents_shape
                if parents_shape is not None:
                    # Uncomment to leave broadcasting completely up to NumPy's random functions
                    # if bindshape[-np.alen(parents_shape):]!=parents_shape:
                    # Uncomment to limit broadcasting flexibility to what the
                    # Fortran likelihoods can handle.
                    if bindshape < parents_shape:
                        shape_error()

                if random is not None:
                    random = bind_size(random, bindshape)

            elif 'size' in kwds.keys():
                raise ValueError(
                    'No size argument allowed for multivariate stochastic variables.')

            # Call base class initialization method
            if arg_dict_out.pop('debug'):
                logp = debug_wrapper(logp)
                random = debug_wrapper(random)
            else:
                Stochastic.__init__(
                    self,
                    logp=logp,
                    random=random,
                    logp_partial_gradients=logp_partial_gradients,
                    dtype=dtype,
                    **arg_dict_out)

    new_class.__name__ = name
    new_class.parent_names = parent_names
    new_class.parents_default = parents_default
    new_class.dtype = dtype
    new_class.mv = mv
    new_class.raw_fns = {'logp': logp, 'random': random}

    return new_class


def stochastic_from_dist(
        name, logp, random=None, logp_partial_gradients={}, dtype=np.float, mv=False):
    """
    Return a Stochastic subclass made from a particular distribution.

    :Parameters:
      name : string
        The name of the new class.
      logp : function
        The log-probability function.
      random : function
        The random function
      dtype : numpy dtype
        The dtype of values of instances.
      mv : boolean
        A flag indicating whether this class represents
        array-valued variables.

    :Example:
      >>> Exponential = stochastic_from_dist('exponential',
                                              logp=exponential_like,
                                              random=rexponential,
                                              dtype=np.float,
                                              mv=False)
      >>> A = Exponential(self_name, value, beta)

    .. note::
      new_dist_class is a more flexible class factory. Also consider
      subclassing Stochastic directly.
      stochastic_from_data is suited for non-parametric distributions.

    :SeeAlso:
      new_dist_class, stochastic_from_data
    """

    (args, defaults) = utils.get_signature(logp)
    parent_names = args[1:]
    try:
        parents_default = dict(zip(args[-len(defaults):], defaults))
    except TypeError:  # No parents at all.
        parents_default = {}

    name = capitalize(name)

    # Build docstring from distribution
    parents_str = ''
    if parent_names:
        parents_str = ', '.join(parent_names) + ', '
    docstr = name[0] + ' = ' + name + \
        '(name, ' + parents_str + 'value=None, observed=False,'
    if not mv:
        docstr += ' size=1,'
    docstr += ' trace=True, rseed=True, doc=None, verbose=-1, debug=False)\n\n'
    docstr += 'Stochastic variable with ' + name + \
        ' distribution.\nParents are: ' + ', '.join(parent_names) + '.\n\n'
    docstr += 'Docstring of log-probability function:\n'
    try:
        docstr += logp.__doc__
    except TypeError:
        pass  # This will happen when logp doesn't have a docstring

    logp = valuewrapper(logp)
    distribution_arguments = logp.__dict__

    wrapped_logp_partial_gradients = {}

    for parameter, func in six.iteritems(logp_partial_gradients):
        wrapped_logp_partial_gradients[parameter] = valuewrapper(
            logp_partial_gradients[parameter],
            arguments=distribution_arguments)

    return new_dist_class(dtype, name, parent_names, parents_default, docstr,
                          logp, random, mv, wrapped_logp_partial_gradients)


def stochastic_from_data(name, data, lower=-np.inf, upper=np.inf,
                         value=None, observed=False, trace=True, verbose=-1, debug=False):
    """
    Return a Stochastic subclass made from arbitrary data.

    The histogram for the data is fitted with Kernel Density Estimation.

    :Parameters:
      - `data`  : An array with samples (e.g. trace[:])
      - `lower` : Lower bound on possible outcomes
      - `upper` : Upper bound on possible outcomes

    :Example:
       >>> from pymc import stochastic_from_data
       >>> pos = stochastic_from_data('posterior', posterior_samples)
       >>> prior = pos # update the prior with arbitrary distributions

    :Alias:
      Histogram
    """
    pdf = gaussian_kde(data)  # automatic bandwidth selection

    # account for tail contribution
    lower_tail = upper_tail = 0.
    if lower > -np.inf:
        lower_tail = pdf.integrate_box(-np.inf, lower)
    if upper < np.inf:
        upper_tail = pdf.integrate_box(upper, np.inf)
    factor = 1. / (1. - (lower_tail + upper_tail))

    def logp(value):
        prob = factor * pdf(value)
        if value < lower or value > upper:
            return -np.inf
        elif prob <= 0.:
            return -np.inf
        else:
            return np.log(prob)

    def random():
        res = pdf.resample(1)[0][0]
        while res < lower or res > upper:
            res = pdf.resample(1)[0][0]
        return res

    if value is None:
        value = random()

    return Stochastic(logp=logp,
                      doc='Non-parametric density with Gaussian Kernels.',
                      name=name,
                      parents={},
                      random=random,
                      trace=trace,
                      value=value,
                      dtype=float,
                      observed=observed,
                      verbose=verbose)


# Alias following Stochastics naming convention
Histogram = stochastic_from_data

#-------------------------------------------------------------
# Light decorators
#-------------------------------------------------------------


def randomwrap(func):
    """
    Decorator for random value generators

    Allows passing of sequence of parameters, as well as a size argument.


    Convention:

      - If size=1 and the parameters are all scalars, return a scalar.
      - If size=1, the random variates are 1D.
      - If the parameters are scalars and size > 1, the random variates are 1D.
      - If size > 1 and the parameters are sequences, the random variates are
        aligned as (size, max(length)), where length is the parameters size.


    :Example:
      >>> rbernoulli(.1)
      0
      >>> rbernoulli([.1,.9])
      np.asarray([0, 1])
      >>> rbernoulli(.9, size=2)
      np.asarray([1, 1])
      >>> rbernoulli([.1,.9], 2)
      np.asarray([[0, 1],
             [0, 1]])
    """

    # Find the order of the arguments.
    refargs, defaults = utils.get_signature(func)
    # vfunc = np.vectorize(self.func)
    npos = len(refargs) - len(defaults)  # Number of pos. arg.
    nkwds = len(defaults)  # Number of kwds args.
    mv = func.__name__[
        1:] in mv_continuous_distributions + mv_discrete_distributions

    # Use the NumPy random function directly if this is not a multivariate
    # distribution
    if not mv:
        return func

    def wrapper(*args, **kwds):
        # First transform keyword arguments into positional arguments.
        n = len(args)
        if nkwds > 0:
            args = list(args)
            for i, k in enumerate(refargs[n:]):
                if k in kwds.keys():
                    args.append(kwds[k])
                else:
                    args.append(defaults[n - npos + i])

        r = []
        s = []
        largs = []
        nr = args[-1]
        length = [np.atleast_1d(a).shape[0] for a in args]
        dimension = [np.atleast_1d(a).ndim for a in args]
        N = max(length)
        if len(set(dimension)) > 2:
            raise('Dimensions do not agree.')
        # Make sure all elements are iterable and have consistent lengths, ie
        # 1 or n, but not m and n.

        for arg, s in zip(args, length):
            t = type(arg)
            arr = np.empty(N, type)
            if s == 1:
                arr.fill(arg)
            elif s == N:
                arr = np.asarray(arg)
            else:
                raise RuntimeError('Arguments size not allowed: %s.' % s)
            largs.append(arr)

        if mv and N > 1 and max(dimension) > 1 and nr > 1:
            raise ValueError(
                'Multivariate distributions cannot take s>1 and multiple values.')

        if mv:
            for i, arg in enumerate(largs[:-1]):
                largs[0] = np.atleast_2d(arg)

        for arg in zip(*largs):
            r.append(func(*arg))

        size = arg[-1]
        vec_stochastics = len(r) > 1
        if mv:
            if nr == 1:
                return r[0]
            else:
                return np.vstack(r)
        else:
            if size > 1 and vec_stochastics:
                return np.atleast_2d(r).T
            elif vec_stochastics or size > 1:
                return np.concatenate(r)
            else:  # Scalar case
                return r[0][0]

    wrapper.__doc__ = func.__doc__
    wrapper.__name__ = func.__name__
    return wrapper


def debug_wrapper(func, name):
    # Wrapper to debug distributions

    import pdb

    def wrapper(*args, **kwargs):

        print_('Debugging inside %s:' % name)
        print_('\tPress \'s\' to step into function for debugging')
        print_('\tCall \'args\' to list function arguments')

        # Set debugging trace
        pdb.set_trace()

        # Call function
        return func(*args, **kwargs)

    return wrapper


#-------------------------------------------------------------
# Utility functions
#-------------------------------------------------------------

def constrain(value, lower=-np.Inf, upper=np.Inf, allow_equal=False):
    """
    Apply interval constraint on stochastic value.
    """

    ok = flib.constrain(value, lower, upper, allow_equal)
    if ok == 0:
        raise ZeroProbability


def standardize(x, loc=0, scale=1):
    """
    Standardize x

    Return (x-loc)/scale
    """

    return flib.standardize(x, loc, scale)

# ==================================
# = vectorize causes memory leaks. =
# ==================================
# @Vectorize


def gammaln(x):
    """
    Logarithm of the Gamma function
    """

    return flib.gamfun(x)


def expand_triangular(X, k):
    """
    Expand flattened triangular matrix.
    """

    X = X.tolist()
    # Unflatten matrix
    Y = np.asarray(
        [[0] * i + X[i * k - (i * (i - 1)) / 2: i * k + (k - i)] for i in range(k)])
    # Loop over rows
    for i in range(k):
        # Loop over columns
        for j in range(k):
            Y[j, i] = Y[i, j]
    return Y


# Loss functions

absolute_loss = lambda o, e: absolute(o - e)

squared_loss = lambda o, e: (o - e) ** 2

chi_square_loss = lambda o, e: (1. * (o - e) ** 2) / e

loss_functions = {
    'absolute': absolute_loss,
    'squared': squared_loss,
    'chi_square': chi_square_loss}


def GOFpoints(x, y, expval, loss):
    # Return pairs of points for GOF calculation
    return np.sum(np.transpose([loss(x, expval), loss(y, expval)]), 0)


def gofwrapper(f, loss_function='squared'):
    """
    Goodness-of-fit decorator function for likelihoods
    ==================================================
    Generates goodness-of-fit points for data likelihoods.

    Wrap function f(*args, **kwds) where f is a likelihood.

    Assume args = (x, parameter1, parameter2, ...)
    Before passing the arguments to the function, the wrapper makes sure that
    the parameters have the same shape as x.
    """

    name = f.__name__[:-5]
    # Take a snapshot of the main namespace.

    # Find the functions needed to compute the gof points.
    expval_func = eval(name + '_expval')
    random_func = eval('r' + name)

    def wrapper(*args, **kwds):
        """
        This wraps a likelihood.
        """

        """Return gof points."""

        # Calculate loss
        loss = kwds.pop('gof', loss_functions[loss_function])

        # Expected value, given parameters
        expval = expval_func(*args[1:], **kwds)
        y = random_func(size=len(args[0]), *args[1:], **kwds)
        f.gof_points = GOFpoints(args[0], y, expval, loss)

        """Return likelihood."""

        return f(*args, **kwds)

    # Assign function attributes to wrapper.
    wrapper.__doc__ = f.__doc__
    wrapper.__name__ = f.__name__
    wrapper.name = name

    return wrapper

#--------------------------------------------------------
# Statistical distributions
# random generator, expval, log-likelihood
#--------------------------------------------------------

# Autoregressive lognormal


def rarlognormal(a, sigma, rho, size=1):
    R"""
    Autoregressive normal random variates.

    If a is a scalar, generates one series of length size.
    If a is a sequence, generates size series of the same length
    as a.
    """
    f = utils.ar1
    if np.isscalar(a):
        r = f(rho, 0, sigma, size)
    else:
        n = len(a)
        r = [f(rho, 0, sigma, n) for i in range(size)]
        if size == 1:
            r = r[0]
    return a * np.exp(r)


def arlognormal_like(x, a, sigma, rho):
    R"""
    Autoregressive lognormal log-likelihood.

    .. math::
        x_i & = a_i \exp(e_i) \\
        e_i & = \rho e_{i-1} + \epsilon_i

    where :math:`\epsilon_i \sim N(0,\sigma)`.
    """
    return flib.arlognormal(x, np.log(a), sigma, rho, beta=1)


# Bernoulli----------------------------------------------
@randomwrap
def rbernoulli(p, size=None):
    """
    Random Bernoulli variates.
    """

    return np.random.random(size) < p


def bernoulli_expval(p):
    """
    Expected value of bernoulli distribution.
    """

    return p


[docs]def bernoulli_like(x, p): R"""Bernoulli log-likelihood The Bernoulli distribution describes the probability of successes (x=1) and failures (x=0). .. math:: f(x \mid p) = p^{x} (1-p)^{1-x} :Parameters: - `x` : Series of successes (1) and failures (0). :math:`x=0,1` - `p` : Probability of success. :math:`0 < p < 1`. :Example: >>> from pymc import bernoulli_like >>> bernoulli_like([0,1,0,1], .4) -2.854232711280291 .. note:: - :math:`E(x)= p` - :math:`Var(x)= p(1-p)` """ return flib.bernoulli(x, p)
bernoulli_grad_like = {'p': flib.bern_grad_p} # Beta---------------------------------------------- @randomwrap def rbeta(alpha, beta, size=None): """ Random beta variates. """ from scipy.stats.distributions import beta as sbeta return sbeta.ppf(np.random.random(size), alpha, beta) # return np.random.beta(alpha, beta, size) def beta_expval(alpha, beta): """ Expected value of beta distribution. """ return 1.0 * alpha / (alpha + beta)
[docs]def beta_like(x, alpha, beta): R""" Beta log-likelihood. The conjugate prior for the parameter :math:`p` of the binomial distribution. .. math:: f(x \mid \alpha, \beta) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)} x^{\alpha - 1} (1 - x)^{\beta - 1} :Parameters: - `x` : 0 < x < 1 - `alpha` : alpha > 0 - `beta` : beta > 0 :Example: >>> from pymc import beta_like >>> beta_like(.4,1,2) 0.182321556793954 .. note:: - :math:`E(X)=\frac{\alpha}{\alpha+\beta}` - :math:`Var(X)=\frac{\alpha \beta}{(\alpha+\beta)^2(\alpha+\beta+1)}` """ # try: # constrain(alpha, lower=0, allow_equal=True) # constrain(beta, lower=0, allow_equal=True) # constrain(x, 0, 1, allow_equal=True) # except ZeroProbability: # return -np.Inf return flib.beta_like(x, alpha, beta)
beta_grad_like = {'value': flib.beta_grad_x, 'alpha': flib.beta_grad_a, 'beta': flib.beta_grad_b} # Binomial---------------------------------------------- @randomwrap def rbinomial(n, p, size=None): """ Random binomial variates. """ # return np.random.binomial(n,p,size) return np.random.binomial(np.ravel(n), np.ravel(p), size) def binomial_expval(n, p): """ Expected value of binomial distribution. """ return p * n
[docs]def binomial_like(x, n, p): R""" Binomial log-likelihood. The discrete probability distribution of the number of successes in a sequence of n independent yes/no experiments, each of which yields success with probability p. .. math:: f(x \mid n, p) = \frac{n!}{x!(n-x)!} p^x (1-p)^{n-x} :Parameters: - `x` : [int] Number of successes, > 0. - `n` : [int] Number of Bernoulli trials, > x. - `p` : Probability of success in each trial, :math:`p \in [0,1]`. .. note:: - :math:`E(X)=np` - :math:`Var(X)=np(1-p)` """ # Temporary hack to avoid issue #614 return flib.binomial(x, n, p)
binomial_grad_like = {'p': flib.binomial_gp} # Beta---------------------------------------------- @randomwrap def rbetabin(alpha, beta, n, size=None): """ Random beta-binomial variates. """ phi = np.random.beta(alpha, beta, size) return np.random.binomial(n, phi) def betabin_expval(alpha, beta, n): """ Expected value of beta-binomial distribution. """ return n * alpha / (alpha + beta) def betabin_like(x, alpha, beta, n): R""" Beta-binomial log-likelihood. Equivalent to binomial random variables with probabilities drawn from a :math:`\texttt{Beta}(\alpha,\beta)` distribution. .. math:: f(x \mid \alpha, \beta, n) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \frac{\Gamma(n+1)}{\Gamma(x+1)\Gamma(n-x+1)} \frac{\Gamma(\alpha + x)\Gamma(n+\beta-x)}{\Gamma(\alpha+\beta+n)} :Parameters: - `x` : x=0,1,\ldots,n - `alpha` : alpha > 0 - `beta` : beta > 0 - `n` : n=x,x+1,\ldots :Example: >>> betabin_like(3,1,1,10) -2.3978952727989 .. note:: - :math:`E(X)=n\frac{\alpha}{\alpha+\beta}` - :math:`Var(X)=n\frac{\alpha \beta}{(\alpha+\beta)^2(\alpha+\beta+1)}` """ return flib.betabin_like(x, alpha, beta, n) betabin_grad_like = {'alpha': flib.betabin_ga, 'beta': flib.betabin_gb} # Categorical---------------------------------------------- # Note that because categorical elements are not ordinal, there # is no expected value. #@randomwrap def rcategorical(p, size=None): """ Categorical random variates. """ out = flib.rcat(p, np.random.random(size=size)) if sum(out.shape) == 1: return out.squeeze() else: return out
[docs]def categorical_like(x, p): R""" Categorical log-likelihood. The most general discrete distribution. .. math:: f(x=i \mid p) = p_i for :math:`i \in 0 \ldots k-1`. :Parameters: - `x` : [int] :math:`x \in 0\ldots k-1` - `p` : [float] :math:`p > 0`, :math:`\sum p = 1` """ p = np.atleast_2d(p) if np.any(abs(np.sum(p, 1) - 1) > 0.0001): print_("Probabilities in categorical_like sum to", np.sum(p, 1)) return flib.categorical(np.array(x).astype(int), p)
# Cauchy---------------------------------------------- @randomwrap def rcauchy(alpha, beta, size=None): """ Returns Cauchy random variates. """ return alpha + beta * np.tan(pi * random_number(size) - pi / 2.0) def cauchy_expval(alpha, beta): """ Expected value of cauchy distribution. """ return alpha # In wikipedia, the arguments name are k, x0.
[docs]def cauchy_like(x, alpha, beta): R""" Cauchy log-likelihood. The Cauchy distribution is also known as the Lorentz or the Breit-Wigner distribution. .. math:: f(x \mid \alpha, \beta) = \frac{1}{\pi \beta [1 + (\frac{x-\alpha}{\beta})^2]} :Parameters: - `alpha` : Location parameter. - `beta` : Scale parameter > 0. .. note:: - Mode and median are at alpha. """ return flib.cauchy(x, alpha, beta)
cauchy_grad_like = {'value': flib.cauchy_grad_x, 'alpha': flib.cauchy_grad_a, 'beta': flib.cauchy_grad_b} # Chi square---------------------------------------------- @randomwrap def rchi2(nu, size=None): """ Random :math:`\chi^2` variates. """ return np.random.chisquare(nu, size) def chi2_expval(nu): """ Expected value of Chi-squared distribution. """ return nu
[docs]def chi2_like(x, nu): R""" Chi-squared :math:`\chi^2` log-likelihood. .. math:: f(x \mid \nu) = \frac{x^{(\nu-2)/2}e^{-x/2}}{2^{\nu/2}\Gamma(\nu/2)} :Parameters: - `x` : > 0 - `nu` : [int] Degrees of freedom ( nu > 0 ) .. note:: - :math:`E(X)=\nu` - :math:`Var(X)=2\nu` """ return flib.gamma(x, 0.5 * nu, 1. / 2)
chi2_grad_like = {'value': lambda x, nu: flib.gamma_grad_x(x, 0.5 * nu, 1. / 2), 'nu': lambda x, nu: flib.gamma_grad_alpha(x, 0.5 * nu, 1. / 2) * .5} # chi2_grad_like = {'x' : lambda x, nu : (nu / 2 - 1) / x -.5, # 'nu' : flib.chi2_grad_nu } # Degenerate--------------------------------------------- @randomwrap def rdegenerate(k, size=1): """ Random degenerate variates. """ return np.ones(size) * k def degenerate_expval(k): """ Expected value of degenerate distribution. """ return k
[docs]def degenerate_like(x, k): R""" Degenerate log-likelihood. .. math:: f(x \mid k) = \left\{ \begin{matrix} 1 \text{ if } x = k \\ 0 \text{ if } x \ne k\end{matrix} \right. :Parameters: - `x` : Input value. - `k` : Degenerate value. """ x = np.atleast_1d(x) return sum(np.log([i == k for i in x]))
# def degenerate_grad_like(x, k): # R""" # degenerate_grad_like(x, k) # # Degenerate gradient log-likelihood. # # .. math:: # f(x \mid k) = \left\{ \begin{matrix} 1 \text{ if } x = k \\ 0 \text{ if } x \ne k\end{matrix} \right. # # :Parameters: # - `x` : Input value. # - `k` : Degenerate value. # """ # return np.zeros(np.size(x))*k # Dirichlet---------------------------------------------- @randomwrap def rdirichlet(theta, size=1): """ Dirichlet random variates. """ gammas = np.vstack([rgamma(theta, 1) for i in xrange(size)]) if size > 1 and np.size(theta) > 1: return (gammas.T / gammas.sum(1))[:-1].T elif np.size(theta) > 1: return (gammas[0] / gammas[0].sum())[:-1] else: return 1. def dirichlet_expval(theta): """ Expected value of Dirichlet distribution. """ return theta / np.sum(theta).astype(float)
[docs]def dirichlet_like(x, theta): R""" Dirichlet log-likelihood. This is a multivariate continuous distribution. .. math:: f(\mathbf{x}) = \frac{\Gamma(\sum_{i=1}^k \theta_i)}{\prod \Gamma(\theta_i)}\prod_{i=1}^{k-1} x_i^{\theta_i - 1} \cdot\left(1-\sum_{i=1}^{k-1}x_i\right)^\theta_k :Parameters: x : (n, k-1) array Array of shape (n, k-1) where `n` is the number of samples and `k` the dimension. :math:`0 < x_i < 1`, :math:`\sum_{i=1}^{k-1} x_i < 1` theta : array An (n,k) or (1,k) array > 0. .. note:: Only the first `k-1` elements of `x` are expected. Can be used as a parent of Multinomial and Categorical nevertheless. """ x = np.atleast_2d(x) theta = np.atleast_2d(theta) if (np.shape(x)[-1] + 1) != np.shape(theta)[-1]: raise ValueError('The dimension of x in dirichlet_like must be k-1.') return flib.dirichlet(x, theta)
# Exponential---------------------------------------------- @randomwrap def rexponential(beta, size=None): """ Exponential random variates. """ return np.random.exponential(1. / beta, size) def exponential_expval(beta): """ Expected value of exponential distribution. """ return 1. / beta
[docs]def exponential_like(x, beta): R""" Exponential log-likelihood. The exponential distribution is a special case of the gamma distribution with alpha=1. It often describes the time until an event. .. math:: f(x \mid \beta) = \beta e^{-\beta x} :Parameters: - `x` : x > 0 - `beta` : Rate parameter (beta > 0). .. note:: - :math:`E(X) = 1/\beta` - :math:`Var(X) = 1/\beta^2` - PyMC's beta is named 'lambda' by Wikipedia, SciPy, Wolfram MathWorld and other sources. """ return flib.gamma(x, 1, beta)
exponential_grad_like = {'value': lambda x, beta: flib.gamma_grad_x(x, 1.0, beta), 'beta': lambda x, beta: flib.gamma_grad_beta(x, 1.0, beta)} # Exponentiated Weibull----------------------------------- @randomwrap def rexponweib(alpha, k, loc=0, scale=1, size=None): """ Random exponentiated Weibull variates. """ q = np.random.uniform(size=size) r = flib.exponweib_ppf(q, alpha, k) return loc + r * scale def exponweib_expval(alpha, k, loc, scale): # Not sure how we can do this, since the first moment is only # tractable at particular values of k raise NotImplementedError('exponweib_expval has not been implemented yet.')
[docs]def exponweib_like(x, alpha, k, loc=0, scale=1): R""" Exponentiated Weibull log-likelihood. The exponentiated Weibull distribution is a generalization of the Weibull family. Its value lies in being able to model monotone and non-monotone failure rates. .. math:: f(x \mid \alpha,k,loc,scale) & = \frac{\alpha k}{scale} (1-e^{-z^k})^{\alpha-1} e^{-z^k} z^{k-1} \\ z & = \frac{x-loc}{scale} :Parameters: - `x` : x > 0 - `alpha` : Shape parameter - `k` : k > 0 - `loc` : Location parameter - `scale` : Scale parameter (scale > 0). """ return flib.exponweib(x, alpha, k, loc, scale)
""" commented out because tests fail exponweib_grad_like = {'value' : flib.exponweib_gx, 'alpha' : flib.exponweib_ga, 'k' : flib.exponweib_gk, 'loc' : flib.exponweib_gl, 'scale' : flib.exponweib_gs} """ # Gamma---------------------------------------------- @randomwrap def rgamma(alpha, beta, size=None): """ Random gamma variates. """ return np.random.gamma(shape=alpha, scale=1. / beta, size=size) def gamma_expval(alpha, beta): """ Expected value of gamma distribution. """ return 1. * np.asarray(alpha) / beta
[docs]def gamma_like(x, alpha, beta): R""" Gamma log-likelihood. Represents the sum of alpha exponentially distributed random variables, each of which has rate parameter beta. .. math:: f(x \mid \alpha, \beta) = \frac{\beta^{\alpha}x^{\alpha-1}e^{-\beta x}}{\Gamma(\alpha)} :Parameters: - `x` : math:`x \ge 0` - `alpha` : Shape parameter (alpha > 0). - `beta` : Rate parameter (beta > 0). .. note:: - :math:`E(X) = \frac{\alpha}{\beta}` - :math:`Var(X) = \frac{\alpha}{\beta^2}` """ return flib.gamma(x, alpha, beta)
gamma_grad_like = {'value': flib.gamma_grad_x, 'alpha': flib.gamma_grad_alpha, 'beta': flib.gamma_grad_beta} # GEV Generalized Extreme Value ------------------------ # Modify parameterization -> Hosking (kappa, xi, alpha) @randomwrap def rgev(xi, mu=0, sigma=1, size=None): """ Random generalized extreme value (GEV) variates. """ q = np.random.uniform(size=size) z = flib.gev_ppf(q, xi) return z * sigma + mu def gev_expval(xi, mu=0, sigma=1): """ Expected value of generalized extreme value distribution. """ return mu - (sigma / xi) + (sigma / xi) * flib.gamfun(1 - xi) def gev_like(x, xi, mu=0, sigma=1): R""" Generalized Extreme Value log-likelihood .. math:: pdf(x \mid \xi,\mu,\sigma) = \frac{1}{\sigma}(1 + \xi \left[\frac{x-\mu}{\sigma}\right])^{-1/\xi-1}\exp{-(1+\xi \left[\frac{x-\mu}{\sigma}\right])^{-1/\xi}} .. math:: \sigma & > 0,\\ x & > \mu-\sigma/\xi \text{ if } \xi > 0,\\ x & < \mu-\sigma/\xi \text{ if } \xi < 0\\ x & \in [-\infty,\infty] \text{ if } \xi = 0 """ return flib.gev(x, xi, mu, sigma) # Geometric---------------------------------------------- # Changed the return value @randomwrap def rgeometric(p, size=None): """ Random geometric variates. """ return np.random.geometric(p, size) def geometric_expval(p): """ Expected value of geometric distribution. """ return 1. / p
[docs]def geometric_like(x, p): R""" Geometric log-likelihood. The probability that the first success in a sequence of Bernoulli trials occurs on the x'th trial. .. math:: f(x \mid p) = p(1-p)^{x-1} :Parameters: - `x` : [int] Number of trials before first success (x > 0). - `p` : Probability of success on an individual trial, :math:`p \in [0,1]` .. note:: - :math:`E(X)=1/p` - :math:`Var(X)=\frac{1-p}{p^2}` """ return flib.geometric(x, p)
geometric_grad_like = {'p': flib.geometric_gp} # Half Cauchy---------------------------------------------- @randomwrap def rhalf_cauchy(alpha, beta, size=None): """ Returns half-Cauchy random variates. """ return abs(alpha + beta * np.tan(pi * random_number(size) - pi / 2.0)) def half_cauchy_expval(alpha, beta): """ Expected value of cauchy distribution is undefined. """ return inf # In wikipedia, the arguments name are k, x0.
[docs]def half_cauchy_like(x, alpha, beta): R""" Half-Cauchy log-likelihood. Simply the absolute value of Cauchy. .. math:: f(x \mid \alpha, \beta) = \frac{2}{\pi \beta [1 + (\frac{x-\alpha}{\beta})^2]} :Parameters: - `alpha` : Location parameter. - `beta` : Scale parameter (beta > 0). .. note:: - x must be non-negative. """ x = np.atleast_1d(x) if sum(x.ravel() < 0): return -inf return flib.cauchy(x, alpha, beta) + len(x) * np.log(2)
# Half-normal---------------------------------------------- @randomwrap def rhalf_normal(tau, size=None): """ Random half-normal variates. """ return abs(np.random.normal(0, np.sqrt(1 / tau), size)) def half_normal_expval(tau): """ Expected value of half normal distribution. """ return np.sqrt(2. * pi / np.asarray(tau))
[docs]def half_normal_like(x, tau): R""" Half-normal log-likelihood, a normal distribution with mean 0 limited to the domain :math:`x \in [0, \infty)`. .. math:: f(x \mid \tau) = \sqrt{\frac{2\tau}{\pi}}\exp\left\{ {\frac{-x^2 \tau}{2}}\right\} :Parameters: - `x` : :math:`x \ge 0` - `tau` : tau > 0 """ return flib.hnormal(x, tau)
half_normal_grad_like = {'value': flib.hnormal_gradx, 'tau': flib.hnormal_gradtau} # Hypergeometric---------------------------------------------- def rhypergeometric(n, m, N, size=None): """ Returns hypergeometric random variates. """ if n == 0: return np.zeros(size, dtype=int) elif n == N: out = np.empty(size, dtype=int) out.fill(m) return out return np.random.hypergeometric(n, N - n, m, size) def hypergeometric_expval(n, m, N): """ Expected value of hypergeometric distribution. """ return 1. * n * m / N
[docs]def hypergeometric_like(x, n, m, N): R""" Hypergeometric log-likelihood. Discrete probability distribution that describes the number of successes in a sequence of draws from a finite population without replacement. .. math:: f(x \mid n, m, N) = \frac{\left({ \begin{array}{c} {m} \\ {x} \\ \end{array} }\right)\left({ \begin{array}{c} {N-m} \\ {n-x} \\ \end{array}}\right)}{\left({ \begin{array}{c} {N} \\ {n} \\ \end{array}}\right)} :Parameters: - `x` : [int] Number of successes in a sample drawn from a population. - `n` : [int] Size of sample drawn from the population. - `m` : [int] Number of successes in the population. - `N` : [int] Total number of units in the population. .. note:: :math:`E(X) = \frac{n n}{N}` """ return flib.hyperg(x, n, m, N)
# Inverse gamma---------------------------------------------- @randomwrap def rinverse_gamma(alpha, beta, size=None): """ Random inverse gamma variates. """ return 1. / np.random.gamma(shape=alpha, scale=1. / beta, size=size) def inverse_gamma_expval(alpha, beta): """ Expected value of inverse gamma distribution. """ return 1. * np.asarray(beta) / (alpha - 1.)
[docs]def inverse_gamma_like(x, alpha, beta): R""" Inverse gamma log-likelihood, the reciprocal of the gamma distribution. .. math:: f(x \mid \alpha, \beta) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{-\alpha - 1} \exp\left(\frac{-\beta}{x}\right) :Parameters: - `x` : x > 0 - `alpha` : Shape parameter (alpha > 0). - `beta` : Scale parameter (beta > 0). .. note:: :math:`E(X)=\frac{\beta}{\alpha-1}` for :math:`\alpha > 1` :math:`Var(X)=\frac{\beta^2}{(\alpha-1)^2(\alpha-2)}` for :math:`\alpha > 2` """ return flib.igamma(x, alpha, beta)
inverse_gamma_grad_like = {'value': flib.igamma_grad_x, 'alpha': flib.igamma_grad_alpha, 'beta': flib.igamma_grad_beta} # Inverse Wishart--------------------------------------------------- # def rinverse_wishart(n, C): # """ # Return an inverse Wishart random matrix. # # :Parameters: # - `n` : [int] Degrees of freedom (n > 0). # - `C` : Symmetric and positive definite scale matrix # """ # wi = rwishart(n, np.asmatrix(C).I).I # flib.symmetrize(wi) # return wi # # def inverse_wishart_expval(n, C): # """ # Expected value of inverse Wishart distribution. # # :Parameters: # - `n` : [int] Degrees of freedom (n > 0). # - `C` : Symmetric and positive definite scale matrix # # """ # return np.asarray(C)/(n-len(C)-1) # # def inverse_wishart_like(X, n, C): # R""" # Inverse Wishart log-likelihood. The inverse Wishart distribution # is the conjugate prior for the covariance matrix of a multivariate # normal distribution. # # .. math:: # f(X \mid n, T) = \frac{{\mid T \mid}^{n/2}{\mid X # \mid}^{-(n+k+1)/2} \exp\left\{ -\frac{1}{2} Tr(TX^{-1}) # \right\}}{2^{nk/2} \Gamma_p(n/2)} # # where :math:`k` is the rank of X. # # :Parameters: # - `X` : Symmetric, positive definite matrix. # - `n` : [int] Degrees of freedom (n > 0). # - `C` : Symmetric and positive definite scale matrix # # .. note:: # Step method MatrixMetropolis will preserve the symmetry of # Wishart variables. # # """ # return flib.blas_wishart(X.I, n, C.I) # Double exponential (Laplace)-------------------------------------------- @randomwrap def rlaplace(mu, tau, size=None): """ Laplace (double exponential) random variates. """ u = np.random.uniform(-0.5, 0.5, size) return mu - np.sign(u) * np.log(1 - 2 * np.abs(u)) / tau rdexponential = rlaplace def laplace_expval(mu, tau): """ Expected value of Laplace (double exponential) distribution. """ return mu dexponential_expval = laplace_expval
[docs]def laplace_like(x, mu, tau): R""" Laplace (double exponential) log-likelihood. The Laplace (or double exponential) distribution describes the difference between two independent, identically distributed exponential events. It is often used as a heavier-tailed alternative to the normal. .. math:: f(x \mid \mu, \tau) = \frac{\tau}{2}e^{-\tau |x-\mu|} :Parameters: - `x` : :math:`-\infty < x < \infty` - `mu` : Location parameter :math:`-\infty < mu < \infty` - `tau` : Scale parameter :math:`\tau > 0` .. note:: - :math:`E(X) = \mu` - :math:`Var(X) = \frac{2}{\tau^2}` """ return flib.gamma(np.abs(np.array(x) - mu), 1, tau) - \ np.size(x) * np.log(2)
laplace_grad_like = {'value': lambda x, mu, tau: flib.gamma_grad_x(np.abs(x - mu), 1, tau) * np.sign(x - mu), 'mu': lambda x, mu, tau: -flib.gamma_grad_x(np.abs(x - mu), 1, tau) * np.sign(x - mu), 'tau': lambda x, mu, tau: flib.gamma_grad_beta(np.abs(x - mu), 1, tau)} dexponential_like = laplace_like dexponential_grad_like = laplace_grad_like # Logistic----------------------------------- @randomwrap def rlogistic(mu, tau, size=None): """ Logistic random variates. """ u = np.random.random(size) return mu + np.log(u / (1 - u)) / tau def logistic_expval(mu, tau): """ Expected value of logistic distribution. """ return mu
[docs]def logistic_like(x, mu, tau): R""" Logistic log-likelihood. The logistic distribution is often used as a growth model; for example, populations, markets. Resembles a heavy-tailed normal distribution. .. math:: f(x \mid \mu, tau) = \frac{\tau \exp(-\tau[x-\mu])}{[1 + \exp(-\tau[x-\mu])]^2} :Parameters: - `x` : :math:`-\infty < x < \infty` - `mu` : Location parameter :math:`-\infty < mu < \infty` - `tau` : Scale parameter (tau > 0) .. note:: - :math:`E(X) = \mu` - :math:`Var(X) = \frac{\pi^2}{3\tau^2}` """ return flib.logistic(x, mu, tau)
# Lognormal---------------------------------------------- @randomwrap def rlognormal(mu, tau, size=None): """ Return random lognormal variates. """ return np.random.lognormal(mu, np.sqrt(1. / tau), size) def lognormal_expval(mu, tau): """ Expected value of log-normal distribution. """ return np.exp(mu + 1. / 2 / tau)
[docs]def lognormal_like(x, mu, tau): R""" Log-normal log-likelihood. Distribution of any random variable whose logarithm is normally distributed. A variable might be modeled as log-normal if it can be thought of as the multiplicative product of many small independent factors. .. math:: f(x \mid \mu, \tau) = \sqrt{\frac{\tau}{2\pi}}\frac{ \exp\left\{ -\frac{\tau}{2} (\ln(x)-\mu)^2 \right\}}{x} :Parameters: - `x` : x > 0 - `mu` : Location parameter. - `tau` : Scale parameter (tau > 0). .. note:: :math:`E(X)=e^{\mu+\frac{1}{2\tau}}` :math:`Var(X)=(e^{1/\tau}-1)e^{2\mu+\frac{1}{\tau}}` """ return flib.lognormal(x, mu, tau)
lognormal_grad_like = {'value': flib.lognormal_gradx, 'mu': flib.lognormal_gradmu, 'tau': flib.lognormal_gradtau} # Multinomial---------------------------------------------- #@randomwrap def rmultinomial(n, p, size=None): """ Random multinomial variates. """ # Leaving size=None as the default means return value is 1d array # if not specified-- nicer. # Single value for p: if len(np.shape(p)) == 1: return np.random.multinomial(n, p, size) # Multiple values for p: if np.isscalar(n): n = n * np.ones(np.shape(p)[0], dtype=np.int) out = np.empty(np.shape(p)) for i in xrange(np.shape(p)[0]): out[i, :] = np.random.multinomial(n[i], p[i,:], size) return out def multinomial_expval(n, p): """ Expected value of multinomial distribution. """ return np.asarray([pr * n for pr in p])
[docs]def multinomial_like(x, n, p): R""" Multinomial log-likelihood. Generalization of the binomial distribution, but instead of each trial resulting in "success" or "failure", each one results in exactly one of some fixed finite number k of possible outcomes over n independent trials. 'x[i]' indicates the number of times outcome number i was observed over the n trials. .. math:: f(x \mid n, p) = \frac{n!}{\prod_{i=1}^k x_i!} \prod_{i=1}^k p_i^{x_i} :Parameters: x : (ns, k) int Random variable indicating the number of time outcome i is observed. :math:`\sum_{i=1}^k x_i=n`, :math:`x_i \ge 0`. n : int Number of trials. p : (k,) Probability of each one of the different outcomes. :math:`\sum_{i=1}^k p_i = 1)`, :math:`p_i \ge 0`. .. note:: - :math:`E(X_i)=n p_i` - :math:`Var(X_i)=n p_i(1-p_i)` - :math:`Cov(X_i,X_j) = -n p_i p_j` - If :math:`\sum_i p_i < 0.999999` a log-likelihood value of -inf will be returned. """ # flib expects 2d arguments. Do we still want to support multiple p # values along realizations ? x = np.atleast_2d(x) p = np.atleast_2d(p) return flib.multinomial(x, n, p)
# Multivariate hypergeometric------------------------------ def rmultivariate_hypergeometric(n, m, size=None): """ Random multivariate hypergeometric variates. Parameters: - `n` : Number of draws. - `m` : Number of items in each categoy. """ N = len(m) urn = np.repeat(np.arange(N), m) if size: draw = np.array([[urn[i] for i in np.random.permutation(len(urn))[:n]] for j in range(size)]) r = [[np.sum(draw[j] == i) for i in range(len(m))] for j in range(size)] else: draw = np.array([urn[i] for i in np.random.permutation(len(urn))[:n]]) r = [np.sum(draw == i) for i in range(len(m))] return np.asarray(r) def multivariate_hypergeometric_expval(n, m): """ Expected value of multivariate hypergeometric distribution. Parameters: - `n` : Number of draws. - `m` : Number of items in each categoy. """ m = np.asarray(m, float) return n * (m / m.sum())
[docs]def multivariate_hypergeometric_like(x, m): R""" Multivariate hypergeometric log-likelihood Describes the probability of drawing x[i] elements of the ith category, when the number of items in each category is given by m. .. math:: \frac{\prod_i \left({ \begin{array}{c} {m_i} \\ {x_i} \\ \end{array}}\right)}{\left({ \begin{array}{c} {N} \\ {n} \\ \end{array}}\right)} where :math:`N = \sum_i m_i` and :math:`n = \sum_i x_i`. :Parameters: - `x` : [int sequence] Number of draws from each category, (x < m). - `m` : [int sequence] Number of items in each categoy. """ return flib.mvhyperg(x, m)
# Multivariate normal-------------------------------------- def rmv_normal(mu, tau, size=1): """ Random multivariate normal variates. """ sig = np.linalg.cholesky(tau) mu_size = np.shape(mu) if size == 1: out = np.random.normal(size=mu_size) try: flib.dtrsm_wrap(sig, out, 'L', 'T', 'L', 1.) except: out = np.linalg.solve(sig, out) out += mu return out else: if not hasattr(size, '__iter__'): size = (size,) tot_size = np.prod(size) out = np.random.normal(size=(tot_size,) + mu_size) for i in xrange(tot_size): try: flib.dtrsm_wrap(sig, out[i, :], 'L', 'T', 'L', 1.) except: out[i, :] = np.linalg.solve(sig, out[i,:]) out[i, :] += mu return out.reshape(size + mu_size) def mv_normal_expval(mu, tau): """ Expected value of multivariate normal distribution. """ return mu
[docs]def mv_normal_like(x, mu, tau): R""" Multivariate normal log-likelihood .. math:: f(x \mid \pi, T) = \frac{|T|^{1/2}}{(2\pi)^{1/2}} \exp\left\{ -\frac{1}{2} (x-\mu)^{\prime}T(x-\mu) \right\} :Parameters: - `x` : (n,k) - `mu` : (k) Location parameter sequence. - `Tau` : (k,k) Positive definite precision matrix. .. seealso:: :func:`mv_normal_chol_like`, :func:`mv_normal_cov_like` """ # TODO: Vectorize in Fortran if len(np.shape(x)) > 1: return np.sum([flib.prec_mvnorm(r, mu, tau) for r in x]) else: return flib.prec_mvnorm(x, mu, tau)
# Multivariate normal, parametrized with covariance--------------------------- def rmv_normal_cov(mu, C, size=1): """ Random multivariate normal variates. """ mu_size = np.shape(mu) if size == 1: return np.random.multivariate_normal(mu, C, size).reshape(mu_size) else: return np.random.multivariate_normal( mu, C, size).reshape((size,) + mu_size) def mv_normal_cov_expval(mu, C): """ Expected value of multivariate normal distribution. """ return mu
[docs]def mv_normal_cov_like(x, mu, C): R""" Multivariate normal log-likelihood parameterized by a covariance matrix. .. math:: f(x \mid \pi, C) = \frac{1}{(2\pi|C|)^{1/2}} \exp\left\{ -\frac{1}{2} (x-\mu)^{\prime}C^{-1}(x-\mu) \right\} :Parameters: - `x` : (n,k) - `mu` : (k) Location parameter. - `C` : (k,k) Positive definite covariance matrix. .. seealso:: :func:`mv_normal_like`, :func:`mv_normal_chol_like` """ # TODO: Vectorize in Fortran if len(np.shape(x)) > 1: return np.sum([flib.cov_mvnorm(r, mu, C) for r in x]) else: return flib.cov_mvnorm(x, mu, C)
# Multivariate normal, parametrized with Cholesky factorization.---------- def rmv_normal_chol(mu, sig, size=1): """ Random multivariate normal variates. """ mu_size = np.shape(mu) if size == 1: out = np.random.normal(size=mu_size) try: flib.dtrmm_wrap(sig, out, 'L', 'N', 'L', 1.) except: out = np.dot(sig, out) out += mu return out else: if not hasattr(size, '__iter__'): size = (size,) tot_size = np.prod(size) out = np.random.normal(size=(tot_size,) + mu_size) for i in xrange(tot_size): try: flib.dtrmm_wrap(sig, out[i, :], 'L', 'N', 'L', 1.) except: out[i, :] = np.dot(sig, out[i,:]) out[i, :] += mu return out.reshape(size + mu_size) def mv_normal_chol_expval(mu, sig): """ Expected value of multivariate normal distribution. """ return mu
[docs]def mv_normal_chol_like(x, mu, sig): R""" Multivariate normal log-likelihood. .. math:: f(x \mid \pi, \sigma) = \frac{1}{(2\pi)^{1/2}|\sigma|)} \exp\left\{ -\frac{1}{2} (x-\mu)^{\prime}(\sigma \sigma^{\prime})^{-1}(x-\mu) \right\} :Parameters: - `x` : (n,k) - `mu` : (k) Location parameter. - `sigma` : (k,k) Lower triangular matrix. .. seealso:: :func:`mv_normal_like`, :func:`mv_normal_cov_like` """ # TODO: Vectorize in Fortran if len(np.shape(x)) > 1: return np.sum([flib.chol_mvnorm(r, mu, sig) for r in x]) else: return flib.chol_mvnorm(x, mu, sig)
# Negative binomial---------------------------------------- @randomwrap def rnegative_binomial(mu, alpha, size=None): """ Random negative binomial variates. """ # Using gamma-poisson mixture rather than numpy directly # because numpy apparently rounds mu = np.asarray(mu, dtype=float) pois_mu = np.random.gamma(alpha, mu / alpha, size) return np.random.poisson(pois_mu, size) # return np.random.negative_binomial(alpha, alpha / (mu + alpha), size) def negative_binomial_expval(mu, alpha): """ Expected value of negative binomial distribution. """ return mu
[docs]def negative_binomial_like(x, mu, alpha): R""" Negative binomial log-likelihood. The negative binomial distribution describes a Poisson random variable whose rate parameter is gamma distributed. PyMC's chosen parameterization is based on this mixture interpretation. .. math:: f(x \mid \mu, \alpha) = \frac{\Gamma(x+\alpha)}{x! \Gamma(\alpha)} (\alpha/(\mu+\alpha))^\alpha (\mu/(\mu+\alpha))^x :Parameters: - `x` : x = 0,1,2,... - `mu` : mu > 0 - `alpha` : alpha > 0 .. note:: - :math:`E[x]=\mu` - In Wikipedia's parameterization, :math:`r=\alpha` :math:`p=\alpha/(\mu+\alpha)` :math:`\mu=rp/(1-p)` """ alpha = np.array(alpha) if (alpha > 1e10).any(): if (alpha > 1e10).all(): # Return Poisson when alpha gets very large return flib.poisson(x, mu) # Split big and small dispersion values big = alpha > 1e10 return flib.poisson(x[big], mu[big]) + flib.negbin2(x[big - True], mu[big - True], alpha[big - True]) return flib.negbin2(x, mu, alpha)
negative_binomial_grad_like = {'mu': flib.negbin2_gmu, 'alpha': flib.negbin2_ga} # Normal--------------------------------------------------- @randomwrap def rnormal(mu, tau, size=None): """ Random normal variates. """ return np.random.normal(mu, 1. / np.sqrt(tau), size) def normal_expval(mu, tau): """ Expected value of normal distribution. """ return mu
[docs]def normal_like(x, mu, tau): R""" Normal log-likelihood. .. math:: f(x \mid \mu, \tau) = \sqrt{\frac{\tau}{2\pi}} \exp\left\{ -\frac{\tau}{2} (x-\mu)^2 \right\} :Parameters: - `x` : Input data. - `mu` : Mean of the distribution. - `tau` : Precision of the distribution, which corresponds to :math:`1/\sigma^2` (tau > 0). .. note:: - :math:`E(X) = \mu` - :math:`Var(X) = 1/\tau` """ # try: # constrain(tau, lower=0) # except ZeroProbability: # return -np.Inf return flib.normal(x, mu, tau)
def t_normal_grad_x(x, mu, tau): return flib.normal_grad_x(x, mu, tau) def t_normal_grad_mu(x, mu, tau): return flib.normal_grad_mu(x, mu, tau) def t_normal_grad_tau(x, mu, tau): return flib.normal_grad_tau(x, mu, tau) normal_grad_like = {'value': t_normal_grad_x, 'mu': t_normal_grad_mu, 'tau': t_normal_grad_tau} # normal_grad_like = {'x' : flib.normal_grad_x, # 'mu' : flib.normal_grad_mu, # 'tau' : flib.normal_grad_tau} # von Mises-------------------------------------------------- @randomwrap def rvon_mises(mu, kappa, size=None): """ Random von Mises variates. """ # TODO: Just return straight from numpy after release 1.3 return (np.random.mtrand.vonmises( mu, kappa, size) + np.pi) % (2. * np.pi) - np.pi def von_mises_expval(mu, kappa): """ Expected value of von Mises distribution. """ return mu
[docs]def von_mises_like(x, mu, kappa): R""" von Mises log-likelihood. .. math:: f(x \mid \mu, k) = \frac{e^{k \cos(x - \mu)}}{2 \pi I_0(k)} where `I_0` is the modified Bessel function of order 0. :Parameters: - `x` : Input data. - `mu` : Mean of the distribution. - `kappa` : Dispersion of the distribution .. note:: - :math:`E(X) = \mu` """ return flib.vonmises(x, mu, kappa)
# Pareto--------------------------------------------------- @randomwrap def rpareto(alpha, m, size=None): """ Random Pareto variates. """ return m / (random_number(size) ** (1. / alpha)) def pareto_expval(alpha, m): """ Expected value of Pareto distribution. """ if alpha <= 1: return inf return alpha * m / (alpha - 1)
[docs]def pareto_like(x, alpha, m): R""" Pareto log-likelihood. The Pareto is a continuous, positive probability distribution with two parameters. It is often used to characterize wealth distribution, or other examples of the 80/20 rule. .. math:: f(x \mid \alpha, m) = \frac{\alpha m^{\alpha}}{x^{\alpha+1}} :Parameters: - `x` : Input data (x > m) - `alpha` : Shape parameter (alpha>0) - `m` : Scale parameter (m>0) .. note:: - :math:`E(x)=\frac{\alpha m}{\alpha-1} if \alpha > 1` - :math:`Var(x)=\frac{m^2 \alpha}{(\alpha-1)^2(\alpha-2)} if \alpha > 2` """ return flib.pareto(x, alpha, m)
# Truncated Pareto--------------------------------------------------- @randomwrap def rtruncated_pareto(alpha, m, b, size=None): """ Random bounded Pareto variates. """ u = random_number(size) return (-(u * b ** alpha - u * m ** alpha - b ** alpha) / (b ** alpha * m ** alpha)) ** (-1. / alpha) def truncated_pareto_expval(alpha, m, b): """ Expected value of truncated Pareto distribution. """ if alpha <= 1: return inf part1 = (m ** alpha) / (1. - (m / b) ** alpha) part2 = 1. * alpha / (alpha - 1) part3 = (1. / (m ** (alpha - 1)) - 1. / (b ** (alpha - 1.))) return part1 * part2 * part3
[docs]def truncated_pareto_like(x, alpha, m, b): R""" Truncated Pareto log-likelihood. The Pareto is a continuous, positive probability distribution with two parameters. It is often used to characterize wealth distribution, or other examples of the 80/20 rule. .. math:: f(x \mid \alpha, m, b) = \frac{\alpha m^{\alpha} x^{-\alpha}}{1-(m/b)**{\alpha}} :Parameters: - `x` : Input data (x > m) - `alpha` : Shape parameter (alpha>0) - `m` : Scale parameter (m>0) - `b` : Upper bound (b>m) """ return flib.truncated_pareto(x, alpha, m, b)
# Poisson-------------------------------------------------- @randomwrap def rpoisson(mu, size=None): """ Random poisson variates. """ return np.random.poisson(mu, size) def poisson_expval(mu): """ Expected value of Poisson distribution. """ return mu
[docs]def poisson_like(x, mu): R""" Poisson log-likelihood. The Poisson is a discrete probability distribution. It is often used to model the number of events occurring in a fixed period of time when the times at which events occur are independent. The Poisson distribution can be derived as a limiting case of the binomial distribution. .. math:: f(x \mid \mu) = \frac{e^{-\mu}\mu^x}{x!} :Parameters: - `x` : [int] :math:`x \in {0,1,2,...}` - `mu` : Expected number of occurrences during the given interval, :math:`\mu \geq 0`. .. note:: - :math:`E(x)=\mu` - :math:`Var(x)=\mu` """ return flib.poisson(x, mu)
poisson_grad_like = {'mu': flib.poisson_gmu} # Truncated Poisson-------------------------------------------------- @randomwrap def rtruncated_poisson(mu, k, size=None): """ Random truncated Poisson variates with minimum value k, generated using rejection sampling. """ # Calculate m try: m = max(0, np.floor(k - mu)) except (TypeError, ValueError): # More than one mu return np.array([rtruncated_poisson(x, i) for x, i in zip(mu, np.resize(k, np.size(mu)))]).squeeze() k -= 1 # Calculate constant for acceptance probability C = np.exp(flib.factln(k + 1) - flib.factln(k + 1 - m)) # Empty array to hold random variates rvs = np.empty(0, int) total_size = np.prod(size or 1) while(len(rvs) < total_size): # Propose values by sampling from untruncated Poisson with mean mu + m proposals = np.random.poisson( mu + m, (total_size * 4, np.size(m))).squeeze() # Acceptance probability a = C * np.array([np.exp(flib.factln(y - m) - flib.factln(y)) for y in proposals]) a *= proposals > k # Uniform random variates u = np.random.random(total_size * 4) rvs = np.append(rvs, proposals[u < a]) return np.reshape(rvs[:total_size], size) def truncated_poisson_expval(mu, k): """ Expected value of Poisson distribution truncated to be no smaller than k. """ return mu / (1. - poiscdf(k, mu))
[docs]def truncated_poisson_like(x, mu, k): R""" Truncated Poisson log-likelihood. The Truncated Poisson is a discrete probability distribution that is arbitrarily truncated to be greater than some minimum value k. For example, zero-truncated Poisson distributions can be used to model counts that are constrained to be non-negative. .. math:: f(x \mid \mu, k) = \frac{e^{-\mu}\mu^x}{x!(1-F(k|\mu))} :Parameters: - `x` : [int] :math:`x \in {0,1,2,...}` - `mu` : Expected number of occurrences during the given interval, :math:`\mu \geq 0`. - `k` : Truncation point representing the minimum allowable value. .. note:: - :math:`E(x)=\frac{\mu}{1-F(k|\mu)}` - :math:`Var(x)=\frac{\mu}{1-F(k|\mu)}` """ return flib.trpoisson(x, mu, k)
truncated_poisson_grad_like = {'mu': flib.trpoisson_gmu} # Truncated normal distribution-------------------------- @randomwrap def rtruncated_normal(mu, tau, a=-np.inf, b=np.inf, size=None): """ Random truncated normal variates. """ sigma = 1. / np.sqrt(tau) na = utils.normcdf((a - mu) / sigma) nb = utils.normcdf((b - mu) / sigma) # Use the inverse CDF generation method. U = np.random.mtrand.uniform(size=size) q = U * nb + (1 - U) * na R = utils.invcdf(q) # Unnormalize return R * sigma + mu rtruncnorm = rtruncated_normal def truncated_normal_expval(mu, tau, a, b): """Expected value of the truncated normal distribution. .. math:: E(X) =\mu + \frac{\sigma(\varphi_1-\varphi_2)}{T} where .. math:: T & =\Phi\left(\frac{B-\mu}{\sigma}\right)-\Phi \left(\frac{A-\mu}{\sigma}\right)\text \\ \varphi_1 &= \varphi\left(\frac{A-\mu}{\sigma}\right) \\ \varphi_2 &= \varphi\left(\frac{B-\mu}{\sigma}\right) \\ and :math:`\varphi = N(0,1)` and :math:`tau & 1/sigma**2`. :Parameters: - `mu` : Mean of the distribution. - `tau` : Precision of the distribution, which corresponds to 1/sigma**2 (tau > 0). - `a` : Left bound of the distribution. - `b` : Right bound of the distribution. """ phia = np.exp(normal_like(a, mu, tau)) phib = np.exp(normal_like(b, mu, tau)) sigma = 1. / np.sqrt(tau) Phia = utils.normcdf((a - mu) / sigma) if b == np.inf: Phib = 1.0 else: Phib = utils.normcdf((b - mu) / sigma) return (mu + (phia - phib) / (Phib - Phia))[0] truncnorm_expval = truncated_normal_expval
[docs]def truncated_normal_like(x, mu, tau, a=None, b=None): R""" Truncated normal log-likelihood. .. math:: f(x \mid \mu, \tau, a, b) = \frac{\phi(\frac{x-\mu}{\sigma})} {\Phi(\frac{b-\mu}{\sigma}) - \Phi(\frac{a-\mu}{\sigma})}, where :math:`\sigma^2=1/\tau`, `\phi` is the standard normal PDF and `\Phi` is the standard normal CDF. :Parameters: - `x` : Input data. - `mu` : Mean of the distribution. - `tau` : Precision of the distribution, which corresponds to 1/sigma**2 (tau > 0). - `a` : Left bound of the distribution. - `b` : Right bound of the distribution. """ x = np.atleast_1d(x) if a is None: a = -np.inf a = np.atleast_1d(a) if b is None: b = np.inf b = np.atleast_1d(b) mu = np.atleast_1d(mu) sigma = (1. / np.atleast_1d(np.sqrt(tau))) if (x < a).any() or (x > b).any(): return -np.inf else: n = len(x) phi = normal_like(x, mu, tau) lPhia = utils.normcdf((a - mu) / sigma, log=True) lPhib = utils.normcdf((b - mu) / sigma, log=True) try: d = utils.log_difference(lPhib, lPhia) except ValueError: return -np.inf # d = np.log(Phib-Phia) if len(d) == n: Phi = d.sum() else: Phi = n * d if np.isnan(Phi) or np.isinf(Phi): return -np.inf return phi - Phi
truncnorm_like = truncated_normal_like # Azzalini's skew-normal----------------------------------- @randomwrap def rskew_normal(mu, tau, alpha, size=()): """ Skew-normal random variates. """ size_ = size or (1,) len_ = np.prod(size_) return flib.rskewnorm( len_, mu, tau, alpha, np.random.normal(size=2 * len_)).reshape(size) def skew_normal_expval(mu, tau, alpha): """ Expectation of skew-normal random variables. """ delta = alpha / np.sqrt(1. + alpha ** 2) return mu + np.sqrt(2 / pi / tau) * delta
[docs]def skew_normal_like(x, mu, tau, alpha): R""" Azzalini's skew-normal log-likelihood .. math:: f(x \mid \mu, \tau, \alpha) = 2 \Phi((x-\mu)\sqrt{\tau}\alpha) \phi(x,\mu,\tau) where :math:\Phi is the normal CDF and :math: \phi is the normal PDF. :Parameters: - `x` : Input data. - `mu` : Mean of the distribution. - `tau` : Precision of the distribution (> 0). - `alpha` : Shape parameter of the distribution. .. note:: See http://azzalini.stat.unipd.it/SN/ """ return flib.sn_like(x, mu, tau, alpha)
# Student's t----------------------------------- @randomwrap def rt(nu, size=None): """ Student's t random variates. """ return rnormal(0, 1, size) / np.sqrt(rchi2(nu, size) / nu) def t_expval(nu): """ Expectation of Student's t random variables. """ return 0
[docs]def t_like(x, nu): R""" Student's T log-likelihood. Describes a zero-mean normal variable whose precision is gamma distributed. Alternatively, describes the mean of several zero-mean normal random variables divided by their sample standard deviation. .. math:: f(x \mid \nu) = \frac{\Gamma(\frac{\nu+1}{2})}{\Gamma(\frac{\nu}{2}) \sqrt{\nu\pi}} \left( 1 + \frac{x^2}{\nu} \right)^{-\frac{\nu+1}{2}} :Parameters: - `x` : Input data. - `nu` : Degrees of freedom. """ nu = np.asarray(nu) return flib.t(x, nu)
# Non-central Student's t----------------------------------- @randomwrap def rnoncentral_t(mu, lam, nu, size=None): """ Non-central Student's t random variates. """ tau = rgamma(nu / 2., nu / (2. * lam), size) return rnormal(mu, tau) def noncentral_t_expval(mu, lam, nu): """noncentral_t_expval(mu, lam, nu) Expectation of non-central Student's t random variables. Only defined for nu>1. """ if nu > 1: return mu return inf
[docs]def noncentral_t_like(x, mu, lam, nu): R""" Non-central Student's T log-likelihood. Describes a normal variable whose precision is gamma distributed. .. math:: f(x|\mu,\lambda,\nu) = \frac{\Gamma(\frac{\nu + 1}{2})}{\Gamma(\frac{\nu}{2})} \left(\frac{\lambda}{\pi\nu}\right)^{\frac{1}{2}} \left[1+\frac{\lambda(x-\mu)^2}{\nu}\right]^{-\frac{\nu+1}{2}} :Parameters: - `x` : Input data. - `mu` : Location parameter. - `lambda` : Scale parameter. - `nu` : Degrees of freedom. """ mu = np.asarray(mu) lam = np.asarray(lam) nu = np.asarray(nu) return flib.nct(x, mu, lam, nu)
def t_grad_setup(x, nu, f): nu = np.asarray(nu) return f(x, nu) t_grad_like = {'value': lambda x, nu: t_grad_setup(x, nu, flib.t_grad_x), 'nu': lambda x, nu: t_grad_setup(x, nu, flib.t_grad_nu)} # DiscreteUniform-------------------------------------------------- @randomwrap def rdiscrete_uniform(lower, upper, size=None): """ Random discrete_uniform variates. """ return np.random.randint(lower, upper + 1, size) def discrete_uniform_expval(lower, upper): """ Expected value of discrete_uniform distribution. """ return (upper - lower) / 2.
[docs]def discrete_uniform_like(x, lower, upper): R""" Discrete uniform log-likelihood. .. math:: f(x \mid lower, upper) = \frac{1}{upper-lower+1} :Parameters: - `x` : [int] :math:`x \in \{lower, lower+1, \ldots, upper-1, upper\}` - `lower` : Lower limit. - `upper` : Upper limit (upper > lower). """ return flib.duniform_like(x, lower, upper)
# Uniform-------------------------------------------------- @randomwrap def runiform(lower, upper, size=None): """ Random uniform variates. """ return np.random.uniform(lower, upper, size) def uniform_expval(lower, upper): """ Expected value of uniform distribution. """ return (upper - lower) / 2.
[docs]def uniform_like(x, lower, upper): R""" Uniform log-likelihood. .. math:: f(x \mid lower, upper) = \frac{1}{upper-lower} :Parameters: - `x` : :math:`lower \leq x \leq upper` - `lower` : Lower limit. - `upper` : Upper limit (upper > lower). """ return flib.uniform_like(x, lower, upper)
uniform_grad_like = {'value': flib.uniform_grad_x, 'lower': flib.uniform_grad_l, 'upper': flib.uniform_grad_u} # Weibull-------------------------------------------------- @randomwrap def rweibull(alpha, beta, size=None): """ Weibull random variates. """ tmp = -np.log(runiform(0, 1, size)) return beta * (tmp ** (1. / alpha)) def weibull_expval(alpha, beta): """ Expected value of weibull distribution. """ return beta * np.exp(gammaln((alpha + 1.) / alpha))
[docs]def weibull_like(x, alpha, beta): R""" Weibull log-likelihood .. math:: f(x \mid \alpha, \beta) = \frac{\alpha x^{\alpha - 1} \exp(-(\frac{x}{\beta})^{\alpha})}{\beta^\alpha} :Parameters: - `x` : :math:`x \ge 0` - `alpha` : alpha > 0 - `beta` : beta > 0 .. note:: - :math:`E(x)=\beta \Gamma(1+\frac{1}{\alpha})` - :math:`Var(x)=\beta^2 \Gamma(1+\frac{2}{\alpha} - \mu^2)` """ return flib.weibull(x, alpha, beta)
weibull_grad_like = {'value': flib.weibull_gx, 'alpha': flib.weibull_ga, 'beta': flib.weibull_gb} # Wishart--------------------------------------------------- def rwishart(n, Tau): """ Return a Wishart random matrix. Tau is the inverse of the 'covariance' matrix :math:`C`. """ p = np.shape(Tau)[0] sig = np.linalg.cholesky(Tau) if n <= (p-1): raise ValueError('Wishart parameter n must be greater ' 'than size of matrix.') norms = np.random.normal(size=(p * (p - 1)) // 2) chi_sqs = np.sqrt(np.random.chisquare(df=np.arange(n, n - p, -1))) A = flib.expand_triangular(chi_sqs, norms) flib.dtrsm_wrap(sig, A, side='L', uplo='L', transa='T', alpha=1.) w = np.asmatrix(np.dot(A, A.T)) flib.symmetrize(w) return w def wishart_expval(n, Tau): """ Expected value of wishart distribution. """ return n * np.asarray(Tau.I)
[docs]def wishart_like(X, n, Tau): R""" Wishart log-likelihood. The Wishart distribution is the probability distribution of the maximum-likelihood estimator (MLE) of the precision matrix of a multivariate normal distribution. If Tau=1, the distribution is identical to the chi-square distribution with n degrees of freedom. For an alternative parameterization based on :math:`C=T{-1}`, see `wishart_cov_like`. .. math:: f(X \mid n, T) = \frac{{\mid T \mid}^{n/2}{\mid X \mid}^{(n-k-1)/2}}{2^{nk/2} \Gamma_p(n/2)} \exp\left\{ -\frac{1}{2} Tr(TX) \right\} where :math:`k` is the rank of X. :Parameters: X : matrix Symmetric, positive definite. n : int Degrees of freedom, > 0. Tau : matrix Symmetric and positive definite .. note:: Step method MatrixMetropolis will preserve the symmetry of Wishart variables. """ return flib.blas_wishart(X, n, Tau)
# Wishart, parametrized by covariance ------------------------------------ def rwishart_cov(n, C): """ Return a Wishart random matrix. :Parameters: n : int Degrees of freedom, > 0. C : matrix Symmetric and positive definite """ # return rwishart(n, np.linalg.inv(C)) p = np.shape(C)[0] # Need cholesky decomposition of precision matrix C^-1? sig = np.linalg.cholesky(C) if n <= (p-1): raise ValueError('Wishart parameter n must be greater ' 'than size of matrix.') norms = np.random.normal(size=(p * (p - 1)) // 2) chi_sqs = np.sqrt(np.random.chisquare(df=np.arange(n, n - p, -1))) A = flib.expand_triangular(chi_sqs, norms) flib.dtrmm_wrap(sig, A, side='L', uplo='L', transa='N', alpha=1.) w = np.asmatrix(np.dot(A, A.T)) flib.symmetrize(w) return w def wishart_cov_expval(n, C): """ Expected value of wishart distribution. :Parameters: n : int Degrees of freedom, > 0. C : matrix Symmetric and positive definite """ return n * np.asarray(C)
[docs]def wishart_cov_like(X, n, C): R""" wishart_like(X, n, C) Wishart log-likelihood. The Wishart distribution is the probability distribution of the maximum-likelihood estimator (MLE) of the covariance matrix of a multivariate normal distribution. If C=1, the distribution is identical to the chi-square distribution with n degrees of freedom. For an alternative parameterization based on :math:`T=C^{-1}`, see `wishart_like`. .. math:: f(X \mid n, C) = {\mid C^{-1} \mid}^{n/2}{\mid X \mid}^{(n-k-1)/2} \exp\left\{ -\frac{1}{2} Tr(C^{-1}X) \right\} where :math:`k` is the rank of X. :Parameters: X : matrix Symmetric, positive definite. n : int Degrees of freedom, > 0. C : matrix Symmetric and positive definite """ return flib.blas_wishart_cov(X, n, C)
# ----------------------------------------------------------- # DECORATORS # ----------------------------------------------------------- def name_to_funcs(name, module): if isinstance(module, types.ModuleType): module = copy(module.__dict__) elif isinstance(module, dict): module = copy(module) else: raise AttributeError try: logp = module[name + "_like"] except: raise KeyError("No likelihood found with this name ", name + "_like") try: random = module['r' + name] except: random = None try: grad_logp = module[name + "_grad_like"] except: grad_logp = {} return logp, random, grad_logp def valuewrapper(f, arguments=None): """Return a likelihood accepting value instead of x as a keyword argument. This is specifically intended for the instantiator above. """ def wrapper(**kwds): value = kwds.pop('value') return f(value, **kwds) if arguments is None: wrapper.__dict__.update(f.__dict__) else: wrapper.__dict__.update(arguments) return wrapper """ Decorate the likelihoods """ snapshot = locals().copy() likelihoods = {} for name, obj in six.iteritems(snapshot): if name[-5:] == '_like' and name[:-5] in availabledistributions: likelihoods[name[:-5]] = snapshot[name] def local_decorated_likelihoods(obj): """ New interface likelihoods """ for name, like in six.iteritems(likelihoods): obj[name + '_like'] = gofwrapper(like, snapshot) # local_decorated_likelihoods(locals()) # Decorating the likelihoods breaks the creation of distribution # instantiators -DH. # Create Stochastic instantiators def _inject_dist(distname, kwargs={}, ns=locals()): """ Reusable function to inject Stochastic subclasses into module namespace """ dist_logp, dist_random, grad_logp = name_to_funcs(distname, ns) classname = capitalize(distname) ns[classname] = stochastic_from_dist(distname, dist_logp, dist_random, grad_logp, **kwargs) for dist in sc_continuous_distributions: _inject_dist(dist) for dist in mv_continuous_distributions: _inject_dist(dist, kwargs={'mv': True}) for dist in sc_bool_distributions: _inject_dist(dist, kwargs={'dtype': np.bool}) for dist in sc_discrete_distributions: _inject_dist(dist, kwargs={'dtype': np.int}) for dist in mv_discrete_distributions: _inject_dist(dist, kwargs={'dtype': np.int, 'mv': True}) def uninformative_like(x): """ Uninformative log-likelihood. Returns 0 regardless of the value of x. """ return 0. def one_over_x_like(x): """ returns -np.Inf if x<0, -np.log(x) otherwise. """ if np.any(x < 0): return -np.Inf else: return -np.sum(np.log(x)) Uninformative = stochastic_from_dist('uninformative', logp=uninformative_like) DiscreteUninformative = stochastic_from_dist( 'uninformative', logp=uninformative_like, dtype=np.int) DiscreteUninformative.__name__ = 'DiscreteUninformative' OneOverX = stochastic_from_dist('one_over_x', logp=one_over_x_like) # Conjugates of Dirichlet get special treatment, can be parametrized # by first k-1 'p' values def extend_dirichlet(p): """ extend_dirichlet(p) Concatenates 1-sum(p) to the end of p and returns. """ if len(np.shape(p)) > 1: return np.hstack((p, np.atleast_2d(1. - np.sum(p)))) else: return np.hstack((p, 1. - np.sum(p))) def mod_categorical_like(x, p): """ Categorical log-likelihood with parent p of length k-1. An implicit k'th category is assumed to exist with associated probability 1-sum(p). ..math:: f(x=i \mid p, m, s) = p_i, ..math:: i \in 0\ldots k-1 :Parameters: x : integer :math: `x \in 0\ldots k-1` p : (k-1) float :math: `p > 0` :math: `\sum p < 1` minval : integer step : integer :math: `s \ge 1` """ return categorical_like(x, extend_dirichlet(p)) def mod_categorical_expval(p): """ Expected value of categorical distribution with parent p of length k-1. An implicit k'th category is assumed to exist with associated probability 1-sum(p). """ p = extend_dirichlet(p) return np.sum([p * i for i, p in enumerate(p)]) def rmod_categor(p, size=None): """ Categorical random variates with parent p of length k-1. An implicit k'th category is assumed to exist with associated probability 1-sum(p). """ return rcategorical(extend_dirichlet(p), size) class Categorical(Stochastic): __doc__ = """ C = Categorical(name, p, value=None, dtype=np.int, observed=False, size=1, trace=True, rseed=False, cache_depth=2, plot=None) Stochastic variable with Categorical distribution. Parent is: p If parent p is Dirichlet and has length k-1, an implicit k'th category is assumed to exist with associated probability 1-sum(p.value). Otherwise parent p's value should sum to 1. Docstring of categorical_like (case where P is not a Dirichlet): """\ + categorical_like.__doc__ +\ """ Docstring of categorical_like (case where P is a Dirichlet): """\ + categorical_like.__doc__ parent_names = ['p', 'minval', 'step'] def __init__(self, name, p, value=None, dtype=np.int, observed=False, size=None, trace=True, rseed=False, cache_depth=2, plot=None, verbose=-1, **kwds): if value is not None: if np.isscalar(value): self.size = None else: self.size = len(value) else: self.size = size if isinstance(p, Dirichlet): Stochastic.__init__(self, logp=valuewrapper(mod_categorical_like), doc='A Categorical random variable', name=name, parents={'p': p}, random=bind_size(rmod_categor, self.size), trace=trace, value=value, dtype=dtype, rseed=rseed, observed=observed, cache_depth=cache_depth, plot=plot, verbose=verbose, **kwds) else: Stochastic.__init__(self, logp=valuewrapper(categorical_like), doc='A Categorical random variable', name=name, parents={'p': p}, random=bind_size(rcategorical, self.size), trace=trace, value=value, dtype=dtype, rseed=rseed, observed=observed, cache_depth=cache_depth, plot=plot, verbose=verbose, **kwds) # class ModCategorical(Stochastic): # __doc__ = """ # C = ModCategorical(name, p, minval, step[, trace=True, value=None, # rseed=False, observed=False, cache_depth=2, plot=None, verbose=0]) # # Stochastic variable with ModCategorical distribution. # Parents are: p, minval, step. # # If parent p is Dirichlet and has length k-1, an implicit k'th # category is assumed to exist with associated probability 1-sum(p.value). # # Otherwise parent p's value should sum to 1. # # Docstring of mod_categorical_like (case where P is not a Dirichlet): # """\ # + mod_categorical_like.__doc__ +\ # """ # Docstring of mod_categorical_like (case where P is a Dirichlet): # """\ # + mod_categorical_like.__doc__ # # # parent_names = ['p', 'minval', 'step'] # # def __init__(self, name, p, minval=0, step=1, value=None, dtype=np.float, observed=False, size=1, trace=True, rseed=False, cache_depth=2, plot=None, verbose=0, **kwds): # # if value is not None: # if np.isscalar(value): # self.size = 1 # else: # self.size = len(value) # else: # self.size = size # # if isinstance(p, Dirichlet): # Stochastic.__init__(self, logp=valuewrapper(mod_categorical_like), doc='A ModCategorical random variable', name=name, # parents={'p':p,'minval':minval,'step':step}, random=bind_size(rmod_categor, self.size), trace=trace, value=value, dtype=dtype, # rseed=rseed, observed=observed, cache_depth=cache_depth, plot=plot, verbose=verbose, **kwds) # else: # Stochastic.__init__(self, logp=valuewrapper(mod_categorical_like), doc='A ModCategorical random variable', name=name, # parents={'p':p,'minval':minval,'step':step}, random=bind_size(rmod_categorical, self.size), trace=trace, value=value, dtype=dtype, # rseed=rseed, observed=observed, cache_depth=cache_depth, plot=plot, # verbose=verbose, **kwds) def mod_rmultinom(n, p): return rmultinomial(n, extend_dirichlet(p)) def mod_multinom_like(x, n, p): return multinomial_like(x, n, extend_dirichlet(p)) class Multinomial(Stochastic): """ M = Multinomial(name, n, p, trace=True, value=None, rseed=False, observed=False, cache_depth=2, plot=None]) A multinomial random variable. Parents are p, minval, step. If parent p is Dirichlet and has length k-1, an implicit k'th category is assumed to exist with associated probability 1-sum(p.value). Otherwise parent p's value should sum to 1. """ parent_names = ['n', 'p'] def __init__(self, name, n, p, trace=True, value=None, rseed=False, observed=False, cache_depth=2, plot=None, verbose=-1, **kwds): if isinstance(p, Dirichlet): Stochastic.__init__(self, logp=valuewrapper(mod_multinom_like), doc='A Multinomial random variable', name=name, parents={'n': n, 'p': p}, random=mod_rmultinom, trace=trace, value=value, dtype=np.int, rseed=rseed, observed=observed, cache_depth=cache_depth, plot=plot, verbose=verbose, **kwds) else: Stochastic.__init__(self, logp=valuewrapper(multinomial_like), doc='A Multinomial random variable', name=name, parents={'n': n, 'p': p}, random=rmultinomial, trace=trace, value=value, dtype=np.int, rseed=rseed, observed=observed, cache_depth=cache_depth, plot=plot, verbose=verbose, **kwds) def Impute(name, dist_class, imputable, **parents): """ This function accomodates missing elements for the data of simple Stochastic distribution subclasses. The masked_values argument is an object of type numpy.ma.MaskedArray, which contains the raw data and a boolean mask indicating missing values. The resulting list contains a list of stochastics of type dist_class, with the extant values as data stochastics and the missing values as variable stochastics. :Arguments: - name : string Name of the data stochastic - dist_class : Stochastic Stochastic subclass such as Poisson, Normal, etc. - imputable : numpy.ma.core.MaskedArray or iterable A masked array with missing elements (where mask=True, value is assumed missing), or any iterable that contains None elements that will be imputed. - parents (optional): dict Arbitrary keyword arguments. """ dims = np.shape(imputable) masked_values = np.ravel(imputable) if not isinstance(masked_values, np.ma.core.MaskedArray): # Generate mask mask = [v is None or np.isnan(v) for v in masked_values] # Generate masked array masked_values = np.ma.masked_array(masked_values, mask) # Initialise list vars = [] for i in xrange(len(masked_values)): # Name of element this_name = name + '[%i]' % i # Dictionary to hold parents these_parents = {} # Parse parents for key, parent in six.iteritems(parents): try: # If parent is a PyMCObject shape = np.shape(parent.value) except AttributeError: shape = np.shape(parent) if shape == dims: these_parents[key] = Lambda(key + '[%i]' % i, lambda p=np.ravel(parent), i=i: p[i]) elif shape == np.shape(masked_values): these_parents[key] = Lambda(key + '[%i]' % i, lambda p=parent, i=i: p[i]) else: these_parents[key] = parent if masked_values.mask[i]: # Missing values vars.append(dist_class(this_name, **these_parents)) else: # Observed values vars.append(dist_class(this_name, value=masked_values[i], observed=True, **these_parents)) return np.reshape(vars, dims) ImputeMissing = Impute if __name__ == "__main__": import doctest doctest.testmod()