Creates the structure of the model (model matrices, etc) for
        a mini-batch Normal family ARIMAX model.

        Here the structure is the same as for _normal_model() but we are going to
        sample a random choice of data points (of length mini_batch).

        Parameters
        ----------
        beta : np.ndarray
            Contains untransformed starting values for the latent variables

        mini_batch : int
            Mini batch size for the data sampling

        Returns
        ----------
        mu : np.ndarray
            Contains the predicted values (location) for the time series

        Y : np.ndarray
            Contains the length-adjusted time series (accounting for lags)

def _mb_normal_model(self, beta, mini_batch):
        """ Creates the structure of the model (model matrices, etc) for
        a mini-batch Normal family ARIMAX model.

        Here the structure is the same as for _normal_model() but we are going to
        sample a random choice of data points (of length mini_batch).

        Parameters
        ----------
        beta : np.ndarray
            Contains untransformed starting values for the latent variables

        mini_batch : int
            Mini batch size for the data sampling

        Returns
        ----------
        mu : np.ndarray
            Contains the predicted values (location) for the time series

        Y : np.ndarray
            Contains the length-adjusted time series (accounting for lags)
        """     

        rand_int =  np.random.randint(low=0, high=self.data_length-mini_batch-self.max_lag+1)
        sample = np.arange(start=rand_int, stop=rand_int+mini_batch)

        data = self.y[sample]
        X = self.X[sample, :]
        Y = data[self.max_lag:]

        if self.ar != 0:
            ar_matrix = data[(self.max_lag-1):-1]
            for i in range(1, self.ar):
                ar_matrix = np.vstack((ar_matrix, data[(self.max_lag-i-1):-i-1]))
        else:
            ar_matrix = np.zeros(data.shape[0]-self.max_lag)

        # Transform latent variables
        z = np.array([self.latent_variables.z_list[k].prior.transform(beta[k]) for k in range(beta.shape[0])])

        # Constant and AR terms
        if self.ar == 0:
            mu = np.transpose(ar_matrix)
        elif self.ar == 1:
            mu = np.transpose(ar_matrix)*z[:-self.family_z_no-self.ma-len(self.X_names)][0]
        else:
            mu = np.matmul(np.transpose(ar_matrix),z[:-self.family_z_no-self.ma-len(self.X_names)])

        # X terms
        mu = mu + np.matmul(X[self.integ+self.max_lag:],z[self.ma+self.ar:(self.ma+self.ar+len(self.X_names))])

        # MA terms
        if self.ma != 0:
            mu = arimax_recursion(z, mu, Y, self.max_lag, Y.shape[0], self.ar, self.ma)

        return mu, Y