<jupyter_start><jupyter_text>`mixture.py` This notebook tests the `mixture` module.This module contains the classes to reconstruct a probability density and some utility functions.  Utilities  `numba_gammaln` This function is a Numba wrapper for the Cython implementation of `scipy.special.gammaln`, takes a double as input.<jupyter_code>import numpy as np
from figaro.mixture import _numba_gammaln
from scipy.special import gammaln

x = 0.5

# Numba first run is significantly slower due to compilation
_ = _numba_gammaln(x)

print("Scipy:")
%timeit gammaln(x)
print("Numba:")
%timeit _numba_gammaln(x)

y = np.arange(0.1,5, 1000)

numba_val = np.array([_numba_gammaln(yi) for yi in y])
scipy_val = gammaln(y)

np.alltrue(numba_val == scipy_val)<jupyter_output><empty_output><jupyter_text>`student_t` [<NAME>'s implementation](http://gregorygundersen.com/blog/2020/01/20/multivariate-t/) of multivariate Student-t distribution (pulled to Scipy). This function takes 2d arrays as inputs (as FIGARO).<jupyter_code>from figaro.mixture import _student_t
from scipy.stats import multivariate_t

# Parameters
df    = 5
dim   = 3
t     = np.ones(dim)*0.5
mu    = np.ones(dim)
sigma = np.identity(dim)*0.7**2

# Instance multivariate_t
scipy_t = multivariate_t(loc = mu, shape = sigma, df = df)

# Numba first run is significantly slower due to compilation
_ = _student_t(df = df, t = np.atleast_2d(t), mu = mu, sigma = sigma, dim = dim)

print("Scipy:")
%timeit scipy_t.logpdf(t)
# FIGARO handles 2d arrays
t  = np.atleast_2d(t)
mu = np.atleast_2d(mu)
print("Numba:")
%timeit _student_t(df = df, t = t, mu = mu, sigma = sigma, dim = dim)

T = scipy_t.rvs(1000)
scipy_val = np.array([scipy_t.logpdf(ti) for ti in T])
T = np.array([np.atleast_2d(ti) for ti in T])
numba_val = np.array([_student_t(df = df, t = ti, mu = mu, sigma = sigma, dim = dim) for ti in T])

np.allclose(numba_val, scipy_val, atol = 1e-15)<jupyter_output><empty_output><jupyter_text>`update_alpha` MH sampler for concentration parameter. Implements likelihood in Eq. (3.6) of [Rinaldi & <NAME> (2022a)](https://arxiv.org/pdf/2109.05960.pdf):\$p(\alpha|K,n) \propto \frac{\Gamma(\alpha)}{\Gamma(n+\alpha)} \alpha^K e^{-1/\alpha}\,.$\`alpha` is the initial value for concentration parameter, `n` the number of samples inclued in the mixture and `K` is the number of active components. Since for the inference algorithm we need only one `alpha` sample per iteration, this method returns only one draw from the distribution.<jupyter_code>from figaro.mixture import update_alpha
from scipy.special import gammaln, logsumexp

def prob_alpha(a, k, n):
    return gammaln(a) - gammaln(a+n) + k*np.log(a) - 1/a

alpha0 = 1
n = 300
K = 20
n_draws = 10000

A = [alpha0]
for i in range(n_draws):
    A.append(update_alpha(A[-1], n, K))
A = np.array(A)

import matplotlib.pyplot as plt
alphas = np.linspace(0,np.max(A),1001)[1:]
plt.hist(A, bins = int(np.sqrt(n_draws+1)), histtype = 'step', density = True, label = '$\\mathrm{samples}$')

p = prob_alpha(alphas, K, n)

plt.plot(alphas, np.exp(p - logsumexp(p) - np.log(alphas[1]-alphas[0])), lw = 0.8, c = 'r', label = '$K = {0},\ n={1}$'.format(K,n))
plt.xlabel('$\\alpha$')
plt.ylabel('$p(\\alpha|K,n)$')
plt.legend(loc = 0, frameon = False)
plt.grid(alpha = 0.5)<jupyter_output><empty_output><jupyter_text>The samples are uncorrelated:<jupyter_code>def autocorrelation(v, tau, l, m, s):
    return np.mean([(v[i]-m)*(v[(i+tau)%l]-m) for i in range(l)])/s

m = np.mean(A)
s = np.var(A-m)
C = []

max_tau = 100

taus = np.arange(0, np.min([max_tau, n_draws//2]))
for t in taus:
    C.append(autocorrelation(A, t, n_draws, m, s))
C = np.array(C)
    
plt.axhline(0, lw = 0.5, ls = '--')
plt.plot(taus, C, lw = 0.7)
plt.xlabel('$\\tau$')
plt.ylabel('$C(\\tau)$')
plt.grid(alpha = 0.5)<jupyter_output><empty_output><jupyter_text>`compute_component_suffstats` This method iteratively updates mean and scatter matrix of the samples assigned to a cluster and returns them along with the expected values for mean and covariance matrix of the cluster given a NIW prior. Here we will check that the updated mean and covariance are the same as they were computed directly with all the samples.<jupyter_code>from figaro.mixture import compute_component_suffstats

dim = 3
n_draws = 100

samples = np.random.uniform(size = (n_draws, dim))
p_mu = np.zeros(dim)
p_L  = np.identity(dim)
p_k  = 1e-2
p_nu = 5

means      = [samples[0]]
covs       = [np.zeros(shape = (dim, dim))]
iter_means = [samples[0]]
iter_covs  = [np.zeros(shape = (dim, dim))]

for i in range(1,n_draws):
    # Mean
    means.append(np.mean(samples[:i+1], axis = 0))
    # Scatter matrix
    c = np.zeros(shape = (dim,dim))
    for vi in samples[:i+1]:
        c += (vi - means[-1]).T@(vi - means[-1])
    covs.append(c)
    # FIGARO
    n_mean, n_cov, N, n_mu, n_sigma = compute_component_suffstats(samples[i], iter_means[-1], iter_covs[-1], i, p_mu, p_k, p_nu, p_L)
    iter_means.append(n_mean)
    iter_covs.append(n_cov)


print('Mean: {0}'.format(np.allclose(means, iter_means, atol = 1e-15)))
print('Scatter matrix: {0}'.format(np.allclose(covs, iter_covs, atol = 1e-15)))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15167343780228285, 'fraction_numerical': 0.023408783130199264, 'mean_word_length': 4.619565217391305, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 23, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '640681', 'license': 'MIT', 'n_tokens_mistral': 2182, 'n_tokens_neox': 2046, 'n_words': 525, 'path': 'notebooks/tests/mixture.ipynb', 'repo_name': 'sterinaldi/figaro', 'size': '5178'}