Python numpy.random.gamma() Examples

def corrupt(Y, p=0.05):
    """Corrupt a count tensor with anomalies.

    The corruption noise model is:

        corrupt(y) = y * g, where g ~ Gamma(10, 2)

    PARAMS:
    p -- (float) proportion of tensor entries to corrupt

    RETURNS:
    out -- (np.ndarray) corrupted count tensor
    mask -- (np.ndarray) boolean array, same shape as count tensor
                         True means that entry was corrupted.
    """
    out = Y.copy()
    mask = (rn.random(size=out.shape) < p).astype(bool)
    out[mask] = rn.poisson(out[mask] * rn.gamma(10., 2., size=out[mask].shape))
    return out, mask 

def generate(shp=(30, 30, 20, 10), K=5, alpha=0.1, beta=0.1):
    """Generate a count tensor from the BPTF model.

    PARAMS:
    shp -- (tuple) shape of the generated count tensor
    K -- (int) number of latent components
    alpha -- (float) shape parameter of gamma prior over factors
    beta -- (float) rate parameter of gamma prior over factors

    RETURNS:
    Mu -- (np.ndarray) true Poisson rates
    Y -- (np.ndarray) generated count tensor
    """
    Theta_DK_M = [rn.gamma(alpha, 1./beta, size=(D, K)) for D in shp]
    Mu = parafac(Theta_DK_M)
    assert Mu.shape == shp
    Y = rn.poisson(Mu)
    return Mu, Y 

def _init_component(self, m, dim):
        assert self.mode_dims[m] == dim
        K = self.n_components
        s = self.smoothness
        if not self.debug:
            gamma_DK = s * rn.gamma(s, 1. / s, size=(dim, K))
            delta_DK = s * rn.gamma(s, 1. / s, size=(dim, K))
        else:
            gamma_DK = s * np.ones((dim, K))
            delta_DK = s * np.ones((dim, K))
        self.gamma_DK_M[m] = gamma_DK
        self.delta_DK_M[m] = delta_DK
        self.E_DK_M[m] = gamma_DK / delta_DK
        self.sumE_MK[m, :] = self.E_DK_M[m].sum(axis=0)
        self.G_DK_M[m] = np.exp(sp.psi(gamma_DK) - np.log(delta_DK))
        if m == 0 or not self.debug:
            self.beta_M[m] = 1. / self.E_DK_M[m].mean() 

def rand_mask(im_size):
    mask = np.zeros(shape=(im_size, im_size), dtype=np.uint8);
    cx = (im_size - 1) / 2.0;
    cy = (im_size - 1) / 2.0;
    a = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
    b = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
    ratio = float(max(a, b)) / min(a, b);
    mask = add_ellipse(mask, cx, cy, uniform() * math.pi, a, b);
    for i in range(np.random.randint(2, 5)):
        x = cx;
        y = cy;
        while ((x - cx)**2 + (y - cy)**2)**0.5 < im_size * 0.3:
            x = np.random.randint(0, im_size);
            y = np.random.randint(0, im_size);

        a = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
        b = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
        mask = add_ellipse(mask, x, y, uniform() * math.pi, a, b);

    return (mask, ratio); 

def rand_mask(im_size):
    mask = np.zeros(shape=(im_size, im_size), dtype=np.uint8);
    cx = (im_size - 1) / 2.0;
    cy = (im_size - 1) / 2.0;
    a = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
    b = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
    ratio = float(max(a, b)) / min(a, b);
    mask = add_ellipse(mask, cx, cy, uniform() * math.pi, a, b);
    for i in range(np.random.randint(0, 4)):
        x = cx;
        y = cy;
        while ((x - cx)**2 + (y - cy)**2)**0.5 < im_size * 0.3:
            x = np.random.randint(0, im_size);
            y = np.random.randint(0, im_size);

        a = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
        b = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
        mask = add_ellipse(mask, x, y, uniform() * math.pi, a, b);

    return (mask, ratio); 

def rand_mask(im_size):
    mask = np.zeros(shape=(im_size, im_size), dtype=np.uint8);
    cx = (im_size - 1) / 2.0;
    cy = (im_size - 1) / 2.0;
    a = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
    b = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
    ratio = float(max(a, b)) / min(a, b);
    mask = add_ellipse(mask, cx, cy, uniform() * math.pi, a, b);
    for i in range(np.random.randint(2, 5)):
        x = cx;
        y = cy;
        while ((x - cx)**2 + (y - cy)**2)**0.5 < im_size * 0.3:
            x = np.random.randint(0, im_size);
            y = np.random.randint(0, im_size);

        a = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
        b = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
        mask = add_ellipse(mask, x, y, uniform() * math.pi, a, b);

    return (mask, ratio); 

def rand_mask(im_size):
    mask = np.zeros(shape=(im_size, im_size), dtype=np.uint8);
    cx = (im_size - 1) / 2.0;
    cy = (im_size - 1) / 2.0;
    a = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
    b = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
    ratio = float(max(a, b)) / min(a, b);
    mask = add_ellipse(mask, cx, cy, uniform() * math.pi, a, b);
    for i in range(np.random.randint(2, 5)):
        x = cx;
        y = cy;
        while ((x - cx)**2 + (y - cy)**2)**0.5 < im_size * 0.3:
            x = np.random.randint(0, im_size);
            y = np.random.randint(0, im_size);

        a = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
        b = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
        mask = add_ellipse(mask, x, y, uniform() * math.pi, a, b);

    return (mask, ratio); 

def rand_mask(im_size):
    mask = np.zeros(shape=(im_size, im_size), dtype=np.uint8);
    cx = (im_size - 1) / 2.0;
    cy = (im_size - 1) / 2.0;
    a = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
    b = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
    ratio = float(max(a, b)) / min(a, b);
    mask = add_ellipse(mask, cx, cy, uniform() * math.pi, a, b);
    for i in range(np.random.randint(2, 5)):
        x = cx;
        y = cy;
        while ((x - cx)**2 + (y - cy)**2)**0.5 < im_size * 0.3:
            x = np.random.randint(0, im_size);
            y = np.random.randint(0, im_size);

        a = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
        b = im_size * (gamma(2.2, 0.6) / 10.0 + 0.04);
        mask = add_ellipse(mask, x, y, uniform() * math.pi, a, b);

    return (mask, ratio); 

def detect(Y, K=5, alpha=0.1, thresh=1e-5):
    """Detect anomalies using BPTF.

    This method fits BPTF to Y and obtains Mu, which is the model's
    reconstruction of Y (computed from the inferred latent factors).
    Anomalies are then all entries of Y whose probability given Mu
    is less than a given threshold.

        If P(y | mu) < thresh ==> y is  anomaly!

        Here P(y | mu) = Pois(y; mu), the PMF of the Poisson distribution.

    PARAMS:
    Y -- (np.ndarray) data count tensor
    K -- (int) number of latent components
    alpha -- (float) shape parameter of gamma prior over factors
    thresh -- (float) anomaly threshold (between 0 and 1).
    """
    bptf = BPTF(n_modes=Y.ndim,
                n_components=K,
                max_iter=100,
                tol=1e-4,
                smoothness=100,
                verbose=False,
                alpha=alpha,
                debug=False)
    bptf.fit(Y)
    Mu = bptf.reconstruct()

    return st.poisson.pmf(Y, Mu) < thresh 

def _create_prices(t):
    last_average = 100 if t==0 else source.data['average'][-1]
    returns = asarray(lognormal(mean.value, stddev.value, 1))
    average =  last_average * cumprod(returns)
    high = average * exp(abs(gamma(1, 0.03, size=1)))
    low = average / exp(abs(gamma(1, 0.03, size=1)))
    delta = high - low
    open = low + delta * uniform(0.05, 0.95, size=1)
    close = low + delta * uniform(0.05, 0.95, size=1)
    return open[0], high[0], low[0], close[0], average[0] 

def ppf(self, u):
        return stats.gamma.ppf(u, self.a, scale=self.scale) 

def logpdf(self, x):
        return stats.gamma.logpdf(x, self.a, scale=self.scale) 

def rvs(self, size=None):
        return random.gamma(self.a, scale=self.scale, size=size) 

def __init__(self, *args, **kwargs):
        """
        :name: Name of the parameter
        :shape: The shape of the gamma distribution.
        :scale: The scale of the gamme distribution
        :step:     (optional) number for step size required for some algorithms,
                eg. mcmc need a parameter of the variance for the next step
                default is median of rndfunc(*rndargs, size=1000)
        :optguess: (optional) number for start point of parameter
                default is quantile(0.5) - quantile(0.4) of
                rndfunc(*rndargs, size=1000)
        """

        super(Gamma, self).__init__(rnd.gamma, 'Gamma', *args, **kwargs) 

def _create_prices(t):
    last_average = 100 if t==0 else source.data['average'][-1]
    returns = asarray(lognormal(mean.value, stddev.value, 1))
    average =  last_average * cumprod(returns)
    high = average * exp(abs(gamma(1, 0.03, size=1)))
    low = average / exp(abs(gamma(1, 0.03, size=1)))
    delta = high - low
    open = low + delta * uniform(0.05, 0.95, size=1)
    close = low + delta * uniform(0.05, 0.95, size=1)
    return open[0], high[0], low[0], close[0], average[0] 

def _create_prices(t):
    last_average = 100 if t==0 else source.data['average'][-1]
    returns = asarray(lognormal(mean.value, stddev.value, 1))
    average =  last_average * cumprod(returns)
    high = average * exp(abs(gamma(1, 0.03, size=1)))
    low = average / exp(abs(gamma(1, 0.03, size=1)))
    delta = high - low
    open = low + delta * uniform(0.05, 0.95, size=1)
    close = low + delta * uniform(0.05, 0.95, size=1)
    return open[0], high[0], low[0], close[0], average[0] 

def _create_prices(t):
    global last_average
    returns = asarray(lognormal(mean, stddev, 1))
    average =  last_average * cumprod(returns)
    last_average = average

    high = average * exp(abs(gamma(1, 0.03, size=1)))
    low = average / exp(abs(gamma(1, 0.03, size=1)))
    delta = high - low
    open = low + delta * uniform(0.05, 0.95, size=1)
    close = low + delta * uniform(0.05, 0.95, size=1)
    return open[0], high[0], low[0], close[0], average[0] 

def beta_geometric_beta_binom_model(N, alpha, beta, gamma, delta, size=1):
    """
    Generate artificial data according to the Beta-Geometric/Beta-Binomial
    Model.

    You may wonder why we can have frequency = n_periods, when frequency excludes their
    first order. When a customer purchases something, they are born, _and in the next
    period_ we start asking questions about their alive-ness. So really they customer has
    bought frequency + 1, and been observed for n_periods + 1

    Parameters
    ----------
    N: array_like
        Number of transaction opportunities for new customers.
    alpha, beta, gamma, delta: float
        Parameters in the model. See [1]_
    size: int, optional
        The number of customers to generate

    Returns
    -------
    DataFrame
        with index as customer_ids and the following columns:
        'frequency', 'recency', 'n_periods', 'lambda', 'p', 'alive', 'customer_id'

    References
    ----------
    .. [1] Fader, Peter S., Bruce G.S. Hardie, and Jen Shang (2010),
       "Customer-Base Analysis in a Discrete-Time Noncontractual Setting,"
       Marketing Science, 29 (6), 1086-1108.

    """

    if type(N) in [float, int, np.int64]:
        N = N * np.ones(size)
    else:
        N = np.asarray(N)

    probability_of_post_purchase_death = random.beta(a=alpha, b=beta, size=size)
    thetas = random.beta(a=gamma, b=delta, size=size)

    columns = ["frequency", "recency", "n_periods", "p", "theta", "alive", "customer_id"]
    df = pd.DataFrame(np.zeros((size, len(columns))), columns=columns)
    for i in range(size):
        p = probability_of_post_purchase_death[i]
        theta = thetas[i]

        # hacky until I can find something better
        current_t = 0
        alive = True
        times = []
        while current_t < N[i] and alive:
            alive = random.binomial(1, theta) == 0
            if alive and random.binomial(1, p) == 1:
                times.append(current_t)
            current_t += 1
        # adding in final death opportunity to agree with [1]
        if alive:
            alive = random.binomial(1, theta) == 0
        df.iloc[i] = len(times), times[-1] + 1 if len(times) != 0 else 0, N[i], p, theta, alive, i
    return df