Published: December 27, 2020
Author: Xinyu Chen [GitHub homepage]
Download: This Jupyter notebook is at our GitHub repository. If you want to evaluate the code, please download the notebook from the transdim repository.
This notebook shows how to implement the Bayesian Temporal Tensor Factorization (BTTF), a fully Bayesian matrix factorization model, on some real-world data sets. To overcome the missing data problem in multivariate time series, BTTF takes into account both low-rank matrix structure and time series autoregression. For an in-depth discussion of BTTF, please see [1].
import numpy as np
from numpy.linalg import inv as inv
from numpy.random import normal as normrnd
from numpy.random import multivariate_normal as mvnrnd
from scipy.linalg import khatri_rao as kr_prod
from scipy.stats import wishart
from scipy.stats import invwishart
from numpy.linalg import solve as solve
from numpy.linalg import cholesky as cholesky_lower
from scipy.linalg import cholesky as cholesky_upper
from scipy.linalg import solve_triangular as solve_ut
import matplotlib.pyplot as plt
%matplotlib inline
def mvnrnd_pre(mu, Lambda):
src = normrnd(size = (mu.shape[0],))
return solve_ut(cholesky_upper(Lambda, overwrite_a = True, check_finite = False),
src, lower = False, check_finite = False, overwrite_b = True) + mu
def cov_mat(mat, mat_bar):
mat = mat - mat_bar
return mat.T @ mat
def ten2mat(tensor, mode):
return np.reshape(np.moveaxis(tensor, mode, 0), (tensor.shape[mode], -1), order = 'F')
def sample_factor_u(tau_sparse_tensor, tau_ind, U, V, X, beta0 = 1):
"""Sampling M-by-R factor matrix U and its hyperparameters (mu_u, Lambda_u)."""
dim1, rank = U.shape
U_bar = np.mean(U, axis = 0)
temp = dim1 / (dim1 + beta0)
var_mu_hyper = temp * U_bar
var_U_hyper = inv(np.eye(rank) + cov_mat(U, U_bar) + temp * beta0 * np.outer(U_bar, U_bar))
var_Lambda_hyper = wishart.rvs(df = dim1 + rank, scale = var_U_hyper)
var_mu_hyper = mvnrnd_pre(var_mu_hyper, (dim1 + beta0) * var_Lambda_hyper)
var1 = kr_prod(X, V).T
var2 = kr_prod(var1, var1)
var3 = (var2 @ ten2mat(tau_ind, 0).T).reshape([rank, rank, dim1]) + var_Lambda_hyper[:, :, None]
var4 = var1 @ ten2mat(tau_sparse_tensor, 0).T + (var_Lambda_hyper @ var_mu_hyper)[:, None]
for i in range(dim1):
U[i, :] = mvnrnd_pre(solve(var3[:, :, i], var4[:, i]), var3[:, :, i])
return U
def sample_factor_v(tau_sparse_tensor, tau_ind, U, V, X, beta0 = 1):
"""Sampling N-by-R factor matrix V and its hyperparameters (mu_v, Lambda_v)."""
dim2, rank = V.shape
V_bar = np.mean(V, axis = 0)
temp = dim2 / (dim2 + beta0)
var_mu_hyper = temp * V_bar
var_V_hyper = inv(np.eye(rank) + cov_mat(V, V_bar) + temp * beta0 * np.outer(V_bar, V_bar))
var_Lambda_hyper = wishart.rvs(df = dim2 + rank, scale = var_V_hyper)
var_mu_hyper = mvnrnd_pre(var_mu_hyper, (dim2 + beta0) * var_Lambda_hyper)
var1 = kr_prod(X, U).T
var2 = kr_prod(var1, var1)
var3 = (var2 @ ten2mat(tau_ind, 1).T).reshape([rank, rank, dim2]) + var_Lambda_hyper[:, :, None]
var4 = var1 @ ten2mat(tau_sparse_tensor, 1).T + (var_Lambda_hyper @ var_mu_hyper)[:, None]
for j in range(dim2):
V[j, :] = mvnrnd_pre(solve(var3[:, :, j], var4[:, j]), var3[:, :, j])
return V
def mnrnd(M, U, V):
"""
Generate matrix normal distributed random matrix.
M is a m-by-n matrix, U is a m-by-m matrix, and V is a n-by-n matrix.
"""
dim1, dim2 = M.shape
X0 = np.random.randn(dim1, dim2)
P = cholesky_lower(U)
Q = cholesky_lower(V)
return M + P @ X0 @ Q.T
def sample_var_coefficient(X, time_lags):
dim, rank = X.shape
d = time_lags.shape[0]
tmax = np.max(time_lags)
Z_mat = X[tmax : dim, :]
Q_mat = np.zeros((dim - tmax, rank * d))
for k in range(d):
Q_mat[:, k * rank : (k + 1) * rank] = X[tmax - time_lags[k] : dim - time_lags[k], :]
var_Psi0 = np.eye(rank * d) + Q_mat.T @ Q_mat
var_Psi = inv(var_Psi0)
var_M = var_Psi @ Q_mat.T @ Z_mat
var_S = np.eye(rank) + Z_mat.T @ Z_mat - var_M.T @ var_Psi0 @ var_M
Sigma = invwishart.rvs(df = rank + dim - tmax, scale = var_S)
return mnrnd(var_M, var_Psi, Sigma), Sigma
def sample_factor_x(tau_sparse_tensor, tau_ind, time_lags, U, V, X, A, Lambda_x):
"""Sampling T-by-R factor matrix X."""
dim3, rank = X.shape
tmax = np.max(time_lags)
tmin = np.min(time_lags)
d = time_lags.shape[0]
A0 = np.dstack([A] * d)
for k in range(d):
A0[k * rank : (k + 1) * rank, :, k] = 0
mat0 = Lambda_x @ A.T
mat1 = np.einsum('kij, jt -> kit', A.reshape([d, rank, rank]), Lambda_x)
mat2 = np.einsum('kit, kjt -> ij', mat1, A.reshape([d, rank, rank]))
var1 = kr_prod(V, U).T
var2 = kr_prod(var1, var1)
var3 = (var2 @ ten2mat(tau_ind, 2).T).reshape([rank, rank, dim3]) + Lambda_x[:, :, None]
var4 = var1 @ ten2mat(tau_sparse_tensor, 2).T
for t in range(dim3):
Mt = np.zeros((rank, rank))
Nt = np.zeros(rank)
Qt = mat0 @ X[t - time_lags, :].reshape(rank * d)
index = list(range(0, d))
if t >= dim3 - tmax and t < dim3 - tmin:
index = list(np.where(t + time_lags < dim3))[0]
elif t < tmax:
Qt = np.zeros(rank)
index = list(np.where(t + time_lags >= tmax))[0]
if t < dim3 - tmin:
Mt = mat2.copy()
temp = np.zeros((rank * d, len(index)))
n = 0
for k in index:
temp[:, n] = X[t + time_lags[k] - time_lags, :].reshape(rank * d)
n += 1
temp0 = X[t + time_lags[index], :].T - np.einsum('ijk, ik -> jk', A0[:, :, index], temp)
Nt = np.einsum('kij, jk -> i', mat1[index, :, :], temp0)
var3[:, :, t] = var3[:, :, t] + Mt
if t < tmax:
var3[:, :, t] = var3[:, :, t] - Lambda_x + np.eye(rank)
X[t, :] = mvnrnd_pre(solve(var3[:, :, t], var4[:, t] + Nt + Qt), var3[:, :, t])
return X
def compute_mape(var, var_hat):
return np.sum(np.abs(var - var_hat) / var) / var.shape[0]
def compute_rmse(var, var_hat):
return np.sqrt(np.sum((var - var_hat) ** 2) / var.shape[0])
def ar4cast(A, X, Sigma, time_lags, multi_step):
dim, rank = X.shape
d = time_lags.shape[0]
X_new = np.append(X, np.zeros((multi_step, rank)), axis = 0)
for t in range(multi_step):
var = A.T @ X_new[dim + t - time_lags, :].reshape(rank * d)
X_new[dim + t, :] = mvnrnd(var, Sigma)
return X_new
def BTTF(dense_tensor, sparse_tensor, init, rank, time_lags, burn_iter, gibbs_iter, multi_step = 1):
"""Bayesian Temporal Tensor Factorization, BTTF."""
dim1, dim2, dim3 = sparse_tensor.shape
d = time_lags.shape[0]
U = init["U"]
V = init["V"]
X = init["X"]
if np.isnan(sparse_tensor).any() == False:
ind = sparse_tensor != 0
pos_obs = np.where(ind)
pos_test = np.where((dense_tensor != 0) & (sparse_tensor == 0))
elif np.isnan(sparse_tensor).any() == True:
pos_test = np.where((dense_tensor != 0) & (np.isnan(sparse_tensor)))
ind = ~np.isnan(sparse_tensor)
pos_obs = np.where(ind)
sparse_tensor[np.isnan(sparse_tensor)] = 0
dense_test = dense_tensor[pos_test]
del dense_tensor
U_plus = np.zeros((dim1, rank, gibbs_iter))
V_plus = np.zeros((dim2, rank, gibbs_iter))
X_plus = np.zeros((dim3 + multi_step, rank, gibbs_iter))
A_plus = np.zeros((rank * d, rank, gibbs_iter))
tau_plus = np.zeros(gibbs_iter)
Sigma_plus = np.zeros((rank, rank, gibbs_iter))
temp_hat = np.zeros(len(pos_test[0]))
show_iter = 500
tau = 1
tensor_hat_plus = np.zeros(sparse_tensor.shape)
tensor_new_plus = np.zeros((dim1, dim2, multi_step))
for it in range(burn_iter + gibbs_iter):
tau_ind = tau * ind
tau_sparse_tensor = tau * sparse_tensor
U = sample_factor_u(tau_sparse_tensor, tau_ind, U, V, X)
V = sample_factor_v(tau_sparse_tensor, tau_ind, U, V, X)
A, Sigma = sample_var_coefficient(X, time_lags)
X = sample_factor_x(tau_sparse_tensor, tau_ind, time_lags, U, V, X, A, inv(Sigma))
tensor_hat = np.einsum('is, js, ts -> ijt', U, V, X)
tau = np.random.gamma(1e-6 + 0.5 * np.sum(ind),
1 / (1e-6 + 0.5 * np.sum(((sparse_tensor - tensor_hat) ** 2) * ind)))
temp_hat += tensor_hat[pos_test]
if (it + 1) % show_iter == 0 and it < burn_iter:
temp_hat = temp_hat / show_iter
print('Iter: {}'.format(it + 1))
print('MAPE: {:.6}'.format(compute_mape(dense_test, temp_hat)))
print('RMSE: {:.6}'.format(compute_rmse(dense_test, temp_hat)))
temp_hat = np.zeros(len(pos_test[0]))
print()
if it + 1 > burn_iter:
U_plus[:, :, it - burn_iter] = U
V_plus[:, :, it - burn_iter] = V
A_plus[:, :, it - burn_iter] = A
Sigma_plus[:, :, it - burn_iter] = Sigma
tau_plus[it - burn_iter] = tau
tensor_hat_plus += tensor_hat
X0 = ar4cast(A, X, Sigma, time_lags, multi_step)
X_plus[:, :, it - burn_iter] = X0
tensor_new_plus += np.einsum('is, js, ts -> ijt', U, V, X0[- multi_step :, :])
tensor_hat = tensor_hat_plus / gibbs_iter
print('Imputation MAPE: {:.6}'.format(compute_mape(dense_test, tensor_hat[:, :, : dim3][pos_test])))
print('Imputation RMSE: {:.6}'.format(compute_rmse(dense_test, tensor_hat[:, :, : dim3][pos_test])))
print()
tensor_hat = np.append(tensor_hat, tensor_new_plus / gibbs_iter, axis = 2)
tensor_hat[tensor_hat < 0] = 0
return tensor_hat, U_plus, V_plus, X_plus, A_plus, Sigma_plus, tau_plus
def sample_factor_x_partial(tau_sparse_tensor, tau_ind, time_lags, U, V, X, A, Lambda_x, back_step):
"""Sampling T-by-R factor matrix X."""
dim3, rank = X.shape
tmax = np.max(time_lags)
tmin = np.min(time_lags)
d = time_lags.shape[0]
A0 = np.dstack([A] * d)
for k in range(d):
A0[k * rank : (k + 1) * rank, :, k] = 0
mat0 = Lambda_x @ A.T
mat1 = np.einsum('kij, jt -> kit', A.reshape([d, rank, rank]), Lambda_x)
mat2 = np.einsum('kit, kjt -> ij', mat1, A.reshape([d, rank, rank]))
var1 = kr_prod(V, U).T
var2 = kr_prod(var1, var1)
var3 = (var2 @ ten2mat(tau_ind[:, :, - back_step :], 2).T).reshape([rank, rank, back_step]) + Lambda_x[:, :, None]
var4 = var1 @ ten2mat(tau_sparse_tensor[:, :, - back_step :], 2).T
for t in range(dim3 - back_step, dim3):
Mt = np.zeros((rank, rank))
Nt = np.zeros(rank)
Qt = mat0 @ X[t - time_lags, :].reshape(rank * d)
index = list(range(0, d))
if t >= dim3 - tmax and t < dim3 - tmin:
index = list(np.where(t + time_lags < dim3))[0]
if t < dim3 - tmin:
Mt = mat2.copy()
temp = np.zeros((rank * d, len(index)))
n = 0
for k in index:
temp[:, n] = X[t + time_lags[k] - time_lags, :].reshape(rank * d)
n += 1
temp0 = X[t + time_lags[index], :].T - np.einsum('ijk, ik -> jk', A0[:, :, index], temp)
Nt = np.einsum('kij, jk -> i', mat1[index, :, :], temp0)
var3[:, :, t + back_step - dim3] = var3[:, :, t + back_step - dim3] + Mt
X[t, :] = mvnrnd_pre(solve(var3[:, :, t + back_step - dim3],
var4[:, t + back_step - dim3] + Nt + Qt), var3[:, :, t + back_step - dim3])
return X
def BTTF_partial(dense_tensor, sparse_tensor, init, rank, time_lags, burn_iter, gibbs_iter, multi_step = 1, gamma = 10):
"""Bayesian Temporal Tensor Factorization, BTTF."""
dim1, dim2, dim3 = sparse_tensor.shape
U_plus = init["U_plus"]
V_plus = init["V_plus"]
X_plus = init["X_plus"]
A_plus = init["A_plus"]
Sigma_plus = init["Sigma_plus"]
tau_plus = init["tau_plus"]
if np.isnan(sparse_tensor).any() == False:
ind = sparse_tensor != 0
pos_obs = np.where(ind)
elif np.isnan(sparse_tensor).any() == True:
ind = ~np.isnan(sparse_tensor)
pos_obs = np.where(ind)
sparse_tensor[np.isnan(sparse_tensor)] = 0
X_new_plus = np.zeros((dim3 + multi_step, rank, gibbs_iter))
tensor_new_plus = np.zeros((dim1, dim2, multi_step))
back_step = gamma * multi_step
for it in range(gibbs_iter):
tau_ind = tau_plus[it] * ind
tau_sparse_tensor = tau_plus[it] * sparse_tensor
X = sample_factor_x_partial(tau_sparse_tensor, tau_ind, time_lags, U_plus[:, :, it], V_plus[:, :, it],
X_plus[:, :, it], A_plus[:, :, it], inv(Sigma_plus[:, :, it]), back_step)
X0 = ar4cast(A_plus[:, :, it], X, Sigma_plus[:, :, it], time_lags, multi_step)
X_new_plus[:, :, it] = X0
tensor_new_plus += np.einsum('is, js, ts -> ijt', U_plus[:, :, it], V_plus[:, :, it], X0[- multi_step :, :])
tensor_hat = tensor_new_plus / gibbs_iter
tensor_hat[tensor_hat < 0] = 0
return tensor_hat, U_plus, V_plus, X_new_plus, A_plus, Sigma_plus, tau_plus
from ipywidgets import IntProgress
from IPython.display import display
def BTTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 10):
dim1, dim2, T = dense_tensor.shape
start_time = T - pred_step
max_count = int(np.ceil(pred_step / multi_step))
tensor_hat = np.zeros((dim1, dim2, max_count * multi_step))
f = IntProgress(min = 0, max = max_count) # instantiate the bar
display(f) # display the bar
for t in range(max_count):
if t == 0:
init = {"U": 0.1 * np.random.randn(dim1, rank),
"V": 0.1 * np.random.randn(dim2, rank),
"X": 0.1 * np.random.randn(start_time, rank)}
tensor, U, V, X_new, A, Sigma, tau = BTTF(dense_tensor[:, :, : start_time],
sparse_tensor[:, :, : start_time], init, rank, time_lags, burn_iter, gibbs_iter, multi_step)
else:
init = {"U_plus": U, "V_plus": V, "X_plus": X_new, "A_plus": A, "Sigma_plus": Sigma, "tau_plus": tau}
tensor, U, V, X_new, A, Sigma, tau = BTTF_partial(dense_tensor[:, :, : start_time + t * multi_step],
sparse_tensor[:, :, : start_time + t * multi_step], init,
rank, time_lags, burn_iter, gibbs_iter, multi_step, gamma)
tensor_hat[:, :, t * multi_step : (t + 1) * multi_step] = tensor[:, :, - multi_step :]
f.value = t
small_dense_tensor = dense_tensor[:, :, start_time : T]
pos = np.where(small_dense_tensor != 0)
print('Prediction MAPE: {:.6}'.format(compute_mape(small_dense_tensor[pos], tensor_hat[pos])))
print('Prediction RMSE: {:.6}'.format(compute_rmse(small_dense_tensor[pos], tensor_hat[pos])))
print()
return tensor_hat
Scenario setting:
import scipy.io
import warnings
warnings.simplefilter('ignore')
dense_tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')['tensor'].astype(np.float32)
sparse_tensor = dense_tensor.copy()
Model setting:
import time
rank = 30
pred_step = 7 * 24
time_lags = np.array([1, 2, 3, 24, 25, 26, 7 * 24, 7 * 24 + 1, 7 * 24 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [2, 4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=84)
Iter: 500 MAPE: nan RMSE: nan Iter: 1000 MAPE: nan RMSE: nan Imputation MAPE: nan Imputation RMSE: nan Prediction MAPE: 0.539154 Prediction RMSE: 5.10179 Running time: 1525 seconds Prediction time horizon (delta) = 4.
IntProgress(value=0, max=42)
Iter: 500 MAPE: nan RMSE: nan Iter: 1000 MAPE: nan RMSE: nan Imputation MAPE: nan Imputation RMSE: nan Prediction MAPE: 0.546757 Prediction RMSE: 5.17932 Running time: 1441 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=28)
Iter: 500 MAPE: nan RMSE: nan Iter: 1000 MAPE: nan RMSE: nan Imputation MAPE: nan Imputation RMSE: nan Prediction MAPE: 0.56155 Prediction RMSE: 5.25294 Running time: 1312 seconds
Scenario setting:
import scipy.io
import numpy as np
np.random.seed(1000)
dense_tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')['tensor'].astype(np.float32)
dim = dense_tensor.shape
missing_rate = 0.4 # Random missing (RM)
sparse_tensor = dense_tensor * np.round(np.random.rand(dim[0], dim[1], dim[2]) + 0.5 - missing_rate)
Model setting:
import time
rank = 30
pred_step = 7 * 24
time_lags = np.array([1, 2, 3, 24, 25, 26, 7 * 24, 7 * 24 + 1, 7 * 24 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [2, 4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=84)
Iter: 500 MAPE: 0.476403 RMSE: 5.1438 Iter: 1000 MAPE: 0.489101 RMSE: 5.30296 Imputation MAPE: 0.490007 Imputation RMSE: 5.30489 Prediction MAPE: 0.482274 Prediction RMSE: 5.95222 Running time: 1410 seconds Prediction time horizon (delta) = 4.
IntProgress(value=0, max=42)
Iter: 500 MAPE: 0.503474 RMSE: 5.15323 Iter: 1000 MAPE: 0.505224 RMSE: 5.30834 Imputation MAPE: 0.505934 Imputation RMSE: 5.33052 Prediction MAPE: 0.537489 Prediction RMSE: 5.32858 Running time: 1362 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=28)
Iter: 500 MAPE: 0.50377 RMSE: 5.14515 Iter: 1000 MAPE: 0.521959 RMSE: 5.29029 Imputation MAPE: 0.523315 Imputation RMSE: 5.30101 Prediction MAPE: 0.561496 Prediction RMSE: 5.4019 Running time: 1328 seconds
Scenario setting:
import scipy.io
import numpy as np
np.random.seed(1000)
dense_tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')['tensor'].astype(np.float32)
dim = dense_tensor.shape
missing_rate = 0.6 # Random missing (RM)
sparse_tensor = dense_tensor * np.round(np.random.rand(dim[0], dim[1], dim[2]) + 0.5 - missing_rate)
Model setting:
import time
rank = 30
pred_step = 7 * 24
time_lags = np.array([1, 2, 3, 24, 25, 26, 7 * 24, 7 * 24 + 1, 7 * 24 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [2, 4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=84)
Iter: 500 MAPE: 0.485907 RMSE: 5.31316 Iter: 1000 MAPE: 0.494965 RMSE: 5.57816 Imputation MAPE: 0.495821 Imputation RMSE: 5.60975 Prediction MAPE: 0.488657 Prediction RMSE: 6.72835 Running time: 1554 seconds Prediction time horizon (delta) = 4.
IntProgress(value=0, max=42)
Iter: 500 MAPE: 0.508972 RMSE: 5.29687 Iter: 1000 MAPE: 0.523079 RMSE: 5.52528 Imputation MAPE: 0.510566 Imputation RMSE: 5.54642 Prediction MAPE: 0.529534 Prediction RMSE: 5.43116 Running time: 1382 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=28)
Iter: 500 MAPE: 0.513433 RMSE: 5.32022 Iter: 1000 MAPE: 0.535674 RMSE: 5.49721 Imputation MAPE: 0.535587 Imputation RMSE: 5.50608 Prediction MAPE: 0.552921 Prediction RMSE: 5.46179 Running time: 1375 seconds
Scenario setting:
import scipy.io
import numpy as np
np.random.seed(1000)
dense_tensor = scipy.io.loadmat('../datasets/NYC-data-set/tensor.mat')['tensor']
dim = dense_tensor.shape
nm_tensor = np.random.rand(dim[0], dim[1], dim[2])
missing_rate = 0.4 # Non-random missing (NM)
binary_tensor = np.zeros(dense_tensor.shape)
for i1 in range(dim[0]):
for i2 in range(dim[1]):
for i3 in range(61):
binary_tensor[i1, i2, i3 * 24 : (i3 + 1) * 24] = np.round(nm_tensor[i1, i2, i3] + 0.5 - missing_rate)
sparse_tensor = dense_tensor * binary_tensor
Model setting:
import time
rank = 30
pred_step = 7 * 24
time_lags = np.array([1, 2, 3, 24, 25, 26, 7 * 24, 7 * 24 + 1, 7 * 24 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [2, 4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=84)
Iter: 500 MAPE: 0.48005 RMSE: 5.33599 Iter: 1000 MAPE: 0.483847 RMSE: 5.50755 Imputation MAPE: 0.483285 Imputation RMSE: 5.50448 Prediction MAPE: 0.481164 Prediction RMSE: 6.14635 Running time: 1492 seconds Prediction time horizon (delta) = 4.
IntProgress(value=0, max=42)
Iter: 500 MAPE: 0.515712 RMSE: 5.27328 Iter: 1000 MAPE: 0.528275 RMSE: 5.40139 Imputation MAPE: 0.532998 Imputation RMSE: 5.44158 Prediction MAPE: 0.548801 Prediction RMSE: 5.31528 Running time: 1446 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=28)
Iter: 500 MAPE: 0.508245 RMSE: 5.30643 Iter: 1000 MAPE: 0.522899 RMSE: 5.47282 Imputation MAPE: 0.525004 Imputation RMSE: 5.43607 Prediction MAPE: 0.548346 Prediction RMSE: 5.35651 Running time: 1308 seconds
Scenario setting:
import numpy as np
import warnings
warnings.simplefilter('ignore')
dense_tensor = np.load('../datasets/Temperature-data-set/tensor.npy').astype(np.float32)
pos = np.where(dense_tensor[:, 0, :] > 50)
dense_tensor[pos[0], :, pos[1]] = 0
sparse_tensor = dense_tensor.copy()
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [2, 4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=60)
Iter: 500 MAPE: nan RMSE: nan Iter: 1000 MAPE: nan RMSE: nan Imputation MAPE: nan Imputation RMSE: nan Prediction MAPE: 0.0285055 Prediction RMSE: 0.924666 Running time: 1222 seconds Prediction time horizon (delta) = 4.
IntProgress(value=0, max=30)
Iter: 500 MAPE: nan RMSE: nan Iter: 1000 MAPE: nan RMSE: nan Imputation MAPE: nan Imputation RMSE: nan Prediction MAPE: 0.0258453 Prediction RMSE: 0.83646 Running time: 1127 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=20)
Iter: 500 MAPE: nan RMSE: nan Iter: 1000 MAPE: nan RMSE: nan Imputation MAPE: nan Imputation RMSE: nan Prediction MAPE: 0.0281043 Prediction RMSE: 0.908133 Running time: 1100 seconds
Scenario setting:
import numpy as np
np.random.seed(1000)
dense_tensor = np.load('../datasets/Temperature-data-set/tensor.npy').astype(np.float32)
pos = np.where(dense_tensor[:, 0, :] > 50)
dense_tensor[pos[0], :, pos[1]] = 0
random_tensor = np.random.rand(dense_tensor.shape[0], dense_tensor.shape[1], dense_tensor.shape[2])
missing_rate = 0.4
## Random missing (RM)
binary_tensor = np.round(random_tensor + 0.5 - missing_rate)
sparse_tensor = dense_tensor.copy()
sparse_tensor[binary_tensor == 0] = np.nan
sparse_tensor[sparse_tensor == 0] = np.nan
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
multi_step = 2
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=60)
Iter: 500 MAPE: 0.0144821 RMSE: 0.482159 Iter: 1000 MAPE: 0.0143581 RMSE: 0.476049 Imputation MAPE: 0.0142973 Imputation RMSE: 0.474065 Prediction MAPE: 0.0231417 Prediction RMSE: 0.758236 Running time: 1349 seconds
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 4.
IntProgress(value=0, max=30)
Iter: 500 MAPE: 0.014985 RMSE: 0.498448 Iter: 1000 MAPE: 0.0147059 RMSE: 0.487322 Imputation MAPE: 0.0145875 Imputation RMSE: 0.483447 Prediction MAPE: 0.0249583 Prediction RMSE: 0.809355 Running time: 1179 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=20)
Iter: 500 MAPE: 0.0150258 RMSE: 0.499392 Iter: 1000 MAPE: 0.0148409 RMSE: 0.491953 Imputation MAPE: 0.0147268 Imputation RMSE: 0.488867 Prediction MAPE: 0.0257521 Prediction RMSE: 0.840748 Running time: 1240 seconds
Scenario setting:
import numpy as np
np.random.seed(1000)
dense_tensor = np.load('../datasets/Temperature-data-set/tensor.npy').astype(np.float32)
pos = np.where(dense_tensor[:, 0, :] > 50)
dense_tensor[pos[0], :, pos[1]] = 0
random_tensor = np.random.rand(dense_tensor.shape[0], dense_tensor.shape[1], dense_tensor.shape[2])
missing_rate = 0.6
## Random missing (RM)
binary_tensor = np.round(random_tensor + 0.5 - missing_rate)
sparse_tensor = dense_tensor.copy()
sparse_tensor[binary_tensor == 0] = np.nan
sparse_tensor[sparse_tensor == 0] = np.nan
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
multi_step = 2
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=60)
Iter: 500 MAPE: 0.0145609 RMSE: 0.485104 Iter: 1000 MAPE: 0.0144643 RMSE: 0.4799 Imputation MAPE: 0.0144504 Imputation RMSE: 0.479404 Prediction MAPE: 0.461062 Prediction RMSE: 12.2373 Running time: 1320 seconds
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
multi_step = 2
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 40)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=60)
Iter: 500 MAPE: 0.0147971 RMSE: 0.493339 Iter: 1000 MAPE: 0.0146462 RMSE: 0.487619 Imputation MAPE: 0.0145599 Imputation RMSE: 0.484679 Prediction MAPE: 0.0235214 Prediction RMSE: 0.76378 Running time: 1645 seconds
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 4.
IntProgress(value=0, max=30)
Iter: 500 MAPE: 0.0147971 RMSE: 0.493339 Iter: 1000 MAPE: 0.0146462 RMSE: 0.487619 Imputation MAPE: 0.0145574 Imputation RMSE: 0.484643 Prediction MAPE: 0.0243273 Prediction RMSE: 0.787939 Running time: 1199 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=20)
Iter: 500 MAPE: 0.0149478 RMSE: 0.499435 Iter: 1000 MAPE: 0.0146285 RMSE: 0.484588 Imputation MAPE: 0.0145979 Imputation RMSE: 0.48322 Prediction MAPE: 0.0322035 Prediction RMSE: 1.03053 Running time: 1261 seconds
Scenario setting:
import numpy as np
np.random.seed(1000)
dense_tensor = np.load('../datasets/Temperature-data-set/tensor.npy').astype(np.float32)
pos = np.where(dense_tensor[:, 0, :] > 50)
dense_tensor[pos[0], :, pos[1]] = 0
random_tensor = np.random.rand(dense_tensor.shape[0], dense_tensor.shape[1], int(dense_tensor.shape[2] / 3))
missing_rate = 0.4
## Non-random missing (NM)
binary_tensor = np.zeros(dense_tensor.shape)
for i1 in range(dense_tensor.shape[0]):
for i2 in range(dense_tensor.shape[1]):
for i3 in range(int(dense_tensor.shape[2] / 3)):
binary_tensor[i1, i2, i3 * 3 : (i3 + 1) * 3] = np.round(random_tensor[i1, i2, i3] + 0.5 - missing_rate)
sparse_tensor = dense_tensor.copy()
sparse_tensor[binary_tensor == 0] = np.nan
sparse_tensor[sparse_tensor == 0] = np.nan
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
multi_step = 2
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 40)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()s
Prediction time horizon (delta) = 2.
IntProgress(value=0, max=60)
Iter: 500 MAPE: 0.0148901 RMSE: 0.494736 Iter: 1000 MAPE: 0.0146504 RMSE: 0.485204 Imputation MAPE: 0.0145857 Imputation RMSE: 0.482443 Prediction MAPE: 0.0231671 Prediction RMSE: 0.753352 Running time: 1435 seconds
import time
rank = 30
pred_step = 10 * 12
time_lags = np.array([1, 2, 3, 12, 13, 14, 2 * 12, 2 * 12 + 1, 2 * 12 + 2])
burn_iter = 1000
gibbs_iter = 200
for multi_step in [4, 6]:
start = time.time()
print('Prediction time horizon (delta) = {}.'.format(multi_step))
tensor_hat = BTTF_forecast(dense_tensor, sparse_tensor, pred_step, multi_step,
rank, time_lags, burn_iter, gibbs_iter, gamma = 30)
end = time.time()
print('Running time: %d seconds'%(end - start))
print()
Prediction time horizon (delta) = 4.
IntProgress(value=0, max=30)
Iter: 500 MAPE: 0.0149142 RMSE: 0.495803 Iter: 1000 MAPE: 0.0145458 RMSE: 0.482243 Imputation MAPE: 0.0144844 Imputation RMSE: 0.47983 Prediction MAPE: 0.0241001 Prediction RMSE: 0.788512 Running time: 1308 seconds Prediction time horizon (delta) = 6.
IntProgress(value=0, max=20)
Iter: 500 MAPE: 0.014927 RMSE: 0.497119 Iter: 1000 MAPE: 0.0144465 RMSE: 0.479791 Imputation MAPE: 0.0144483 Imputation RMSE: 0.478898 Prediction MAPE: 0.0244661 Prediction RMSE: 0.799278 Running time: 1346 seconds