# Path: bmlingam/prob.py
# def ll_laplace(e):
#     """Return of log likelihood of the standard Laplace distribution. 
# 
#     :param e: Sample values.
#     :type e: ndarray, shape=(n_theta_sampling, n_samples)
#     :return: pdf values at given samples. 
#     :rtype: ndarray, shape=(n_theta_sampling, n_samples)
#     """
#     b = 1 / sqrt(2)
#     return -log(2 * b) - abs(e) / b
# 
# def ll_gg(e, beta):
#     """Return of log likelihood of generalized Gaussian distributions. 
#     """
#     beta = float(beta)
#     m = gamma(0.5 / beta) / ((2**(1 / beta)) * gamma(3 / (2 * beta)))
#     
#     return - 0.5 * power((e**2) / m, beta) + log(beta) \
#            - gammaln(0.5 / beta) - (0.5 / beta) * log(2) - 0.5 * log(m)
# 
# def sample_gg(scov, beta, n_samples, rng, dim, normalize):
#     """Draw samples from GG distribution. 
# 
#     See: Gómez, E., Gomez-Viilegas, M. A., & Marin, J. M. (1998). 
#     A multivariate generalization of the power exponential family of distributions. 
#     Communications in Statistics-Theory and Methods, 27(3), 589-600.
# 
#     Parameters
#     ----------
#     scov: numpy.ndarray, shape=(dim, dim)
#         A scaled covariance matrix. For beta=1, it becomes the covariance matrix. 
# 
#     beta: float
#         Shape parameter, a positive value. 
# 
#     n_samples: int
#         Number of samples. 
# 
#     rng: numpy.random.RandomState
#         Random number generator. 
# 
#     dim: int 
#         The dimension of the distribution. 
# 
#     normalize: bool
#         If true, each dimension is normalized such that its variance is 1. 
#     """
#     # Draw samples from unit sphere
#     gs = rng.normal(0, 1, size=(n_samples, dim))
#     ss = gs + sign(sign(gs) + 1e-10) * 1e-10
#     ns = np.linalg.norm(ss, ord=2, axis=1)[:, np.newaxis]
#     us = ss / ns
# 
#     # Draw samples from GG
#     S = np.linalg.cholesky(scov).T
#     ts = rng.gamma(shape=dim / (2.0 * beta), scale=2.0, 
#                    size=n_samples)[:, np.newaxis]
#     xs = (ts**(1.0 / (2.0 * beta))) * us.dot(S)
# 
#     if normalize:
#         v = (2**(1.0 / beta) * gamma((dim + 2.0) / (2.0 * beta))) / \
#             (dim * gamma(dim / (2.0 * beta)))
#         xs = xs / np.sqrt(v * np.diag(scov))
# 
#     return xs


import matplotlib.pyplot as plt
import numpy as np
from numpy.testing import assert_allclose
from scipy import stats
from scipy.special import gamma
from bmlingam.prob import ll_laplace, ll_gg, sample_gg

    return k

def test_sample_gg(n_samples=1000000, plot=False):
    """Tests for generalized Gaussian. 
    """
    rng = np.random.RandomState(0)

    # Test 1
    print('Test1')
    dim = 2
    scov = np.eye(dim)
    beta = 1.0
    xs = sample_gg(scov, beta, n_samples, rng, dim, normalize=True)
    _describe_and_check('xs[:, 0]', xs[:, 0], [0, 1, 0, 0])
    _describe_and_check('xs[:, 1]', xs[:, 1], [0, 1, 0, 0])

    # Test 2
    print('\nTest2')
    dim = 2
    scov = np.array([[1.0, 0.5], [0.5, 1.0]])
    beta = 1.0
    xs = sample_gg(scov, beta, n_samples, rng, dim, normalize=True)
    _describe_and_check('xs[:, 0]', xs[:, 0], [0, 1, 0, 0])
    _describe_and_check('xs[:, 1]', xs[:, 1], [0, 1, 0, 0])

    # Test 3
    print('\nTest3')
    dim = 1
    scov = np.eye(dim)
    beta = 0.5

# Path: bmlingam/prob.py
# def ll_laplace(e):
#     """Return of log likelihood of the standard Laplace distribution. 
# 
#     :param e: Sample values.
#     :type e: ndarray, shape=(n_theta_sampling, n_samples)
#     :return: pdf values at given samples. 
#     :rtype: ndarray, shape=(n_theta_sampling, n_samples)
#     """
#     b = 1 / sqrt(2)
#     return -log(2 * b) - abs(e) / b
# 
# def ll_gg(e, beta):
#     """Return of log likelihood of generalized Gaussian distributions. 
#     """
#     beta = float(beta)
#     m = gamma(0.5 / beta) / ((2**(1 / beta)) * gamma(3 / (2 * beta)))
#     
#     return - 0.5 * power((e**2) / m, beta) + log(beta) \
#            - gammaln(0.5 / beta) - (0.5 / beta) * log(2) - 0.5 * log(m)
# 
# def sample_gg(scov, beta, n_samples, rng, dim, normalize):
#     """Draw samples from GG distribution. 
# 
#     See: Gómez, E., Gomez-Viilegas, M. A., & Marin, J. M. (1998). 
#     A multivariate generalization of the power exponential family of distributions. 
#     Communications in Statistics-Theory and Methods, 27(3), 589-600.
# 
#     Parameters
#     ----------
#     scov: numpy.ndarray, shape=(dim, dim)
#         A scaled covariance matrix. For beta=1, it becomes the covariance matrix. 
# 
#     beta: float
#         Shape parameter, a positive value. 
# 
#     n_samples: int
#         Number of samples. 
# 
#     rng: numpy.random.RandomState
#         Random number generator. 
# 
#     dim: int 
#         The dimension of the distribution. 
# 
#     normalize: bool
#         If true, each dimension is normalized such that its variance is 1. 
#     """
#     # Draw samples from unit sphere
#     gs = rng.normal(0, 1, size=(n_samples, dim))
#     ss = gs + sign(sign(gs) + 1e-10) * 1e-10
#     ns = np.linalg.norm(ss, ord=2, axis=1)[:, np.newaxis]
#     us = ss / ns
# 
#     # Draw samples from GG
#     S = np.linalg.cholesky(scov).T
#     ts = rng.gamma(shape=dim / (2.0 * beta), scale=2.0, 
#                    size=n_samples)[:, np.newaxis]
#     xs = (ts**(1.0 / (2.0 * beta))) * us.dot(S)
# 
#     if normalize:
#         v = (2**(1.0 / beta) * gamma((dim + 2.0) / (2.0 * beta))) / \
#             (dim * gamma(dim / (2.0 * beta)))
#         xs = xs / np.sqrt(v * np.diag(scov))
# 
#     return xs


# Path: bmlingam/tests/test_prob.py
import matplotlib.pyplot as plt
import numpy as np
from numpy.testing import assert_allclose
from scipy import stats
from scipy.special import gamma
from bmlingam.prob import ll_laplace, ll_gg, sample_gg

    return k

def test_sample_gg(n_samples=1000000, plot=False):
    """Tests for generalized Gaussian. 
    """
    rng = np.random.RandomState(0)

    # Test 1
    print('Test1')
    dim = 2
    scov = np.eye(dim)
    beta = 1.0
    xs = sample_gg(scov, beta, n_samples, rng, dim, normalize=True)
    _describe_and_check('xs[:, 0]', xs[:, 0], [0, 1, 0, 0])
    _describe_and_check('xs[:, 1]', xs[:, 1], [0, 1, 0, 0])

    # Test 2
    print('\nTest2')
    dim = 2
    scov = np.array([[1.0, 0.5], [0.5, 1.0]])
    beta = 1.0
    xs = sample_gg(scov, beta, n_samples, rng, dim, normalize=True)
    _describe_and_check('xs[:, 0]', xs[:, 0], [0, 1, 0, 0])
    _describe_and_check('xs[:, 1]', xs[:, 1], [0, 1, 0, 0])

    # Test 3
    print('\nTest3')
    dim = 1
    scov = np.eye(dim)
    beta = 0.5

    xs = sample_gg(scov, beta, n_samples, rng, dim, normalize=True)
