added random effect rho for all DGLMs

pybats/analysis.py

@@ -11,7 +11,8 @@ def analysis(Y, X, k, forecast_start, forecast_end,
              model_prior = None, prior_length=20, ntrend=1,
              dates = None, holidays = [],
              seasPeriods = [], seasHarmComponents = [],
-             ret=['model', 'forecast'], mean_only = False,
+             ret=['model', 'forecast'],
+             mean_only = False, forecast_path = False,
              **kwargs):
     """
     This is a helpful function to run a standard analysis. The function will:
@@ -53,13 +54,12 @@ def analysis(Y, X, k, forecast_start, forecast_end,
     :param mean_only: True/False - return the mean only when forecasting, instead of samples?
     :param kwargs: Further key word arguments to define the model prior. Common arguments include discount factors deltrend, delregn, delseas, and delhol.
 
-    :return: Default is a list with [model, forecast_samples]. forecast_samples will have dimension (nsamps by forecast_length by k), where forecast_length is the number of time steps between forecast_start and forecast_end. The output can be customized with the 'ret' parameter, which is a list of values to return.
+    :return: Default is a list with [model, forecast_samples]. forecast_samples will have dimension (nsamps by forecast_length by k), where forecast_length is the number of time steps between forecast_start and forecast_end. The output can be customized with the 'ret' parameter, which is a list of objects to return.
     """
 
     # Add the holiday indicator variables to the regression matrix
     nhol = len(holidays)
-    if nhol > 0:
-        X = define_holiday_regressors(X, dates, holidays)
+    X = define_holiday_regressors(X, dates, holidays)
 
     if model_prior is None:
         mod = define_dglm(Y, X, family=family, n=n, prior_length=prior_length, ntrend=ntrend, nhol=nhol,
@@ -109,17 +109,20 @@ def analysis(Y, X, k, forecast_start, forecast_end,
                 print('beginning forecasting')
 
             if ret.__contains__('forecast'):
-                if family == "binomial":
-                    forecast[:, t - forecast_start, :] = np.array(list(map(
-                        lambda k, n, x:
-                        mod.forecast_marginal(k=k, n=n, X=x, nsamps=nsamps, mean_only=mean_only),
-                        horizons, n[t + horizons - 1], X[t + horizons - 1, :]))).squeeze().T.reshape(-1, k)  # .reshape(-1, 1)
+                if forecast_path:
+                    forecast[:, t - forecast_start, :] = mod.forecast_path(k=k, X = X[t + horizons - 1, :], nsamps=nsamps)
                 else:
-                    # Get the forecast samples for all the items over the 1:k step ahead marginal forecast distributions
-                    forecast[:, t - forecast_start, :] = np.array(list(map(
-                        lambda k, x:
-                        mod.forecast_marginal(k=k, X=x, nsamps=nsamps, mean_only=mean_only),
-                        horizons, X[t + horizons - 1, :]))).squeeze().T.reshape(-1, k)#.reshape(-1, 1)
+                    if family == "binomial":
+                        forecast[:, t - forecast_start, :] = np.array(list(map(
+                            lambda k, n, x:
+                            mod.forecast_marginal(k=k, n=n, X=x, nsamps=nsamps, mean_only=mean_only),
+                            horizons, n[t + horizons - 1], X[t + horizons - 1, :]))).squeeze().T.reshape(-1, k)  # .reshape(-1, 1)
+                    else:
+                        # Get the forecast samples for all the items over the 1:k step ahead marginal forecast distributions
+                        forecast[:, t - forecast_start, :] = np.array(list(map(
+                            lambda k, x:
+                            mod.forecast_marginal(k=k, X=x, nsamps=nsamps, mean_only=mean_only),
+                            horizons, X[t + horizons - 1, :]))).squeeze().T.reshape(-1, k)#.reshape(-1, 1)
 
 
         # Now observe the true y value, and update:

pybats/dglm.py

@@ -1,3 +1,9 @@
+"""
+Core class of dynamic generalized linear models (DGLMs).
+
+The module also includes definitions for sub classes for poisson, bernoulli, and binomial DGLMs and standard dynamic linear models (DLMs).
+"""
+
 import numpy as np
 import scipy as sc
 from collections.abc import Iterable
@@ -19,7 +25,8 @@ class dglm:
 
     Children include Poisson, Bernoulli, and Binomial DGLMs, as well as the Normal DLM.
     """
-
+
+
     def __init__(self,
                  a0=None,
                  R0=None,
@@ -30,6 +37,7 @@ def __init__(self,
                  seasHarmComponents=[],
                  deltrend=1, delregn=1,
                  delhol=1, delseas=1,
+                 rho=1,
                  interpolate=True,
                  adapt_discount=False,
                  adapt_factor=0.5,
@@ -49,9 +57,10 @@ def __init__(self,
         :param delhol: Discount factor on holiday components (currently deprecated)
         :param delseas: Discount factor on seasonal components
         :param interpolate: Whether to use interpolation for conjugate parameters (provides a computational speedup)
-        :param adapt_discount: What method of discount adaption. False = None, 'positive-regn' = only discount if regression information is available, 'info' = information based,\
+        :param adapt_discount: What method of discount adaption. False = None, 'positive_regn' = only discount if regression information is available, 'info' = information based,\
         :param adapt_factor: If adapt_discount='info', then a higher value adapt_factor leads to a quicker adaptation (with less discounting) on overly uncertain parameters
         :param discount_forecast: Whether to use discounting when forecasting
+        :return An object of class dglm
         """
 
         # Setting up trend F, G matrices
@@ -72,27 +81,28 @@ def __init__(self,
             Ftrend = np.array([1, 0]).reshape(-1, 1)
 
         # Setting up regression F, G matrices
+        self.iregn = list(range(i, i + nregn))
         if nregn == 0:
             Fregn = np.empty([0]).reshape(-1, 1)
             Gregn = np.empty([0, 0])
         else:
             Gregn = np.identity(nregn)
             Fregn = np.ones([nregn]).reshape(-1, 1)
-            self.iregn = list(range(i, i + nregn))
             i += nregn
 
         # Setting up holiday F, G matrices (additional regression indicators components)
+        self.ihol = list(range(i, i + nhol))
+        self.iregn.extend(self.ihol)  # Adding on to the self.iregn
         if nhol == 0:
             Fhol = np.empty([0]).reshape(-1, 1)
             Ghol = np.empty([0, 0])
         else:
             Ghol = np.identity(nhol)
             Fhol = np.ones([nhol]).reshape(-1, 1)
-            self.ihol = list(range(i, i + nhol))
-            self.iregn.extend(self.ihol)  # Adding on to the self.iregn
             i += nhol
 
         # Setting up seasonal F, G matrices
+        self.iseas = []
         if len(seasPeriods) == 0:
             Fseas = np.empty([0]).reshape(-1, 1)
             Gseas = np.empty([0, 0])
@@ -106,7 +116,6 @@ def __init__(self,
             Fseas = np.zeros([nseas]).reshape(-1, 1)
             Gseas = np.zeros([nseas, nseas])
             idx = 0
-            self.iseas = []
             for harmComponents in seasHarmComponents:
                 self.iseas.append(list(range(i, i + 2 * len(harmComponents))))
                 i += 2 * len(harmComponents)
@@ -126,15 +135,22 @@ def __init__(self,
         self.delhol = delhol
         self.delseas = delseas
 
+        # Random effect to inflate variance (if rho < 1)
+        self.rho = rho
+
         self.ntrend = ntrend
         self.nregn = nregn + nhol  # Adding on nhol
         self.nregn_exhol = nregn
         self.nhol = nhol
         self.nseas = nseas
 
+        self.adapt_discount = adapt_discount
+        self.k = adapt_factor
+
         # Set up discount matrix
         self.discount_forecast = discount_forecast
         Discount = self.build_discount_matrix()
+        self.Discount = Discount
 
         self.param1 = 2  # Random initial guess
         self.param2 = 2  # Random initial guess
@@ -143,13 +159,10 @@ def __init__(self,
         self.seasHarmComponents = seasHarmComponents
         self.F = F
         self.G = G
-        self.Discount = Discount
         self.a = a0.reshape(-1, 1)
         self.R = R0
         self.t = 0
         self.interpolate = interpolate
-        self.adapt_discount = adapt_discount
-        self.k = adapt_factor
         self.W = self.get_W()
 
     def build_discount_matrix(self, X=None):
@@ -162,7 +175,7 @@ def build_discount_matrix(self, X=None):
         component_discounts = np.ones([p, p])
         i = 0 # this will be the offset of the current block
         for discount_pair, n in zip([('std', self.deltrend), ('regn', self.delregn),
-                                     ('regn', self.delhol),('std', self.delseas)],
+                                     ('hol', self.delhol),('std', self.delseas)],
                                     [self.ntrend, self.nregn_exhol, self.nhol, self.nseas]):
             discount_type, discount = discount_pair
             if n > 0:
@@ -177,9 +190,15 @@ def build_discount_matrix(self, X=None):
                     component_discounts[i:(i+n), i:(i+n)] = discount
 
                 # overwrite with ones if doing the positive logic
-                if X is not None and self.adapt_discount == 'positive_regn' and discount_type == 'regn':
-                    # offset of the regression params
-                    regn_i = 0
+                if X is not None and self.adapt_discount == 'positive_regn' and (discount_type == 'regn' or discount_type == 'hol'):
+                    if not isinstance(X, Iterable):
+                        X = [X]
+
+                    if discount_type == 'regn':
+                        # offset of the regression params
+                        regn_i = 0
+                    elif discount_type == 'hol':
+                        regn_i = self.nregn_exhol
                     # look through the regression params and set that slice on the
                     # discount to 1 if 0
                     for j in range(n):
@@ -198,12 +217,14 @@ def update(self, y=None, X=None, **kwargs):
         """
         Update the DGLM state vector mean and covariance after observing 'y', with covariates 'X'.
 
+        Example usage:
+
         >>> mod.update(y[t], X[t])
 
         Posterior mean and covariance:
 
         >>> [mod.m, mod.C]
-        
+
         You can also access the state vector *prior* mean and variance for the next time step.
         The state vector prior mean will be the same as the posterior mean. The variance will be larger, due to discounting.
 
@@ -215,7 +236,7 @@ def update(self, y=None, X=None, **kwargs):
         :return: No output; DGLM state vector is updated.
         """
 
-        update(self, y, X, **kwargs)
+        update(self, y, X)
 
     def forecast_marginal(self, k, X=None, nsamps=1, mean_only=False, state_mean_var=False):
         """
@@ -252,7 +273,7 @@ def forecast_state_mean_and_var(self, k, X = None):
         return forecast_state_mean_and_var(self, k, X)
 
     def get_mean_and_var(self, F, a, R):
-        mean, var = F.T @ a, F.T @ R @ F
+        mean, var = F.T @ a, F.T @ R @ F / self.rho
         return np.ravel(mean)[0], np.ravel(var)[0]
 
     def get_W(self, X=None):
@@ -271,11 +292,6 @@ def get_W(self, X=None):
 
 class bern_dglm(dglm):
 
-    def get_mean_and_var(self, F, a, R):
-        ft, qt = F.T @ a, F.T @ R @ F
-        ft, qt = ft.flatten()[0], qt.flatten()[0]
-        return ft, qt
-
     def get_conjugate_params(self, ft, qt, alpha_init, beta_init):
         # Choose conjugate prior, beta, and match mean & variance
         return bern_conjugate_params(ft, qt, alpha_init, beta_init, interp=self.interpolate)
@@ -329,13 +345,6 @@ def get_prior_var(self, alpha, beta):
 
 class pois_dglm(dglm):
 
-    def __init__(self, *args, rho=1, **kwargs):
-        self.rho = rho
-        super().__init__(*args, **kwargs)
-
-    def get_mean_and_var(self, F, a, R):
-        return F.T @ a, F.T @ R @ F / self.rho
-
     def get_conjugate_params(self, ft, qt, alpha_init, beta_init):
         # Choose conjugate prior, gamma, and match mean & variance
         return pois_conjugate_params(ft, qt, alpha_init, beta_init, interp=self.interpolate)
@@ -406,12 +415,12 @@ def simulate(self, mean, var, nsamps):
         return mean + np.sqrt(var) * np.random.standard_t(self.n, size=[nsamps])
 
     def simulate_from_sampling_model(self, mean, var, nsamps):
-        return np.random.normal(mean, var, nsamps)
+        return np.random.normal(mean, np.sqrt(var), nsamps)
 
-    def update(self, y=None, X=None):
+    def update(self, y=None, X=None, **kwargs):
         update_dlm(self, y, X)
 
-    def forecast_path(self, k, X=None, nsamps=1):
+    def forecast_path(self, k, X=None, nsamps=1, **kwargs):
         return forecast_path_dlm(self, k, X, nsamps)
 
 class bin_dglm(dglm):

pybats/loss_functions.py

@@ -64,7 +64,6 @@ def MAPE(y, f):
     :param f: Point forecast vector
     :return: Mean absolute percent error (MAPE)
     """
-
     y = np.ravel(y)
     f = np.ravel(f)
     return 100*np.mean(np.abs((y - f)) / y)
@@ -89,6 +88,7 @@ def WAPE(y, f):
     :param f: Point forecast vector
     :return: Weighted absolute percent error (WAPE)
     """
+
     y = np.ravel(y)
     f = np.ravel(f)
     return 100*np.sum(np.abs(y-f)) / np.sum(y)
@@ -112,6 +112,7 @@ def WAFE(y, f):
     :param f: Point forecast vector
     :return: Weighted absolute forecast error (WAFE)
     """
+
     y = np.ravel(y)
     f = np.ravel(f)
     return 100*np.sum(np.abs(y-f)) / ((np.sum(y) + np.sum(f))/2)
@@ -136,6 +137,7 @@ def ZAPE(y, f):
     :param f: Point forecast vector
     :return: The mean Zero-Adjusted absolute percent error (ZAPE)
     """
+
     y = np.ravel(y)
     f = np.ravel(f)
     nonzeros = y.nonzero()[0]

pybats/plot.py

@@ -30,6 +30,8 @@ def plot_data_forecast(fig, ax, y, f, samples, dates, linewidth=1, linecolor='b'
     lower = np.percentile(samples, [alpha], axis=0).reshape(-1)
     ax.fill_between(dates, upper, lower, alpha=.3, color=linecolor)
 
+    if kwargs.get('xlim') is None:
+        kwargs.update({'xlim':[dates[0], dates[-1]]})
     ax = ax_style(ax, **kwargs)
 
     # If dates are actually dates, then format the dates on the x-axis

pybats/point_forecast.py

@@ -40,7 +40,7 @@ def weighted_quantile(samp, weights, quantile=0.5):
     return np.round((upper + lower) / 2)
 
 
-# Optimal for APE. Always less than the median. Returns nan if some samples are 0.
+# Optimal for APE. Always less than the median. Ignores samples that are 0.
 def m_one_median(samps):
     """
     Return the (-1)-median point forecasts, given samples from the analysis function.
@@ -55,6 +55,8 @@ def m_one_median(samp):
         weights = 1/samp[nz]
         norm = np.sum(weights)
         weights = weights/norm
+        if len(nz) < 5:
+            print('hello')
         return weighted_quantile(samp[nz], weights)
 
     forecast = np.apply_along_axis(m_one_median, 0, samps)
@@ -63,7 +65,9 @@ def m_one_median(samp):
 
 
 # Here we get the joint one_median, where the rows are forecast samples
-# Assume that the forecast is 'joint' across the last dimension
+# Assume that the forecast is 'joint' across the second dimension
+# This is optimal for the WAPE loss, where the denominator in the WAPE score is the sum over the second dimension
+# If the forecast samples are from a standard analysis function, that will be the sum over all forecast dates
 def joint_m_one_median(samps):
 
     def joint_m_one_median(samp):
@@ -188,3 +192,36 @@ def constrained_joint_m_one_median(samp, F):
                              samps.transpose([1, 0, 2]),
                              F)))[:,0,:]
 
+
+# Optimal for ZAPE. Always less than the (-1)-median.
+def zape_point_estimate(samps):
+    """
+    Return the optimal point forecast for ZAPE loss, given samples from the analysis function.
+
+    This forecast is theoretically optimal for minimizing ZAPE loss, which is defined as:
+
+    .. math:: ZAPE(y, f) = \\frac{1}{n} \sum_{i=1:n} I(y_i = 0) * f_i + I(y_i = 1) * |y_i-f_i| / y_i
+
+    :param samps: Forecast samples, returned from the analysis function. Will have 3-dimensions (nsamps * time * forecast horizon)
+    :return: Array of (-1)-median forecasts. Will have dimension (time * forecast horizon)
+    """
+    def est_c_hat(samp):
+        nz = samp.nonzero()[0]
+        weights = 1/samp[nz]
+        c_hat = 1 / (1/len(nz) * np.sum(weights))
+        return c_hat
+
+    def zape_point_est(samp):
+        nz = samp.nonzero()[0]
+        pi_0 = len(nz) / len(samp) # probability of 0
+        weights = 1 / samp[nz]
+        norm = np.sum(weights)
+        weights = weights / norm
+        c_hat = est_c_hat(samp)
+        quantile = (1 - c_hat*pi_0)/2
+
+        return weighted_quantile(samp[nz], weights, quantile)
+
+    forecast = np.apply_along_axis(m_one_median, 0, samps)
+
+    return forecast

pybats/seasonal.py

@@ -1,5 +1,4 @@
 import numpy as np
-
 from .forecast import forecast_aR, forecast_R_cov