| \begin{comment} |
| |
| Licensed to the Apache Software Foundation (ASF) under one |
| or more contributor license agreements. See the NOTICE file |
| distributed with this work for additional information |
| regarding copyright ownership. The ASF licenses this file |
| to you under the Apache License, Version 2.0 (the |
| "License"); you may not use this file except in compliance |
| with the License. You may obtain a copy of the License at |
| |
| http://www.apache.org/licenses/LICENSE-2.0 |
| |
| Unless required by applicable law or agreed to in writing, |
| software distributed under the License is distributed on an |
| "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY |
| KIND, either express or implied. See the License for the |
| specific language governing permissions and limitations |
| under the License. |
| |
| \end{comment} |
| |
| \subsection{Generalized Linear Models (GLM)} |
| \label{sec:GLM} |
| |
| \noindent{\bf Description} |
| \smallskip |
| |
| Generalized Linear Models~\cite{Gill2000:GLM,McCullagh1989:GLM,Nelder1972:GLM} |
| extend the methodology of linear and logistic regression to a variety of |
| distributions commonly assumed as noise effects in the response variable. |
| As before, we are given a collection |
| of records $(x_1, y_1)$, \ldots, $(x_n, y_n)$ where $x_i$ is a numerical vector of |
| explanatory (feature) variables of size~\mbox{$\dim x_i = m$}, and $y_i$ is the |
| response (dependent) variable observed for this vector. GLMs assume that some |
| linear combination of the features in~$x_i$ determines the \emph{mean}~$\mu_i$ |
| of~$y_i$, while the observed $y_i$ is a random outcome of a noise distribution |
| $\Prob[y\mid \mu_i]\,$\footnote{$\Prob[y\mid \mu_i]$ is given by a density function |
| if $y$ is continuous.} |
| with that mean~$\mu_i$: |
| \begin{equation*} |
| x_i \,\,\,\,\mapsto\,\,\,\, \eta_i = \beta_0 + \sum\nolimits_{j=1}^m \beta_j x_{i,j} |
| \,\,\,\,\mapsto\,\,\,\, \mu_i \,\,\,\,\mapsto \,\,\,\, y_i \sim \Prob[y\mid \mu_i] |
| \end{equation*} |
| |
| In linear regression the response mean $\mu_i$ \emph{equals} some linear combination |
| over~$x_i$, denoted above by~$\eta_i$. |
| In logistic regression with $y\in\{0, 1\}$ (Bernoulli) the mean of~$y$ is the same |
| as $\Prob[y=1]$ and equals $1/(1+e^{-\eta_i})$, the logistic function of~$\eta_i$. |
| In GLM, $\mu_i$ and $\eta_i$ can be related via any given smooth monotone function |
| called the \emph{link function}: $\eta_i = g(\mu_i)$. The unknown linear combination |
| parameters $\beta_j$ are assumed to be the same for all records. |
| |
| The goal of the regression is to estimate the parameters~$\beta_j$ from the observed |
| data. Once the~$\beta_j$'s are accurately estimated, we can make predictions |
| about~$y$ for a new feature vector~$x$. To do so, compute $\eta$ from~$x$ and use |
| the inverted link function $\mu = g^{-1}(\eta)$ to compute the mean $\mu$ of~$y$; |
| then use the distribution $\Prob[y\mid \mu]$ to make predictions about~$y$. |
| Both $g(\mu)$ and $\Prob[y\mid \mu]$ are user-provided. Our GLM script supports |
| a standard set of distributions and link functions, see below for details. |
| |
| \smallskip |
| \noindent{\bf Usage} |
| \smallskip |
| |
| {\hangindent=\parindent\noindent\it% |
| {\tt{}-f }path/\/{\tt{}GLM.dml} |
| {\tt{} -nvargs} |
| {\tt{} X=}path/file |
| {\tt{} Y=}path/file |
| {\tt{} B=}path/file |
| {\tt{} fmt=}format |
| {\tt{} O=}path/file |
| {\tt{} Log=}path/file |
| {\tt{} dfam=}int |
| {\tt{} vpow=}double |
| {\tt{} link=}int |
| {\tt{} lpow=}double |
| {\tt{} yneg=}double |
| {\tt{} icpt=}int |
| {\tt{} reg=}double |
| {\tt{} tol=}double |
| {\tt{} disp=}double |
| {\tt{} moi=}int |
| {\tt{} mii=}int |
| |
| } |
| |
| \smallskip |
| \noindent{\bf Arguments} |
| \begin{Description} |
| \item[{\tt X}:] |
| Location (on HDFS) to read the matrix of feature vectors; each row constitutes |
| an example. |
| \item[{\tt Y}:] |
| Location to read the response matrix, which may have 1 or 2 columns |
| \item[{\tt B}:] |
| Location to store the estimated regression parameters (the $\beta_j$'s), with the |
| intercept parameter~$\beta_0$ at position {\tt B[}$m\,{+}\,1$, {\tt 1]} if available |
| \item[{\tt fmt}:] (default:\mbox{ }{\tt "text"}) |
| Matrix file output format, such as {\tt text}, {\tt mm}, or {\tt csv}; |
| see read/write functions in SystemML Language Reference for details. |
| \item[{\tt O}:] (default:\mbox{ }{\tt " "}) |
| Location to write certain summary statistics described in Table~\ref{table:GLM:stats}, |
| by default it is standard output. |
| \item[{\tt Log}:] (default:\mbox{ }{\tt " "}) |
| Location to store iteration-specific variables for monitoring and debugging purposes, |
| see Table~\ref{table:GLM:log} for details. |
| \item[{\tt dfam}:] (default:\mbox{ }{\tt 1}) |
| Distribution family code to specify $\Prob[y\mid \mu]$, see Table~\ref{table:commonGLMs}:\\ |
| {\tt 1} = power distributions with $\Var(y) = \mu^{\alpha}$; |
| {\tt 2} = binomial or Bernoulli |
| \item[{\tt vpow}:] (default:\mbox{ }{\tt 0.0}) |
| When {\tt dfam=1}, this provides the~$q$ in $\Var(y) = a\mu^q$, the power |
| dependence of the variance of~$y$ on its mean. In particular, use:\\ |
| {\tt 0.0} = Gaussian, |
| {\tt 1.0} = Poisson, |
| {\tt 2.0} = Gamma, |
| {\tt 3.0} = inverse Gaussian |
| \item[{\tt link}:] (default:\mbox{ }{\tt 0}) |
| Link function code to determine the link function~$\eta = g(\mu)$:\\ |
| {\tt 0} = canonical link (depends on the distribution family), see Table~\ref{table:commonGLMs};\\ |
| {\tt 1} = power functions, |
| {\tt 2} = logit, |
| {\tt 3} = probit, |
| {\tt 4} = cloglog, |
| {\tt 5} = cauchit |
| \item[{\tt lpow}:] (default:\mbox{ }{\tt 1.0}) |
| When {\tt link=1}, this provides the~$s$ in $\eta = \mu^s$, the power link |
| function; {\tt lpow=0.0} gives the log link $\eta = \log\mu$. Common power links:\\ |
| {\tt -2.0} = $1/\mu^2$, |
| {\tt -1.0} = reciprocal, |
| {\tt 0.0} = log, |
| {\tt 0.5} = sqrt, |
| {\tt 1.0} = identity |
| \item[{\tt yneg}:] (default:\mbox{ }{\tt 0.0}) |
| When {\tt dfam=2} and the response matrix $Y$ has 1~column, |
| this specifies the $y$-value used for Bernoulli ``No'' label. |
| All other $y$-values are treated as the ``Yes'' label. |
| For example, {\tt yneg=-1.0} may be used when $y\in\{-1, 1\}$; |
| either {\tt yneg=1.0} or {\tt yneg=2.0} may be used when $y\in\{1, 2\}$. |
| \item[{\tt icpt}:] (default:\mbox{ }{\tt 0}) |
| Intercept and shifting/rescaling of the features in~$X$:\\ |
| {\tt 0} = no intercept (hence no~$\beta_0$), no shifting/rescaling of the features;\\ |
| {\tt 1} = add intercept, but do not shift/rescale the features in~$X$;\\ |
| {\tt 2} = add intercept, shift/rescale the features in~$X$ to mean~0, variance~1 |
| \item[{\tt reg}:] (default:\mbox{ }{\tt 0.0}) |
| L2-regularization parameter (lambda) |
| \item[{\tt tol}:] (default:\mbox{ }{\tt 0.000001}) |
| Tolerance (epsilon) used in the convergence criterion: we terminate the outer iterations |
| when the deviance changes by less than this factor; see below for details |
| \item[{\tt disp}:] (default:\mbox{ }{\tt 0.0}) |
| Dispersion parameter, or {\tt 0.0} to estimate it from data |
| \item[{\tt moi}:] (default:\mbox{ }{\tt 200}) |
| Maximum number of outer (Fisher scoring) iterations |
| \item[{\tt mii}:] (default:\mbox{ }{\tt 0}) |
| Maximum number of inner (conjugate gradient) iterations, or~0 if no maximum |
| limit provided |
| \end{Description} |
| |
| |
| \begin{table}[t]\small\centerline{% |
| \begin{tabular}{|ll|} |
| \hline |
| Name & Meaning \\ |
| \hline |
| {\tt TERMINATION\_CODE} & A positive integer indicating success/failure as follows: \\ |
| & $1 = {}$Converged successfully; |
| $2 = {}$Maximum \# of iterations reached; \\ |
| & $3 = {}$Input ({\tt X}, {\tt Y}) out of range; |
| $4 = {}$Distribution/link not supported \\ |
| {\tt BETA\_MIN} & Smallest beta value (regression coefficient), excluding the intercept \\ |
| {\tt BETA\_MIN\_INDEX} & Column index for the smallest beta value \\ |
| {\tt BETA\_MAX} & Largest beta value (regression coefficient), excluding the intercept \\ |
| {\tt BETA\_MAX\_INDEX} & Column index for the largest beta value \\ |
| {\tt INTERCEPT} & Intercept value, or NaN if there is no intercept (if {\tt icpt=0}) \\ |
| {\tt DISPERSION} & Dispersion used to scale deviance, provided in {\tt disp} input argument \\ |
| & or estimated (same as {\tt DISPERSION\_EST}) if {\tt disp} argument is${} \leq 0$ \\ |
| {\tt DISPERSION\_EST} & Dispersion estimated from the dataset \\ |
| {\tt DEVIANCE\_UNSCALED} & Deviance from the saturated model, assuming dispersion${} = 1.0$ \\ |
| {\tt DEVIANCE\_SCALED} & Deviance from the saturated model, scaled by {\tt DISPERSION} value \\ |
| \hline |
| \end{tabular}} |
| \caption{Besides~$\beta$, GLM regression script computes a few summary statistics listed above. |
| They are provided in CSV format, one comma-separated name-value pair per each line.} |
| \label{table:GLM:stats} |
| \end{table} |
| |
| |
| |
| |
| |
| |
| \begin{table}[t]\small\centerline{% |
| \begin{tabular}{|ll|} |
| \hline |
| Name & Meaning \\ |
| \hline |
| {\tt NUM\_CG\_ITERS} & Number of inner (Conj.\ Gradient) iterations in this outer iteration \\ |
| {\tt IS\_TRUST\_REACHED} & $1 = {}$trust region boundary was reached, $0 = {}$otherwise \\ |
| {\tt POINT\_STEP\_NORM} & L2-norm of iteration step from old point ($\beta$-vector) to new point \\ |
| {\tt OBJECTIVE} & The loss function we minimize (negative partial log-likelihood) \\ |
| {\tt OBJ\_DROP\_REAL} & Reduction in the objective during this iteration, actual value \\ |
| {\tt OBJ\_DROP\_PRED} & Reduction in the objective predicted by a quadratic approximation \\ |
| {\tt OBJ\_DROP\_RATIO} & Actual-to-predicted reduction ratio, used to update the trust region \\ |
| {\tt GRADIENT\_NORM} & L2-norm of the loss function gradient (omitted if point is rejected) \\ |
| {\tt LINEAR\_TERM\_MIN} & The minimum value of $X \pxp \beta$, used to check for overflows \\ |
| {\tt LINEAR\_TERM\_MAX} & The maximum value of $X \pxp \beta$, used to check for overflows \\ |
| {\tt IS\_POINT\_UPDATED} & $1 = {}$new point accepted; $0 = {}$new point rejected, old point restored \\ |
| {\tt TRUST\_DELTA} & Updated trust region size, the ``delta'' \\ |
| \hline |
| \end{tabular}} |
| \caption{ |
| The {\tt Log} file for GLM regression contains the above \mbox{per-}iteration |
| variables in CSV format, each line containing triple (Name, Iteration\#, Value) with Iteration\# |
| being~0 for initial values.} |
| \label{table:GLM:log} |
| \end{table} |
| |
| \begin{table}[t]\hfil |
| \begin{tabular}{|ccccccc|} |
| \hline |
| \multicolumn{4}{|c}{INPUT PARAMETERS} & Distribution & Link & Cano- \\ |
| {\tt dfam} & {\tt vpow} & {\tt link} & {\tt\ lpow} & family & function & nical?\\ |
| \hline |
| {\tt 1} & {\tt 0.0} & {\tt 1} & {\tt -1.0} & Gaussian & inverse & \\ |
| {\tt 1} & {\tt 0.0} & {\tt 1} & {\tt\ 0.0} & Gaussian & log & \\ |
| {\tt 1} & {\tt 0.0} & {\tt 1} & {\tt\ 1.0} & Gaussian & identity & Yes \\ |
| {\tt 1} & {\tt 1.0} & {\tt 1} & {\tt\ 0.0} & Poisson & log & Yes \\ |
| {\tt 1} & {\tt 1.0} & {\tt 1} & {\tt\ 0.5} & Poisson & sq.root & \\ |
| {\tt 1} & {\tt 1.0} & {\tt 1} & {\tt\ 1.0} & Poisson & identity & \\ |
| {\tt 1} & {\tt 2.0} & {\tt 1} & {\tt -1.0} & Gamma & inverse & Yes \\ |
| {\tt 1} & {\tt 2.0} & {\tt 1} & {\tt\ 0.0} & Gamma & log & \\ |
| {\tt 1} & {\tt 2.0} & {\tt 1} & {\tt\ 1.0} & Gamma & identity & \\ |
| {\tt 1} & {\tt 3.0} & {\tt 1} & {\tt -2.0} & Inverse Gauss & $1/\mu^2$ & Yes \\ |
| {\tt 1} & {\tt 3.0} & {\tt 1} & {\tt -1.0} & Inverse Gauss & inverse & \\ |
| {\tt 1} & {\tt 3.0} & {\tt 1} & {\tt\ 0.0} & Inverse Gauss & log & \\ |
| {\tt 1} & {\tt 3.0} & {\tt 1} & {\tt\ 1.0} & Inverse Gauss & identity & \\ |
| \hline |
| {\tt 2} & {\tt *} & {\tt 1} & {\tt\ 0.0} & Binomial & log & \\ |
| {\tt 2} & {\tt *} & {\tt 1} & {\tt\ 0.5} & Binomial & sq.root & \\ |
| {\tt 2} & {\tt *} & {\tt 2} & {\tt\ *} & Binomial & logit & Yes \\ |
| {\tt 2} & {\tt *} & {\tt 3} & {\tt\ *} & Binomial & probit & \\ |
| {\tt 2} & {\tt *} & {\tt 4} & {\tt\ *} & Binomial & cloglog & \\ |
| {\tt 2} & {\tt *} & {\tt 5} & {\tt\ *} & Binomial & cauchit & \\ |
| \hline |
| \end{tabular}\hfil |
| \caption{Common GLM distribution families and link functions. |
| (Here ``{\tt *}'' stands for ``any value.'')} |
| \label{table:commonGLMs} |
| \end{table} |
| |
| \noindent{\bf Details} |
| \smallskip |
| |
| In GLM, the noise distribution $\Prob[y\mid \mu]$ of the response variable~$y$ |
| given its mean~$\mu$ is restricted to have the \emph{exponential family} form |
| \begin{equation} |
| Y \sim\, \Prob[y\mid \mu] \,=\, \exp\left(\frac{y\theta - b(\theta)}{a} |
| + c(y, a)\right),\,\,\textrm{where}\,\,\,\mu = \E(Y) = b'(\theta). |
| \label{eqn:GLM} |
| \end{equation} |
| Changing the mean in such a distribution simply multiplies all \mbox{$\Prob[y\mid \mu]$} |
| by~$e^{\,y\hspace{0.2pt}\theta/a}$ and rescales them so that they again integrate to~1. |
| Parameter $\theta$ is called \emph{canonical}, and the function $\theta = b'^{\,-1}(\mu)$ |
| that relates it to the mean is called the~\emph{canonical link}; constant~$a$ is called |
| \emph{dispersion} and rescales the variance of~$y$. Many common distributions can be put |
| into this form, see Table~\ref{table:commonGLMs}. The canonical parameter~$\theta$ |
| is often chosen to coincide with~$\eta$, the linear combination of the regression features; |
| other choices for~$\eta$ are possible too. |
| |
| Rather than specifying the canonical link, GLM distributions are commonly defined |
| by their variance $\Var(y)$ as the function of the mean~$\mu$. It can be shown |
| from Eq.~(\ref{eqn:GLM}) that $\Var(y) = a\,b''(\theta) = a\,b''(b'^{\,-1}(\mu))$. |
| For example, for the Bernoulli distribution $\Var(y) = \mu(1-\mu)$, for the Poisson |
| distribution \mbox{$\Var(y) = \mu$}, and for the Gaussian distribution |
| $\Var(y) = a\cdot 1 = \sigma^2$. |
| It turns out that for many common distributions $\Var(y) = a\mu^q$, a power function. |
| We support all distributions where $\Var(y) = a\mu^q$, as well as the Bernoulli and |
| the binomial distributions. |
| |
| For distributions with $\Var(y) = a\mu^q$ the canonical link is also a power function, |
| namely $\theta = \mu^{1-q}/(1-q)$, except for the Poisson ($q = 1$) whose canonical link is |
| $\theta = \log\mu$. We support all power link functions in the form $\eta = \mu^s$, |
| dropping any constant factor, with $\eta = \log\mu$ for $s=0$. The binomial distribution |
| has its own family of link functions, which includes logit (the canonical link), |
| probit, cloglog, and cauchit (see Table~\ref{table:binomial_links}); we support these |
| only for the binomial and Bernoulli distributions. Links and distributions are specified |
| via four input parameters: {\tt dfam}, {\tt vpow}, {\tt link}, and {\tt lpow} (see |
| Table~\ref{table:commonGLMs}). |
| |
| \begin{table}[t]\hfil |
| \begin{tabular}{|cc|cc|} |
| \hline |
| Name & Link function & Name & Link function \\ |
| \hline |
| Logit & $\displaystyle \eta = 1 / \big(1 + e^{-\mu}\big)^{\mathstrut}$ & |
| Cloglog & $\displaystyle \eta = \log \big(\!- \log(1 - \mu)\big)^{\mathstrut}$ \\ |
| Probit & $\displaystyle \mu = \frac{1}{\sqrt{2\pi}}\int\nolimits_{-\infty_{\mathstrut}}^{\,\eta\mathstrut} |
| \!\!\!\!\! e^{-\frac{t^2}{2}} dt$ & |
| Cauchit & $\displaystyle \eta = \tan\pi(\mu - 1/2)$ \\ |
| \hline |
| \end{tabular}\hfil |
| \caption{The supported non-power link functions for the Bernoulli and the binomial |
| distributions. (Here $\mu$~is the Bernoulli mean.)} |
| \label{table:binomial_links} |
| \end{table} |
| |
| The observed response values are provided to the regression script as matrix~$Y$ |
| having 1 or 2 columns. If a power distribution family is selected ({\tt dfam=1}), |
| matrix $Y$ must have 1~column that provides $y_i$ for each~$x_i$ in the corresponding |
| row of matrix~$X$. When {\tt dfam=2} and $Y$ has 1~column, we assume the Bernoulli |
| distribution for $y_i\in\{y_{\mathrm{neg}}, y_{\mathrm{pos}}\}$ with $y_{\mathrm{neg}}$ |
| from the input parameter {\tt yneg} and with $y_{\mathrm{pos}} \neq y_{\mathrm{neg}}$. |
| When {\tt dfam=2} and $Y$ has 2~columns, we assume the |
| binomial distribution; for each row~$i$ in~$X$, cells $Y[i, 1]$ and $Y[i, 2]$ provide |
| the positive and the negative binomial counts respectively. Internally we convert |
| the 1-column Bernoulli into the 2-column binomial with 0-versus-1 counts. |
| |
| We estimate the regression parameters via L2-regularized negative log-likelihood |
| minimization: |
| \begin{equation*} |
| f(\beta; X, Y) \,\,=\,\, -\sum\nolimits_{i=1}^n \big(y_i\theta_i - b(\theta_i)\big) |
| \,+\,(\lambda/2) \sum\nolimits_{j=1}^m \beta_j^2\,\,\to\,\,\min |
| \end{equation*} |
| where $\theta_i$ and $b(\theta_i)$ are from~(\ref{eqn:GLM}); note that $a$ |
| and $c(y, a)$ are constant w.r.t.~$\beta$ and can be ignored here. |
| The canonical parameter $\theta_i$ depends on both $\beta$ and~$x_i$: |
| \begin{equation*} |
| \theta_i \,\,=\,\, b'^{\,-1}(\mu_i) \,\,=\,\, b'^{\,-1}\big(g^{-1}(\eta_i)\big) \,\,=\,\, |
| \big(b'^{\,-1}\circ g^{-1}\big)\left(\beta_0 + \sum\nolimits_{j=1}^m \beta_j x_{i,j}\right) |
| \end{equation*} |
| The user-provided (via {\tt reg}) regularization coefficient $\lambda\geq 0$ can be used |
| to mitigate overfitting and degeneracy in the data. Note that the intercept is never |
| regularized. |
| |
| Our iterative minimizer for $f(\beta; X, Y)$ uses the Fisher scoring approximation |
| to the difference $\varDelta f(z; \beta) = f(\beta + z; X, Y) \,-\, f(\beta; X, Y)$, |
| recomputed at each iteration: |
| \begin{gather*} |
| \varDelta f(z; \beta) \,\,\,\approx\,\,\, 1/2 \cdot z^T A z \,+\, G^T z, |
| \,\,\,\,\textrm{where}\,\,\,\, A \,=\, X^T\!\diag(w) X \,+\, \lambda I\\ |
| \textrm{and}\,\,\,\,G \,=\, - X^T u \,+\, \lambda\beta, |
| \,\,\,\textrm{with $n\,{\times}\,1$ vectors $w$ and $u$ given by}\\ |
| \forall\,i = 1\ldots n: \,\,\,\, |
| w_i = \big[v(\mu_i)\,g'(\mu_i)^2\big]^{-1} |
| \!\!\!\!\!\!,\,\,\,\,\,\,\,\,\, |
| u_i = (y_i - \mu_i)\big[v(\mu_i)\,g'(\mu_i)\big]^{-1} |
| \!\!\!\!\!\!.\,\,\,\, |
| \end{gather*} |
| Here $v(\mu_i)=\Var(y_i)/a$, the variance of $y_i$ as the function of the mean, and |
| $g'(\mu_i) = d \eta_i/d \mu_i$ is the link function derivative. The Fisher scoring |
| approximation is minimized by trust-region conjugate gradient iterations (called the |
| \emph{inner} iterations, with the Fisher scoring iterations as the \emph{outer} |
| iterations), which approximately solve the following problem: |
| \begin{equation*} |
| 1/2 \cdot z^T A z \,+\, G^T z \,\,\to\,\,\min\,\,\,\,\textrm{subject to}\,\,\,\, |
| \|z\|_2 \leq \delta |
| \end{equation*} |
| The conjugate gradient algorithm closely follows Algorithm~7.2 on page~171 |
| of~\cite{Nocedal2006:Optimization}. |
| The trust region size $\delta$ is initialized as $0.5\sqrt{m}\,/ \max\nolimits_i \|x_i\|_2$ |
| and updated as described in~\cite{Nocedal2006:Optimization}. |
| The user can specify the maximum number of the outer and the inner iterations with |
| input parameters {\tt moi} and {\tt mii}, respectively. The Fisher scoring algorithm |
| terminates successfully if $2|\varDelta f(z; \beta)| < (D_1(\beta) + 0.1)\hspace{0.5pt}\eps$ |
| where $\eps > 0$ is a tolerance supplied by the user via {\tt tol}, and $D_1(\beta)$ is |
| the unit-dispersion deviance estimated as |
| \begin{equation*} |
| D_1(\beta) \,\,=\,\, 2 \cdot \big(\Prob[Y \mid \! |
| \begin{smallmatrix}\textrm{saturated}\\\textrm{model}\end{smallmatrix}, a\,{=}\,1] |
| \,\,-\,\,\Prob[Y \mid X, \beta, a\,{=}\,1]\,\big) |
| \end{equation*} |
| The deviance estimate is also produced as part of the output. Once the Fisher scoring |
| algorithm terminates, if requested by the user, we estimate the dispersion~$a$ from |
| Eq.~\ref{eqn:GLM} using Pearson residuals |
| \begin{equation} |
| \hat{a} \,\,=\,\, \frac{1}{n-m}\cdot \sum_{i=1}^n \frac{(y_i - \mu_i)^2}{v(\mu_i)} |
| \label{eqn:dispersion} |
| \end{equation} |
| and use it to adjust our deviance estimate: $D_{\hat{a}}(\beta) = D_1(\beta)/\hat{a}$. |
| If input argument {\tt disp} is {\tt 0.0} we estimate $\hat{a}$, otherwise we use its |
| value as~$a$. Note that in~(\ref{eqn:dispersion}) $m$~counts the intercept |
| ($m \leftarrow m+1$) if it is present. |
| |
| \smallskip |
| \noindent{\bf Returns} |
| \smallskip |
| |
| The estimated regression parameters (the $\hat{\beta}_j$'s) are populated into |
| a matrix and written to an HDFS file whose path/name was provided as the ``{\tt B}'' |
| input argument. What this matrix contains, and its size, depends on the input |
| argument {\tt icpt}, which specifies the user's intercept and rescaling choice: |
| \begin{Description} |
| \item[{\tt icpt=0}:] No intercept, matrix~$B$ has size $m\,{\times}\,1$, with |
| $B[j, 1] = \hat{\beta}_j$ for each $j$ from 1 to~$m = {}$ncol$(X)$. |
| \item[{\tt icpt=1}:] There is intercept, but no shifting/rescaling of~$X$; matrix~$B$ |
| has size $(m\,{+}\,1) \times 1$, with $B[j, 1] = \hat{\beta}_j$ for $j$ from 1 to~$m$, |
| and $B[m\,{+}\,1, 1] = \hat{\beta}_0$, the estimated intercept coefficient. |
| \item[{\tt icpt=2}:] There is intercept, and the features in~$X$ are shifted to |
| mean${} = 0$ and rescaled to variance${} = 1$; then there are two versions of |
| the~$\hat{\beta}_j$'s, one for the original features and another for the |
| shifted/rescaled features. Now matrix~$B$ has size $(m\,{+}\,1) \times 2$, with |
| $B[\cdot, 1]$ for the original features and $B[\cdot, 2]$ for the shifted/rescaled |
| features, in the above format. Note that $B[\cdot, 2]$ are iteratively estimated |
| and $B[\cdot, 1]$ are obtained from $B[\cdot, 2]$ by complementary shifting and |
| rescaling. |
| \end{Description} |
| Our script also estimates the dispersion $\hat{a}$ (or takes it from the user's input) |
| and the deviances $D_1(\hat{\beta})$ and $D_{\hat{a}}(\hat{\beta})$, see |
| Table~\ref{table:GLM:stats} for details. A log file with variables monitoring |
| progress through the iterations can also be made available, see Table~\ref{table:GLM:log}. |
| |
| \smallskip |
| \noindent{\bf Examples} |
| \smallskip |
| |
| {\hangindent=\parindent\noindent\tt |
| \hml -f GLM.dml -nvargs X=/user/biadmin/X.mtx Y=/user/biadmin/Y.mtx |
| B=/user/biadmin/B.mtx fmt=csv dfam=2 link=2 yneg=-1.0 icpt=2 reg=0.01 tol=0.00000001 |
| disp=1.0 moi=100 mii=10 O=/user/biadmin/stats.csv Log=/user/biadmin/log.csv |
| |
| } |
| |
| \smallskip |
| \noindent{\bf See Also} |
| \smallskip |
| |
| In case of binary classification problems, consider using L2-SVM or binary logistic |
| regression; for multiclass classification, use multiclass~SVM or multinomial logistic |
| regression. For the special cases of linear regression and logistic regression, it |
| may be more efficient to use the corresponding specialized scripts instead of~GLM. |
