real comparison code

pybasicbayes/distributions/dynglm_optimisation_test.py

+"""
+Need to find out whether loglikelihood is computed correctly.
+
+Or whether a bug here allows states to invade each other more easily.
+We'll test this by comparing to pymc results.
+"""
+import numpy as np
+from dynamic_glm import Dynamic_GLM
+import pickle
+import matplotlib.pyplot as plt
+
+# Testing of Dynamic_GLM implementation
+np.set_printoptions(suppress=True)
+
+seed = np.random.randint(10000)  # 215
+print(seed)
+seed = 5489
+np.random.seed(seed)
+
+
+T = 16
+n_regressors = 3
+step_size = 0.2
+Q = np.tile(np.eye(n_regressors), (T, 1, 1))
+test = Dynamic_GLM(n_regressors=n_regressors, T=T, P_0=4 * np.eye(n_regressors), Q=Q * step_size, prior_mean=np.zeros(n_regressors))
+w = np.zeros(n_regressors)
+
+test_points = [0]
+predictors = []
+a, b = np.zeros(1000), np.zeros(1000)
+a[250:500] = 1
+a[750:1000] = 1
+b[500:] = 1
+for _ in range(T):
+    if _ in [0, 1, 8, 9, 10, 11, 12, 14, 15]:
+        predictors.append(np.empty((0, n_regressors)))
+        continue
+    t = 1000
+    pred = np.empty((t, n_regressors))
+    pred[:, 0] = a
+    pred[:, 1] = b
+    pred[:, 2] = 1
+    predictors.append(pred)
+sample = test.rvs(predictors, list(range(T)))
+pickle.dump(sample, open('test_data', 'wb'))
+pickle.dump(test, open('truth', 'wb'))
+
+plt.plot(test.weights)
+plt.show()