carl.distributions module

def cdf(self, X, **kwargs):
    """Cumulative distribution function.
    This function returns the value of the cumulative distribution function
    `F(x_i) = p(x <= x_i)`, at all `x_i` in `X`.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `cdf` [array, shape=(n_samples,)]:
        `cdf(x_i)` for all `x_i` in `X`.
    """
    raise NotImplementedError

def fit(self, X, **kwargs):
    """Fit the distribution parameters to data.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `self` [object]:
        `self`.
    """
    return self

def nll(self, X, **kwargs):
    """Negative log-likelihood.
    This function returns the value of the negative log-likelihood
    `- log(p(x_i))`, at all `x_i` in `X`.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `nll` [array, shape=(n_samples,)]:
        `-log(p(x_i))` for all `x_i` in `X`.
    """
    raise NotImplementedError

def pdf(self, X, **kwargs):
    """Probability density function.
    This function returns the value of the probability density function
    `p(x_i)`, at all `x_i` in `X`.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `pdf` [array, shape=(n_samples,)]:
        `p(x_i)` for all `x_i` in `X`.
    """
    raise NotImplementedError

def ppf(self, X, **kwargs):
    """Percent point function (inverse of cdf).
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `ppf` [array, shape=(n_samples,)]:
        Percent point function for all `x_i` in `X`.
    """
    raise NotImplementedError

def rvs(self, n_samples, random_state=None, **kwargs):
    """Draw samples.
    Parameters
    ----------
    * `n_samples` [int]:
        The number of samples to draw.
    * `random_state` [integer or RandomState object]:
        The random seed.
    Returns
    -------
    * `X` [array, shape=(n_samples, n_features)]:
        The generated samples.
    """
    raise NotImplementedError

def score(self, X, **kwargs):
    """Evaluate the goodness of fit of the distribution w.r.t. `X`.
    The higher, the better.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `score` [float]:
        The goodness of fit.
    """
    return NotImplementedError

def __init__(self, inverse_scale=1.0):
    """Constructor.
    Parameters
    ----------
    * `inverse_scale` [float]:
        The inverse scale.
    """
    super(Exponential, self).__init__(inverse_scale=inverse_scale)
    # pdf
    self.pdf_ = T.switch(
        T.lt(self.X, 0.),
        0.,
        self.inverse_scale * T.exp(-self.inverse_scale * self.X)).ravel()
    self._make(self.pdf_, "pdf")
    # -log pdf
    self.nll_ = bound(
        -T.log(self.inverse_scale) + self.inverse_scale * self.X,
        np.inf,
        self.inverse_scale > 0.).ravel()
    self._make(self.nll_, "nll")
    # cdf
    self.cdf_ = (1. - T.exp(-self.inverse_scale * self.X)).ravel()
    self._make(self.cdf_, "cdf")
    # ppf
    self.ppf_ = -T.log(1. - self.p) / self.inverse_scale
    self._make(self.ppf_, "ppf", args=[self.p])

def cdf(self, X, **kwargs):
    """Cumulative distribution function.
    This function returns the value of the cumulative distribution function
    `F(x_i) = p(x <= x_i)`, at all `x_i` in `X`.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `cdf` [array, shape=(n_samples,)]:
        `cdf(x_i)` for all `x_i` in `X`.
    """
    raise NotImplementedError

def fit(self, X, bounds=None, constraints=None, use_gradient=True,
        optimizer=None, **kwargs):
    """Fit the distribution parameters to data by minimizing the negative
    log-likelihood of the data.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    * `bounds` [list of (parameter, (low, high))]:
        The parameter bounds.
    * `constraints`:
        The constraints on the parameters.
    * `use_gradient` [boolean, default=True]:
        Whether to use exact gradients (if `True`) or numerical gradients
        (if `False`).
    * `optimizer` [string]:
        The optimization method.
    Returns
    -------
    * `self` [object]:
        `self`.
    """
    # Map parameters to placeholders
    param_to_placeholder = []
    param_to_index = {}
    for i, v in enumerate(self.parameters_):
        w = T.TensorVariable(v.type)
        param_to_placeholder.append((v, w))
        param_to_index[v] = i
    # Build bounds
    mapped_bounds = None
    if bounds is not None:
        mapped_bounds = [(None, None) for v in param_to_placeholder]
        for b in bounds:
            mapped_bounds[param_to_index[b["param"]]] = b["bounds"]
    # Build constraints
    mapped_constraints = None
    if constraints is not None:
        mapped_constraints = []
        for c in constraints:
            args = c["param"]
            if isinstance(args, SharedVariable):
                args = (args, )
            m_c = {
                "type": c["type"],
                "fun": lambda x: c["fun"](*[x[param_to_index[a]]
                                            for a in args])
            }
            if "jac" in c:
                m_c["jac"] = lambda x: c["jac"](*[x[param_to_index[a]]
                                                  for a in args])
            mapped_constraints.append(m_c)
    # Derive objective and gradient
    objective_ = theano.function(
        [self.X] + [w for _, w in param_to_placeholder] +
        [theano.In(v, name=v.name) for v in self.observeds_],
        T.sum(self.nll_),
        givens=param_to_placeholder,
        allow_input_downcast=True)
    def objective(x):
        return objective_(X, *x, **kwargs) / len(X)
    if use_gradient:
        gradient_ = theano.function(
            [self.X] + [w for _, w in param_to_placeholder] +
            [theano.In(v, name=v.name) for v in self.observeds_],
            theano.grad(T.sum(self.nll_),
                        [v for v, _ in param_to_placeholder]),
            givens=param_to_placeholder,
            allow_input_downcast=True)
        def gradient(x):
            return np.array(gradient_(X, *x, **kwargs)) / len(X)
    # Solve!
    x0 = np.array([v.get_value() for v, _ in param_to_placeholder])
    r = minimize(objective,
                 jac=gradient if use_gradient else None,
                 x0=x0,
                 method=optimizer,
                 bounds=mapped_bounds,
                 constraints=mapped_constraints)
    if r.success:
        # Assign the solution
        for i, value in enumerate(r.x):
            param_to_placeholder[i][0].set_value(value)
    else:
        print("Parameter fitting failed!")
        print(r)
    return self

def nll(self, X, **kwargs):
    """Negative log-likelihood.
    This function returns the value of the negative log-likelihood
    `- log(p(x_i))`, at all `x_i` in `X`.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `nll` [array, shape=(n_samples,)]:
        `-log(p(x_i))` for all `x_i` in `X`.
    """
    raise NotImplementedError

def pdf(self, X, **kwargs):
    """Probability density function.
    This function returns the value of the probability density function
    `p(x_i)`, at all `x_i` in `X`.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `pdf` [array, shape=(n_samples,)]:
        `p(x_i)` for all `x_i` in `X`.
    """
    raise NotImplementedError

def ppf(self, X, **kwargs):
    """Percent point function (inverse of cdf).
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `ppf` [array, shape=(n_samples,)]:
        Percent point function for all `x_i` in `X`.
    """
    raise NotImplementedError

def rvs(self, n_samples, random_state=None, **kwargs):
    rng = check_random_state(random_state)
    p = rng.rand(n_samples, 1)
    return self.ppf(p, **kwargs)

def score(self, X, **kwargs):
    """Evaluate the log-likelihood `-self.nll(X).sum()` of `X`.
    The higher, the better.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `score` [float]:
        The log-likelihood of `X`.
    """
    return -self.nll(X, **kwargs).sum()

def set_params(self, **params):
    for name, value in params.items():
        var = getattr(self, name, None)
        if isinstance(var, SharedVariable):
            var.set_value(value)
        elif (isinstance(var, T.TensorVariable) or
              isinstance(var, T.TensorConstant)):
            raise ValueError("Only shared variables can be updated.")
        else:
            super(DistributionMixin, self).set_params(**{name: value})

def __init__(self, bins=10, range=None, interpolation=None,
             variable_width=False):
    """Constructor.
    Parameters
    ----------
    * `bins` [int or string]:
        The number of bins or `"blocks"` for bayesian blocks.
    * `range` [list of bounds (low, high)]:
        The boundaries. If `None`, bounds are inferred from data.
    * `interpolation` [string, optional]
        Specifies the kind of interpolation between bins as a string
        (`"linear"`, `"nearest"`, `"zero"`, `"slinear"`, `"quadratic"`,
        `"cubic"`).
    * `variable_width` [boolean, optional]
        If True use equal probability variable length bins.
    """
    self.bins = bins
    self.range = range
    self.interpolation = interpolation
    self.variable_width = variable_width

def cdf(self, X, **kwargs):
    """Not supported."""
    raise NotImplementedError

def fit(self, X, sample_weight=None, **kwargs):
    # Checks
    X = check_array(X)
    if sample_weight is not None and len(sample_weight) != len(X):
        raise ValueError
    # Compute histogram and edges
    if self.bins == "blocks":
        bins = bayesian_blocks(X.ravel(), fitness="events", p0=0.0001)
        range_ = self.range[0] if self.range else None
        h, e = np.histogram(X.ravel(), bins=bins, range=range_,
                            weights=sample_weight, normed=False)
        e = [e]
    elif self.variable_width:
        ticks = [np.percentile(X.ravel(), 100. * k / self.bins) for k
                 in range(self.bins + 1)]
        ticks[-1] += 1e-5
        range_ = self.range[0] if self.range else None
        h, e = np.histogram(X.ravel(), bins=ticks, range=range_,
                            normed=False, weights=sample_weight)
        h, e = h.astype(float), e.astype(float)
        widths = e[1:] - e[:-1]
        h = h / widths / h.sum()
        e = [e]
    else:
        bins = self.bins
        h, e = np.histogramdd(X, bins=bins, range=self.range,
                              weights=sample_weight, normed=True)
    # Add empty bins for out of bound samples
    for j in range(X.shape[1]):
        h = np.insert(h, 0, 0., axis=j)
        h = np.insert(h, h.shape[j], 0., axis=j)
        e[j] = np.insert(e[j], 0, -np.inf)
        e[j] = np.insert(e[j], len(e[j]), np.inf)
    if X.shape[1] == 1 and self.interpolation:
        inputs = e[0][2:-1] - (e[0][2] - e[0][1]) / 2.
        inputs[0] = e[0][1]
        inputs[-1] = e[0][-2]
        outputs = h[1:-1]
        self.interpolation_ = interp1d(inputs, outputs,
                                       kind=self.interpolation,
                                       bounds_error=False,
                                       fill_value=0.)
    self.histogram_ = h
    self.edges_ = e
    self.ndim_ = X.shape[1]
    return self

def nll(self, X, **kwargs):
    return -np.log(self.pdf(X, **kwargs))

def pdf(self, X, **kwargs):
    X = check_array(X)
    if self.ndim_ == 1 and self.interpolation:
        return self.interpolation_(X[:, 0])
    all_indices = []
    for j in range(X.shape[1]):
        indices = np.searchsorted(self.edges_[j],
                                  X[:, j],
                                  side="right") - 1
        # For the last bin, the upper is inclusive
        indices[X[:, j] == self.edges_[j][-2]] -= 1
        all_indices.append(indices)
    return self.histogram_[all_indices]

def ppf(self, X, **kwargs):
    """Not supported."""
    raise NotImplementedError

def rvs(self, n_samples, random_state=None, **kwargs):
    rng = check_random_state(random_state)
    # Draw random bins with respect to their densities
    h = (self.histogram_ / self.histogram_.sum()).ravel()
    flat_indices = np.searchsorted(np.cumsum(h),
                                   rng.rand(n_samples),
                                   side="right")
    # Build bin corners
    indices = np.unravel_index(flat_indices, self.histogram_.shape)
    indices_end = [a + 1 for a in indices]
    shape = [len(d) for d in self.edges_] + [len(self.edges_)]
    corners = np.array(list(product(*self.edges_))).reshape(shape)
    # Draw uniformly within bins
    low = corners[indices]
    high = corners[indices_end]
    u = rng.rand(*low.shape)
    return low + u * (high - low)

def score(self, X, **kwargs):
    """Evaluate the goodness of fit of the distribution w.r.t. `X`.
    The higher, the better.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `score` [float]:
        The goodness of fit.
    """
    return NotImplementedError

def __init__(self, components):
    """Constructor.
    Parameters
    ----------
    * `components` [list of `DistributionMixin`]:
        The components to join together.
    """
    super(Join, self).__init__()
    self.components = components
    for i, component in enumerate(components):
        # Add component parameters, constants and observeds
        if isinstance(component, TheanoDistribution):
            for p_i in component.parameters_:
                self.parameters_.add(p_i)
            for c_i in component.constants_:
                self.constants_.add(c_i)
            for o_i in component.observeds_:
                self.observeds_.add(o_i)
    # Derive and overide pdf and nll analytically if possible
    if all([hasattr(c, "pdf_") for c in self.components]):
        # pdf
        c0 = self.components[0]
        self.pdf_ = theano.clone(c0.pdf_, {c0.X: self.X[:, 0:c0.ndim]})
        start = c0.ndim
        for c in self.components[1:]:
            self.pdf_ *= theano.clone(
                c.pdf_, {c.X: self.X[:, start:start+c.ndim]})
            start += c.ndim
        self._make(self.pdf_, "pdf")
    if all([hasattr(c, "nll_") for c in self.components]):
        # nll
        c0 = self.components[0]
        self.nll_ = theano.clone(c0.nll_, {c0.X: self.X[:, 0:c0.ndim]})
        start = c0.ndim
        for c in self.components[1:]:
            self.nll_ += theano.clone(
                c.nll_, {c.X: self.X[:, start:start+c.ndim]})
            start += c.ndim
        self._make(self.nll_, "nll")

def cdf(self, X, **kwargs):
    """Not supported."""
    raise NotImplementedError

def fit(self, X, **kwargs):
    if hasattr(self, "nll_"):
        return super(Join, self).fit(X, **kwargs)
    else:
        raise NotImplementedError

def nll(self, X, **kwargs):
    out = np.zeros(len(X))
    start = 0
    for i, component in enumerate(self.components):
        out += component.nll(X[:, start:start+component.ndim], **kwargs)
        start += component.ndim
    return out

def pdf(self, X, **kwargs):
    out = np.ones(len(X))
    start = 0
    for i, component in enumerate(self.components):
        out *= component.pdf(X[:, start:start+component.ndim], **kwargs)
        start += component.ndim
    return out

def ppf(self, X, **kwargs):
    """Not supported."""
    raise NotImplementedError

def rvs(self, n_samples, random_state=None, **kwargs):
    rng = check_random_state(random_state)
    out = np.zeros((n_samples, self.ndim))
    start = 0
    for i, component in enumerate(self.components):
        out[:, start:start+component.ndim] = component.rvs(
            n_samples, random_state=rng, **kwargs)
        start += component.ndim
    return out

def score(self, X, **kwargs):
    """Evaluate the log-likelihood `-self.nll(X).sum()` of `X`.
    The higher, the better.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `score` [float]:
        The log-likelihood of `X`.
    """
    return -self.nll(X, **kwargs).sum()

def set_params(self, **params):
    for name, value in params.items():
        var = getattr(self, name, None)
        if isinstance(var, SharedVariable):
            var.set_value(value)
        elif (isinstance(var, T.TensorVariable) or
              isinstance(var, T.TensorConstant)):
            raise ValueError("Only shared variables can be updated.")
        else:
            super(DistributionMixin, self).set_params(**{name: value})

def __init__(self, bandwidth=None):
    """Constructor.
    Parameters
    ----------
    * `bandwidth` [string or float, optional]:
        The method used to calculate the estimator bandwidth.
    """
    self.bandwidth = bandwidth

def cdf(self, X, **kwargs):
    """Not supported."""
    raise NotImplementedError

def fit(self, X, **kwargs):
    """Fit the KDE estimator to data.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `self` [object]:
        `self`.
    """
    X = check_array(X).T
    self.kde_ = gaussian_kde(X, bw_method=self.bandwidth)
    return self

def nll(self, X, **kwargs):
    X = check_array(X)
    return -self.kde_.logpdf(X.T)

def pdf(self, X, **kwargs):
    X = check_array(X)
    return self.kde_.pdf(X.T)

def ppf(self, X, **kwargs):
    """Not supported."""
    raise NotImplementedError

def rvs(self, n_samples, random_state=None, **kwargs):
    # XXX gaussian_kde uses Numpy global random state...
    return self.kde_.resample(n_samples).T

def score(self, X, **kwargs):
    """Evaluate the goodness of fit of the distribution w.r.t. `X`.
    The higher, the better.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `score` [float]:
        The goodness of fit.
    """
    return NotImplementedError

def __init__(self, p, A):
    """Constructor.
    Parameters
    ----------
    * `p` [`DistributionMixin`]:
        The base distribution.
    * `A` [array, shape=(p.ndim, p.ndim)]:
        The linear operator.
    """
    super(LinearTransform, self).__init__()
    self.p = p
    self.A = A
    self.inv_A = np.linalg.inv(A)
    if isinstance(p, TheanoDistribution):
        for p_i in p.parameters_:
            self.parameters_.add(p_i)
        for c_i in p.constants_:
            self.constants_.add(c_i)
        for o_i in p.observeds_:
            self.observeds_.add(o_i)

def cdf(self, X, **kwargs):
    """Not supported."""
    raise NotImplementedError

def fit(self, X, bounds=None, constraints=None, use_gradient=True,
        optimizer=None, **kwargs):
    """Fit the distribution parameters to data by minimizing the negative
    log-likelihood of the data.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    * `bounds` [list of (parameter, (low, high))]:
        The parameter bounds.
    * `constraints`:
        The constraints on the parameters.
    * `use_gradient` [boolean, default=True]:
        Whether to use exact gradients (if `True`) or numerical gradients
        (if `False`).
    * `optimizer` [string]:
        The optimization method.
    Returns
    -------
    * `self` [object]:
        `self`.
    """
    # Map parameters to placeholders
    param_to_placeholder = []
    param_to_index = {}
    for i, v in enumerate(self.parameters_):
        w = T.TensorVariable(v.type)
        param_to_placeholder.append((v, w))
        param_to_index[v] = i
    # Build bounds
    mapped_bounds = None
    if bounds is not None:
        mapped_bounds = [(None, None) for v in param_to_placeholder]
        for b in bounds:
            mapped_bounds[param_to_index[b["param"]]] = b["bounds"]
    # Build constraints
    mapped_constraints = None
    if constraints is not None:
        mapped_constraints = []
        for c in constraints:
            args = c["param"]
            if isinstance(args, SharedVariable):
                args = (args, )
            m_c = {
                "type": c["type"],
                "fun": lambda x: c["fun"](*[x[param_to_index[a]]
                                            for a in args])
            }
            if "jac" in c:
                m_c["jac"] = lambda x: c["jac"](*[x[param_to_index[a]]
                                                  for a in args])
            mapped_constraints.append(m_c)
    # Derive objective and gradient
    objective_ = theano.function(
        [self.X] + [w for _, w in param_to_placeholder] +
        [theano.In(v, name=v.name) for v in self.observeds_],
        T.sum(self.nll_),
        givens=param_to_placeholder,
        allow_input_downcast=True)
    def objective(x):
        return objective_(X, *x, **kwargs) / len(X)
    if use_gradient:
        gradient_ = theano.function(
            [self.X] + [w for _, w in param_to_placeholder] +
            [theano.In(v, name=v.name) for v in self.observeds_],
            theano.grad(T.sum(self.nll_),
                        [v for v, _ in param_to_placeholder]),
            givens=param_to_placeholder,
            allow_input_downcast=True)
        def gradient(x):
            return np.array(gradient_(X, *x, **kwargs)) / len(X)
    # Solve!
    x0 = np.array([v.get_value() for v, _ in param_to_placeholder])
    r = minimize(objective,
                 jac=gradient if use_gradient else None,
                 x0=x0,
                 method=optimizer,
                 bounds=mapped_bounds,
                 constraints=mapped_constraints)
    if r.success:
        # Assign the solution
        for i, value in enumerate(r.x):
            param_to_placeholder[i][0].set_value(value)
    else:
        print("Parameter fitting failed!")
        print(r)
    return self

def nll(self, X, **kwargs):
    return self.p.nll(np.dot(self.inv_A, X.T).T, **kwargs)

def pdf(self, X, **kwargs):
    return self.p.pdf(np.dot(self.inv_A, X.T).T, **kwargs)

def ppf(self, X, **kwargs):
    """Not supported."""
    raise NotImplementedError

def rvs(self, n_samples, random_state=None, **kwargs):
    rng = check_random_state(random_state)
    out = self.p.rvs(n_samples, random_state=rng, **kwargs)
    return np.dot(self.A, out.T).T

def score(self, X, **kwargs):
    """Evaluate the log-likelihood `-self.nll(X).sum()` of `X`.
    The higher, the better.
    Parameters
    ----------
    * `X` [array-like, shape=(n_samples, n_features)]:
        The samples.
    Returns
    -------
    * `score` [float]:
        The log-likelihood of `X`.
    """
    return -self.nll(X, **kwargs).sum()

def set_params(self, **params):
    for name, value in params.items():
        var = getattr(self, name, None)
        if isinstance(var, SharedVariable):
            var.set_value(value)
        elif (isinstance(var, T.TensorVariable) or
              isinstance(var, T.TensorConstant)):
            raise ValueError("Only shared variables can be updated.")
        else:
            super(DistributionMixin, self).set_params(**{name: value})