text:
\documentclass{article}

\usepackage{hyperref}    % for making the table of contents clickable
\hypersetup{linktoc=all}

\usepackage{titlesec}     % for changing the format of my titles
\usepackage[titletoc]{appendix}    % for adding nice appendix title to the table of contents

\usepackage{graphicx}    % for handling figures in the document

\usepackage{tabularx}    % for handling tables with variable width columns and text wrapping
\usepackage{tabulary}     % for handling tables with variable width columns and text wrapping

\usepackage{array}
\newcolumntype{K}[1]{>{\raggedright\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}    %for having a single variable width column

\usepackage{booktabs}    % for making tables pretty
\usepackage{caption}

\usepackage[font=normalsize, skip=5pt]{caption}    % for customizing captions

\usepackage{amsmath}   % for math environments
\usepackage{amssymb}   % for special math symbols
\DeclareMathOperator*{\argmax}{arg\,max} %create an argmax operator

\usepackage{color, soul}    % for highlighting

\usepackage[round,
				    authoryear,
				    colon]{natbib}    %for the bibliography
\bibliographystyle{plainnat}

\usepackage{url}    %for handling of urls in bibtex file
\def\UrlBreaks{\do\/\do-}    %use to break urls at hyphens

\usepackage{pdflscape}    % for landscape orientation of pages

\usepackage[margin=1in]{geometry}    % to get 1-inch margins

\title{PyLogit: Maximum Likelihood Estimation for Conditional ``Logit-type'' Models in Python}
\author{Timothy Brathwaite}
\date{}

\begin{document}
\maketitle

\section{Introduction}
In many fields---including transportation, marketing, economics, finance, medicine, public health, and operations research---it is necessary to model the occurrence of unordered, discrete outcomes. Often, such outcomes represent an individual's choice from a set of alternatives. A simple way to model such an outcome is with what is ambiguously known as a ``logit model.'' Across the many fields mentioned above, ``logit models'' are also referred to as:
\begin{enumerate}
\item \textit{\textbf{binary logistic regression}} or simply \textit{\textbf{logistic regression}} when there are only two alternatives and all covariates are individual specific. These names are often used in statistics, computer science, medicine, and public health.

\item \textit{\textbf{binary logit}} when there are only two alternatives and covariates can vary across alternatives as well as across individuals. This name is often used in transportation, marketing, and economics.

\item \textit{\textbf{polychotomous regression}} or \textit{\textbf{multinomial logistic regression}} when there are two or more (typically more) alternatives and all covariates are individual specific. These names are often used in medicine and public health.

\item \textit{\textbf{softmax regression}} when there are two or more (typically more) alternatives and all covariates are individual specific. This name is often used in computer science and machine learning.

\item \textit{\textbf{conditional logit}} or \textit{\textbf{multinomial logit}} when there are two or more (typically more) alternatives and covariates may vary across alternatives as well as across individuals. These names are often used in transportation, marketing, and economics.
\end{enumerate}

Of the various logit model varieties listed above, conditional logit models are the most general. Polychotomous/softmax regression is a special case of conditional logit where the covariates only vary across individuals. Binary logit is a special case of conditional logit where there are only two alternatives. Lastly, (binary) logistic regression is a special case of conditional logit where there are only two alternatives and when covariates only vary across individuals.

In Python's scientific-computing ecosystem today, logit models are capable of being estimated by packages such as statsmodels and sci-kit-learn. However, these packages can only estimate binary logistic regression and polychotomous/softmax/multinomial regression where the covariates are individual specific. There are currently no Python packages on the Python Packaging Index (PyPi) that allow one to estimate conditional logit models. Note, the Python Biogeme software is able to estimate conditional logit models. However, Python Biogeme is not available on PyPi, so it cannot easily be installed. Moreover, Python Biogeme is largely a wrapper for routines written in C, and it does not fully accept all Python syntax.

Beyond Python, R is also heavily used for free and open-source scientific-computing. In R, the packages mlogit and mnlogit are capable of estimating conditional logit models. However, these packages cannot estimate conditional logit models when the choice set differs across individual. Additionally, these packages do not directly support the estimation of models where variables are associated with coefficients that vary only across subsets of alternatives. For instance, in a travel-mode choice context, mlogit and mnlogit do not directly support models where the travel time variable has one coefficient in the index of a drive-alone alternative and one coefficient for both the bus and train alternatives.

The pylogit package that is being introduced is capable of
\begin{itemize}
\item estimating conditional logit models (as well as other conditional ``logit-type'' models to be defined later)
\item estimating models using datasets where the choice set differs between individuals
\item estimating models with very general model specifications where the coefficients for a given variable may be
\begin{itemize}
\item completely alternative-specific (i.e. one coefficient per alternative, subject to identification of the coefficients),
\item subset-specific (i.e. one coefficient per subset of alternatives, where each alternative belongs to only one subset, and there are more than 1 but less than $J$ subsets, where $J$ is the maximum number of available alternatives in the dataset),
\item completely generic (i.e. one coefficient across all alternatives).
\end{itemize}
\end{itemize}

\section{Logit-type Models}
This section describes the type of models that are capable of being estimated with PyLogit, and in the process, it introduces some of the notation that will be used later in describing how the maximum likelihood estimation (MLE) is performed for this class of models.

``Logit-type'' models are models where the probability of individual $i$ choosing alternative $j$ (or being associated with outcome $j$) is given by:
\begin{equation}
\label{eqn:logit_type}
\begin{aligned}
P \left( y_{ij} = 1 | x_{ij}, \beta, \tau, \gamma \right) &= \frac{\exp \left[ S \left( V_{ij} \mid \tau _j, \gamma _j \right) \right]}{ \sum _{\ell \in C_i} \exp \left[ S \left( V_{i \ell} \mid \tau _{\ell}, \gamma _{\ell} \right) \right]}  = \frac{\exp \left( S_{ij} \right)}{ \sum _{\ell \in C_i} \exp \left( S_{i \ell} \right)}\\
\textrm{where } y_{ij} &= \textrm{a binary (0 or 1) indicator of whether individual $i$} \\
&\quad \  \textrm{is associated with outcome $j$.} \\
x_{ij} &= \textrm{a row vector} = r \left( z_j, \zeta _i \right)  \\
r \left( \cdot  \right) &= \textrm{a function that returns a row vector.} \\
z _j &= \textrm{attributes of alternative $j$ for individual $i$.} \\
\zeta _i &= \textrm{characteristics of individual $i$.} \\
\beta &= \textrm{a column vector of unknown population parameters.} \\
\tau &= \textrm{a 1-dimensional vector of constants, with one value}\\
&\quad \ \textrm{for each alternative in the dataset.}\\
\gamma &= \textrm{a 1-dimensional vector of ``shape parameters'', with}\\
&\quad \ \textrm{one value for each alternative in the dataset.}\\
\tau _j &= \textrm{a constant associated with alternative $j$.}\\
\gamma _j &= \textrm{a ``shape parameter'' associated with alternative $j$.} \\
V_{ij} &= x_{ij} \beta = \textrm{the \textit{index} for alternative $j$}\\
S \left( \cdot \right) &= \textrm{a model-specific function of $V_{ij}$, $\tau_j$ and $\gamma_j$.} \\
&\quad \ \textrm{It is monotonically increasing in $V_{ij}$.}\\
C_i &= \textrm{the choice set for individual $i$.}
\end{aligned}
\end{equation}
All models estimated by PyLogit have this form. In a standard conditional logit model, $S \left( V_{ij} \mid \tau _j, \gamma _j \right) = \tau _j + V_{ij}$ where $\tau_j$ is the alternative specific constant or intercept term for alternative $j$ and $V_{ij} $ does not contain a constant of its own.

\section{Computation}
In performing the MLE of logit-type models, PyLogit evaluates as many operations as it can through the use of matrix dot products. In particular, the calculation of the log-likelihood, the gradient of the log-likelihood with respect to the estimated parameters, and the hessian of the log-likelihood of log-likelihood with respect to the estimated parameters are all calculated through the use of one or more matrix dot products. Often, to be able to use dot products to perform the desired calculations, we require the use of ``mapping matrices'': matrices of zeros and ones that map one quantity to another related quantity. More will be said about the mapping matrices later.

\subsection{The log-likelihood}

Maximum likelihood estimation works as follows to find a desired vector of parameters $\hat{\theta}$:
\begin{equation}
\label{eqn:log_likelihood}
\begin{aligned}
\hat{\theta} &= \argmax_{\theta} \prod _{i = 1} ^N \prod _{j \in C_i} P_{ij} ^{y_{ij}} \\
&= \argmax_{\theta} \sum _i ^N \sum _{j \in C_i} y_{ij} \ln \left( P_{ij} \right) \\
&= \argmax_{\theta} Y^T \ln \left( P \right) \\
&= \argmax_{\theta} \mathcal{L} \left( \theta \right) \\
\textrm{where } \theta &= \left[ \hat{\gamma} ^T \mid \hat{\tau} ^T \mid \hat{\beta} ^T \right]^T = \textrm{a column vector created by stacking}\\
&\quad \ \textrm{ $\hat{\gamma}$ on top of $\hat{\tau}$ on top of $\hat{\beta}$. Note that terms with ``hats''} \\
&\quad \ \textrm{on them are estimates of their corresponding, unknown} \\
&\quad \ \textrm{quantities.} \\
P_{ij} &= \textrm{the probability of individual $i$ choosing alternative $j$} \\
&\quad \ \textrm{or being associated with outcome $j$. This probability is}\\
&\quad \ \textrm{given in Eqution \ref{eqn:logit_type}.}\\
Y &= \textrm{a column vector in $\mathbb{R}^{N_r}$, made by vertically stacking the} \\
&\quad \ \textrm{$y_{ij}$ for every available alternative for each individual.}\\
N_r &= \sum _i \parallel C_i \parallel \\
\parallel C_i \parallel &= \textrm{the size of the choice set for individual $i$.}\\
P &= \textrm{a column vector in $\mathbb{R}^{N_r}$, made by vertically stacking the} \\
&\quad \ \textrm{ $P_{ij}$ for every available alternative for every individual.}\\
\ln \left( P \right) &= \textrm{a column vector in $\mathbb{R}^{N_r}$, made by taking the natural} \\
&\quad \ \textrm{logarithm, element-wise, of all the elements in  $P$.}\\
\mathcal{L} \left( \theta \right) &= \textrm{the log-likelihood at the current value of $\theta$.} = Y^T \ln \left( P \right)
\end{aligned}
\end{equation}

To compute the maximum-likelihood estimate, $\hat{\theta}$, we therefore need at minimum the vectors $Y$ and $P$. While introduced above through the notion of vertically stacking each of the individual elements, $Y$ and $P$ are not typically obtained this way by PyLogit. First, PyLogit requires that all data used to estimate its logit-type models be in so-called ``long-format.'' Long-format data has one row for every available alternative for every individual. Moreover, the column used to indicate the choice or outcome for an individual is required to be a binary column---i.e. filled with zeros and ones. A ``$1$'' indicates that the alternative on the given row was chosen or associated with the corresponding individual for that row. In other words, $Y$ is available directly from the dataset provided to PyLogit.

A series of internal steps is used to efficiently calculate $P$. First, from the specification of the desired index coefficients, $\beta$, the indices those coefficients belong to, and the variables being multiplied by those coefficients, a design matrix $X$ is formed. Second, the vector $V$ of desired index values for each available alternative for each individual is formed through the dot product of $X$ and $\hat{\beta}$. An example of a set of desired index values is given in the \href{https://github.com/timothyb0912/pylogit/blob/master/examples/Main\%20PyLogit\%20Example.ipynb}{``Main PyLogit Example''} ($\leftarrow$ click me). Thirdly, the necessary transformations $S \left( V_{ij} \mid \tau _j, \gamma _j \right)$ are computed. As much as possible, such transformations are vectorized for speed using numpy and scipy. Finally, $P$ is formed through a set of elementwise operations and matrix dot products. Formally,

\begin{equation}
\label{eqn:index_and_probs}
\begin{aligned}
V &= X \hat{\beta} \quad \textrm{and } P = \frac{\exp \left( \vec{S} \right)}{\lambda \lambda^T \exp \left( \vec{S} \right)}\\
\textrm{where }X &\in \mathbb{R}^{N_r \times E}\\
X &= \textrm{The design matrix, formed by vertically stacking all} \\
&\quad \ \textrm{of the $X_{ij}$ on top of each other.} \\
X_{ij} &\in \mathbb{R}^{1 \times E} \\
N_r &= \textrm{the number of rows in X. Defined computationally in} \\
&\quad \ \textrm{Equation \ref{eqn:log_likelihood}.} \\
E &= \textrm{the number of explanatory variables in the model.}\\
&\quad \ \textrm{I.e. the number of index coefficients being estimated.} \\
\vec{S} &\in \mathbb{R}^{N_r \times 1}\\
\vec{S} &= \textrm{a column vector formed by stacking all of the} \\
&\quad \ \textrm{$S \left( V_{ij} \mid \hat{\tau _j}, \hat{\gamma _j} \right) = S_{ij}$ on top of each other.}\\
\exp \left( \vec{S} \right) &= \textrm{a column vector in $\mathbb{R}^{N_r}$, made by raising} \\
&\quad \ \textrm{the natural number $e$, element-wise, to the power of the}\\
&\quad \ \textrm{elements in  $\vec{S}$.} \\
\lambda &\in \mathbb{R}^{N_r \times N} \\
\lambda &= \textrm{a mapping matrix of zeros and ones that indicates the} \\
&\quad \ \textrm{observation (given on the columns) that each row of the} \\
&\quad \ \textrm{design matrix corresponds to.} \\
N &= \textrm{the number of observations in one's dataset.}
\end{aligned}
\end{equation}
Note that in the formula given above for $P$, the fraction indicates element-wise division.

\subsection{The Gradient, ${\partial \mathcal{L} \left( \hat{\theta} \right)} / {\partial \hat{\theta}}$}
The gradient of the log-likelihood with respect to our estimates $\hat{\theta}$ is used by many of the optimization algorithms in Scipy to numerically find a local optima of our log-likelihood function. Therefore, we need an efficient way of calculating the gradient.

The gradient of the log-likelihood function with respect to $\hat{\theta}$ is defined as follows:
\begin{equation*}
\begin{aligned}
\nabla \mathcal{L} \left( \hat{\theta} \right) &= \frac{\partial \mathcal{L} \left( \hat{\theta} \right)}{\partial \hat{\theta}} \\
&= \frac{\partial \left[ Y^T \ln \left( P \right) \right] }{\partial \hat{\theta}} \\
&= \frac{\partial}{\partial \hat{\theta}} \left[ \sum _i \sum _{j \in C_i} y_{ij} \ln \left( P_{ij} \right) \right], \quad i \in \left\lbrace 1, 2, ..., N \right\rbrace \\
&= \frac{\partial}{\partial \hat{\theta}} \left[ \sum _i \sum _{j \in C_i} y_{ij} \left( \ln \left[ \exp \left( S_{ij} \right) \right] - \ln \left[ \sum _{\ell \in C_i} \exp \left( S_{i \ell} \right) \right] \right) \right] \\
&= \frac{\partial}{\partial \hat{\theta}} \left[ \sum _i \sum _{j \in C_i} y_{ij} \left(  S_{ij} - \ln \left[ \sum _{\ell \in C_i} \exp \left( S_{i \ell} \right) \right] \right) \right] \\
&= \frac{\partial}{\partial \hat{\theta}} \left[ \sum _i \sum _{j \in C_i} y_{ij} S_{ij} - y_{ij} \ln \left( \sum _{\ell \in C_i} \exp \left[ S_{i \ell} \right] \right) \right] \\
&= \frac{\partial}{\partial \hat{\theta}} \left[ \sum _i \sum _{j \in C_i} y_{ij} S_{ij} -  \sum _i \sum _{j \in C_i} y_{ij} \ln \left( \sum _{\ell \in C_i} \exp \left[ S_{i \ell} \right] \right) \right] \\
&= \frac{\partial}{\partial \hat{\theta}} \left[ \sum _i \sum _{j \in C_i} y_{ij} S_{ij} -  \sum _i \ln \left( \sum _{\ell \in C_i} \exp \left[ S_{i \ell} \right] \right) \sum _{j \in C_i} y_{ij} \right] \\
&= \frac{\partial}{\partial \hat{\theta}} \left[ \sum _i \sum _{j \in C_i} y_{ij} S_{ij} -  \sum _i \ln \left( \sum _{\ell \in C_i} \exp \left[ S_{i \ell} \right] \right) \right] \\
&= \frac{\partial}{\partial \hat{\theta}} \left[ \sum _i  \left( \sum _{j \in C_i} y_{ij} S_{ij} \right) -  \ln \left( \sum _{\ell \in C_i} \exp \left[ S_{i \ell} \right] \right) \right] \\
&= \sum _i  \left[  \frac{\partial}{\partial \hat{\theta}} \left( \sum _{j \in C_i} y_{ij} S_{ij} \right) - \frac{\partial}{\partial \hat{\theta}} \ln \left( \sum _{\ell \in C_i} \exp \left[ S_{i \ell} \right] \right) \right] \\
&= \sum _i  \left[ \left( \sum _{j \in C_i} \frac{\partial}{\partial \hat{\theta}} \left[ y_{ij} S_{ij} \right] \right)  - \frac{\partial}{\partial \hat{\theta}} \ln \left( \sum _{\ell \in C_i} \exp \left[ S_{i \ell} \right] \right) \right] \\
&= \sum _i  \left[ \left( \sum _{j \in C_i} \frac{\partial}{\partial S_{ij}} \left[ y_{ij} S_{ij} \right] \frac{\partial S_{ij}}{\partial \hat{\theta}} \right)  - \frac{\partial}{\partial \hat{\theta}} \ln \left( \sum _{\ell \in C_i} \exp \left[ S_{i \ell} \right] \right) \right] \\
&= \sum _i  \left[ \left( \sum _{j \in C_i}  y_{ij} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right)  - \frac{\partial}{\partial \hat{\theta}} \ln \left( \sum _{\ell \in C_i} \exp \left[ S_{i \ell} \right] \right) \right] \\
&= \sum _i  \left[ \left( \sum _{j \in C_i}  y_{ij} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right)  - \frac{\partial}{\partial \varphi} \left[ \ln \left( \varphi \right) \right] \frac{\partial \varphi}{\partial \hat{\theta}} \right]  \quad \textrm{where $\varphi =  \sum _{\ell \in C_i} \exp \left[ S_{i \ell} \right] $}\\
&= \sum _i  \left[ \left( \sum _{j \in C_i}  y_{ij} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right)  - \frac{1}{\varphi} \frac{\partial \varphi}{\partial \hat{\theta}} \right]
\end{aligned}
\end{equation*}

\begin{equation}
\label{eqn:gradient_log_likelihood}
\begin{aligned}
\nabla \mathcal{L} \left( \hat{\theta} \right) &= \sum _i  \left[ \left( \sum _{j \in C_i}  y_{ij} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right)  - \frac{1}{\varphi} \frac{\partial}{\partial \hat{\theta}} \left( \sum _{j \in C_i} \exp \left[ S_{ij} \right] \right) \right] \\
&= \sum _i  \left[ \left( \sum _{j \in C_i}  y_{ij} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right)  - \frac{1}{\varphi} \left( \sum _{j \in C_i} \frac{\partial}{\partial \hat{\theta}} \exp \left[ S_{i j} \right] \right) \right] \\
&= \sum _i  \left[ \left( \sum _{j \in C_i}  y_{ij} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right)  - \frac{1}{\varphi} \left( \sum _{j \in C_i} \frac{\partial}{\partial S_{ij}} \exp \left[ S_{i j} \right] \frac{\partial S_{ij}}{\partial \hat{\theta}} \right) \right] \\
&= \sum _i  \left[ \left( \sum _{j \in C_i}  y_{ij} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right)  - \frac{1}{\varphi} \left( \sum _{j \in C_i} \exp \left[ S_{i j} \right] \frac{\partial S_{ij}}{\partial \hat{\theta}} \right) \right] \\
&= \sum _i  \left[ \left( \sum _{j \in C_i}  y_{ij} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right)  - \left( \sum _{j \in C_i} \frac{\exp \left[ S_{i j} \right]}{\varphi} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right) \right] \\
&= \sum _i  \left[ \left( \sum _{j \in C_i}  y_{ij} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right)  - \left( \sum _{j \in C_i} \frac{\exp \left[ S_{i j} \right]}{\sum _{\ell \in C_i} \exp \left[ S_{i \ell} \right]} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right) \right] \\
&= \sum _i  \left[ \left( \sum _{j \in C_i}  y_{ij} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right)  - \left( \sum _{j \in C_i} P_{ij} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right) \right] \\
&= \sum _i  \left[ \left( \sum _{j \in C_i}  y_{ij} \frac{\partial S_{ij}}{\partial \hat{\theta}} - P_{ij} \frac{\partial S_{ij}}{\partial \hat{\theta}} \right) \right] \\
&= \sum _i  \sum _{j \in C_i}  \left[ \left( y_{ij}  - P_{ij} \right) \frac{\partial S_{ij}}{\partial \hat{\theta}} \right] \\
&= \left( Y - P \right)^T \frac{\partial \vec{S}}{\partial \hat{\theta}} \\
&= \left( Y - P \right)^T \left[ \frac{\partial \vec{S}}{\partial \hat{\gamma}} \mid \frac{\partial \vec{S}}{\partial \hat{\tau}} \mid \frac{\partial \vec{S}}{\partial \hat{\beta}} \right] \\
&= \left( Y - P \right)^T \left[ \frac{\partial \vec{S}}{\partial \hat{\gamma}} \mid \frac{\partial \vec{S}}{\partial \hat{\tau}} \mid \frac{\partial \vec{S}}{\partial V} \frac{\partial V}{\partial \hat{\beta}} \right] \\
&= \left( Y - P \right)^T \left[ \frac{\partial \vec{S}}{\partial \hat{\gamma}} \mid \frac{\partial \vec{S}}{\partial \hat{\tau}} \mid \frac{\partial \vec{S}}{\partial V} X \right] \\
&= \left( Y - P \right)^T \left[ \frac{\partial \vec{S}}{\partial \hat{\gamma}} \mid \xi^{(-1)} \mid \frac{\partial \vec{S}}{\partial V} X \right] \\
\textrm{where } \xi &= \textrm{a mapping matrix of zeros and ones that indicates the} \\
&\quad \ \textrm{alternative (given on the columns) that each row of the} \\
&\quad \ \textrm{design matrix corresponds to.} \\
\xi^{(-1)} &= \textrm{the mapping matrix $\xi$, without the column that corresponds}\\
&\quad \ \textrm{to the alternative whose intercept, $\tau _j$, is not being estimated.}
\end{aligned}
\end{equation}
Note that by the chain rule, $\dfrac{\partial \mathcal{L} \left( \hat{\theta} \right)}{\partial \hat{\theta}} = \dfrac{\partial \mathcal{L} \left( \hat{\theta} \right)}{\partial \vec{S}} \dfrac{\partial \vec{S}}{\partial \hat{\theta}}$. From above, this implies that $ \dfrac{\partial \mathcal{L} \left( \hat{\theta} \right)}{\partial \vec{S}} = \left( Y - P \right)^T$. \newpage

In PyLogit, the calculation of the gradient is based on Equation \ref{eqn:gradient_log_likelihood}. The vectors $Y$ and $P$ are readily available based on the provided data and Equation \ref{eqn:index_and_probs}. The matrices $\dfrac{\partial \vec{S}}{\partial \hat{\gamma}}$ and $\dfrac{\partial \vec{S}}{\partial V}$ differ from model to model, and custom functions to calculate them must\footnote{Automatic Differentiation or Finite Differences could be used to avoid having to manually write such functions, but the estimation times would increase greatly as a result.} be written for each model that PyLogit supports. $\xi^{(-1)}$ differs from specification to specification, but it can always be easily calculated once per model estimation. Note that in PyLogit, all of the mapping matrices are implemented as Scipy sparse matrices for fast matrix dot products.

\subsection{The Hessian, $\dfrac{\partial ^2 \mathcal{L} \left( \hat{\theta} \right) }{{\partial \hat{\theta}}^2}$}
Although it is not required by all of the optimizers in scipy.optimize, being able to calculate the Hessian in closed-form is useful for optimizing non-concave log-likelihood functions. Below, the form of the Hessian used by PyLogit is derived for logit-type models.

The hessian is given by:
\begin{equation}
\label{eqn:partitioned_hessian}
\begin{aligned}
\textrm{Hessian} = \frac{\partial ^2 \mathcal{L} \left( \hat{\theta} \right) }{ {\partial \hat{\theta}}^2 } &= \frac{\partial}{\partial \hat{\theta}} \left[ \frac{\partial \left( Y^T \ln \left[ P \right] \right) }{\partial \hat{\theta}} \right] \\
&= \frac{\partial}{\partial \hat{\theta}} \left[ \frac{\partial}{\partial \hat{\gamma}}  \left( Y^T \ln \left[ P \right] \right)   \mid \frac{\partial}{\partial \hat{\tau}}  \left( Y^T \ln \left[ P \right] \right) \mid \frac{\partial}{\partial \hat{\beta}}  \left( Y^T \ln \left[ P \right] \right) \right]
\\
&= \left[ \begin{array}{c}
\frac{\partial}{\partial \hat{\theta}} \left( \frac{\partial}{\partial \hat{\gamma}}  \left[ Y^T \ln \left( P \right) \right] \right)
\\[1.2ex]
\frac{\partial}{\partial \hat{\theta}} \left( \frac{\partial}{\partial \hat{\tau}}  \left[ Y^T \ln \left( P \right) \right] \right)
\\[1.2ex]
\frac{\partial}{\partial \hat{\theta}} \left( \frac{\partial}{\partial \hat{\beta}} \left[ Y^T \ln \left( P \right) \right] \right)
\end{array}  \right]
\\
&= \left[ \begin{array}{c|c|c}
\frac{\partial}{\partial \hat{\gamma}}  \left( \frac{\partial}{\partial \hat{\gamma}}  \left[ Y^T \ln \left( P \right) \right] \right)  &

\frac{\partial}{\partial \hat{\tau}}  \left( \frac{\partial}{\partial \hat{\gamma}}  \left[ Y^T \ln \left( P \right) \right] \right) &

\frac{\partial}{\partial \hat{\beta}}  \left( \frac{\partial}{\partial \hat{\gamma}}  \left[ Y^T \ln \left( P \right) \right] \right) \\[1.2ex]

 \frac{\partial}{\partial \hat{\gamma}} \left( \frac{\partial}{\partial \hat{\tau}}  \left[ Y^T \ln \left( P \right) \right] \right) &

\frac{\partial}{\partial \hat{\tau}} \left( \frac{ \partial }{ \partial \hat{\tau} }  \left[ Y^T \ln \left( P \right) \right] \right) &

 \frac{\partial}{\partial \hat{\beta}} \left( \frac{\partial}{\partial \hat{\tau}}  \left[ Y^T \ln \left( P \right) \right] \right) \\[1.2ex]

 \frac{\partial}{\partial \hat{\gamma}} \left( \frac{\partial}{\partial \hat{\beta}} \left[ Y^T \ln \left( P \right) \right] \right) &

 \frac{\partial}{\partial \hat{\tau}} \left( \frac{\partial}{\partial \hat{\beta}} \left[ Y^T \ln \left( P \right) \right] \right) &

\frac{\partial}{\partial \hat{\beta}} \left( \frac{\partial}{\partial \hat{\beta}} \left[ Y^T \ln \left( P \right) \right] \right)
\end{array}  \right] \\
&= \left[ \begin{array}{c|c|c}
H_{11} & H_{12} & H_{13} \\
H_{21} & H_{22} & H_{23} \\
H_{31} & H_{32} & H_{33}
\end{array}  \right]
\end{aligned}
\end{equation}
Based on Equation \ref{eqn:partitioned_hessian}, we can see that the Hessian is actually a partitioned matrix with nine sub-matrices. In the cases where either $\gamma$ and/or $\tau$ do not exist in one's model, the Hessian matrix can be derived from Equation \ref{eqn:partitioned_hessian} by deleting the rows and columns that contain partial derivatives with respect to the parameters that are not present in $\theta$. In the simplest case, where there are no shape parameters ($\gamma$) and no intercept parameters ($\tau$) outside of the index, the Hessian reduces to the usual expression for conditional logit models: $\dfrac{\partial}{\partial \hat{\beta}} \left( \dfrac{\partial}{\partial \hat{\beta}} \left[ Y^T \ln \left( P \right) \right] \right)$.

Below, we will derive expressions for the nine sub-matrices $H_{mn}$ where $m, n \in \left\lbrace 1, 2, 3 \right\rbrace$. Note however, that we do not actually need to derive all nine sub-matrices because mixed partial derivatives of vectors are transposes of each other. Stated more precisely, $\dfrac{\partial}{\partial A} \left[ \dfrac{\partial \alpha \left(A, B \right)}{\partial B} \right] = \left(\dfrac{\partial}{\partial B} \left[ \dfrac{\partial \alpha \left(A, B \right)}{\partial A} \right] \right)^T $ for a function $\alpha \left( \cdot \right)$ that maps input vectors $A$ and $B$ to scalar outputs.

\subsubsection{$H_{11}$}
\begin{equation}
\begin{aligned}
H_{11} &= \frac{\partial}{\partial \hat{\gamma}}  \left( \frac{\partial}{\partial \hat{\gamma}}  \left[ Y^T \ln \left( P \right) \right] \right)
\\
&= \frac{\partial}{\partial \hat{\gamma}}  \left( \frac{\partial}{\partial \vec{S}}  \left[ Y^T \ln \left( P \right) \right] \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right)
\\
&= \frac{\partial}{\partial \hat{\gamma}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right)
\\
&= \frac{\partial}{\partial \vec{S}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right) \frac{\partial \vec{S}}{\partial \hat{\gamma}}
\\
&= \frac{\partial}{\partial \left( Y - P \right)}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right) \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\gamma}}
\\
&= \left( \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right)^T \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\gamma}}
\\
&= - \left( \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right)^T \frac{\partial P}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\gamma}}
\end{aligned}
\end{equation}

\subsubsection{$H_{12}$}
\begin{equation}
\begin{aligned}
H_{12} &= \frac{\partial}{\partial \hat{\tau}}  \left( \frac{\partial}{\partial \hat{\gamma}}  \left[ Y^T \ln \left( P \right) \right] \right)
\\
&= \frac{\partial}{\partial \hat{\tau}}  \left( \frac{\partial}{\partial \vec{S}}  \left[ Y^T \ln \left( P \right) \right] \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right)
\\
&= \frac{\partial}{\partial \hat{\tau}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right)
\\
&= \frac{\partial}{\partial \vec{S}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right) \frac{\partial \vec{S}}{\partial \hat{\tau}}
\\
&= \frac{\partial}{\partial \left( Y - P \right)}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right) \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\tau}}
\\
&= \left( \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right)^T \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\tau}}
\\
&= - \left( \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right)^T \frac{\partial P}{\partial \vec{S}} \xi^{(-1)}
\end{aligned}
\end{equation}

\subsubsection{$H_{13}$}
\begin{equation}
\begin{aligned}
H_{13} &= \frac{\partial}{\partial \hat{\beta}}  \left( \frac{\partial}{\partial \hat{\gamma}}  \left[ Y^T \ln \left( P \right) \right] \right)
\\
&= \frac{\partial}{\partial \hat{\beta}}  \left( \frac{\partial}{\partial \vec{S}}  \left[ Y^T \ln \left( P \right) \right] \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right)
\\
&= \frac{\partial}{\partial \hat{\beta}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right)
\\
&= \frac{\partial}{\partial \vec{S}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right) \frac{\partial \vec{S}}{\partial \hat{\beta}}
\\
&= \frac{\partial}{\partial \left( Y - P \right)}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right) \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\beta}}
\\
&= \left( \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right)^T \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\beta}}
\\
&= - \left( \frac{\partial \vec{S}}{\partial \hat{\gamma}} \right)^T \frac{\partial P}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial V} X
\end{aligned}
\end{equation}

\subsubsection{$H_{21}$}
\begin{equation}
\begin{aligned}
H_{21} &= \frac{\partial}{\partial \hat{\gamma}}  \left( \frac{\partial}{\partial \hat{\tau}}  \left[ Y^T \ln \left( P \right) \right] \right)
\\
&= \frac{\partial}{\partial \hat{\gamma}}  \left( \frac{\partial}{\partial \vec{S}}  \left[ Y^T \ln \left( P \right) \right] \frac{\partial \vec{S}}{\partial \hat{\tau}} \right)
\\
&= \frac{\partial}{\partial \hat{\gamma}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\tau}} \right)
\\
&= \frac{\partial}{\partial \vec{S}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\tau}} \right) \frac{\partial \vec{S}}{\partial \hat{\gamma}}
\\
&= \frac{\partial}{\partial \left( Y - P \right)}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\tau}} \right) \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\gamma}}
\\
&= \left( \frac{\partial \vec{S}}{\partial \hat{\tau}} \right)^T \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\gamma}}
\\
&= - \left( \xi^{(-1)} \right)^T \frac{\partial P}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\gamma}}
\end{aligned}
\end{equation}

\subsubsection{$H_{22}$}
\begin{equation}
\begin{aligned}
H_{22} &= \frac{\partial}{\partial \hat{\tau}}  \left( \frac{\partial}{\partial \hat{\tau}}  \left[ Y^T \ln \left( P \right) \right] \right)
\\
&= \frac{\partial}{\partial \hat{\tau}}  \left( \frac{\partial}{\partial \vec{S}}  \left[ Y^T \ln \left( P \right) \right] \frac{\partial \vec{S}}{\partial \hat{\tau}} \right)
\\
&= \frac{\partial}{\partial \hat{\tau}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\tau}} \right)
\\
&= \frac{\partial}{\partial \vec{S}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\tau}} \right) \frac{\partial \vec{S}}{\partial \hat{\tau}}
\\
&= \frac{\partial}{\partial \left( Y - P \right)}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\tau}} \right) \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\tau}}
\\
&= \left( \frac{\partial \vec{S}}{\partial \hat{\tau}} \right)^T \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\tau}}
\\
&= - \left( \xi^{(-1)} \right)^T \frac{\partial P}{\partial \vec{S}} \xi^{(-1)}
\end{aligned}
\end{equation}

\subsubsection{$H_{23}$}
\begin{equation}
\begin{aligned}
H_{23} &= \frac{\partial}{\partial \hat{\beta}}  \left( \frac{\partial}{\partial \hat{\tau}}  \left[ Y^T \ln \left( P \right) \right] \right)
\\
&= \frac{\partial}{\partial \hat{\beta}}  \left( \frac{\partial}{\partial \vec{S}}  \left[ Y^T \ln \left( P \right) \right] \frac{\partial \vec{S}}{\partial \hat{\tau}} \right)
\\
&= \frac{\partial}{\partial \hat{\beta}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\tau}} \right)
\\
&= \frac{\partial}{\partial \vec{S}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\tau}} \right) \frac{\partial \vec{S}}{\partial \hat{\beta}}
\\
&= \frac{\partial}{\partial \left( Y - P \right)}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\tau}} \right) \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\beta}}
\\
&= \left( \frac{\partial \vec{S}}{\partial \hat{\tau}} \right)^T \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\beta}}
\\
&= - \left( \xi^{(-1)} \right)^T \frac{\partial P}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\beta}} \\
&= - \left( \xi^{(-1)} \right)^T \frac{\partial P}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial V} X
\end{aligned}
\end{equation}

\subsubsection{$H_{31}$}
\begin{equation}
\begin{aligned}
H_{31} &= \frac{\partial}{\partial \hat{\gamma}}  \left( \frac{\partial}{\partial \hat{\beta}}  \left[ Y^T \ln \left( P \right) \right] \right)
\\
&= \frac{\partial}{\partial \hat{\gamma}}  \left( \frac{\partial}{\partial \vec{S}}  \left[ Y^T \ln \left( P \right) \right] \frac{\partial \vec{S}}{\partial \hat{\beta}} \right)
\\
&= \frac{\partial}{\partial \hat{\gamma}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\beta}} \right)
\\
&= \frac{\partial}{\partial \vec{S}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\beta}} \right) \frac{\partial \vec{S}}{\partial \hat{\gamma}}
\\
&= \frac{\partial}{\partial \left( Y - P \right)}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\beta}} \right) \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\gamma}}
\\
&= \left( \frac{\partial \vec{S}}{\partial \hat{\beta}} \right)^T \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\gamma}}
\\
&= - \left( \frac{\partial \vec{S}}{\partial \hat{\beta}} \right)^T \frac{\partial P}{\partial \vec{S}}  \frac{\partial \vec{S}}{\partial \hat{\gamma}}
\\
&= - \left( \frac{\partial \vec{S}}{\partial V} X \right)^T \frac{\partial P}{\partial \vec{S}}  \frac{\partial \vec{S}}{\partial \hat{\gamma}}
 \\
&= - X^T  \left( \frac{\partial \vec{S}}{\partial V} \right)^T \frac{\partial P}{\partial \vec{S}}  \frac{\partial \vec{S}}{\partial \hat{\gamma}}
\end{aligned}
\end{equation}

\subsubsection{$H_{32}$}
\begin{equation}
\begin{aligned}
H_{32} &= \frac{\partial}{\partial \hat{\tau}}  \left( \frac{\partial}{\partial \hat{\beta}}  \left[ Y^T \ln \left( P \right) \right] \right)
\\
&= \frac{\partial}{\partial \hat{\tau}}  \left( \frac{\partial}{\partial \vec{S}}  \left[ Y^T \ln \left( P \right) \right] \frac{\partial \vec{S}}{\partial \hat{\beta}} \right)
\\
&= \frac{\partial}{\partial \hat{\tau}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\beta}} \right)
\\
&= \frac{\partial}{\partial \vec{S}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\beta}} \right) \frac{\partial \vec{S}}{\partial \hat{\tau}}
\\
&= \frac{\partial}{\partial \left( Y - P \right)}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\beta}} \right) \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\tau}}
\\
&= \left( \frac{\partial \vec{S}}{\partial \hat{\beta}} \right)^T \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\tau}}
\\
&= - \left( \frac{\partial \vec{S}}{\partial \hat{\beta}} \right)^T \frac{\partial P}{\partial \vec{S}} \xi^{(-1)}  \\
&= - \left( \frac{\partial \vec{S}}{\partial V} X  \right)^T \frac{\partial P}{\partial \vec{S}} \xi^{(-1)} \\
&= - X^T \left( \frac{\partial \vec{S}}{\partial V} \right)^T \frac{\partial P}{\partial \vec{S}} \xi^{(-1)}
\end{aligned}
\end{equation}

\subsubsection{$H_{33}$}
\begin{equation}
\begin{aligned}
H_{33} &= \frac{\partial}{\partial \hat{\beta}}  \left( \frac{\partial}{\partial \hat{\beta}}  \left[ Y^T \ln \left( P \right) \right] \right)
\\
&= \frac{\partial}{\partial \hat{\beta}}  \left( \frac{\partial}{\partial \vec{S}}  \left[ Y^T \ln \left( P \right) \right] \frac{\partial \vec{S}}{\partial \hat{\beta}} \right)
\\
&= \frac{\partial}{\partial \hat{\beta}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\beta}} \right)
\\
&= \frac{\partial}{\partial \vec{S}}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\beta}} \right) \frac{\partial \vec{S}}{\partial \hat{\beta}}
\\
&= \frac{\partial}{\partial \left( Y - P \right)}  \left( \left[ Y - P \right]^T \frac{\partial \vec{S}}{\partial \hat{\beta}} \right) \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\beta}}
\\
&= \left( \frac{\partial \vec{S}}{\partial \hat{\beta}} \right)^T \frac{\partial \left( Y - P \right)}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial \hat{\beta}}
\\
&= - \left( \frac{\partial \vec{S}}{\partial V} X \right)^T \frac{\partial P}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial V} X \\
&= - X^T \left( \frac{\partial \vec{S}}{\partial V} \right)^T \frac{\partial P}{\partial \vec{S}} \frac{\partial \vec{S}}{\partial V} X
\end{aligned}
\end{equation}

\subsubsection{$\dfrac{\partial P}{\partial \vec{S}}$}
The last nine sub-sections showed that it is possible to express all of the sub-matrices of the Hessian in terms of dot products of matrices that were used to form the gradient and one unknown matrix $\dfrac{\partial P}{\partial \vec{S}}$. Below, we derive an expression for $\dfrac{\partial P}{\partial \vec{S}}$, thereby allowing us to compute a closed-form solution for the Hessian of the logit-type models used in PyLogit.

First, note that we can write:
\begin{equation*}
\begin{aligned}
P &= \left[ \begin{array}{c}
P_1 \\
P_2 \\
 ... \\
 P_N
\end{array} \right] \\
\textrm{where } P_i &= \left[ \begin{array}{c}
P_{i1} \\
P_{i2} \\
 ... \\
 P_{ij}
\end{array} \right], \quad j \in C_i
\end{aligned}
\end{equation*}
It then follows that
\begin{equation}
\frac{\partial P}{\partial \vec{S}} = \left[ \begin{array}{c}
\frac{\partial P_1}{\partial \vec{S}} \\[1.2ex]
\frac{\partial P_2}{\partial \vec{S}} \\
... \\
\frac{\partial P_N}{\partial \vec{S}}
\end{array} \right]
\end{equation}
A generic partition, $\dfrac{\partial P_i}{\partial \vec{S}}$ of $\dfrac{\partial P}{\partial \vec{S}}$ can further be decomposed as follows:
\begin{equation}
\begin{aligned}
\frac{\partial P_i}{\partial \vec{S}} &= \left[ \frac{\partial P_i}{\partial \vec{S}_1} \mid \frac{\partial P_i}{\partial \vec{S}_2} \mid ... \mid \frac{\partial P_i}{\partial \vec{S}_N} \right] \\
\textrm{where } \frac{\partial P_i}{\partial \vec{S}_{i'}} &= \vec{0}_{i, i'}, \quad \forall i' \neq i
\\
\vec{0}_{i, i'} &\in \mathbb{R}^{ \parallel C_i \parallel \times \parallel C_{i'} \parallel}
\\
i, i' &\in \left\lbrace 1, 2, ... N \right\rbrace
\\
\vec{S}_i &= \left[ S_{i1} \mid S_{i2} \mid ... \mid S_{ij} \right]^T, \quad j \in C_i
\end{aligned}
\end{equation}
Because $\dfrac{\partial P_i}{\partial \vec{S}_{i'}} = \vec{0}_{i, i'}, \  \forall i' \neq i$, we can focus on computing $\dfrac{\partial P_i}{\partial \vec{S}_{i}}$.
\begin{equation}
\frac{\partial P_i}{\partial \vec{S}_{i}} = \left[ \begin{array}{cccc}
\frac{\partial P_{i1}}{\partial S_{i1}} & \frac{\partial P_{i1}}{\partial S_{i2}} & ... & \frac{\partial P_{i1}}{\partial S_{ij}}
\\[1.2ex]
\frac{\partial P_{i2}}{\partial S_{i1}} & \frac{\partial P_{i2}}{\partial S_{i2}} & ... & \frac{\partial P_{i2}}{\partial S_{ij}}
\\
... & ... & ... & ...
\\[1.2ex]
\frac{\partial P_{ij}}{\partial S_{i1}} & \frac{\partial P_{ij}}{\partial S_{i2}} & ... & \frac{\partial P_{ij}}{\partial S_{ij}}
\end{array} \right] \quad j \in C_i
\end{equation}
The individual elements of $\dfrac{\partial P_i}{\partial \vec{S}_{i}}$ can be computed as follows:
\begin{equation*}
\begin{aligned}
\frac{\partial P_{i \ell}}{\partial S_{i \ell}} &= \frac{\partial}{\partial S_{i \ell}} \left[ \frac{ \exp \left( S_{i \ell} \right) }{\sum _{k \in C_i} \exp \left( S_{i k} \right) } \right]
\\
&= \frac{ \left[ \sum _{k \in C_i} \exp \left( S_{i k} \right) \right] \exp \left( S_{i \ell} \right) - \exp \left( S_{i \ell} \right) \exp \left( S_{i \ell} \right) }{ \left[ \sum _{k \in C_i} \exp \left( S_{i k} \right) \right]^2 }
 \\
&= P_{i \ell} - \left( P_{i \ell} \right)^2
\\
\frac{\partial P_{i \ell}}{\partial S_{i \ell'}} &= \frac{\partial}{\partial S_{i \ell'}} \left[ \frac{ \exp \left( S_{i \ell} \right) }{\sum _{k \in C_i} \exp \left( S_{i k} \right) } \right] \quad \forall \ell' \neq \ell
\\
&=  \frac{ \left[ \sum _{k \in C_i} \exp \left( S_{i k} \right) \right] * 0 - \exp \left( S_{i \ell} \right) \exp \left( S_{i \ell'} \right) }{ \left[ \sum _{k \in C_i} \exp \left( S_{i k} \right) \right]^2 }
\\
&= -P_{i \ell} P_{i \ell'}
\end{aligned}
\end{equation*}
With these results in hand, $\dfrac{\partial P_i}{\partial \vec{S}_{i}}$ can be rewritten as follows:
\begin{equation}
\label{eqn:dp_ds_submatrices}
\begin{aligned}
\frac{\partial P_i}{\partial \vec{S}_{i}} &= \left[ \begin{array}{cccc}
\frac{\partial P_{i1}}{\partial S_{i1}} & \frac{\partial P_{i1}}{\partial S_{i2}} & ... & \frac{\partial P_{i1}}{\partial S_{ij}}
\\[1.2ex]
\frac{\partial P_{i2}}{\partial S_{i1}} & \frac{\partial P_{i2}}{\partial S_{i2}} & ... & \frac{\partial P_{i2}}{\partial S_{ij}}
\\
... & ... & ... & ...
\\[1.2ex]
\frac{\partial P_{ij}}{\partial S_{i1}} & \frac{\partial P_{ij}}{\partial S_{i2}} & ... & \frac{\partial P_{ij}}{\partial S_{ij}}
\end{array} \right], \quad j \in C_i
\\
&= \left[ \begin{array}{cccc}
P_{i 1} - \left( P_{i 1} \right)^2 &-P_{i 1} P_{i 2} & ... & -P_{i 1} P_{i j}
\\[1.2ex]
-P_{i 1} P_{i 2}& P_{i 2} - \left( P_{i 2} \right)^2 & ... & -P_{i 2} P_{i j}
\\
... & ... & ... & ...
\\[1.2ex]
-P_{i 1} P_{i j} & -P_{i 2} P_{i j} & ... & P_{i j} - \left( P_{i j} \right)^2
\end{array} \right]
\\
&= \textrm{diag} \left( P_i \right) - P_i {P_i}^T
\\
\textrm{where diag} \left( \cdot \right) &= \textrm{a diagonal matrix with the entries of the argument}
\\
&\quad \ \textrm{on the main diagonal.}
\end{aligned}
\end{equation}
Returning finally to $\dfrac{\partial P}{\partial \vec{S}}$ we can convince ourselves that $\dfrac{\partial P}{\partial \vec{S}}$ is a block diagonal matrix formed by the submatrices $\dfrac{\partial P_i}{\partial \vec{S}_i}$ for $i \in \left\lbrace 1, 2, ..., N \right\rbrace$.
\begin{equation}
\label{eqn:dp_ds}
\begin{aligned}
\frac{\partial P}{\partial \vec{S}} &= \left[ \begin{array}{c}
\frac{\partial P_1}{\partial \vec{S}} \\[1.2ex]
\frac{\partial P_2}{\partial \vec{S}} \\
... \\
\frac{\partial P_N}{\partial \vec{S}}
\end{array} \right]
\\
&= \left[ \begin{array}{cccc}
\frac{\partial P_1}{\partial \vec{S}_1} & \vec{0}_{1, 2}  & \cdots &  \vec{0}_{1, N}\\[1.2ex]
\vec{0}_{2, 1} & \frac{\partial P_2}{\partial \vec{S}_2}  & \cdots &  \vec{0}_{2, N}\\
\cdots & \cdots & \ddots & \cdots \\
\vec{0}_{N, 1} & \vec{0}_{N, 1} & \cdots & \frac{\partial P_N}{\partial \vec{S}_N}
\end{array} \right]
\end{aligned}
\end{equation}

As written above, we only have a vectorized implementation of $\dfrac{\partial P_i}{\partial \vec{S}_{i}}$. To compute $\dfrac{\partial P}{\partial \vec{S}}$, a for-loop would need to be employed to create each $\dfrac{\partial P_i}{\partial \vec{S}_{i}}$ and to assemble them block-diagonally. Aside from being computationally inefficient, this would also be memory inefficient since $\dfrac{\partial P}{\partial \vec{S}}$ will have shape $N_r \times N_r$.

To avoid these difficulties, we'll first extend the vectorized implementation of $\dfrac{\partial P_i}{\partial \vec{S}_{i}}$ from individuals to the entire dataset. In particular, to form $\dfrac{\partial P}{\partial \vec{S}}$, we use
\begin{equation}
\label{eq:dp-ds-vectorized}
\begin{aligned}
\frac{\partial P}{\partial \vec{S}} &= \textrm{diag} \left( P \right) - \left[ \textrm{diag} \left( P \right) \lambda \right] \left[ \textrm{diag} \left( P \right) \lambda \right]^T\\
\textrm{where } \lambda &\in \mathbb{R}^{N_r \times N} \\
\lambda &= \textrm{a mapping matrix of zeros and ones that indicates the} \\
&\quad \ \textrm{observation (given on the columns) that each row of the} \\
&\quad \ \textrm{design matrix corresponds to.} \\
N_r &= \textrm{The number of rows in $X$.}\\
N &= \textrm{The number of decision makers in the dataset.}
\end{aligned}
\end{equation}

Next, we note that matrix multiplication is associative, so we can calculate the hessian without explicitly calculating $\dfrac{\partial P}{\partial \vec{S}}$. To do so, we distribute all matrix multiplication operations across the terms in Equation \ref{eq:dp-ds-vectorized}: $\textrm{diag} \left( P \right)$, $\left[ \textrm{diag} \left( P \right) \lambda \right]$, and $\left[ \textrm{diag} \left( P \right) \lambda \right]^T$.
\end{document}

max_stars_repo_name: mathijsvdv/pylogit

max_issues_repo_name: mathijsvdv/pylogit

max_forks_repo_name: mathijsvdv/pylogit


text:
% --- [ Evaluation Process ] ---------------------------------------------------

\subsection{Evaluation Process}
\label{sec:evaluation_process}

% TODO: rephrase

A preliminary metric for the effectiveness of control flow recovery is how many of the original high-level control flow primitives that are recovered and how many are omitted. The metric for effectiveness should also includes a notion of false positives for high-level primitives that are recovered but are not part of the original source code.

For each test program, the evaluation process consists of the following steps.

\begin{enumerate}
	\item \label{itm:step_c_to_json} Parse the original C source code of the test program to determine the occurrence of each high-level control flow primitive, and store this truth table in JSON format (which serves as a \textit{ground truth} during evaluation).
	\item \label{itm:step_c_to_ll} Convert C source code to LLVM IR using Clang. For test programs consisting of several source files, the LLVM IR is produced using the WLLVM\footnote{Whole Program LLVM: \url{https://github.com/travitch/whole-program-llvm}} project which provides tools for building whole-program LLVM bitcode files from unmodified C or C++ source code.
	\item \label{itm:step_ll_to_dot} Generate a control flow graph for each function defined in the LLVM IR, and store the CFGs in Graphviz DOT format.
	\item \label{itm:step_dot_to_json} Apply each control flow recovery method on the generated CFGs, and compare the recovered high-level primitives against the truth table produced in step 1.
\end{enumerate}

Special consideration has been taken to ensure that the C source code is parsed identically in step 1 and 2, since pre-processor directives and compiler specific macros may otherwise introduce uncertainty into the results. The Python script which generates truth tables in step 1 parses C source files using the official bindings for the Clang library, and the C source files are converted to LLVM IR in step 2 using Clang (also when invoked through WLLVM).

ext: tex

max_stars_repo_name: decomp/doc

max_issues_repo_name: decomp/doc

max_forks_repo_name: decomp/doc


text:
\section{Data Analysis and Results}
% description of any analytically steps: parameter estimation, error estimation, model-fitting...

\subsection{Detecting Molecular Gas Features}

We compare the final sky-subtracted and wavelength calibrated spectrum of the
Earthshine to predictions and existing measurements. First, we
detect the molecular bands of O$_2$, O$_3$, and H$_2$O by comparing the
 absorption lines obtained from \cite{nuclear}.  We verify that the abundance of
those components are compatible to the
literature (see  Fig. \ref{mol}).

\begin{figure}[htb]
\begin{center}
\includegraphics[scale=0.5]{plots/mo1.png}
  \caption{Atmospheric molecular bands detected in the Earthshine spectrum.}
\label{mol}
\end{center}
\end{figure}

The presence of strong
absorption lines in the near infrared when comparing to  shorter regions (blue)
of the spectra indicates that if human vision were sensitive a little further
toward the red, the natural world would be very red and exceeding bright. In the
next session we also see that the blue spectrum from the subsurface ocean water
has a very small contribution ranging for wavelengths around 500 to 600 nm.

\subsection{Detecting the Vegetation Red Edge}

The main contributors to the optical spectrum of Earthshine \cite{writeup}
\cite{USGS} are
\begin{enumerate}
 \item The {\bf neutral reflectivity from high clouds} (same as the Sun
(blackbody) with $T \sim 5700$K and independent of the wavelength).
\item The blue  spectrum of subsurface {\bf ocean water} (see Fig. \ref{veg} (right))
.
\item The {\bf transmission of Earth's  clear atmosphere}.
\item The {\bf Rayleigh scattering} (proportional to $1/\lambda^4$).
\item The {\bf vegetation} reflection spectrum from land chlorophyll
plants (see  Fig. \ref{veg} (left)).
\end{enumerate}

\begin{figure}[htb]
\begin{center}
\includegraphics[scale=0.22]{figs/green.png}
\includegraphics[scale=0.22]{figs/sea.png}
\caption{Some of the models we fitted our spectra: (left) Vegetation reflection
and (right) Ocean Surface.}
\label{veg}
\end{center}
\end{figure}

 We fit together each of the previous five models to our spectra of the
Eartshine, inferring the partial contributions to
the combined spectrum (see Fig. \ref{models}), \ie our model has 5 parameters. We calculate the best values for these parameters
to minimize $\chi^2_{point}$ on each point then integrate to an acceptable
final value of  $\chi^2=7.8$ (see fitting code in ROOT in the appendix). We
assume the data are normally distributed with a variance equal to the mean and
the data points are independent from each other, so we can use the level of
significance of $\alpha >0.1$ for the fit. For five parameters, the $\chi^2_5$
should be less than 9, so we see that our fit is indeed significant: $\chi^2$ is
not too large for a good fit and it is large enough to reject the {\it null
hypothesis} \cite{writeup}.

\begin{figure}[htb]
\begin{center}
\includegraphics[scale=0.6]{plots/fit.png}
\caption{The main models to the optical spectrum of Earthshine which we fit
to our data. They are scaled in the plot for a better display.}
\label{models}
\end{center}
\end{figure}

The most important components in the fitting are the Rayleigh scattering
(high cloud continuous) in the beginning of the spectra and the clear atmosphere
(molecular spectrum) in the whole range.

We calculate the vegetation edge from equation
\ref{vge}. With mean reflectances in the [600, 670 nm] and [740, 800 nm]
windows in the spectrum, we obtain
$VGE =4 \pm 5\%$, which was small but compatible with
the literature \cite{USGS}.
The errors were estimating by looking at the standard deviation during
a ``flat'' part of the spectra and dividing by the mean in that region.


text:
\section{Introduction}
\label{sec:introduction_chip}

The cost of \gls{SNP} arrays have declined over the years. This has made it possible to increase study sizes at a fixed budget. However, most current \gls{SNP} arrays do not capture variation in African populations well. As shown in table \ref{tab:chip_SNPs} on page \pageref{tab:chip_SNPs}, figure \ref{fig:SN02f2} and elsewhere\cite{Gurdasani2015} approximately 20\% of the autosomal SNPs on the Illumina HumanOmni2.5 SNP array are monomorphic in African populations. And as shown in figure \ref{fig:SN10f2} on page \pageref{fig:SN10f2} \gls{SNP} density also alters imputation accuracy quite significantly. And SN2 figure 2 on page 58 of the supplementary material to the AGV project publication\cite{Gurdasani2015} shows that the site frequency spectrum on the Omni2.5 array in 3 African populations has ascertainment bias by being skewed towards common variants. Therefore it is desirable to explore whether a continent specific \gls{SNP} array can capture variation in African populations better. I do so here. This work should be seen in relation to the results on page \pageref{sec:rp_results} showing that common variation is captured well relative to low coverage sequencing, when continent specific reference panels are used. The development of continent specific reference panels and \gls{SNP} arrays could possibly be a great combination that could compete with study designs involving sequencing. I refer to the \gls{SNP} array as hypothetical here, because no manufacturer has confirmed, that the selected tag \glspl{SNP} are viable to go on a \gls{SNP} array. But by using nearby \glspl{SNP} that can be assayed it should be possible to achieve results very comparable to those presented here.

\input{fig/SN02f2}


text:
\documentclass[11pt]{article}

%%% PAGE DIMENSIONS
\usepackage{geometry}
\geometry{a4paper}

%%% PACKAGES
\usepackage[english]{babel}
\usepackage[utf8]{inputenc}
\usepackage{graphicx}
\usepackage[parfill]{parskip}
\usepackage{booktabs}
\usepackage{array}
\usepackage{paralist}
\usepackage{verbatim}
\usepackage{subfig}
\usepackage{cite}
\usepackage{amsmath}
\usepackage[colorinlistoftodos]{todonotes}

%%% IMAGE FILES PATH
\graphicspath{{Images/}}

%%% HEADERS & FOOTERS
\usepackage{fancyhdr}
\pagestyle{fancy}
\renewcommand{\headrulewidth}{0pt}
\lhead{}\chead{}\rhead{}
\lfoot{}\cfoot{\thepage}\rfoot{}

%%% SECTION TITLE APPEARANCE
\usepackage{sectsty}
\allsectionsfont{\sffamily\mdseries\upshape}

%%% TABLE OF CONTENTS APPEARANCE
\usepackage[nottoc,notlof,notlot]{tocbibind}
\usepackage[titles,subfigure]{tocloft}
\renewcommand{\cftsecfont}{\rmfamily\mdseries\upshape}
\renewcommand{\cftsecpagefont}{\rmfamily\mdseries\upshape}

%%% DOCUMENT
\title{Mapping the Brain: An Introduction to Connectomics\\Streamlining Analysis for Automated Neuronal Membrane Segmentation Pipelines}

\author{Thomas Keady, Albert Lee, Augusto Ramirez}

\date{\today}

\begin{document}
\maketitle

\section{Abstract}

Creating a complete wiring diagram of neuronal connectivity is one of the most ambitious goals of connectomics and contemporary neuroscience \cite{roncal1}. Thanks to technical and computational advances that automate the collection of high resolution electron-microscopy data, complex networks of neurons and their synaptic connections can be directly observed and eventually fully mapped. An important step in this process is the annotation of membranes, which allows for deeper levels of analysis. However, the most reliable method of membrane detection to date is performed manually \cite{plaza}. Although automated membrane segmentation tools exist, the landscape of such algorithms consists mainly of poor documentation and untested parameters/optimization. The neuroXplorer team ran tested popular image processing functions such as watershed and gala in different testing environments (Matlab, Python) \cite{nguyen}.

By comparing a range of parameters on a set of images, relevant statistical data can be inferred and utilized to select the arguments needed to optimize the tools mentioned earlier. Reducing the average error rates for metrics such as split\_VI, adjusted rand index, and Fowlkes-Mallows index can help us indicate optimized parameters for each image.

\section{Results}

Three parameters within our program were varied to find their optimal values. These were the size of the closing footprint, the size of the median filter, and the size of the watershed footprint. To obtain the ideal parameters, we tested a sparse but evenly distributed sampling of the possible combinations of values on our training dataset. We then found the parameters that led to the most similar segmentation to ground truth according to four metrics: the rate of under segmentation (merge errors), the rate of over segmentation (split errors), the adjusted rand index, and the Fowlkes-Mallows index. For each image, we recorded the parameters that led to the best segmentation according to each of the four metrics. We then used the resulting ranges of parameters to do a narrow, more denser sweep to more finely tune the parameters. After again identifying the best segmentations and their corresponding parameters, we averaged each parameter and rounded to the nearest integer to get our final parameters. The size of the closing footprint became 24, the size of the median filter was 7, and the size of the watershed footprint was 1. Across all metrics the optimized parameters were very close to each other, +/- 1. There was one exception: the undersegmentation metric was best optimized when the watershed footprint was 6 (6.07 averaged), while all three other metrics optimized it to 1.00 exactly. We decided to use the optimized parameter as 1 instead of averaging three 1's with a 6, and applied these parameters to our testing dataset. All of the metrics appeared similar between the training and testing datasets. We had slightly better segmentations in the training set, according to the graphs on our poster. We speculate that this was due to underfitting, where the algorithm was not subjected to enough possible cases to optimize its parameters as best as possible, rather than overfitting where its optimization became so tailored to the training set such that it did not extrapolate well. We primarily believe this because of the relatively small size of our training set (50 sections). Also, in both the training and testing datasets, the algorithm was prone to oversegmentation as can be seen in our graphs and images. However, we point out that our algorithm does considerably better with merge errors, and believe this is valuable as merge errors are more difficult to rectify than split errors.

\begin{figure}
	\includegraphics{truth_anno1}
	\caption{Fully annotated ground truth.}
\end{figure}

\begin{figure}
	\includegraphics{fig_0_24_6_29}
	\caption{Watershed optimized membrane segmentation.}
\end{figure}

\begin{thebibliography}{1}
\bibitem{roncal1} W. G. Roncal, D. M. Kleissas, J. T. Vogelstein, "An automated images-to-graphs pipeline for high resolution connectomics." \textit{Frontiers in Neuroinformatics}, vol. 9, no. 20, Aug. 2015.
\bibitem{plaza} S. Plaza, "Focused Proofreading: Efficiently Extracting Connectomes from Segmented EM Images." \textit{Cornell University Library}, Sep. 2014.
\bibitem{nguyen} Q. Nguyen,  "Parallel and scalable neural image segmentation for connectome graph extraction.” \textit{DSpace at MIT}. Aug. 2015.
\end{thebibliography}

\end{document}

max_stars_repo_name: Connectomics-Classes/Neuro-x-polorer

max_issues_repo_name: Connectomics-Classes/Neuro-x-polorer

max_forks_repo_name: Connectomics-Classes/Neuro-x-polorer