import math
import numpy as np
import scipy
import scipy.special
import scipy.stats
rng = np.random.RandomState(42)
def gaussian_predictive_theoretical_prob_slow(x_1n, x_1n1, mu, cov, var):
n = len(x_1n)
cov_matrix_x1n = np.ones((n, n)) * cov
cov_matrix_x1n[np.diag_indices(n)] = var
pdf_x_1n = scipy.stats.multivariate_normal.pdf(x_1n, mean=np.ones(n) * mu, cov=cov_matrix_x1n)
cov_matrix_x1n1 = np.ones((n + 1, n + 1)) * cov
cov_matrix_x1n1[np.diag_indices(n + 1)] = var
pdf_x_1n1 = scipy.stats.multivariate_normal.pdf(x_1n1, mean=np.ones(n + 1) * mu, cov=cov_matrix_x1n1)
return pdf_x_1n1 / pdf_x_1n
def predictive_mdim_gaussian_dist(X, cov_u, cov_v, var_u, var_v, mu, recursive=False):
n = X.shape[0]
p = X.shape[1]
mu_a = np.ones((p,)) * mu
# U is an (n x n) convariance matrix (intersamples)
# V is a (p x p) convariance matrix (insample)
v = np.ones((p, p)) * cov_v
np.fill_diagonal(v, var_v)
if recursive:
print('using recursive implementation')
mu_a_given_b = mu_a.reshape((1, p))
cov_a_given_b = var_u * v
ident_m = np.eye(p)
for i in range(1, n + 1):
dd_i = cov_u / (var_u + cov_u * (i - 1)) * ident_m
mu_a_given_b = mu_a_given_b.dot(ident_m - dd_i) + X[i - 1, :].dot(dd_i)
cov_a_given_b = cov_a_given_b.dot(ident_m - dd_i) + (var_u - cov_u) * v.dot(dd_i)
mu_a_given_b = mu_a_given_b.flatten()
else:
dd = cov_u / (var_u + cov_u * (n - 1)) * np.eye(p)
mu_a_given_b = mu_a + dd.dot(np.sum(X - mu, axis=0))
cov_a_given_b = var_u * v - n * dd.dot(cov_u * v)
return mu_a_given_b, cov_a_given_b
def predictive_1d_student_dist(X, nu, cov, var, mu, recursive=False):
n = X.shape[0]
p = X.shape[1]
assert p == 1
mu_a = np.ones((p,)) * mu
X_zm = X - mu
if recursive:
# print 'recursive implementation'
if p > 1:
raise NotImplementedError('recursive implementation works only for 1D case')
nu_a_given_b = nu
mu_a_given_b = mu_a.reshape((1, p))
beta = 0
K_aa = var
x_cum_sum = 0
for i in range(1, n + 1):
nu_a_given_b += p
dd_i = cov / (var + cov * (i - 1))
mu_a_given_b = (1. - dd_i) * mu_a_given_b + X[i - 1] * dd_i
a_i = (cov * (i - 2) + var) / ((var - cov) * (cov * (i - 1) + var))
b_i = -1. * cov / ((var - cov) * (cov * (i - 1) + var))
b_i_prev = -1. * cov / ((var - cov) * (cov * (i - 2) + var))
beta += (a_i - b_i) * X_zm[i - 1] ** 2 \
+ b_i * (x_cum_sum + X_zm[i - 1]) ** 2 \
- b_i_prev * x_cum_sum ** 2
x_cum_sum += X_zm[i - 1]
K_aa = (1 - dd_i) * K_aa + (var - cov) * dd_i
cov_a_given_b = (nu + beta - 2.) / (nu_a_given_b - 2.) * K_aa
cov_a_given_b = np.array(cov_a_given_b).reshape((1, 1))
mu_a_given_b = mu_a_given_b.flatten()
else:
v = np.eye(p)
n_b = n * p
dd = cov / (var + cov * (n - 1)) * np.eye(p)
A = (cov * (n - 2) + var) / ((var - cov) * (cov * (n - 1) + var)) * np.linalg.inv(v)
B = -1 * cov / ((var - cov) * (cov * (n - 1) + var)) * np.linalg.inv(v)
nu_a_given_b = nu + n_b
mu_a_given_b = mu_a + dd.dot(np.sum(X - mu, axis=0))
beta = np.sum(np.diagonal(X_zm.dot(A).dot(X_zm.transpose()))) \
+ np.sum(X_zm.dot(B).dot(X_zm.transpose())) - \
np.sum(np.diagonal(X_zm.dot(B).dot(X_zm.transpose())))
cov_a_given_b = (nu + beta - 2.) / (nu + n_b - 2.) * (var * v - n * dd.dot(cov * v))
return nu_a_given_b, mu_a_given_b, cov_a_given_b
def predictive_mdim_student_prob(x_1n, x_next, cov_u, var_u, mu, nu, recursive=False, verbose=False):
p = x_1n.shape[-1]
if p == 1:
nu_a_given_b, mu_a_given_b, cov_a_given_b = predictive_1d_student_dist(x_1n,
nu=nu[0],
cov=cov_u[0],
var=var_u[0],
mu=mu[0], recursive=recursive)
if verbose:
print(nu_a_given_b)
print(mu_a_given_b)
print(cov_a_given_b)
prob = mvt_pdf(x_next, phi=mu_a_given_b, K=cov_a_given_b, nu=nu_a_given_b)
prob_1d = student_pdf_1d(x_next, phi=mu_a_given_b, K=cov_a_given_b, nu=nu_a_given_b)
if verbose:
print('****', prob, prob_1d)
else:
probs_i = []
for i in range(p):
nu_a_given_b, mu_a_given_b, cov_a_given_b = predictive_1d_student_dist(x_1n[:, i:i + 1],
nu=nu[i],
cov=cov_u[i],
var=var_u[i],
mu=mu[i],
recursive=recursive)
if verbose:
print(nu_a_given_b)
print(mu_a_given_b)
print(cov_a_given_b)
probs_i.append(
mvt_pdf(x_next[i:i + 1], phi=mu_a_given_b, K=cov_a_given_b, nu=nu_a_given_b))
prob = np.prod(probs_i)
return prob
def mvt_pdf(X, phi, K, nu):
"""
Multivariate student-t density:
output:
the density of the given element
input:
x = parameter (d dimensional numpy array or scalar)
mu = mean (d dimensional numpy array or scalar)
K = scale matrix (dxd numpy array)
nu = degrees of freedom
"""
d = X.shape[-1]
num = math.gamma((d + nu) / 2.) * pow(
1. + (1. / (nu - 2)) * ((X - phi).dot(np.linalg.inv(K)).dot(np.transpose(X - phi))), -(d + nu) / 2.)
denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi, d / 2.) * pow(np.linalg.det(K), 0.5)
return num / denom
def student_pdf_1d(X, phi, K, nu):
"""
Univariate student-t density:
output:
the density of the given element
input:
x = parameter scalar
mu = mean
var = scale matrix
nu = degrees of freedom
"""
num = math.gamma((1. + nu) / 2.) * pow(
1. + (1. / (nu - 2)) * (1. / K * (X - phi) ** 2), -(1. + nu) / 2.)
denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi * K, 0.5)
return num / denom
