Source code for astroML.stats.random

"""
Statistics for astronomy
"""
import numpy as np
from scipy.stats.distributions import rv_continuous


[docs]def bivariate_normal(mu=[0, 0], sigma_1=1, sigma_2=1, alpha=0, size=None, return_cov=False): """Sample points from a 2D normal distribution Parameters ---------- mu : array-like (length 2) The mean of the distribution sigma_1 : float The unrotated x-axis width sigma_2 : float The unrotated y-axis width alpha : float The rotation counter-clockwise about the origin size : tuple of ints, optional Given a shape of, for example, ``(m,n,k)``, ``m*n*k`` samples are generated, and packed in an `m`-by-`n`-by-`k` arrangement. Because each sample is `N`-dimensional, the output shape is ``(m,n,k,N)``. If no shape is specified, a single (`N`-D) sample is returned. return_cov : boolean, optional If True, return the computed covariance matrix. Returns ------- out : ndarray The drawn samples, of shape *size*, if that was provided. If not, the shape is ``(N,)``. In other words, each entry ``out[i,j,...,:]`` is an N-dimensional value drawn from the distribution. cov : ndarray The 2x2 covariance matrix. Returned only if return_cov == True. Notes ----- This function works by computing a covariance matrix from the inputs, and calling ``np.random.multivariate_normal()``. If the covariance matrix is available, this function can be called directly. """ # compute covariance matrix sigma_xx = ((sigma_1 * np.cos(alpha)) ** 2 + (sigma_2 * np.sin(alpha)) ** 2) sigma_yy = ((sigma_1 * np.sin(alpha)) ** 2 + (sigma_2 * np.cos(alpha)) ** 2) sigma_xy = (sigma_1 ** 2 - sigma_2 ** 2) * np.sin(alpha) * np.cos(alpha) cov = np.array([[sigma_xx, sigma_xy], [sigma_xy, sigma_yy]]) # draw points from the distribution x = np.random.multivariate_normal(mu, cov, size) if return_cov: return x, cov else: return x
#---------------------------------------------------------------------- # Define some new distributions based on rv_continuous class trunc_exp_gen(rv_continuous): """A truncated positive exponential continuous random variable. The probability distribution is:: p(x) ~ exp(k * x) between a and b = 0 otherwise The arguments are (a, b, k) %(before_notes)s %(example)s """ def _argcheck(self, a, b, k): self._const = k / (np.exp(k * b) - np.exp(k * a)) return (a != b) and not np.isinf(k) def _pdf(self, x, a, b, k): pdf = self._const * np.exp(k * x) pdf[(x < a) | (x > b)] = 0 return pdf def _rvs(self, a, b, k): y = np.random.random(self._size) return (1. / k) * np.log(1 + y * k / self._const) trunc_exp = trunc_exp_gen(name="trunc_exp", shapes='a, b, k') class linear_gen(rv_continuous): """A truncated positive exponential continuous random variable. The probability distribution is:: p(x) ~ c * x + d between a and b = 0 otherwise The arguments are (a, b, c). d is set by the normalization %(before_notes)s %(example)s """ def _argcheck(self, a, b, c): return (a != b) and not np.isinf(c) def _pdf(self, x, a, b, c): d = 1. / (b - a) - 0.5 * c * (b + a) pdf = c * x + d pdf[(x < a) | (x > b)] = 0 return pdf def _rvs(self, a, b, c): mu = 0.5 * (a + b) W = (b - a) x0 = 1. / c / W - mu r = np.random.random(self._size) return -x0 + np.sqrt(2. * r / c + a * a + 2. * a * x0 + x0 * x0) linear = linear_gen(name="linear", shapes='a, b, c')