<jupyter_start><jupyter_text>Probabilistic forecasting and multi-layer perceptron.This is the second part of our introduction to forecasting. This notebook starts with a short introduction to __maximum likelihood estimation__ and to __general linear models__. These concepts are used to build a model to generate a __probabilistic forecast of count data__ using a Poisson distribution. At this point you might want to take a quick look at the crash course on probability and statistics in chapter 1. The second part of this notebook is dedicated building a __multi-layer perceptron (MLP) for forecasting__, using the same data and Poisson distribution as in the first part of this notebook.<jupyter_code>from __future__ import print_function
import mxnet as mx
import numpy as np
from mxnet import nd, autograd, ndarray
from mxnet import gluon
from mxnet.gluon.loss import Loss
from math import factorial
import pandas as pd
%matplotlib inline
from matplotlib import pyplot as plt<jupyter_output><empty_output><jupyter_text>Maximum likelihood estimation In the previous notebook we used linear models of the form $Y_t = X_t\beta + \epsilon_t$. The parameters $\beta$ are estimated by minimizing the mean squared error $$\frac{1}{T}\sum_t \hat\epsilon_t^2 = \frac{1}{T}\sum_t (Y_t - X_t\hat\beta)^2.$$ The model as a distributionBy making an extra assumption we can view this model and loss function in a more principled way. Let's assume that __the innovations $\epsilon_t$ are independently drawn from a Gaussian distribution__: $\epsilon_t \sim \mathcal{N}(0,\sigma^2)$. With this extra assumption we can re-write the model as a condition distribution, which is a didtribution that depends on variable $X$ that are determined outside the model. Since we assume the errors $\epsilon$ are Gaussian: $$y_t = X_t\beta + \epsilon_t \implies Y_t \sim \mathcal{N}(X_t\beta,\sigma_2),$$meaning that $y_t$ __follows a Gaussian distribution__ with mean $X_t\beta$ and variance $\sigmaˆ2$. The mean and variance of the distribution are functions of the unknown parameters $\beta$ and $\sigma$. Had we known these parameters we would have known the distribution of $Y_t$ fully, and so we could compute $$Pr(y_t<\mu \mid X_t, \beta, \sigma)=p_\mu,\ \ \mu\in\mathbb{R}.$$$Pr(y_t<\mu \mid X_t, \beta, \sigma)=p_\mu$ is __the probability that we have observed the data $y_t$__ conditional on some value of the parameters ($\beta, \sigma$) and on $X_t$. The problem can be turned around, for what value of the parameters is the probabily that we observe $y_t$ maximized?   Likelihood function and Maximum Likelihood estimationThe __likelihood function__ is defined as:$$ L(\beta; y_t, X_t, \sigma) = Pr(y_t \mid X_t, \beta, \sigma).$$The __maximum likelihood estimator__ (MLE) of our parameters is the maximizer of $L(\beta; y_t, X_t, \sigma)$: $(\hat\beta,\hat\sigma)=arg\,max_{\beta,\sigma}\left(L(\beta; y_t, X_t, \sigma)\right)$, though in practice we generally minimize the opposit of the logarithm of the the likelihood function, the __log likelihood function__ $-\log L(\beta; y_t, X_t, \sigma)=-\log\left(L(\beta)\right)$. I have skipped entirely a lot of extremely important points in this very brief introduction to maximum likelihood, the goal was just to mention the main concepts. There are lots of excellent resourcres out there to read a proper discussion of these topics.   MLE with GluonTo enable MLE with `gluon` we need to define the log likelihood loss function.Let us look at a concrete example. A random variable following a Bernoulli distribution takes value 0 with probability p, or 1 with probability 1-p. For a sample of $N$ independent Bernoulli distributed random variables $y_i$, the log likelihood is $L(p)=\frac{1}{N}\sum_i\left(p y_i + (1-p)(1-y_i)\right)$. Few likelihood functions are included in the current release of `gluon`, but defining a custom loss function is easy:<jupyter_code># Define the Bernoulli loss function

class Bernoulli(Loss):
    """Calculates the Bernoulli loss function
    Output and label must have the same shape. This is a scalar loss function.

    Parameters
    ----------
    weight : float or None
        Global scalar weight for loss.
    sample_weight : Symbol or None
        Per sample weighting. Must be broadcastable to
        the same shape as loss. For example, if loss has
        shape (64, 10) and you want to weight each sample
        in the batch, `sample_weight` should have shape (64, 1).
    batch_axis : int, default 0
        The axis that represents mini-batch.
    """
    def __init__(self, weight=None, batch_axis=0, **kwargs):
        super(Bernoulli, self).__init__(weight, batch_axis, **kwargs)

    def hybrid_forward(self, F, output, label, sample_weight=None):
        label = _reshape_label_as_output(F, output, label)
        loss = -1 * (1 - output) * (1-label) - output * label
        loss = _apply_weighting(F, loss, self._weight, sample_weight)
        return F.mean(loss, axis=self._batch_axis, exclude=True)
    
bernoulli_log_lik = Bernoulli()<jupyter_output><empty_output><jupyter_text>Gaussian MLEWhen we have a Gaussian linear model, it's quite easy to show that the values $\hat\beta$ and $\hat\sigma$, the MLEs, that minimize the log likelihood function are:   1. The value of $\hat\beta$ that minimimze $\frac{1}{T}\sum_t\left( y_t-X_t\beta\right)ˆ2$, that is, __the minimizer of the L2 loss__.   2. $\hat\sigma^2 = \frac{1}{T}\sum_t \epsilon_t^2 = \frac{1}{T}\sum_t (y_t - X_t\hat\beta)^2$, the same estimator we used earlier. In the Gaussian case the maximum likelihood estimator is equivalent to minimizing the L2 loss. This implies that to do Gaussian MLE we do not need to setup a particular function, we can just use L2 loss: `gaussian_log_lik = gluon.loss.L2Loss()`.One advantage of the likelihood framework is that we work with distributions, including in our forecasts. Another advantage is that we can use a wide class of probability distributions to describe the data.    Generalized linear models Count forecastingTo forecast counts we are going to use a discrete probability distribution that takes values ${0,1,2,...}$, the __Poisson distribution__. A random variable following a Poisson distribution with parameter $\lambda_t>0$ takes values $y_t = 0,1,...$ with probability $Pr(y_t=\mu) = \frac{e^{-\lambda_t}\lambda_t^\mu}{\mu!}$. $\lambda$ is both the mean and the variance of the Poisson distribution.Suppose that we would like to model the parameter $\lambda_t$ of a Poisson distributed sample as a function of some explanatory variables $X_t$. We could posit $\lambda_t=X_t\beta$, with $X_t$ a column matrix of features and $\beta$ a column vector of parameters. The specification $\lambda_t=X_t\beta$ makes it possible for the mean and variance of the process to be negative, we can't have that! Instead we assume that $\lambda$ is a function of $X\beta$: $\lambda = f(X\beta)$, $f(.)$ is called the link function.A common choice for the link function is to assume that the parameter of the Poisson distribution is log-linear, $\log\lambda_t = X_t\beta$, implying that the link function is the exponential: $\lambda = f(X\beta)$. This choice can be problematic since we might end up taking the exponential of a large number leading to overflow. Instead we will use a logistic link function: $f(x) = \log(1 + exp(x))$. With link function $f$, we have $f(E(y_t|X_t))=X_t\beta$. The log-likelihood of some parameters $\beta$ at time $t$ is:$$L(\beta | X_t,y_t) =  y_t \log(\lambda_t(\beta)) - \lambda_t(\beta)\\ = y_t \log(f(X_t\beta)) - f(X_t\beta).$$In practice we will minimize $-L(\beta | X_t,y_t) = y_t \log(f(X_t\beta)) - f(X_t\beta)$. We define this log likelihood function below together with some helper functions.<jupyter_code>def logistic(x):
    return nd.log1p(nd.exp(x,dtype='float64'))

def inverse_logistic(x):
    return nd.log(nd.exp(x,dtype='float64')-1)

def _reshape_label_as_output(F, output, label):
    return label.reshape(output.shape) if F is ndarray else label.reshape(())


class Poisson(Loss):
    def __init__(self, weight=None, batch_axis=0, **kwargs):
        super(Poisson, self).__init__(weight, batch_axis, **kwargs)

    def hybrid_forward(self, F, output, label, sample_weight=None):
        label = _reshape_label_as_output(F, output, label)
        loss = logistic(output) - F.log(logistic(output)) * label
        return F.mean(loss, axis=self._batch_axis, exclude=True)<jupyter_output><empty_output><jupyter_text>Data<jupyter_code># Weekly sales for a novelty item (p.37-38: Montgomery) 
# From datamarket.com 
    
df = pd.read_csv("weekly-sales-clean.csv").set_index('week')
print(df.head())
ts = df.values[:,0]
plt.plot(ts);<jupyter_output><empty_output><jupyter_text>FeaturesThe data clearly appears to be trending upward, so the model will include a linear trend. There also might be some autocorrelation so the model will include a lag of the target variable. There does not appear to be any kind of seasonal pattern, so we won't include any seasonal features.<jupyter_code>forecast_length = 10
train_length = len(ts) - forecast_length

tgt_min = ts.min()
ts = ts - tgt_min

# droping the first observation due to lag
target = nd.array(ts[:train_length]).reshape((train_length,1))[1:]

# prediction target
pred_target = nd.array(ts[train_length:]).reshape((forecast_length,1))

# construct lag and trend
trend = nd.arange(train_length).reshape((train_length,1))
lag_sales = nd.array(ts).reshape((train_length,1))

# droping the last observation due to lag
features = nd.concat(trend[:-1],lag_sales[:-1])

# standardize
features_mean = features.mean(axis=0)
features_std = nd.array(features.asnumpy().std(axis=0)).reshape((1,1))
features = (features - features_mean) / features_std

print(features[:5,])
print(target[:5,])
batch_size = 5
train_data = gluon.data.DataLoader(
    gluon.data.ArrayDataset(features, target),
    batch_size=batch_size, shuffle=True)<jupyter_output><empty_output><jupyter_text>ModelThe GLM is still just a linear model, and so we define the network as a single dense layer with a dimension 1 output. The only difference with the Gaussian linear models of the previous notebook is that we now assume the data follows a Poisson distribution and so use the Poisson log likelihood as our loss function.<jupyter_code># Context
ctx = mx.cpu()

# Network
net = gluon.nn.Sequential()
with net.name_scope():
    net.add(gluon.nn.Dense(1))
    
# Parameter initialization    
net.collect_params().initialize(mx.init.Normal(sigma=0.1), ctx=ctx)

# Trainer
trainer = gluon.Trainer(net.collect_params(), 'sgd', {'learning_rate': 0.1})

# Loss
poisson_log_lik = Poisson()<jupyter_output><empty_output><jupyter_text>Training loop<jupyter_code>epochs = 50
smoothing_constant = .05
moving_loss = 0
niter = 0
loss_seq = []

for e in range(epochs):
    for i, (data, label) in enumerate(train_data):
        data = data.as_in_context(ctx)
        label = label.as_in_context(ctx)
        with autograd.record():
            output = net(data)
            loss = poisson_log_lik(output, label)
        loss.backward()
        trainer.step(batch_size)
        
        ##########################
        #  Keep a moving average of the losses
        ##########################
        niter +=1
        curr_loss = nd.mean(loss).asscalar()
        moving_loss = (1 - smoothing_constant) * moving_loss + (smoothing_constant) * curr_loss

        # correct the bias from the moving averages
        est_loss = moving_loss/(1-(1-smoothing_constant)**niter)
        loss_seq.append(est_loss)
    if e % 10 ==0:
        print("Epoch %s. Moving avg of log likelihood: %s" % (e, est_loss)) 
        
        params = net.collect_params() # this returns a ParameterDict

print('The type of "params" is a ',type(params))

# A ParameterDict is a dictionary of Parameter class objects
# therefore, here is how we can read off the parameters from it.

for param in params.values():
    print(param.name,param.data())
# plot the convergence of the estimated loss function
%matplotlib inline

import matplotlib
import matplotlib.pyplot as plt

plt.figure(num=None,figsize=(8, 6),dpi=80, facecolor='w', edgecolor='k')
plt.plot(range(niter),loss_seq, '.')

# adding some additional bells and whistles to the plot
plt.grid(True,which="both")
plt.xlabel('iteration',fontsize=14)
plt.ylabel('est loss',fontsize=14)<jupyter_output><empty_output><jupyter_text>ForecastAs was the case with the Gaussian AR(1) of the previous notebook, the specification of our Poisson model is autoregressive, a feature is the previous value of the target. We are therefore going forecast recursively.The value our model predict is $X_t\beta = f(\lambda_t)$ where $f(.)$ is the link function which in this application is the logistic function. To recover the parameter $\lambda$ of the Poisson distribution, which is what we are ultimately interested in estimating and predicting, we need to transform the output with $fˆ{-1}(.)$. That's what the `inverse_logistic` function is for.<jupyter_code>def forecast_poisson(net, last_obs, features_mean, features_std,
                     forecast_length, train_length):
    forecast = nd.empty((forecast_length,1))
    for t in range(forecast_length):
        if t==0:
            prev_obs = last_obs.reshape((1,1))
        # construct features
        trend = nd.array([train_length + t - 1]).reshape((1,1))
        #fct_feat = trend
        fct_feat = nd.concat(trend,prev_obs,dim=1)
        # normalize features
        fct_feat = (fct_feat - features_mean)/features_std
        # forecast
        fct = net(fct_feat)
        forecast[t,] = fct
        prev_obs = inverse_logistic(fct)
    return forecast
  
fit = net(features)    
fit_trans = inverse_logistic(fit)

fct = forecast_poisson(net,target[train_length-2],
                            features_mean, features_std,
                            forecast_length, train_length)
fct_trans = inverse_logistic(fct)
def plot_forecast(observed, fitted, forecasted):
    plt.plot(fitted.asnumpy(), color="r")
    plt.plot(observed, color="g")
    T = len(fitted)
    plt.plot(np.arange(T, T+len(forecasted)), forecasted.asnumpy(), color="b")
    plt.legend(["Fitted", "Observed", "Forecasted"]);
    
plot_forecast(ts[1:],fit,fct)<jupyter_output><empty_output><jupyter_text>90% prediction intervalsLet's put a 90% prediction interval around our forecast. We proceed as previously, estimating the standard errors of the residuals and relying on the fact that the maximum likelihood estimates are asymptotically Gaussian to construct the intervals. There are several alternative ways to construct intervals relying on weaker or different assumptions, but this is a topic for another day.<jupyter_code>se_fit = (target - fit_trans).asnumpy().std()
interval = np.concatenate((inverse_logistic(fct - 1.65*se_fit).asnumpy(),
                           inverse_logistic(fct + 1.65*se_fit).asnumpy()), axis=1)

def plot_forecast_interval(observed, fitted, forecasted, interval):
    plt.plot(fitted.asnumpy(), color="r")
    T = len(fitted)
    plt.plot(np.arange(T, T+len(forecasted)), forecasted.asnumpy(), color="b")
    plt.fill_between(np.arange(T, T+len(forecasted)), interval[:,0], interval[:,1], alpha=.3)
    plt.plot(observed, color="g")
    plt.legend(["Fitted", "Observed", "Forecasted"]);
    
plot_forecast_interval(ts[1:], fit, fct, interval)<jupyter_output><empty_output><jupyter_text>Mean squared prediction error<jupyter_code>msfe_glm = nd.mean(nd.power(fct - pred_target,2))
print(msfe_glm)<jupyter_output><empty_output><jupyter_text>Forecasting with a multilayer perceptronWe're going to use a multilayer perceptron (MLP) as in chapter 3. The network is composed of three layers, 2 dense layers with 64 hidden units and a _relu_ activation layer, and a one-dimensional output layer.<jupyter_code>num_hidden = 64
num_outputs = 1
mlp_net = gluon.nn.Sequential()
with mlp_net.name_scope():
    mlp_net.add(gluon.nn.Dense(num_hidden, activation="relu"))
    mlp_net.add(gluon.nn.Dense(num_hidden, activation="relu"))
    mlp_net.add(gluon.nn.Dense(num_outputs))<jupyter_output><empty_output><jupyter_text>InitializationHere we initialize the parameters using the Xavier algorithm, set up the trainer, and define the loss function. Since we are modeling and predicting continuous variables, I use a square Loss.<jupyter_code># context
mlp_ctx = mx.cpu()
# Parameters
mlp_net.collect_params().initialize(mx.init.Xavier(magnitude=2.24), ctx=mlp_ctx)
# trainer
trainer = gluon.Trainer(mlp_net.collect_params(), 'sgd', {'learning_rate': .1})
# likelihood
poisson_ll = Poisson()<jupyter_output><empty_output><jupyter_text>Data<jupyter_code>batch_size = 5
train_data = gluon.data.DataLoader(
    gluon.data.ArrayDataset(features, target),
    batch_size=batch_size, shuffle=True)<jupyter_output><empty_output><jupyter_text>Training loop<jupyter_code>epochs = 100
smoothing_constant = .05

for e in range(epochs):
    cumulative_loss = 0
    for i, (data, label) in enumerate(train_data):
        data = data.as_in_context(mlp_ctx)
        label = label.as_in_context(mlp_ctx)
        with autograd.record():
            output = mlp_net(data)
            #print(output)
            #loss = square_loss(output, label)
            loss = poisson_ll(output, label)
            loss.backward()
        trainer.step(data.shape[0])
        cumulative_loss += nd.sum(loss).asscalar()
    if e % 10 ==0:
        print("Epoch %s. Loss: %s" % (e, cumulative_loss))
    
for param in params.values():
    print(param.name,param.data())<jupyter_output><empty_output><jupyter_text>Forecasts and prediction intervals<jupyter_code># fit and residual standard errors
fit = mlp_net(features)
se_fit_mlp = (target - fit).asnumpy().std()

# forecast
fct = forecast_poisson(mlp_net,target[train_length-2],
                            features_mean, features_std,
                            forecast_length, train_length) 
# prediction interval
interval = np.concatenate((inverse_logistic(fct - 1.65*se_fit).asnumpy(),
                           inverse_logistic(fct + 1.65*se_fit).asnumpy()), axis=1)

plot_forecast_interval(ts[1:],inverse_logistic(fit),inverse_logistic(fct),interval)
msfe_mlp = nd.mean(nd.power(fct - pred_target,2))

print("Mean squared forecast error of the GLM: %s, of the MLP: %s" % (msfe_glm, msfe_mlp))<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup