You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@systemml.apache.org by du...@apache.org on 2015/12/02 02:04:41 UTC
[03/47] incubator-systemml git commit: Add Algorithms and Language
References
http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/a891c0a1/Algorithms Reference/DescriptiveUnivarStats.tex
----------------------------------------------------------------------
diff --git a/Algorithms Reference/DescriptiveUnivarStats.tex b/Algorithms Reference/DescriptiveUnivarStats.tex
new file mode 100644
index 0000000..fcea4ec
--- /dev/null
+++ b/Algorithms Reference/DescriptiveUnivarStats.tex
@@ -0,0 +1,582 @@
+\subsection{Univariate Statistics}
+
+\noindent{\bf Description}
+\smallskip
+
+\emph{Univariate statistics} are the simplest form of descriptive statistics in data
+analysis. They are used to quantitatively describe the main characteristics of each
+feature in the data. For a given dataset matrix, script \UnivarScriptName{} computes
+certain univariate statistics for each feature column in the
+matrix. The feature type governs the exact set of statistics computed for that feature.
+For example, the statistic \emph{mean} can only be computed on a quantitative (scale)
+feature like `Height' and `Temperature'. It does not make sense to compute the mean
+of a categorical attribute like `Hair Color'.
+
+
+\smallskip
+\noindent{\bf Usage}
+\smallskip
+
+{\hangindent=\parindent\noindent\it%\tolerance=0
+{\tt{}-f } \UnivarScriptName{}
+{\tt{} -nvargs}
+{\tt{} X=}path/file
+{\tt{} TYPES=}path/file
+{\tt{} STATS=}path/file
+% {\tt{} fmt=}format
+
+}
+
+
+\medskip
+\pagebreak[2]
+\noindent{\bf Arguments}
+\begin{Description}
+\item[{\tt X}:]
+Location (on HDFS) to read the data matrix $X$ whose columns we want to
+analyze as the features.
+\item[{\tt TYPES}:] % (default:\mbox{ }{\tt " "})
+Location (on HDFS) to read the single-row matrix whose $i^{\textrm{th}}$
+column-cell contains the type of the $i^{\textrm{th}}$ feature column
+\texttt{X[,$\,i$]} in the data matrix. Feature types must be encoded by
+integer numbers: $1 = {}$scale, $2 = {}$nominal, $3 = {}$ordinal.
+% The default value means ``treat all $X$-columns as scale.''
+\item[{\tt STATS}:]
+Location (on HDFS) where the output matrix of computed statistics
+will be stored. The format of the output matrix is defined by
+Table~\ref{table:univars}.
+% \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.
+\end{Description}
+
+\begin{table}[t]\hfil
+\begin{tabular}{|rl|c|c|}
+\hline
+\multirow{2}{*}{Row}& \multirow{2}{*}{Name of Statistic} & \multicolumn{2}{c|}{Applies to:} \\
+ & & Scale & Categ.\\
+\hline
+\OutputRowIDMinimum & Minimum & + & \\
+\OutputRowIDMaximum & Maximum & + & \\
+\OutputRowIDRange & Range & + & \\
+\OutputRowIDMean & Mean & + & \\
+\OutputRowIDVariance & Variance & + & \\
+\OutputRowIDStDeviation & Standard deviation & + & \\
+\OutputRowIDStErrorMean & Standard error of mean & + & \\
+\OutputRowIDCoeffVar & Coefficient of variation & + & \\
+\OutputRowIDSkewness & Skewness & + & \\
+\OutputRowIDKurtosis & Kurtosis & + & \\
+\OutputRowIDStErrorSkewness & Standard error of skewness & + & \\
+\OutputRowIDStErrorCurtosis & Standard error of kurtosis & + & \\
+\OutputRowIDMedian & Median & + & \\
+\OutputRowIDIQMean & Inter quartile mean & + & \\
+\OutputRowIDNumCategories & Number of categories & & + \\
+\OutputRowIDMode & Mode & & + \\
+\OutputRowIDNumModes & Number of modes & & + \\
+\hline
+\end{tabular}\hfil
+\caption{The output matrix of \UnivarScriptName{} has one row per each
+univariate statistic and one column per input feature. This table lists
+the meaning of each row. Signs ``+'' show applicability to scale or/and
+to categorical features.}
+\label{table:univars}
+\end{table}
+
+
+\pagebreak[1]
+
+\smallskip
+\noindent{\bf Details}
+\smallskip
+
+Given an input matrix \texttt{X}, this script computes the set of all
+relevant univariate statistics for each feature column \texttt{X[,$\,i$]}
+in~\texttt{X}. The list of statistics to be computed depends on the
+\emph{type}, or \emph{measurement level}, of each column.
+The \textrm{TYPES} command-line argument points to a vector containing
+the types of all columns. The types must be provided as per the following
+convention: $1 = {}$scale, $2 = {}$nominal, $3 = {}$ordinal.
+
+Below we list all univariate statistics computed by script \UnivarScriptName.
+The statistics are collected by relevance into several groups, namely: central
+tendency, dispersion, shape, and categorical measures. The first three groups
+contain statistics computed for a quantitative (also known as: numerical, scale,
+or continuous) feature; the last group contains the statistics for a categorical
+(either nominal or ordinal) feature.
+
+Let~$n$ be the number of data records (rows) with feature values.
+In what follows we fix a column index \texttt{idx} and consider
+sample statistics of feature column \texttt{X[}$\,$\texttt{,}$\,$\texttt{idx]}.
+Let $v = (v_1, v_2, \ldots, v_n)$ be the values of \texttt{X[}$\,$\texttt{,}$\,$\texttt{idx]}
+in their original unsorted order: $v_i = \texttt{X[}i\texttt{,}\,\texttt{idx]}$.
+Let $v^s = (v^s_1, v^s_2, \ldots, v^s_n)$ be the same values in the sorted order,
+preserving duplicates: $v^s_1 \leq v^s_2 \leq \ldots \leq v^s_n$.
+
+\paragraph{Central tendency measures.}
+Sample statistics that describe the location of the quantitative (scale) feature distribution,
+represent it with a single value.
+\begin{Description}
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Mean]
+\OutputRowText{\OutputRowIDMean}
+The arithmetic average over a sample of a quantitative feature.
+Computed as the ratio between the sum of values and the number of values:
+$\left(\sum_{i=1}^n v_i\right)\!/n$.
+Example: the mean of sample $\{$2.2, 3.2, 3.7, 4.4, 5.3, 5.7, 6.1, 6.4, 7.2, 7.8$\}$
+equals~5.2.
+
+Note that the mean is significantly affected by extreme values in the sample
+and may be misleading as a central tendency measure if the feature varies on
+exponential scale. For example, the mean of $\{$0.01, 0.1, 1.0, 10.0, 100.0$\}$
+is 22.222, greater than all the sample values except the~largest.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\begin{figure}[t]
+\setlength{\unitlength}{10pt}
+\begin{picture}(33,12)
+\put( 6.2, 0.0){\small 2.2}
+\put(10.2, 0.0){\small 3.2}
+\put(12.2, 0.0){\small 3.7}
+\put(15.0, 0.0){\small 4.4}
+\put(18.6, 0.0){\small 5.3}
+\put(20.2, 0.0){\small 5.7}
+\put(21.75,0.0){\small 6.1}
+\put(23.05,0.0){\small 6.4}
+\put(26.2, 0.0){\small 7.2}
+\put(28.6, 0.0){\small 7.8}
+\put( 0.5, 0.7){\small 0.0}
+\put( 0.1, 3.2){\small 0.25}
+\put( 0.5, 5.7){\small 0.5}
+\put( 0.1, 8.2){\small 0.75}
+\put( 0.5,10.7){\small 1.0}
+\linethickness{1.5pt}
+\put( 2.0, 1.0){\line(1,0){4.8}}
+\put( 6.8, 1.0){\line(0,1){1.0}}
+\put( 6.8, 2.0){\line(1,0){4.0}}
+\put(10.8, 2.0){\line(0,1){1.0}}
+\put(10.8, 3.0){\line(1,0){2.0}}
+\put(12.8, 3.0){\line(0,1){1.0}}
+\put(12.8, 4.0){\line(1,0){2.8}}
+\put(15.6, 4.0){\line(0,1){1.0}}
+\put(15.6, 5.0){\line(1,0){3.6}}
+\put(19.2, 5.0){\line(0,1){1.0}}
+\put(19.2, 6.0){\line(1,0){1.6}}
+\put(20.8, 6.0){\line(0,1){1.0}}
+\put(20.8, 7.0){\line(1,0){1.6}}
+\put(22.4, 7.0){\line(0,1){1.0}}
+\put(22.4, 8.0){\line(1,0){1.2}}
+\put(23.6, 8.0){\line(0,1){1.0}}
+\put(23.6, 9.0){\line(1,0){3.2}}
+\put(26.8, 9.0){\line(0,1){1.0}}
+\put(26.8,10.0){\line(1,0){2.4}}
+\put(29.2,10.0){\line(0,1){1.0}}
+\put(29.2,11.0){\line(1,0){4.8}}
+\linethickness{0.3pt}
+\put( 6.8, 1.0){\circle*{0.3}}
+\put(10.8, 1.0){\circle*{0.3}}
+\put(12.8, 1.0){\circle*{0.3}}
+\put(15.6, 1.0){\circle*{0.3}}
+\put(19.2, 1.0){\circle*{0.3}}
+\put(20.8, 1.0){\circle*{0.3}}
+\put(22.4, 1.0){\circle*{0.3}}
+\put(23.6, 1.0){\circle*{0.3}}
+\put(26.8, 1.0){\circle*{0.3}}
+\put(29.2, 1.0){\circle*{0.3}}
+\put( 6.8, 1.0){\vector(1,0){27.2}}
+\put( 2.0, 1.0){\vector(0,1){10.8}}
+\put( 2.0, 3.5){\line(1,0){10.8}}
+\put( 2.0, 6.0){\line(1,0){17.2}}
+\put( 2.0, 8.5){\line(1,0){21.6}}
+\put( 2.0,11.0){\line(1,0){27.2}}
+\put(12.8, 1.0){\line(0,1){2.0}}
+\put(19.2, 1.0){\line(0,1){5.0}}
+\put(20.0, 1.0){\line(0,1){5.0}}
+\put(23.6, 1.0){\line(0,1){7.0}}
+\put( 9.0, 4.0){\line(1,0){3.8}}
+\put( 9.2, 2.7){\vector(0,1){0.8}}
+\put( 9.2, 4.8){\vector(0,-1){0.8}}
+\put(19.4, 8.0){\line(1,0){3.0}}
+\put(19.6, 7.2){\vector(0,1){0.8}}
+\put(19.6, 9.3){\vector(0,-1){0.8}}
+\put(13.0, 2.2){\small $q_{25\%}$}
+\put(17.3, 2.2){\small $q_{50\%}$}
+\put(23.8, 2.2){\small $q_{75\%}$}
+\put(20.15,3.5){\small $\mu$}
+\put( 8.0, 3.75){\small $\phi_1$}
+\put(18.35,7.8){\small $\phi_2$}
+\end{picture}
+\label{fig:example_quartiles}
+\caption{The computation of quartiles, median, and interquartile mean from the
+empirical distribution function of the 10-point
+sample $\{$2.2, 3.2, 3.7, 4.4, 5.3, 5.7, 6.1, 6.4, 7.2, 7.8$\}$. Each vertical step in
+the graph has height~$1{/}n = 0.1$. Values $q_{25\%}$, $q_{50\%}$, and $q_{75\%}$ denote
+the $1^{\textrm{st}}$, $2^{\textrm{nd}}$, and $3^{\textrm{rd}}$ quartiles correspondingly;
+value~$\mu$ denotes the median. Values $\phi_1$ and $\phi_2$ show the partial contribution
+of border points (quartiles) $v_3=3.7$ and $v_8=6.4$ into the interquartile mean.}
+\end{figure}
+
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Median]
+\OutputRowText{\OutputRowIDMedian}
+The ``middle'' value that separates the higher half of the sample values
+(in a sorted order) from the lower half.
+To compute the median, we sort the sample in the increasing order, preserving
+duplicates: $v^s_1 \leq v^s_2 \leq \ldots \leq v^s_n$.
+If $n$ is odd, the median equals $v^s_i$ where $i = (n\,{+}\,1)\,{/}\,2$,
+same as the $50^{\textrm{th}}$~percentile of the sample.
+If $n$ is even, there are two ``middle'' values $v^s_{n/2}$ and $v^s_{n/2\,+\,1}$,
+so we compute the median as the mean of these two values.
+(For even~$n$ we compute the $50^{\textrm{th}}$~percentile as~$v^s_{n/2}$,
+not as the median.) Example: the median of sample
+$\{$2.2, 3.2, 3.7, 4.4, 5.3, 5.7, 6.1, 6.4, 7.2, 7.8$\}$
+equals $(5.3\,{+}\,5.7)\,{/}\,2$~${=}$~5.5, see Figure~\ref{fig:example_quartiles}.
+
+Unlike the mean, the median is not sensitive to extreme values in the sample,
+i.e.\ it is robust to outliers. It works better as a measure of central tendency
+for heavy-tailed distributions and features that vary on exponential scale.
+However, the median is sensitive to small sample size.
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Interquartile mean]
+\OutputRowText{\OutputRowIDIQMean}
+For a sample of a quantitative feature, this is
+the mean of the values greater than or equal to the $1^{\textrm{st}}$ quartile
+and less than or equal the $3^{\textrm{rd}}$ quartile.
+In other words, it is a ``truncated mean'' where the lowest 25$\%$ and
+the highest 25$\%$ of the sorted values are omitted in its computation.
+The two ``border values'', i.e.\ the $1^{\textrm{st}}$ and the $3^{\textrm{rd}}$
+quartiles themselves, contribute to this mean only partially.
+This measure is occasionally used as the ``robust'' version of the mean
+that is less sensitive to the extreme values.
+
+To compute the measure, we sort the sample in the increasing order,
+preserving duplicates: $v^s_1 \leq v^s_2 \leq \ldots \leq v^s_n$.
+We set $j = \lceil n{/}4 \rceil$ for the $1^{\textrm{st}}$ quartile index
+and $k = \lceil 3n{/}4 \rceil$ for the $3^{\textrm{rd}}$ quartile index,
+then compute the following weighted mean:
+\begin{equation*}
+\frac{1}{3{/}4 - 1{/}4} \left[
+\left(\frac{j}{n} - \frac{1}{4}\right) v^s_j \,\,+
+\sum_{j<i<k} \left(\frac{i}{n} - \frac{i\,{-}\,1}{n}\right) v^s_i
+\,\,+\,\, \left(\frac{3}{4} - \frac{k\,{-}\,1}{n}\right) v^s_k\right]
+\end{equation*}
+In other words, all sample values between the $1^{\textrm{st}}$ and the $3^{\textrm{rd}}$
+quartile enter the sum with weights $2{/}n$, times their number of duplicates, while the
+two quartiles themselves enter the sum with reduced weights. The weights are proportional
+to the vertical steps in the empirical distribution function of the sample, see
+Figure~\ref{fig:example_quartiles} for an illustration.
+Example: the interquartile mean of sample
+$\{$2.2, 3.2, 3.7, 4.4, 5.3, 5.7, 6.1, 6.4, 7.2, 7.8$\}$ equals the sum
+$0.1 (3.7\,{+}\,6.4) + 0.2 (4.4\,{+}\,5.3\,{+}\,5.7\,{+}\,6.1)$,
+which equals~5.31.
+\end{Description}
+
+
+\paragraph{Dispersion measures.}
+Statistics that describe the amount of variation or spread in a quantitative
+(scale) data feature.
+\begin{Description}
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Variance]
+\OutputRowText{\OutputRowIDVariance}
+A measure of dispersion, or spread-out, of sample values around their mean,
+expressed in units that are the square of those of the feature itself.
+Computed as the sum of squared differences between the values
+in the sample and their mean, divided by one less than the number of
+values: $\sum_{i=1}^n (v_i - \bar{v})^2\,/\,(n\,{-}\,1)$ where
+$\bar{v}=\left(\sum_{i=1}^n v_i\right)\!/n$.
+Example: the variance of sample
+$\{$2.2, 3.2, 3.7, 4.4, 5.3, 5.7, 6.1, 6.4, 7.2, 7.8$\}$ equals~3.24.
+Note that at least two values ($n\geq 2$) are required to avoid division
+by zero. Sample variance is sensitive to outliers, even more than the mean.
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Standard deviation]
+\OutputRowText{\OutputRowIDStDeviation}
+A measure of dispersion around the mean, the square root of variance.
+Computed by taking the square root of the sample variance;
+see \emph{Variance} above on computing the variance.
+Example: the standard deviation of sample
+$\{$2.2, 3.2, 3.7, 4.4, 5.3, 5.7, 6.1, 6.4, 7.2, 7.8$\}$ equals~1.8.
+At least two values are required to avoid division by zero.
+Note that standard deviation is sensitive to outliers.
+
+Standard deviation is used in conjunction with the mean to determine
+an interval containing a given percentage of the feature values,
+assuming the normal distribution. In a large sample from a normal
+distribution, around 68\% of the cases fall within one standard
+deviation and around 95\% of cases fall within two standard deviations
+of the mean. For example, if the mean age is 45 with a standard deviation
+of 10, around 95\% of the cases would be between 25 and 65 in a normal
+distribution.
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Coefficient of variation]
+\OutputRowText{\OutputRowIDCoeffVar}
+The ratio of the standard deviation to the mean, i.e.\ the
+\emph{relative} standard deviation, of a quantitative feature sample.
+Computed by dividing the sample \emph{standard deviation} by the
+sample \emph{mean}, see above for their computation details.
+Example: the coefficient of variation for sample
+$\{$2.2, 3.2, 3.7, 4.4, 5.3, 5.7, 6.1, 6.4, 7.2, 7.8$\}$
+equals 1.8$\,{/}\,$5.2~${\approx}$~0.346.
+
+This metric is used primarily with non-negative features such as
+financial or population data. It is sensitive to outliers.
+Note: zero mean causes division by zero, returning infinity or \texttt{NaN}.
+At least two values (records) are required to compute the standard deviation.
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Minimum]
+\OutputRowText{\OutputRowIDMinimum}
+The smallest value of a quantitative sample, computed as $\min v = v^s_1$.
+Example: the minimum of sample
+$\{$2.2, 3.2, 3.7, 4.4, 5.3, 5.7, 6.1, 6.4, 7.2, 7.8$\}$
+equals~2.2.
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Maximum]
+\OutputRowText{\OutputRowIDMaximum}
+The largest value of a quantitative sample, computed as $\max v = v^s_n$.
+Example: the maximum of sample
+$\{$2.2, 3.2, 3.7, 4.4, 5.3, 5.7, 6.1, 6.4, 7.2, 7.8$\}$
+equals~7.8.
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Range]
+\OutputRowText{\OutputRowIDRange}
+The difference between the largest and the smallest value of a quantitative
+sample, computed as $\max v - \min v = v^s_n - v^s_1$.
+It provides information about the overall spread of the sample values.
+Example: the range of sample
+$\{$2.2, 3.2, 3.7, 4.4, 5.3, 5.7, 6.1, 6.4, 7.2, 7.8$\}$
+equals 7.8$\,{-}\,$2.2~${=}$~5.6.
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Standard error of the mean]
+\OutputRowText{\OutputRowIDStErrorMean}
+A measure of how much the value of the sample mean may vary from sample
+to sample taken from the same (hypothesized) distribution of the feature.
+It helps to roughly bound the distribution mean, i.e.\
+the limit of the sample mean as the sample size tends to infinity.
+Under certain assumptions (e.g.\ normality and large sample), the difference
+between the distribution mean and the sample mean is unlikely to exceed
+2~standard errors.
+
+The measure is computed by dividing the sample standard deviation
+by the square root of the number of values~$n$; see \emph{standard deviation}
+for its computation details. Ensure $n\,{\geq}\,2$ to avoid division by~0.
+Example: for sample
+$\{$2.2, 3.2, 3.7, 4.4, 5.3, 5.7, 6.1, 6.4, 7.2, 7.8$\}$
+with the mean of~5.2 the standard error of the mean
+equals 1.8$\,{/}\sqrt{10}$~${\approx}$~0.569.
+
+Note that the standard error itself is subject to sample randomness.
+Its accuracy as an error estimator may be low if the sample size is small
+or \mbox{non-i.i.d.}, if there are outliers, or if the distribution has
+heavy tails.
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+% \item[\it Quartiles]
+% \OutputRowText{\OutputRowIDQuartiles}
+% %%% dsDefn %%%%
+% The values of a quantitative feature
+% that divide an ordered/sorted set of data records into four equal-size groups.
+% The $1^{\textrm{st}}$ quartile, or the $25^{\textrm{th}}$ percentile, splits
+% the sorted data into the lowest $25\%$ and the highest~$75\%$. In other words,
+% it is the middle value between the minimum and the median. The $2^{\textrm{nd}}$
+% quartile is the median itself, the value that separates the higher half of
+% the data (in the sorted order) from the lower half. Finally, the $3^{\textrm{rd}}$
+% quartile, or the $75^{\textrm{th}}$ percentile, divides the sorted data into
+% lowest $75\%$ and highest~$25\%$.\par
+% %%% dsComp %%%%
+% To compute the quartiles for a data column \texttt{X[,i]} with $n$ numerical values
+% we sort it in the increasing order, preserving duplicates, then return
+% \texttt{X}${}^{\textrm{sort}}$\texttt{[}$k$\texttt{,i]}
+% where $k = \lceil pn \rceil$ for $p = 0.25$, $0.5$, and~$0.75$.
+% When $n$ is even, the $2^{\textrm{nd}}$ quartile (the median) is further adjusted
+% to equal the mean of two middle values
+% $\texttt{X}^{\textrm{sort}}\texttt{[}n{/}2\texttt{,i]}$ and
+% $\texttt{X}^{\textrm{sort}}\texttt{[}n{/}2\,{+}\,1\texttt{,i]}$.
+% %%% dsWarn %%%%
+% We assume that the feature column does not contain \texttt{NaN}s or coded non-numeric values.
+% %%% dsExmpl %%%
+% \textbf{Example(s).}
+\end{Description}
+
+
+\paragraph{Shape measures.}
+Statistics that describe the shape and symmetry of the quantitative (scale)
+feature distribution estimated from a sample of its values.
+\begin{Description}
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Skewness]
+\OutputRowText{\OutputRowIDSkewness}
+It measures how symmetrically the values of a feature are spread out
+around the mean. A significant positive skewness implies a longer (or fatter)
+right tail, i.e. feature values tend to lie farther away from the mean on the
+right side. A significant negative skewness implies a longer (or fatter) left
+tail. The normal distribution is symmetric and has a skewness value of~0;
+however, its sample skewness is likely to be nonzero, just close to zero.
+As a guideline, a skewness value more than twice its standard error is taken
+to indicate a departure from symmetry.
+
+Skewness is computed as the $3^{\textrm{rd}}$~central moment divided by the cube
+of the standard deviation. We estimate the $3^{\textrm{rd}}$~central moment as
+the sum of cubed differences between the values in the feature column and their
+sample mean, divided by the number of values:
+$\sum_{i=1}^n (v_i - \bar{v})^3 / n$
+where $\bar{v}=\left(\sum_{i=1}^n v_i\right)\!/n$.
+The standard deviation is computed
+as described above in \emph{standard deviation}. To avoid division by~0,
+at least two different sample values are required. Example: for sample
+$\{$2.2, 3.2, 3.7, 4.4, 5.3, 5.7, 6.1, 6.4, 7.2, 7.8$\}$
+with the mean of~5.2 and the standard deviation of~1.8
+skewness is estimated as $-1.0728\,{/}\,1.8^3 \approx -0.184$.
+Note: skewness is sensitive to outliers.
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Standard error in skewness]
+\OutputRowText{\OutputRowIDStErrorSkewness}
+A measure of how much the sample skewness may vary from sample to sample,
+assuming that the feature is normally distributed, which makes its
+distribution skewness equal~0.
+Given the number~$n$ of sample values, the standard error is computed as
+\begin{equation*}
+\sqrt{\frac{6n\,(n-1)}{(n-2)(n+1)(n+3)}}
+\end{equation*}
+This measure can tell us, for example:
+\begin{Itemize}
+\item If the sample skewness lands within two standard errors from~0, its
+positive or negative sign is non-significant, may just be accidental.
+\item If the sample skewness lands outside this interval, the feature
+is unlikely to be normally distributed.
+\end{Itemize}
+At least 3~values ($n\geq 3$) are required to avoid arithmetic failure.
+Note that the standard error is inaccurate if the feature distribution is
+far from normal or if the number of samples is small.
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Kurtosis]
+\OutputRowText{\OutputRowIDKurtosis}
+As a distribution parameter, kurtosis is a measure of the extent to which
+feature values cluster around a central point. In other words, it quantifies
+``peakedness'' of the distribution: how tall and sharp the central peak is
+relative to a standard bell curve.
+
+Positive kurtosis (\emph{leptokurtic} distribution) indicates that, relative
+to a normal distribution:
+\begin{Itemize}
+\item observations cluster more about the center (peak-shaped),
+\item the tails are thinner at non-extreme values,
+\item the tails are thicker at extreme values.
+\end{Itemize}
+Negative kurtosis (\emph{platykurtic} distribution) indicates that, relative
+to a normal distribution:
+\begin{Itemize}
+\item observations cluster less about the center (box-shaped),
+\item the tails are thicker at non-extreme values,
+\item the tails are thinner at extreme values.
+\end{Itemize}
+Kurtosis of a normal distribution is zero; however, the sample kurtosis
+(computed here) is likely to deviate from zero.
+
+Sample kurtosis is computed as the $4^{\textrm{th}}$~central moment divided
+by the $4^{\textrm{th}}$~power of the standard deviation, minus~3.
+We estimate the $4^{\textrm{th}}$~central moment as the sum of the
+$4^{\textrm{th}}$~powers of differences between the values in the feature column
+and their sample mean, divided by the number of values:
+$\sum_{i=1}^n (v_i - \bar{v})^4 / n$
+where $\bar{v}=\left(\sum_{i=1}^n v_i\right)\!/n$.
+The standard deviation is computed as described above, see \emph{standard deviation}.
+
+Note that kurtosis is sensitive to outliers, and requires at least two different
+sample values. Example: for sample
+$\{$2.2, 3.2, 3.7, 4.4, 5.3, 5.7, 6.1, 6.4, 7.2, 7.8$\}$
+with the mean of~5.2 and the standard deviation of~1.8,
+sample kurtosis equals $16.6962\,{/}\,1.8^4 - 3 \approx -1.41$.
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Standard error in kurtosis]
+\OutputRowText{\OutputRowIDStErrorCurtosis}
+A measure of how much the sample kurtosis may vary from sample to sample,
+assuming that the feature is normally distributed, which makes its
+distribution kurtosis equal~0.
+Given the number~$n$ of sample values, the standard error is computed as
+\begin{equation*}
+\sqrt{\frac{24n\,(n-1)^2}{(n-3)(n-2)(n+3)(n+5)}}
+\end{equation*}
+This measure can tell us, for example:
+\begin{Itemize}
+\item If the sample kurtosis lands within two standard errors from~0, its
+positive or negative sign is non-significant, may just be accidental.
+\item If the sample kurtosis lands outside this interval, the feature
+is unlikely to be normally distributed.
+\end{Itemize}
+At least 4~values ($n\geq 4$) are required to avoid arithmetic failure.
+Note that the standard error is inaccurate if the feature distribution is
+far from normal or if the number of samples is small.
+\end{Description}
+
+
+\paragraph{Categorical measures.} Statistics that describe the sample of
+a categorical feature, either nominal or ordinal. We represent all
+categories by integers from~1 to the number of categories; we call
+these integers \emph{category~IDs}.
+\begin{Description}
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Number of categories]
+\OutputRowText{\OutputRowIDNumCategories}
+The maximum category~ID that occurs in the sample. Note that some
+categories with~IDs \emph{smaller} than this maximum~ID may have
+no~occurrences in the sample, without reducing the number of categories.
+However, any categories with~IDs \emph{larger} than the maximum~ID with
+no occurrences in the sample will not be counted.
+Example: in sample $\{$1, 3, 3, 3, 3, 4, 4, 5, 7, 7, 7, 7, 8, 8, 8$\}$
+the number of categories is reported as~8. Category~IDs 2 and~6, which have
+zero occurrences, are still counted; but if there is a category with
+ID${}=9$ and zero occurrences, it is not counted.
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Mode]
+\OutputRowText{\OutputRowIDMode}
+The most frequently occurring category value.
+If several values share the greatest frequency of occurrence, then each
+of them is a mode; but here we report only the smallest of these modes.
+Example: in sample $\{$1, 3, 3, 3, 3, 4, 4, 5, 7, 7, 7, 7, 8, 8, 8$\}$
+the modes are 3 and~7, with 3 reported.
+
+Computed by counting the number of occurrences for each category,
+then taking the smallest category~ID that has the maximum count.
+Note that the sample modes may be different from the distribution modes,
+i.e.\ the categories whose (hypothesized) underlying probability is the
+maximum over all categories.
+%%%%%%%%%%%%%%%%%%%% DESCRIPTIVE STATISTIC %%%%%%%%%%%%%%%%%%%%
+\item[\it Number of modes]
+\OutputRowText{\OutputRowIDNumModes}
+The number of category values that each have the largest frequency
+count in the sample.
+Example: in sample $\{$1, 3, 3, 3, 3, 4, 4, 5, 7, 7, 7, 7, 8, 8, 8$\}$
+there are two category IDs (3 and~7) that occur the maximum count of 4~times;
+hence, we return~2.
+
+Computed by counting the number of occurrences for each category,
+then counting how many categories have the maximum count.
+Note that the sample modes may be different from the distribution modes,
+i.e.\ the categories whose (hypothesized) underlying probability is the
+maximum over all categories.
+\end{Description}
+
+
+\smallskip
+\noindent{\bf Returns}
+\smallskip
+
+The output matrix containing all computed statistics is of size $17$~rows and
+as many columns as in the input matrix~\texttt{X}. Each row corresponds to
+a particular statistic, according to the convention specified in
+Table~\ref{table:univars}. The first $14$~statistics are applicable for
+\emph{scale} columns, and the last $3$~statistics are applicable for categorical,
+i.e.\ nominal and ordinal, columns.
+
+
+\pagebreak[2]
+
+\smallskip
+\noindent{\bf Examples}
+\smallskip
+
+{\hangindent=\parindent\noindent\tt
+\hml -f \UnivarScriptName{} -nvargs X=/user/biadmin/X.mtx
+ TYPES=/user/biadmin/types.mtx
+ STATS=/user/biadmin/stats.mtx
+
+}
http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/a891c0a1/Algorithms Reference/GLM.tex
----------------------------------------------------------------------
diff --git a/Algorithms Reference/GLM.tex b/Algorithms Reference/GLM.tex
new file mode 100644
index 0000000..e4d13bf
--- /dev/null
+++ b/Algorithms Reference/GLM.tex
@@ -0,0 +1,409 @@
+\subsection{Generalized Linear Models (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.
http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/a891c0a1/Algorithms Reference/GLMpredict.tex
----------------------------------------------------------------------
diff --git a/Algorithms Reference/GLMpredict.tex b/Algorithms Reference/GLMpredict.tex
new file mode 100644
index 0000000..7fb91c8
--- /dev/null
+++ b/Algorithms Reference/GLMpredict.tex
@@ -0,0 +1,453 @@
+\subsection{Regression Scoring and Prediction}
+
+\noindent{\bf Description}
+\smallskip
+
+Script {\tt GLM-predict.dml} is intended to cover all linear model based regressions,
+including linear regression, binomial and multinomial logistic regression, and GLM
+regressions (Poisson, gamma, binomial with probit link~etc.). Having just one scoring
+script for all these regressions simplifies maintenance and enhancement while ensuring
+compatible interpretations for output statistics.
+
+The script performs two functions, prediction and scoring. To perform prediction,
+the script takes two matrix inputs: a collection of records $X$ (without the response
+attribute) and the estimated regression parameters~$B$, also known as~$\beta$.
+To perform scoring, in addition to $X$ and~$B$, the script takes the matrix of actual
+response values~$Y$ that are compared to the predictions made with $X$ and~$B$. Of course
+there are other, non-matrix, input arguments that specify the model and the output
+format, see below for the full list.
+
+We assume that our test/scoring dataset is given by $n\,{\times}\,m$-matrix $X$ of
+numerical feature vectors, where each row~$x_i$ represents one feature vector of one
+record; we have \mbox{$\dim x_i = m$}. Each record also includes the response
+variable~$y_i$ that may be numerical, single-label categorical, or multi-label categorical.
+A single-label categorical $y_i$ is an integer category label, one label per record;
+a multi-label $y_i$ is a vector of integer counts, one count for each possible label,
+which represents multiple single-label events (observations) for the same~$x_i$. Internally
+we convert single-label categoricals into multi-label categoricals by replacing each
+label~$l$ with an indicator vector~$(0,\ldots,0,1_l,0,\ldots,0)$. In prediction-only
+tasks the actual $y_i$'s are not needed to the script, but they are needed for scoring.
+
+To perform prediction, the script matrix-multiplies $X$ and $B$, adding the intercept
+if available, then applies the inverse of the model's link function.
+All GLMs assume that the linear combination of the features in~$x_i$ and the betas
+in~$B$ determines the means~$\mu_i$ of the~$y_i$'s (in numerical or multi-label
+categorical form) with $\dim \mu_i = \dim y_i$. The observed $y_i$ is assumed to follow
+a specified GLM family distribution $\Prob[y\mid \mu_i]$ with mean(s)~$\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*}
+If $y_i$ is numerical, the predicted mean $\mu_i$ is a real number. Then our script's
+output matrix $M$ is the $n\,{\times}\,1$-vector of these means~$\mu_i$.
+Note that $\mu_i$ predicts the mean of $y_i$, not the actual~$y_i$. For example,
+in Poisson distribution, the mean is usually fractional, but the actual~$y_i$ is
+always integer.
+
+If $y_i$ is categorical, i.e.\ a vector of label counts for record~$i$, then $\mu_i$
+is a vector of non-negative real numbers, one number $\mu_{i,l}$ per each label~$l$.
+In this case we divide the $\mu_{i,l}$ by their sum $\sum_l \mu_{i,l}$ to obtain
+predicted label probabilities~\mbox{$p_{i,l}\in [0, 1]$}. The output matrix $M$ is
+the $n \times (k\,{+}\,1)$-matrix of these probabilities, where $n$ is the number of
+records and $k\,{+}\,1$ is the number of categories\footnote{We use $k+1$ because
+there are $k$ non-baseline categories and one baseline category, with regression
+parameters $B$ having $k$~columns.}. Note again that we do not predict the labels
+themselves, nor their actual counts per record, but we predict the labels' probabilities.
+
+Going from predicted probabilities to predicted labels, in the single-label categorical
+case, requires extra information such as the cost of false positive versus
+false negative errors. For example, if there are 5 categories and we \emph{accurately}
+predicted their probabilities as $(0.1, 0.3, 0.15, 0.2, 0.25)$, just picking the
+highest-probability label would be wrong 70\% of the time, whereas picking the
+lowest-probability label might be right if, say, it represents a diagnosis of cancer
+or another rare and serious outcome. Hence, we keep this step outside the scope of
+{\tt GLM-predict.dml} for now.
+
+\smallskip
+\noindent{\bf Usage}
+\smallskip
+
+{\hangindent=\parindent\noindent\it%
+{\tt{}-f }path/\/{\tt{}GLM-predict.dml}
+{\tt{} -nvargs}
+{\tt{} X=}path/file
+{\tt{} Y=}path/file
+{\tt{} B=}path/file
+{\tt{} M=}path/file
+{\tt{} O=}path/file
+{\tt{} dfam=}int
+{\tt{} vpow=}double
+{\tt{} link=}int
+{\tt{} lpow=}double
+{\tt{} disp=}double
+{\tt{} fmt=}format
+
+}
+
+\smallskip
+\noindent{\bf Arguments}
+\begin{Description}
+\item[{\tt X}:]
+Location (on HDFS) to read the $n\,{\times}\,m$-matrix $X$ of feature vectors, each row
+constitutes one feature vector (one record)
+\item[{\tt Y}:] (default:\mbox{ }{\tt " "})
+Location to read the response matrix $Y$ needed for scoring (but optional for prediction),
+with the following dimensions: \\
+ $n \:{\times}\: 1$: acceptable for all distributions ({\tt dfam=1} or {\tt 2} or {\tt 3}) \\
+ $n \:{\times}\: 2$: for binomial ({\tt dfam=2}) if given by (\#pos, \#neg) counts \\
+ $n \:{\times}\: k\,{+}\,1$: for multinomial ({\tt dfam=3}) if given by category counts
+\item[{\tt M}:] (default:\mbox{ }{\tt " "})
+Location to write, if requested, the matrix of predicted response means (for {\tt dfam=1}) or
+probabilities (for {\tt dfam=2} or {\tt 3}):\\
+ $n \:{\times}\: 1$: for power-type distributions ({\tt dfam=1}) \\
+ $n \:{\times}\: 2$: for binomial distribution ({\tt dfam=2}), col\#~2 is the ``No'' probability \\
+ $n \:{\times}\: k\,{+}\,1$: for multinomial logit ({\tt dfam=3}), col\#~$k\,{+}\,1$ is for the baseline
+\item[{\tt B}:]
+Location to read matrix $B$ of the \mbox{betas}, i.e.\ estimated GLM regression parameters,
+with the intercept at row\#~$m\,{+}\,1$ if available:\\
+ $\dim(B) \,=\, m \:{\times}\: k$: do not add intercept \\
+ $\dim(B) \,=\, m\,{+}\,1 \:{\times}\: k$: add intercept as given by the last $B$-row \\
+ if $k > 1$, use only $B${\tt [, 1]} unless it is Multinomial Logit ({\tt dfam=3})
+\item[{\tt O}:] (default:\mbox{ }{\tt " "})
+Location to store the CSV-file with goodness-of-fit statistics defined in
+Table~\ref{table:GLMpred:stats}, the default is to print them to the standard output
+\item[{\tt dfam}:] (default:\mbox{ }{\tt 1})
+GLM distribution family code to specify the type of distribution $\Prob[y\,|\,\mu]$
+that we assume: \\
+{\tt 1} = power distributions with $\Var(y) = \mu^{\alpha}$, see Table~\ref{table:commonGLMs};\\
+{\tt 2} = binomial;
+{\tt 3} = multinomial logit
+\item[{\tt vpow}:] (default:\mbox{ }{\tt 0.0})
+Power for variance defined as (mean)${}^{\textrm{power}}$ (ignored if {\tt dfam}$\,{\neq}\,1$):
+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)$, ignored for
+multinomial logit ({\tt dfam=3}):\\
+{\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})
+Power for link function defined as (mean)${}^{\textrm{power}}$ (ignored if {\tt link}$\,{\neq}\,1$):
+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 disp}:] (default:\mbox{ }{\tt 1.0})
+Dispersion value, when available; must be positive
+\item[{\tt fmt}:] (default:\mbox{ }{\tt "text"})
+Matrix {\tt M} file output format, such as {\tt text}, {\tt mm}, or {\tt csv};
+see read/write functions in SystemML Language Reference for details.
+\end{Description}
+
+
+\begin{table}[t]\small\centerline{%
+\begin{tabular}{|lccl|}
+\hline
+Name & \hspace{-0.6em}CID\hspace{-0.5em} & \hspace{-0.3em}Disp?\hspace{-0.6em} & Meaning \\
+\hline
+{\tt LOGLHOOD\_Z} & & + & Log-likelihood $Z$-score (in st.\ dev.'s from the mean) \\
+{\tt LOGLHOOD\_Z\_PVAL} & & + & Log-likelihood $Z$-score p-value, two-sided \\
+{\tt PEARSON\_X2} & & + & Pearson residual $X^2$-statistic \\
+{\tt PEARSON\_X2\_BY\_DF} & & + & Pearson $X^2$ divided by degrees of freedom \\
+{\tt PEARSON\_X2\_PVAL} & & + & Pearson $X^2$ p-value \\
+{\tt DEVIANCE\_G2} & & + & Deviance from the saturated model $G^2$-statistic \\
+{\tt DEVIANCE\_G2\_BY\_DF} & & + & Deviance $G^2$ divided by degrees of freedom \\
+{\tt DEVIANCE\_G2\_PVAL} & & + & Deviance $G^2$ p-value \\
+{\tt AVG\_TOT\_Y} & + & & $Y$-column average for an individual response value \\
+{\tt STDEV\_TOT\_Y} & + & & $Y$-column st.\ dev.\ for an individual response value \\
+{\tt AVG\_RES\_Y} & + & & $Y$-column residual average of $Y - \mathop{\mathrm{pred.}\,\mathrm{mean}}(Y|X)$ \\
+{\tt STDEV\_RES\_Y} & + & & $Y$-column residual st.\ dev.\ of $Y - \mathop{\mathrm{pred.}\,\mathrm{mean}}(Y|X)$ \\
+{\tt PRED\_STDEV\_RES} & + & + & Model-predicted $Y$-column residual st.\ deviation\\
+{\tt PLAIN\_R2} & + & & Plain $R^2$ of $Y$-column residual with bias included \\
+{\tt ADJUSTED\_R2} & + & & Adjusted $R^2$ of $Y$-column residual w.\ bias included \\
+{\tt PLAIN\_R2\_NOBIAS} & + & & Plain $R^2$ of $Y$-column residual, bias subtracted \\
+{\tt ADJUSTED\_R2\_NOBIAS} & + & & Adjusted $R^2$ of $Y$-column residual, bias subtracted \\
+\hline
+\end{tabular}}
+\caption{The above goodness-of-fit statistics are provided in CSV format, one per each line, with four
+columns: (Name, [CID], [Disp?], Value). The columns are:
+``Name'' is the string identifier for the statistic, see the table;
+``CID'' is an optional integer value that specifies the $Y$-column index for \mbox{per-}column statistics
+(note that a bi-/multinomial one-column {\tt Y}-input is converted into multi-column);
+``Disp?'' is an optional Boolean value ({\tt TRUE} or {\tt FALSE}) that tells us
+whether or not scaling by the input dispersion parameter {\tt disp} has been applied to this
+statistic;
+``Value'' is the value of the statistic.}
+\label{table:GLMpred:stats}
+\end{table}
+
+\noindent{\bf Details}
+\smallskip
+
+The output matrix $M$ of predicted means (or probabilities) is computed by matrix-multiplying $X$
+with the first column of~$B$ or with the whole~$B$ in the multinomial case, adding the intercept
+if available (conceptually, appending an extra column of ones to~$X$); then applying the inverse
+of the model's link function. The difference between ``means'' and ``probabilities'' in the
+categorical case becomes significant when there are ${\geq}\,2$ observations per record
+(with the multi-label records) or when the labels such as $-1$ and~$1$ are viewed and averaged
+as numerical response values (with the single-label records). To avoid any \mbox{mix-up} or
+information loss, we separately return the predicted probability of each category label for each
+record.
+
+When the ``actual'' response values $Y$ are available, the summary statistics are computed
+and written out as described in Table~\ref{table:GLMpred:stats}. Below we discuss each of
+these statistics in detail. Note that in the categorical case (binomial and multinomial)
+$Y$ is internally represented as the matrix of observation counts for each label in each record,
+rather than just the label~ID for each record. The input~$Y$ may already be a matrix of counts,
+in which case it is used as-is. But if $Y$ is given as a vector of response labels, each
+response label is converted into an indicator vector $(0,\ldots,0,1_l,0,\ldots,0)$ where~$l$
+is the label~ID for this record. All negative (e.g.~$-1$) or zero label~IDs are converted to
+the $1 + {}$maximum label~ID. The largest label~ID is viewed as the ``baseline'' as explained
+in the section on Multinomial Logistic Regression. We assume that there are $k\geq 1$
+non-baseline categories and one (last) baseline category.
+
+We also estimate residual variances for each response value, although we do not output them,
+but use them only inside the summary statistics, scaled and unscaled by the input dispersion
+parameter {\tt disp}, as described below.
+
+\smallskip
+{\tt LOGLHOOD\_Z} and {\tt LOGLHOOD\_Z\_PVAL} statistics measure how far the log-likelihood
+of~$Y$ deviates from its expected value according to the model. The script implements them
+only for the binomial and the multinomial distributions, returning NaN for all other distributions.
+Pearson's~$X^2$ and deviance~$G^2$ often perform poorly with bi- and multinomial distributions
+due to low cell counts, hence we need this extra goodness-of-fit measure. To compute these
+statistics, we use:
+\begin{Itemize}
+\item the $n\times (k\,{+}\,1)$-matrix~$Y$ of multi-label response counts, in which $y_{i,j}$
+is the number of times label~$j$ was observed in record~$i$;
+\item the model-estimated probability matrix~$P$ of the same dimensions that satisfies
+$\sum_{j=1}^{k+1} p_{i,j} = 1$ for all~$i=1,\ldots,n$ and where $p_{i,j}$ is the model
+probability of observing label~$j$ in record~$i$;
+\item the $n\,{\times}\,1$-vector $N$ where $N_i$ is the aggregated count of observations
+in record~$i$ (all $N_i = 1$ if each record has only one response label).
+\end{Itemize}
+We start by computing the multinomial log-likelihood of $Y$ given~$P$ and~$N$, as well as
+the expected log-likelihood given a random~$Y$ and the variance of this log-likelihood if
+$Y$ indeed follows the proposed distribution:
+\begin{align*}
+\ell (Y) \,\,&=\,\, \log \Prob[Y \,|\, P, N] \,\,=\,\, \sum_{i=1}^{n} \,\sum_{j=1}^{k+1} \,y_{i,j}\log p_{i,j} \\
+\E_Y \ell (Y) \,\,&=\,\, \sum_{i=1}^{n}\, \sum_{j=1}^{k+1} \,\mu_{i,j} \log p_{i,j}
+ \,\,=\,\, \sum_{i=1}^{n}\, N_i \,\sum_{j=1}^{k+1} \,p_{i,j} \log p_{i,j} \\
+\Var_Y \ell (Y) \,&=\, \sum_{i=1}^{n} \,N_i \left(\sum_{j=1}^{k+1} \,p_{i,j} \big(\log p_{i,j}\big)^2
+ - \Bigg( \sum_{j=1}^{k+1} \,p_{i,j} \log p_{i,j}\Bigg) ^ {\!\!2\,} \right)
+\end{align*}
+Then we compute the $Z$-score as the difference between the actual and the expected
+log-likelihood $\ell(Y)$ divided by its expected standard deviation, and its two-sided
+p-value in the Normal distribution assumption ($\ell(Y)$ should approach normality due
+to the Central Limit Theorem):
+\begin{equation*}
+Z \,=\, \frac {\ell(Y) - \E_Y \ell(Y)}{\sqrt{\Var_Y \ell(Y)}};\quad
+\mathop{\textrm{p-value}}(Z) \,=\, \Prob \Big[\,\big|\mathop{\textrm{Normal}}(0,1)\big| \, > \, |Z|\,\Big]
+\end{equation*}
+A low p-value would indicate ``underfitting'' if $Z\ll 0$ or ``overfitting'' if $Z\gg 0$. Here
+``overfitting'' means that higher-probability labels occur more often than their probabilities
+suggest.
+
+We also apply the dispersion input ({\tt disp}) to compute the ``scaled'' version of the $Z$-score
+and its p-value. Since $\ell(Y)$ is a linear function of~$Y$, multiplying the GLM-predicted
+variance of~$Y$ by {\tt disp} results in multiplying $\Var_Y \ell(Y)$ by the same {\tt disp}. This, in turn,
+translates into dividing the $Z$-score by the square root of the dispersion:
+\begin{equation*}
+Z_{\texttt{disp}} \,=\, \big(\ell(Y) \,-\, \E_Y \ell(Y)\big) \,\big/\, \sqrt{\texttt{disp}\cdot\Var_Y \ell(Y)}
+\,=\, Z / \sqrt{\texttt{disp}}
+\end{equation*}
+Finally, we recalculate the p-value with this new $Z$-score.
+
+\smallskip
+{\tt PEARSON\_X2}, {\tt PEARSON\_X2\_BY\_DF}, and {\tt PEARSON\_X2\_PVAL}:
+Pearson's residual $X^2$-statistic is a commonly used goodness-of-fit measure for linear models~\cite{McCullagh1989:GLM}.
+The idea is to measure how well the model-predicted means and variances match the actual behavior
+of response values. For each record $i$, we estimate the mean $\mu_i$ and the variance $v_i$
+(or $\texttt{disp}\cdot v_i$) and use them to normalize the residual:
+$r_i = (y_i - \mu_i) / \sqrt{v_i}$. These normalized residuals are then squared, aggregated
+by summation, and tested against an appropriate $\chi^2$~distribution. The computation of~$X^2$
+is slightly different for categorical data (bi- and multinomial) than it is for numerical data,
+since $y_i$ has multiple correlated dimensions~\cite{McCullagh1989:GLM}:
+\begin{equation*}
+X^2\,\textrm{(numer.)} \,=\, \sum_{i=1}^{n}\, \frac{(y_i - \mu_i)^2}{v_i};\quad
+X^2\,\textrm{(categ.)} \,=\, \sum_{i=1}^{n}\, \sum_{j=1}^{k+1} \,\frac{(y_{i,j} - N_i
+\hspace{0.5pt} p_{i,j})^2}{N_i \hspace{0.5pt} p_{i,j}}
+\end{equation*}
+The number of degrees of freedom~\#d.f.\ for the $\chi^2$~distribution is $n - m$ for numerical data and
+$(n - m)k$ for categorical data, where $k = \mathop{\texttt{ncol}}(Y) - 1$. Given the dispersion
+parameter {\tt disp}, the $X^2$ statistic is scaled by division: \mbox{$X^2_{\texttt{disp}} = X^2 / \texttt{disp}$}.
+If the dispersion is accurate, $X^2 / \texttt{disp}$ should be close to~\#d.f. In fact, $X^2 / \textrm{\#d.f.}$
+over the \emph{training} data is the dispersion estimator used in our {\tt GLM.dml} script,
+see~(\ref{eqn:dispersion}). Here we provide $X^2 / \textrm{\#d.f.}$ and $X^2_{\texttt{disp}} / \textrm{\#d.f.}$
+as {\tt PEARSON\_X2\_BY\_DF} to enable dispersion comparison between the training data and
+the test data.
+
+NOTE: For categorical data, both Pearson's $X^2$ and the deviance $G^2$ are unreliable (i.e.\ do not
+approach the $\chi^2$~distribution) unless the predicted means of multi-label counts
+$\mu_{i,j} = N_i \hspace{0.5pt} p_{i,j}$ are fairly large: all ${\geq}\,1$ and 80\% are
+at least~$5$~\cite{Cochran1954:chisq}. They should not be used for ``one label per record'' categoricals.
+
+\smallskip
+{\tt DEVIANCE\_G2}, {\tt DEVIANCE\_G2\_BY\_DF}, and {\tt DEVIANCE\_G2\_PVAL}:
+Deviance $G^2$ is the log of the likelihood ratio between the ``saturated'' model and the
+linear model being tested for the given dataset, multiplied by two:
+\begin{equation}
+G^2 \,=\, 2 \,\log \frac{\Prob[Y \mid \textrm{saturated model}\hspace{0.5pt}]}%
+{\Prob[Y \mid \textrm{tested linear model}\hspace{0.5pt}]}
+\label{eqn:GLMpred:deviance}
+\end{equation}
+The ``saturated'' model sets the mean $\mu_i^{\mathrm{sat}}$ to equal~$y_i$ for every record
+(for categorical data, $p^{\mathrm{sat}}_{i,j} = y_{i,j} / N_i$), which represents the ``perfect fit.''
+For records with $y_{i,j} \in \{0, N_i\}$ or otherwise at a boundary, by continuity we set
+$0 \log 0 = 0$. The GLM~likelihood functions defined in~(\ref{eqn:GLM}) become simplified
+in ratio~(\ref{eqn:GLMpred:deviance}) due to canceling out the term $c(y, a)$ since it is the same
+in both models.
+
+The log of a likelihood ratio between two nested models, times two, is known to approach
+a $\chi^2$ distribution as $n\to\infty$ if both models have fixed parameter spaces.
+But this is not the case for the ``saturated'' model: it adds more parameters with each record.
+In practice, however, $\chi^2$ distributions are used to compute the p-value of~$G^2$~\cite{McCullagh1989:GLM}.
+The number of degrees of freedom~\#d.f.\ and the treatment of dispersion are the same as for
+Pearson's~$X^2$, see above.
+
+\Paragraph{Column-wise statistics.} The rest of the statistics are computed separately
+for each column of~$Y$. As explained above, $Y$~has two or more columns in bi- and multinomial case,
+either at input or after conversion. Moreover, each $y_{i,j}$ in record~$i$ with $N_i \geq 2$ is
+counted as $N_i$ separate observations $y_{i,j,l}$ of 0 or~1 (where $l=1,\ldots,N_i$) with
+$y_{i,j}$~ones and $N_i-y_{i,j}$ zeros.
+For power distributions, including linear regression, $Y$~has only one column and all
+$N_i = 1$, so the statistics are computed for all~$Y$ with each record counted once.
+Below we denote $N = \sum_{i=1}^n N_i \,\geq n$.
+Here is the total average and the residual average (residual bias) of~$y_{i,j,l}$ for each $Y$-column:
+\begin{equation*}
+\texttt{AVG\_TOT\_Y}_j \,=\, \frac{1}{N} \sum_{i=1}^n y_{i,j}; \quad
+\texttt{AVG\_RES\_Y}_j \,=\, \frac{1}{N} \sum_{i=1}^n \, (y_{i,j} - \mu_{i,j})
+\end{equation*}
+Dividing by~$N$ (rather than~$n$) gives the averages for~$y_{i,j,l}$ (rather than~$y_{i,j}$).
+The total variance, and the standard deviation, for individual observations~$y_{i,j,l}$ is
+estimated from the total variance for response values~$y_{i,j}$ using independence assumption:
+$\Var y_{i,j} = \Var \sum_{l=1}^{N_i} y_{i,j,l} = \sum_{l=1}^{N_i} \Var y_{i,j,l}$.
+This allows us to estimate the sum of squares for~$y_{i,j,l}$ via the sum of squares for~$y_{i,j}$:
+\begin{equation*}
+\texttt{STDEV\_TOT\_Y}_j \,=\,
+\Bigg[\frac{1}{N-1} \sum_{i=1}^n \Big( y_{i,j} - \frac{N_i}{N} \sum_{i'=1}^n y_{i'\!,j}\Big)^2\Bigg]^{1/2}
+\end{equation*}
+Analogously, we estimate the standard deviation of the residual $y_{i,j,l} - \mu_{i,j,l}$:
+\begin{equation*}
+\texttt{STDEV\_RES\_Y}_j \,=\,
+\Bigg[\frac{1}{N-m'} \,\sum_{i=1}^n \Big( y_{i,j} - \mu_{i,j} - \frac{N_i}{N} \sum_{i'=1}^n (y_{i'\!,j} - \mu_{i'\!,j})\Big)^2\Bigg]^{1/2}
+\end{equation*}
+Here $m'=m$ if $m$ includes the intercept as a feature and $m'=m+1$ if it does not.
+The estimated standard deviations can be compared to the model-predicted residual standard deviation
+computed from the predicted means by the GLM variance formula and scaled by the dispersion:
+\begin{equation*}
+\texttt{PRED\_STDEV\_RES}_j \,=\, \Big[\frac{\texttt{disp}}{N} \, \sum_{i=1}^n \, v(\mu_{i,j})\Big]^{1/2}
+\end{equation*}
+We also compute the $R^2$ statistics for each column of~$Y$, see Table~\ref{table:GLMpred:R2} for details.
+We compute two versions of~$R^2$: in one version the residual sum-of-squares (RSS) includes any bias in
+the residual that might be present (due to the lack of, or inaccuracy in, the intercept); in the other
+version of~RSS the bias is subtracted by ``centering'' the residual. In both cases we subtract the bias in the total
+sum-of-squares (in the denominator), and $m'$ equals $m$~with the intercept or $m+1$ without the intercept.
+
+\begin{table}[t]\small\centerline{%
+\begin{tabular}{|c|c|}
+\multicolumn{2}{c}{$R^2$ where the residual sum-of-squares includes the bias contribution:} \\
+\hline
+\multicolumn{1}{|l|}{\tt PLAIN\_R2${}_j \,\,= {}$} & \multicolumn{1}{l|}{\tt ADJUSTED\_R2${}_j \,\,= {}$} \\
+$ \displaystyle 1 -
+\frac{\sum\limits_{i=1}^n \,(y_{i,j} - \mu_{i,j})^2}%
+{\sum\limits_{i=1}^n \Big(y_{i,j} - \frac{N_{i\mathstrut}}{N^{\mathstrut}}\!\! \sum\limits_{i'=1}^n \! y_{i'\!,j} \Big)^{\! 2}} $ &
+$ \displaystyle 1 - {\textstyle\frac{N_{\mathstrut} - 1}{N^{\mathstrut} - m}} \,
+\frac{\sum\limits_{i=1}^n \,(y_{i,j} - \mu_{i,j})^2}%
+{\sum\limits_{i=1}^n \Big(y_{i,j} - \frac{N_{i\mathstrut}}{N^{\mathstrut}}\!\! \sum\limits_{i'=1}^n \! y_{i'\!,j} \Big)^{\! 2}} $ \\
+\hline
+\multicolumn{2}{c}{} \\
+\multicolumn{2}{c}{$R^2$ where the residual sum-of-squares is centered so that the bias is subtracted:} \\
+\hline
+\multicolumn{1}{|l|}{\tt PLAIN\_R2\_NOBIAS${}_j \,\,= {}$} & \multicolumn{1}{l|}{\tt ADJUSTED\_R2\_NOBIAS${}_j \,\,= {}$} \\
+$ \displaystyle 1 -
+\frac{\sum\limits_{i=1}^n \Big(y_{i,j} \,{-}\, \mu_{i,j} \,{-}\, \frac{N_{i\mathstrut}}{N^{\mathstrut}}\!\!
+ \sum\limits_{i'=1}^n (y_{i'\!,j} \,{-}\, \mu_{i'\!,j}) \Big)^{\! 2}}%
+{\sum\limits_{i=1}^n \Big(y_{i,j} - \frac{N_{i\mathstrut}}{N^{\mathstrut}}\!\! \sum\limits_{i'=1}^n \! y_{i'\!,j} \Big)^{\! 2}} $ &
+$ \displaystyle 1 - {\textstyle\frac{N_{\mathstrut} - 1}{N^{\mathstrut} - m'}} \,
+\frac{\sum\limits_{i=1}^n \Big(y_{i,j} \,{-}\, \mu_{i,j} \,{-}\, \frac{N_{i\mathstrut}}{N^{\mathstrut}}\!\!
+ \sum\limits_{i'=1}^n (y_{i'\!,j} \,{-}\, \mu_{i'\!,j}) \Big)^{\! 2}}%
+{\sum\limits_{i=1}^n \Big(y_{i,j} - \frac{N_{i\mathstrut}}{N^{\mathstrut}}\!\! \sum\limits_{i'=1}^n \! y_{i'\!,j} \Big)^{\! 2}} $ \\
+\hline
+\end{tabular}}
+\caption{The $R^2$ statistics we compute in {\tt GLM-predict.dml}}
+\label{table:GLMpred:R2}
+\end{table}
+
+\smallskip
+\noindent{\bf Returns}
+\smallskip
+
+The matrix of predicted means (if the response is numerical) or probabilities (if the response
+is categorical), see ``Description'' subsection above for more information. Given {\tt Y}, we
+return some statistics in CSV format as described in Table~\ref{table:GLMpred:stats} and in the
+above text.
+
+\smallskip
+\noindent{\bf Examples}
+\smallskip
+
+Note that in the examples below the value for ``{\tt disp}'' input argument
+is set arbitrarily. The correct dispersion value should be computed from the training
+data during model estimation, or omitted if unknown (which sets it to~{\tt 1.0}).
+
+\smallskip\noindent
+Linear regression example:
+\par\hangindent=\parindent\noindent{\tt
+\hml -f GLM-predict.dml -nvargs dfam=1 vpow=0.0 link=1 lpow=1.0 disp=5.67
+ X=/user/biadmin/X.mtx B=/user/biadmin/B.mtx M=/user/biadmin/Means.mtx fmt=csv
+ Y=/user/biadmin/Y.mtx O=/user/biadmin/stats.csv
+
+}\smallskip\noindent
+Linear regression example, prediction only (no {\tt Y} given):
+\par\hangindent=\parindent\noindent{\tt
+\hml -f GLM-predict.dml -nvargs dfam=1 vpow=0.0 link=1 lpow=1.0
+ X=/user/biadmin/X.mtx B=/user/biadmin/B.mtx M=/user/biadmin/Means.mtx fmt=csv
+
+}\smallskip\noindent
+Binomial logistic regression example:
+\par\hangindent=\parindent\noindent{\tt
+\hml -f GLM-predict.dml -nvargs dfam=2 link=2 disp=3.0004464
+ X=/user/biadmin/X.mtx B=/user/biadmin/B.mtx M=/user/biadmin/Probabilities.mtx fmt=csv
+ Y=/user/biadmin/Y.mtx O=/user/biadmin/stats.csv
+
+}\smallskip\noindent
+Binomial probit regression example:
+\par\hangindent=\parindent\noindent{\tt
+\hml -f GLM-predict.dml -nvargs dfam=2 link=3 disp=3.0004464
+ X=/user/biadmin/X.mtx B=/user/biadmin/B.mtx M=/user/biadmin/Probabilities.mtx fmt=csv
+ Y=/user/biadmin/Y.mtx O=/user/biadmin/stats.csv
+
+}\smallskip\noindent
+Multinomial logistic regression example:
+\par\hangindent=\parindent\noindent{\tt
+\hml -f GLM-predict.dml -nvargs dfam=3
+ X=/user/biadmin/X.mtx B=/user/biadmin/B.mtx M=/user/biadmin/Probabilities.mtx fmt=csv
+ Y=/user/biadmin/Y.mtx O=/user/biadmin/stats.csv
+
+}\smallskip\noindent
+Poisson regression with the log link example:
+\par\hangindent=\parindent\noindent{\tt
+\hml -f GLM-predict.dml -nvargs dfam=1 vpow=1.0 link=1 lpow=0.0 disp=3.45
+ X=/user/biadmin/X.mtx B=/user/biadmin/B.mtx M=/user/biadmin/Means.mtx fmt=csv
+ Y=/user/biadmin/Y.mtx O=/user/biadmin/stats.csv
+
+}\smallskip\noindent
+Gamma regression with the inverse (reciprocal) link example:
+\par\hangindent=\parindent\noindent{\tt
+\hml -f GLM-predict.dml -nvargs dfam=1 vpow=2.0 link=1 lpow=-1.0 disp=1.99118
+ X=/user/biadmin/X.mtx B=/user/biadmin/B.mtx M=/user/biadmin/Means.mtx fmt=csv
+ Y=/user/biadmin/Y.mtx O=/user/biadmin/stats.csv
+
+}