\documentclass[12pt,oneside]{book}
\usepackage[utf8]{inputenc}
\usepackage[a4paper, total={6.4in, 9.7in},top=1in]{geometry}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{hyperref}
\usepackage{parskip}
\usepackage{xcolor}
\usepackage[pdftex]{graphicx}
\graphicspath{{./Figures/}}
\usepackage{enumerate}
\usepackage{flafter}
\usepackage{caption}
\captionsetup{justification   = raggedright, singlelinecheck = false}
\usepackage[section,above,below]{placeins}

\usepackage{listings}
\lstset{ 
  language=R,                     % the language of the code
  basicstyle=\scriptsize\ttfamily, 	% the size of the fonts that are used for the code
  numbers=none,                   % where to put the line-numbers
  backgroundcolor=\color{white},  % choose the background color.
  showspaces=false,               % show spaces adding underscores
  showstringspaces=false,         % underline spaces within strings
  showtabs=false,                 % show tabs within strings
  frame=single,                   % adds a frame around the code
  rulecolor=\color{black},        % if not set, the frame-color may be changed on line-breaks within not-black text
  tabsize=2,                      % sets default tabsize to 2 spaces
  breaklines=true,                % sets automatic line breaking
  breakatwhitespace=false,        % sets if automatic breaks should only happen at whitespace
  keywordstyle=\color{blue!70!black},      % keyword style
  commentstyle=\color{green!50!black},   % comment style
  stringstyle=\color{orange!70!black}      % string literal style
}

\usepackage[backend=biber, style=authoryear, natbib=true, maxbibnames=6, maxcitenames=2, giveninits=true, doi=false]{biblatex}
\addbibresource{SeanPhD.bib}
\setlength{\bibparsep}{6pt}
\DeclareNameAlias{sortname}{last-first}
\renewcommand{\newunitpunct}{,\ }
\renewcommand{\intitlepunct}{}
\DefineBibliographyStrings{english}{in = {{}{}}}

\newcommand{\mytitle}{Topics in the analysis of composition data}
\hypersetup{pdfinfo={
	Title={\mytitle},
	Subject={\mytitle},
	Author={Sean van der Merwe},
	Keywords={Bayes, Dirichlet, Goodness-of-fit, regression, Posterior, Simulation} } }

\renewcommand{\chapterautorefname}{Chapter}
\renewcommand{\sectionautorefname}{Section}
\renewcommand{\subsectionautorefname}{Section}
\renewcommand{\subsubsectionautorefname}{Section}
\newlength{\figwidth} \setlength{\figwidth}{10cm} % I use this to standardise figure sizes

\usepackage{enumerate}
\usepackage{soul}
\newcommand{\com}[1]{{\color{red}#1}} % I'm going to use this command to highlight my changes, remove the red as you see fit. I'll use [] for comments inside it.
\newcommand{\new}[1]{{\color{blue}#1}}   

\newcommand{\bx}{\ensuremath{\mathbf{x}}}
\newcommand{\by}{\ensuremath{\mathbf{y}}}

\renewcommand{\baselinestretch}{1.15}

\usepackage{fancyhdr}
\fancypagestyle{C1}{\fancyhf{} \fancyhead[L]{Chapter 1 --- Introduction} \fancyhead[R]{\thepage} \fancyfoot{} \setlength{\headheight}{15pt}}
\fancypagestyle{C2}{\fancyhf{} \fancyhead[L]{Chapter 2 --- Fitting} \fancyhead[R]{\thepage} \fancyfoot{} \setlength{\headheight}{15pt}}
\fancypagestyle{C3}{\fancyhf{} \fancyhead[L]{Chapter 3 --- Discrimination} \fancyhead[R]{\thepage} \fancyfoot{} \setlength{\headheight}{15pt}}
\fancypagestyle{C4}{\fancyhf{} \fancyhead[L]{Chapter 4 --- Classification} \fancyhead[R]{\thepage} \fancyfoot{} \setlength{\headheight}{15pt}}
\fancypagestyle{C5}{\fancyhf{} \fancyhead[L]{Chapter 5 --- MGAM} \fancyhead[R]{\thepage} \fancyfoot{} \setlength{\headheight}{15pt}}
\fancypagestyle{C6}{\fancyhf{} \fancyhead[L]{Chapter 6 --- Regression} \fancyhead[R]{\thepage} \fancyfoot{} \setlength{\headheight}{15pt}}
\fancypagestyle{C7}{\fancyhf{} \fancyhead[L]{Chapter 7 --- Conclusion} \fancyhead[R]{\thepage} \fancyfoot{} \setlength{\headheight}{15pt}}
\fancypagestyle{plain}{\fancyhf{} \fancyhead{} \fancyfoot[C]{\thepage} \setlength{\headheight}{15pt} \renewcommand{\headrulewidth}{0pt}}
\fancypagestyle{preamble}{\fancyhf{} \fancyhead[L]{Preamble} \fancyhead[R]{\thepage} \fancyfoot{} \setlength{\headheight}{15pt}}
\fancypagestyle{contentsstuff}{\fancyhf{} \fancyhead[L]{\leftmark} \fancyhead[R]{\thepage} \fancyfoot{} \setlength{\headheight}{15pt}}

\usepackage{animate}

\begin{document}
\hfuzz=32pt
\hbadness=4800
\begin{refsection}
\pagenumbering{gobble}
\begin{center}
\includegraphics[width=1.2\figwidth]{UFSmathstatsLogo.png}

\begin{Huge}
\mytitle
\end{Huge}
\vspace{1cm}

\begin{LARGE}
Sean van der Merwe

2001007828
\end{LARGE}
\begin{Large}
\vspace{2cm}

Supervisors: Daniel de Waal, Michael von Maltitz
\end{Large}
\vspace{2cm}

Submitted in fulfilment of the requirements for the degree:\\
Ph.D. (Mathematical Statistics)\\
in the Department of Mathematical Statistics and Actuarial Science in the Faculty of Natural and Agricultural Sciences at the University of the Free State.
\vspace{1cm}

\begin{Large}
\today
\end{Large}

\end{center}

\newpage
\frontmatter
\setcounter{page}{1}
\pagestyle{preamble}
\section*{Declaration}
I, Sean van der Merwe, hereby declare that this work, submitted for the degree Ph.D. (Mathematical Statistics), at the University of the Free State, is original work and has not previously been submitted, for degree purposes or otherwise, to any other institution of higher learning. I further declare that all sources cited or quoted are indicated and acknowledged by means of a comprehensive list of references. Copyright hereby cedes to and is invested in the University of the Free State. I declare that all royalties as regards to intellectual property that was developed during the course of and/or in connection with the study at the University of the Free State, will accrue to the University. 

\includegraphics[width=60mm]{SeanvdMsignature1.png}

\newpage

\section*{Acknowledgement}

I want to thank Jesus first, for keeping my spirits up through failures and rejections, and reminding me that research is a service to the world, not for any personal gain.

I want to thank those that directly helped me with this work, for guidance and direction, even when I stubbornly refused to listen. In particular, I'm referring to Daan de Waal, Roelof Coetzer, Michael von Maltitz, and Robert Schall.

I want to thank my whole department for being such amazingly pleasant people to work with and always supporting me, even when they didn't know what I was doing or weren't supportive of it. They are all role models for me.

I want to thank my friends and family, for helping me maintain balance, and for their interest in my work.

\newpage
\section*{Abstract}
Many scientific and industrial processes produce data that is best analysed as vectors of relative values, often called compositions or proportions. Any data collected as vectors where the proportion measured in each category is of more importance than the quantity surveyed falls into this area of analysis. Such data is multivariate in nature but should not be analysed using standard multivariate analysis methods directly. Techniques specific to composition data are investigated and recommended in this work.

The Dirichlet distribution is a natural distribution for composition or proportion data. It has the advantage of a low number of parameters, making it a parsimonious choice in many cases. This work primarily explores the Dirichlet distribution from various angles. These include discriminating between populations, classifying observations, and regression of compositions on data of other types. Other models for compositions are also discussed. 

Most topics are investigated from the Bayesian perspective. This work particularly contributes to the field of Bayesian analysis using the Dirichlet distribution, expanding the theory of fitting and applying the Dirichlet type I and Dirichlet type II distributions using Bayesian techniques. Prediction and imputation is also explored.

A key contribution of this work is in developing a general approach to building advanced models for composition data, in the Bayesian framework, using the Dirichlet distribution as starting point. These models may include different types of explanatory variables for both expected values and precisions, as well as a mixture of fixed and random effects.



\vspace{4mm}

\noindent%
\textbf{Keywords:}  Composition data, proportions, Dirichlet distribution, Bayes, simulation

\vspace{4mm}

%\noindent%
%{\it Note on media:}  This document is meant to be viewed in Portable Document Format (.pdf). If you did not receive it in this format then you may obtain the latest version here: \url{https://goo.gl/yBhSUu}.

\newpage
\section*{List of Abbreviations}
\begin{tabular}[t]{|l|l|} \hline
\textbf{Abbreviation}&\textbf{Meaning}\\ \hline
D&Dirichlet type I (distribution)\\ \hline
D2&Dirichlet type II (distribution)\\ \hline
LDA&latent Dirichlet allocation (model)\\ \hline
ln&natural logarithm\\ \hline
LN&logistic-normal (distribution)\\ \hline
log&natural logarithm\\ \hline
MD&matrix Dirichlet (distribution)\\ \hline
MDI&maximal data information (prior)\\ \hline
MDirich2&matrix Dirichlet type II (distribution)\\ \hline
MEV&multivariate extreme value (distribution)\\ \hline
MGAM&multivariate gamma (distribution)\\ \hline
ML&(method of) maximum likelihood\\ \hline
SCE&sum of composition errors\\ \hline
\end{tabular}

\newpage
\pagestyle{contentsstuff}
\tableofcontents
\newpage
\listoffigures
\newpage
\listoftables
\newpage

\newpage
\pagestyle{preamble}
\section*{Problem Statement}
In many experiments or data scenarios the situation exists where measurements are relative, or only of interest to researchers when expressed relative to each other in such a way that they sum to one. The result is data that is multivariate, but not directly suitable for analyses using the popular multivariate normal analysis methods. Further, sample sizes may be small, so even if data are transformed to a free scale then the multivariate normal distribution may have too many parameters, leading to over-fitting. Excessive transformations also affect interpretability.

The Dirichlet and related distributions can be seen as a parsimonious solution, but the published literature is not yet extensive. Most notably, there is a need for Bayesian solutions involving the Dirichlet distribution.

This work seeks to address the problems of:
\begin{itemize}
\item Bayesian analysis of the Dirichlet distribution,
\item discrimination and classification with composition data,
\item and regression of composition data using the Dirichlet distribution.
\end{itemize}

\newpage
\section*{Overview}

In \autoref{ch:intro} the main topics are introduced and illustrated, relevant applications are discussed by means of examples from literature, and advances in the theory of the topics are also mentioned.

In \autoref{ch:p1} the Dirichlet type I and type II distributions are discussed in detail. New theory is derived. Objective Bayesian priors are derived. Performance of the derived Bayesian parameter estimation methods are compared to existing methods via a simulation study. Prediction and imputation using these distributions are also explained.

In \autoref{ch:p2} the discriminant function is developed for the Dirichlet and Matrix Dirichlet distributions. Uses are explained both in theory and with examples.

In \autoref{ch:p3} two new approaches to classification of Dirichlet observations are developed and discussed in detail. These are presented in the context of coal properties as inputs to a coal gasification process.

In \autoref{ch:p4} an alternative to the existing distributions for composition data is developed, namely a new class of multivariate gamma distribution. Again, classification in the context of coal gasification is used to explain the use of this approach. This alternative relates to the Dirichlet distribution in that the Dirichlet can be created as a transformation of a special case of this multivariate gamma distribution.

In \autoref{ch:p5} the use of Dirichlet distributions in a regression setting is discussed. A new approach is developed and explained. The advantages of the new approach are illustrated via simulation and example.

In \autoref{ch:conclusion} the advances made through the work above are summarised.

It is worth mentioning that many of the approaches discussed are not compared directly, as they fall in different model classes. Specifically, the Dirichlet models can be seen as the parsimonious class, as the number of parameters increases in line with the dimension of the problem; while the other models discussed (\textit{e.g.}\ multivariate gamma and logistic normal) are referred to as the rich class, since the number of parameters in those models increase in line with the square of the dimension. Models are compared to others in their own class only.

The thread connecting the topics is an expansion of the tools for working with composition data. The focus is on the Dirichlet Type I models. Being in the parsimonious class, these models require less data points for accurate parameter estimation, and they also have the benefit of not requiring transformation of composition data prior to analysis.

\newpage
\section*{Comment on notation}
Chapters 2 to 6 are based on individual papers that were each written separately, by different combinations of authors, from slightly different viewpoints. As such, they originally contained differing definitions and notation. In order to increase the readability of this document, each paper has undergone minor editing. 

In each of these chapters a link is provided to an online version of the relevant paper, in case the reader wishes to view the original. Most changes are related to a consistent notation for the Dirichlet distribution. The main reason for the differences from the original papers is discussed in the first chapter.

\newpage
\section*{Authorship and publication}
Paper 1 (\autoref{ch:p1}) is almost entirely the work of the student, under the direction of Prof.\ Daan de Waal. All derivations are done by the student. All simulation studies were designed and executed by the student. The multiple imputation approach is original to the student, adapted from the multivariate normal imputation approach. It is posted as a technical report as it discusses theory in detail to act as a reference source.

Paper 2 (\autoref{ch:p2}) is mostly the work of Prof.\ Daan de Waal, with the student providing assistance. The theory was derived by the supervisor, along with initial experiments. Results were checked and edited by the student. The student was responsible for typesetting and editing of the publication. This paper is published in the South African Statistics Journal. 

Paper 3 (\autoref{ch:p3}) contains equal contributions from Prof.\ Daan de Waal and the student, with Prof.\ Roelof Coetzer providing leadership and direction. The introduction and interpretations relating to the application were written by Prof.\ Coetzer, who posed the original problem. The first approach to solving the problem was derived by the supervisor. The second approach was derived by the student. The student was responsible for typesetting and editing of the publication, and responding to reviewer feedback. This paper is published in Chemometrics. 

Paper 4 (\autoref{ch:p4}) is mostly the work of Prof.\ Daan de Waal, with the student providing assistance and expanding the theory in places. Prof.\ Roelof Coetzer provided valuable assistance. The new distribution was developed by the supervisor. The student was responsible for the Appendix. The student was responsible for typesetting and editing of the publication, and responding to reviewer feedback. This paper is published in the South African Statistics Journal.

Paper 5 (\autoref{ch:p5}) is entirely the work of the student and is currently under review. The work has been accepted for publication in the South African Statistics Journal. 

The papers are presented in chronological order of being written; however, this order does maintain a logical flow and so the work can be read as a unit.















% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\mainmatter
\chapter{Introduction and Literature} \label{ch:intro}
\pagestyle{C1}
\section{Composition data}

Composition data refers to a set of multivariate observations where each individual observation vector sums to exactly one. Alternatively, it is often defined as a vector summing to a positive real number less than one, with a residual `other' component that is considered non-random making up the difference between one and the sum of the random components. Each of these definitions seem more applicable to some practical situations than others, so examples are discussed in this chapter.

When vectors are measured in whole numbers, for example in a classification process, analysis is often based on the multinomial distribution, with the total count being seen as a nuisance parameter. In cases where the total count is not relevant it seems more natural to work directly with the observed proportions. Sometimes the proportions themselves are observed directly, rather than counts and totals. This class of problem includes any situation where the quantity surveyed does not affect the expectation vector.

A good example of composition data is food composition. Most foods have their compositions analysed and reported for health and safety reasons, and being able to model changes in composition can be valuable to regulators. See the book of \textcite{greenfield} for an explanation of the creation and analysis of this type of data. In \textcite{Davis2004} the composition of food was compared between 1950 and 1999 to investigate possible drops in nutritional value over that period. Astrophysicists have considered using the composition of solar wind to determine its origin \parencite{geiss1995}. Another composition example is given by \textcite{wedepohl2011data}, where the elemental composition of glass was used to determine the glass' region of origin. Compositions are of great interest in industrial applications, where the inputs to an industrial process can affect the financial worth of the output.

\subsection{Note on conflicting definitions}

The most popular approaches to analyse composition data involve designating one component of the vector as a reference component. The other components are considered random, while the reference component is not. Analysis then proceeds relative to the reference component. These approaches were regularly applied out of necessity, due to the absence of alternatives. While there may be practical justifications to these approaches in some cases, it is argued here that applying these approaches without clear justification results in loss of practical interpretability.

In the chapters to come the theory of composition data analysis is expanded in such a way that reliance on a reference category is reduced, freeing the researcher to analyse and interpret all components directly. In accordance with this, definitions and notation have been adapted in places to reflect the idea that all components carry equal weight.

Formally, the following definition is used:

$$
\begin{aligned}
&\text{A random variable }\mathbf{X}\text{, a vector of size }P\text{ by 1, with components }X_1,\ldots,X_P\\
&\text{is called a composition random variable if }\\
&0<X_i<1\ \forall\ i=1,\ldots ,P\text{ and }\sum_{i=1}^P X_i = 1.
\end{aligned}
$$

\section{Composition data analysis methods} \label{sec:OtherMethodsCompos}

Composition data analysis is no different to most types of data analysis and most areas of statistics in that the starting point of analysis is describing the data. Practical statisticians will always advocate `eyeballing' the data before doing any modelling. This means looking at the structure, querying the collection process, looking for anomalies, summarising the data using descriptive measures, and summarising the data visually. 

However, the scale of composition data is not the same as most data types, and so transformations can be useful in the process of summarising and visualising the data. At this stage a researcher may find use in a statistical package such as `compositions' \parencite{van2008compositions} for the R environment \parencite{RCore}. The package authors explain that their package supports multiple scales and moving between them. It allows descriptive and visual analyses on all scales, and can perform principal component analysis or cluster analysis if advanced visualisation would be useful. See \textcite{rencher2002} for an explanation of these techniques.

\textcite{Aitchison1986} advocates the use of the log-ratio transform followed by classical multivariate normal distribution analysis. He calls this the logistic-normal distribution (LN). This allows advanced multivariate analyses to be carried out on composition data. Again, see \textcite{rencher2002} for an introduction to this area of statistics.

As an example, \textcite{lafferty2006} adapt their latent Dirichlet allocation (LDA) model (discussed in \autoref{sec:DirAppliIntro}) to the LN framework and call it Correlated Topic Models.

However, there are two difficulties with the LN approach. The first is interpretability. Since the transformations done are non-linear, all results should be transformed back to the original scale in a multivariate fashion if the researcher wishes to make judgements on the proportions scale. Often this is not possible, or very difficult. This difficulty persists with almost all analyses of composition data, an ever-present concern for researchers in this field.

The second difficulty is parameterisation. If the number of categories or dimensions is small then the multivariate normal approach yields great flexibility and is often useful even when the data does not fit the normal distribution perfectly. However, when the number of dimensions increases then the number of parameters inflates rapidly, and simpler models that still capture the nature of the data are required.

\section{The Dirichlet Distribution} \label{sec:DirIntroCh1}

The origin of the Dirichlet distribution is discussed in detail in the paper of \textcite{Gupta2001}. They explain that the distribution owes its name to Johann Peter Gustav Lejeune Dirichlet, who evaluated the integral that completes the distribution \parencite{Dirichlet1839}.

\subsection{Definition of the Dirichlet Distribution} \label{sec:DirDefCh1}

If $\mathbf{Y}$ (size $P$ by 1) is distributed Dirichlet$(\alpha_1,\ldots,\alpha_{P})$, denoted by $D(\boldsymbol\alpha)$, then the density is given by 
\begin{equation} \label{eq:Dintroduction}
f(\mathbf{y})= \frac{\prod_{i=1}^{P}\Gamma(\alpha_i)}{\Gamma(\alpha_0)}\left\{ \prod_{i=1}^P y_i^{\alpha_i-1}\right\}\ ,\ 0 < y_i < 1,\ \sum_{i=1}^P y_i = 1,\ \alpha_0 = \sum_{i=1}^{P} \alpha_i 
\end{equation}

The above definition is similar to the one given in \textcite{Maier2014}, who uses illustrations to help explain the distribution. Original illustrations are given here for the same purpose. In two dimensions the Dirichlet distribution is the classic beta distribution (\autoref{fig:D2illus}) which is well understood. In three dimensions (\autoref{fig:D3illus}) the flexibility of the beta is extended intuitively. In four dimensions the density becomes difficult to visualise, but a sample from the density can be illustrated as in (\autoref{fig:D4illus}). If one wishes to visualise five or more dimensions, each new dimension must be introduced in a creative way, say by changing the size, colour, or shape of the points of the sample.

\begin{figure}[hbt]
\begin{center}
\includegraphics[width=0.5\figwidth]{Dirichlet2Dg1.pdf}
\includegraphics[width=0.5\figwidth]{Dirichlet2Dg2.pdf}
\end{center}
\caption{Illustrations of the Dirichlet density with dimension two (beta density)}

\label{fig:D2illus}
\end{figure}

\begin{figure}[hbt]
\begin{center}
\includegraphics[width=0.5\figwidth]{DirichletD3g1.png}
\includegraphics[width=0.5\figwidth]{DirichletD3g2.png}
\end{center}
\caption{Illustrations of the Dirichlet density with dimension three}

\label{fig:D3illus}
\end{figure}

\begin{figure}[hbt]
\begin{center}
\animategraphics[controls,autoplay,loop,width=\figwidth ,trim=0 0 0 60]{0.5}{./Figures/DirichletD4ani}{1}{12}
\end{center}
\caption{Random sample from a Dirichlet density with dimension four}

\label{fig:D4illus}
\end{figure}


\subsection{Properties of the Dirichlet Distribution} \label{sec:DirPropsIntro}

Understanding the properties of a distribution helps the researcher to know when to apply the distribution and when not to. Before looking at applications in the next section, the main properties are briefly discussed in an intuitive way. For a more rigorous discussion, see \autoref{sec:intro}.

The Dirichlet distribution is applied to vector observations that sum to one. If one component of the vector increases then at least one other component must decrease, and thus there is a natural negative relationship between components of the vector. The negative relationship is reflected in the moments of the Dirichlet distribution, see \autoref{eq:IntroChapmeancov} below.

\begin{equation} \begin{aligned}\label{eq:IntroChapmeancov}
&E(Y_i)=\frac{\alpha_i}{\alpha_0},\ i = 1,\ldots,P\\
&Cov(Y_i,Y_j)=\left\{\begin{array}{l}
\frac{-\alpha_i \alpha_j}{\alpha_0^2 (1+\alpha_0)},\ i\ne j\\
\frac{\alpha_i (\alpha_0-\alpha_i)}{\alpha_0^2 (1+\alpha_0)},\ i=j
\end{array}\right.
\end{aligned} \end{equation}

Obtaining the mean vector and covariance matrix is a natural step in data analysis, and should the covariances be all negative then that would support the Dirichlet distribution as a candidate for direct modelling. However, in practical situations it is often the case that positive correlations occur, especially as the number of components in the vector (number of dimensions) increase. When clear positive correlations are noted then one should not apply the Dirichlet distribution directly.

There are several ways to address the issue of positive correlations between components without abandoning the Dirichlet distribution. The first is to combine components by summing them. The aggregation property of the Dirichlet distribution \citep[proved in][]{Frigyik2010} says that if two components of a Dirichlet vector are dropped and replaced with their sum then the resulting vector is again Dirichlet with the corresponding parameters replaced in the same way. Aggregating components with positive correlations can only enhance the fit Dirichlet distribution.

The second approach is to consider subsets separately where appropriate. When attempting to discriminate between Dirichlet populations it is important to note that each subset should fit the Dirichlet distribution, not necessarily the sample as a whole. As a crude example, if a researcher is interested in the effect of gender on voting patterns then that researcher can model the male voting patterns separately from the female or other gender voting patterns, and check each model fit on its own merits. See \autoref{ch:p2} for more information.

The last consideration is the case where the vectors are regressed on continuous explanatory variables. Here the requirement of negative correlations generally falls away \citep{Maier2014}. Positive correlations are assumed to be induced via joint dependence on the explanatory variable.

Another important aspect of the Dirichlet distribution is that the components of the observation vectors must be strictly between zero and one. When true zeros are observed then the Dirichlet distribution is not advised, unless the issue can be address by one of the three methods above. When values that are very small and rounded to zero are observed, then adjustments can be made to enable Dirichlet modelling. These include pulling extreme observations in from the edge of the domain \citep{Maier2014}, treating them as censored (implicit imputation) or imputing them explicitly (similar to \autoref{sec:imputation}).

Should none of the adjustments be appropriate then more complex modelling should be considered using alternative distributions (often involving transformations). For more information see \autoref{ch:p4}.

\subsection{Applications of the Dirichlet Distribution} \label{sec:DirAppliIntro}

The most common application of the Dirichlet distribution is as a conjugate prior for Bayesian analyses (\autoref{sec:whatisBayes}) involving the multinomial distribution. The result is known as the Dirichlet-multinomial distribution or Dirichlet compound multinomial distribution, which traces back to \textcite{Mosimann1962}. For a modern example related to text mining see \textcite{Madsen2005}.

The books of \textcite{Ng2011} and \textcite{balakrishnan2006} give some insight into the applications of the Dirichlet and related distributions. Another interesting example is the work by \citet{lange1995} for estimating forensic match probabilities.

The Dirichlet distribution has also been used in economics. Income inequality is a major topic of discussion in economics and it is usually illustrated via the Lorenz curve. \textcite{chotikapanich2002} show that by assuming income proportions are distributed Dirichlet they can improve modelling of such curves.

In medical science, \textcite{briggs2003} use the Dirichlet distribution to develop a fully probabilistic transition matrix for a reliable Markov model for medical interventions.

In computer science, the Latent Dirichlet Allocation model \parencite{blei2003latent} is a very popular model in machine learning, as it allows automatic characterisation of documents according to topic probabilities. The uses of this tool are extensive --- see \textcite{NIPS2010} and \textcite{Lienou2010} as starting points for exploring this interesting research area. Additionally, modern online reputation systems have the power to determine the success of a product or business, so it should be no surprise that their study and optimisation is of great value. \textcite{josang2007} discuss the use of the Dirichlet distribution in this process. Also on the theme of machine learning, \textcite{bouguila2003} develop mixtures based on the Dirichlet Distribution for pattern-recognition in image databases.

\subsection{Dirichlet regression} \label{sec:DirRegIntroCh1}

The Dirichlet distribution has been widely accepted in literature for modelling composition data, subject to the constraint that all the correlations between variables are negative \citep{Maier2014}. A wider class of distributions which allows for positive correlations defined on the same sample space is the logistic-normal (LN) \citep{Aitchison1980}. However, the LN distribution has many parameters to estimate due to the unknown covariance parameter matrix. In contrast, the Dirichlet distribution has only $p$ unknown parameters to estimate for compositions with $p$ components. Further, when the model is expanded to the regression framework then the requirement of negative correlations no longer applies \citep{Maier2014}.

The distinctive characteristic of the regression framework is that it allows for continuous measurements as explanatory variables. The goal is to simultaneously model vector observations in their entirety based on a flexible structure that allows for the analysis of relationships between the observed vectors and the explanatory variables.

As examples of applications, consider the work of \citet{boukal2015analyses} regarding the development of insects, or \citet{EspinGarcia2014} who look at genetic analysis problems, or the work of \citet{smithson2006better} where applications in psychology are discussed. 

Another area in which Dirichlet regression may be useful is in politics, where only the proportion of voters supporting a set of candidates is of interest, possibly dependent on district or demographics. In ecology, preferences of animals for specific types of prey is of interest. In any situation where people choose between a set of options, and the researcher is interested in the choice and not the number of people making it, then the Dirichlet regression model may be of use.

\section{Bayesian analysis} \label{sec:whatisBayes}

Bayesian analysis at its core refers to the idea of starting by understanding one's prior beliefs and then updating those beliefs by observing data \parencite{Efron2013}.

Besides the philosophical differences, there are many practical advantages to the Bayesian approach. It is the experience of the author that the Bayesian approach for capturing model uncertainty is often the biggest advantage as it gives more flexibility in terms of the types of questions that can be answered from a particular analysis. It also simplifies the process of creating uncertain predictions.

Probably the most widely known book on the topic of Bayesian statistics is that of \textcite{gelman2013}, which is a detailed and thorough text. It is mathematically rigorous and can help the reader gain as much insight as desired.

On the applied side, the book of \textcite{Cowles2013} takes the reader quickly from the basic concepts to the practical workings of Bayesian statistics. She explains that people naturally quantify uncertainty in terms of probabilities, but then goes on to explain that the frequentist definition of probability (the proportion of success from an infinite number of trials) is flawed in several ways. Most notably, it relies on the idea that a sequence of repeatable experiments is possible (or at least imaginable) which is often not the case.

In Bayesian statistics a probability is a person's opinion as to how likely a success is. It is unique to the person and can change over time. Sometimes researchers object to the idea of subjectivity in science, saying scientists should strive for objectivity; however, it is impossible to eliminate subjectivity and by attempting to understand and quantify it the researcher can actually increase the trustworthiness of results.

In this work the focus is on `objective' Bayesian analysis, which involves the use of non-informative priors or vague priors, because the solutions discussed are generic. In specific problems where good prior information is available we encourage the researcher to incorporate such knowledge into the analysis.

Prior knowledge is combined with observed data using Bayes' rule. \textcite{Cowles2013} states it in its simplest form as 
\begin{equation} \label{eq:BayesRule}
P(Model|Data)\propto P(Data|Model)\times P(Model)
\end{equation}
where $P(Model)$ is all the knowledge about the model (and in particular parameters) that is available before observing the data, while $P(Data|Model)$ is the likelihood.

The first powerful advantage to the above rule is that it is recursive, allowing knowledge to be updated in a natural way when further data are observed (the posterior of one sample can be used as the prior to a new sample). The second powerful advantage is that extending the definition for making predictions is straightforward:
\begin{equation} \label{eq:BayesPred}
P(New\ Data|Old\ Data)=\int_{Model} P(New\ Data|Model)\times P(Model|Old\ Data)
\end{equation}

Bayesian statistics has risen in popularity \parencite{Efron2013} due to improvements in computing hardware and software over the last few decades, which have made implementation much simpler \parencite{R2OpenBUGS}. The most common way to implement a Bayesian model is to draw random parameter samples from the $Model|Data$ using simulation techniques and then summarise those samples to draw inferences and answer questions. To evaluate the predictive distribution one can simply draw random observations given each of the randomly drawn parameter vectors and summarise those observations in turn.

\section{Outline summary}

As mentioned in the overview, the coming chapters address various topics in the analysis of composition data. The first deals with standard distribution fitting and forms a basis for the topics that follow. The second deals with discrimination between populations. The third deals with classification of observations between populations. The fourth explores a new distribution and its use. The last topic explores the general regression problem and can be seen as the main outcome of the work.

\printbibliography[title=List of References for the introductory chapter]
\end{refsection}























% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{refsection}
\chapter{Paper 1 --- Bayesian Fitting of Dirichlet Type I and II Distributions} \label{ch:p1}
\pagestyle{C2}
\begin{LARGE}
Sean van~der~Merwe and Daan de~Waal
\end{LARGE}
\vspace{10mm}

\textbf{Abstract}

In his 1986 book, Aitchison explains that composition data is regularly mishandled in statistical analyses, a pattern that continues to this day. The Dirichlet Type I distribution is a multivariate distribution commonly used to model a set of proportions that sum to one. Aitchinson goes on to lament the difficulties of Dirichlet modelling and the scarcity of alternatives. While he addresses the second of these issues, this paper address the first. The Dirichlet Type II distribution is a transformation of the Dirichlet Type I distribution and is a multivariate distribution on the positive real numbers with only one more parameter than the number of dimensions. This property of Dirichlet distributions implies advantages over common alternatives as the number of dimensions increase. While not all data are amenable to Dirichlet modelling, there are many cases where the Dirichlet family is the obvious choice. This paper describes the Dirichlet distributions and shows how to fit them using both frequentist and Bayesian methods (the paper derives and applies two objective priors). The beta distribution is discussed as a special case. A small simulation study is reported in order to compare the fitting methods. Conditional distributions and posterior predictive conditional distributions are derived. The flexibility of this distribution family is illustrated via examples, the last of which discusses imputation (using the posterior predictive conditional distributions).

\textit{Keywords}: Bayes, Beta, Dirichlet, Dirichlet Type II, Objective prior, Simulation

\section{Introduction} \label{sec:intro}
The Dirichlet distribution, the multivariate counterpart to the beta distribution, is often used to model a set of proportions that sum to one, or less than one. In the frequentist sense the distribution is well-documented and researched ---~\citet{Minka}, \citet{Ronning1989}, \citet{Wicker2008}, \citet{Narayanan1992}, \textit{etc}.\ all discuss the various approaches to maximising the likelihood, as well as the method of moments; but Bayesian analysis of this distribution appears to focus entirely on its use as a prior and posterior for the multinomial distribution. This paper intends to show that the Dirichlet distribution is useful in its own right and that the Bayesian approach adds to its value.

From \citet{Minka}, if $\mathbf{X}=X_1,\ldots,X_{P} \sim Dirichlet(\alpha_1,\ldots,\alpha_{P})$, then the density is
\begin{equation} \label{eq:density}
f(\mathbf{x}) = \frac{\Gamma(\alpha_0)}{\prod_{i=1}^{P}\Gamma(\alpha_i)}\left(\prod_{i=1}^{P}x_i^{\alpha_i-1}\right)
\end{equation}
where $0 < x_i < 1 , i=1,\ldots,P;\  \sum_{i=1}^{P}x_i =1 $ and $ \alpha_0=\alpha_1+\cdots+\alpha_{P}$.

It can be said that this distribution is $(P-1)$ dimensional because if any $(P-1)$ values of $\mathbf{x}$ are known then so is the last value. However, the practical relevance and implications of this vary between specific applications.

In the case where $P=2$ the density (\autoref{eq:density}) simplifies to the density of the well-known beta distribution:
\begin{equation} 
f(x_1) = \frac{\Gamma(\alpha_1+\alpha_2)}{\Gamma(\alpha_1)\Gamma(\alpha_2)}x_1^{\alpha_1-1}(1-x_1)^{\alpha_2-1}
\end{equation}

In general, the marginal distributions are: $X_i\sim Beta(\alpha_i,\alpha_0-\alpha_i) ,i=1,\ldots,P$, with $E(X_i)=\frac{\alpha_i}{\alpha_0}$. A more interesting property (used to derive the above result on the marginal distributions) is the aggregation property \citep[proved in][]{Frigyik2010}. It says that if two components of a Dirichlet vector are dropped and replaced with their sum then the resulting vector is again Dirichlet with the corresponding parameters replaced in the same way. More explicitly,
\begin{equation}\begin{aligned}
&\text{if }(X_1,\ldots,X_P) \sim D(\alpha_1,\ldots,\alpha_P)\\
&\text{then }(X_1,\ldots,X_i+X_j,\ldots,X_P) \sim D(\alpha_1,\ldots,\alpha_i+\alpha_j,\ldots,\alpha_P ).
\end{aligned}\end{equation}
Another way to interpret this property is that the joint marginal distributions are also Dirichlet, in the sense that if $(X_1,\ldots,X_{P-1},1-\sum_{i=1}^{P-1} X_i) \sim D(\alpha_1,\ldots,\alpha_{P-1},\alpha_{P})$, then $(X_1,\ldots,X_j,1-\sum_{i=1}^j X_i) \sim D(\alpha_1,\ldots,\alpha_j,\sum_{i=j+1}^{P} \alpha_i)$. This result becomes very useful later with regard to the conditional distribution.

The Dirichlet distribution is illustrated in \citet[127--128]{Kotz2000}. Datasets of composition data are given in \citet[354--405]{Aitchison1986} and are used to illustrate the methods in this paper.

\autoref{sec:fitting} discusses methods of estimating the parameters of the distribution. \autoref{sec:type2} discusses the Dirichlet Type II distribution. This is followed by examples (\autoref{sec:examples}) and a simulation study comparing parameter estimation methods (\autoref{sec:simstudy}). \autoref{sec:ppcd} explains the process of prediction using this distribution and \autoref{sec:imputation} shows how to do imputation using this distribution family.

\section{Fitting the Dirichlet distribution} \label{sec:fitting}

Let $\mathbf{x}_1,\ldots,\mathbf{x}_n$ be an observed sample from the Dirichlet Type I distribution. These vectors are combined as rows of a matrix, so that from here on $X_{ij}$ refers to the $j^{th}$ element of the $i^{th}$ observation. In general, $X_j$ where $j \in {1,\ldots,P}$ will refer to a specific margin considered as a random variable.

\subsection{Method of moments} \label{sec:mom}
\citet{Minka} points out that $\alpha_0=\frac{E(X_1 )-E(X_1^2 )}{E(X_1^2 )-E(X_1 )^2 }$. By multiplying each expected value by this expression estimates for each parameter are obtained. Explicitly:
\begin{equation} 
E(X_j )=\frac{\alpha_j}{\alpha_0} \Rightarrow \alpha_j=E(X_j ) \alpha_0=E(X_j )\left[\frac{E(X_1 )-E(X_1^2 )}{E(X_1^2 )-E(X_1 )^2 }\right]
\end{equation}

This method does not correspond exactly to the commonly used estimates for the beta case. An alternative approach is given in \citet{Ronning1989} but is not discussed here as the above method is sufficient for our purpose: using the method of moments as a starting point for the method of maximum likelihood.

\subsection{Method of maximum likelihood} \label{sec:ml}

Assuming an independent sample, the likelihood is obtained by multiplying densities. The resulting likelihood is
\begin{equation}\label{eq:D1likelihoodChp2}
Lik(\mathbf{x}_1,\ldots,\mathbf{x}_n|\boldsymbol\alpha) = \left[\frac{\Gamma(\alpha_0)}{\prod_{j=1}^{P}\Gamma(\alpha_j)}\right]^n\prod_{i=1}^n \left(\prod_{j=1}^{P}x_{ij}^{\alpha_{ij}-1}\right)
\end{equation}

\citet{Minka} suggests the following iteration for maximising the likelihood:
\begin{equation} \label{eq:MLinteration}
\psi(\alpha_j^{new} )=\psi(\alpha_0^{old} )+\frac{1}{n} \sum_{i=1}^n \ln{X_{ij}}  
\end{equation}
$\psi(x)$ refers to the digamma function, which occurs throughout this text as both $\frac{d}{dx}  \ln \Gamma(x)$ and $E[\ln G ]$ where $G\sim Gamma(x,1)$. This equivalence between the derivative of the log gamma function and the expectation of the log of a gamma random variable is useful in derivations that follow. Further, the derivative of the digamma function, denoted $\psi'(x)$, is called the trigamma function and is also used in the coming sections.

A detailed discussion of the methods available to maximise the likelihood is given in \citet{Huang2008}. The method above (\autoref{eq:MLinteration}) was found to converge quickly in all cases. 

\subsection{Posterior distribution} \label{sec:post}

In the Bayesian approach, the parameters of the Dirichlet distribution are considered random. One does not need to estimate specific values for the parameters in order to do Bayesian inference. % However, parameter estimates can be obtained from the posterior distribution under a given loss function (\textit{e.g.\ }the mean of the relevant posterior is the optimal parameter estimate under squared error loss).

First, an appropriate prior distribution must be specified.

\subsubsection{MDI Prior} \label{sec:mdi}

Consider first the maximal data information (MDI) prior \citep[41--53]{zellner1997}. The log of the MDI prior is the expected value of the log density. Since
\begin{equation} 
\ln{f(\mathbf{x})}=\ln {\Gamma(\alpha_0 )}-\sum_{j=1}^{P} \ln\Gamma(\alpha_j ) +\sum_{j=1}^{P}\left[(\alpha_j-1) \ln {x_j} \right]   ,
\end{equation}

the log prior evaluates to
\begin{equation} \label{eq:mdi}
 \begin{aligned}
 E[\ln f(\mathbf{X}) ]&=\ln \Gamma(\alpha_0 )-\sum_{j=1}^{P} \ln \Gamma(\alpha_j ) +\sum_{j=1}^{P}[(\alpha_j-1)E(\ln X_j )] \\ 
&=\ln \Gamma(\alpha_0 )-\sum_{j=1}^{P} \ln \Gamma(\alpha_j ) +\sum_{j=1}^{P} \left[\left(\alpha_j-1\right)\left(\psi(\alpha_j)-\psi(\alpha_0)\right)\right]. 
\end{aligned}
\end{equation}

The log likelihood is
\begin{equation} \label{eq:llike}
 \begin{aligned}
 \ell &=n \ln \Gamma(\alpha_0 )-n \sum_{j=1}^{P} \ln \Gamma(\alpha_j )+\sum_{i=1}^n \sum_{j=1}^{P}\left[\left(\alpha_j-1\right) \ln x_ij \right] \\ 
&=n \ln \Gamma(\alpha_0 )-n \sum_{j=1}^{P} \ln \Gamma(\alpha_j ) +\sum_{j=1}^{P} \left[\left(\alpha_j-1\right) \sum_{i=1}^n \ln x_{ij}\right]. 
\end{aligned}
\end{equation}

Adding \autoref{eq:mdi} and \autoref{eq:llike} then gives the log posterior:
\begin{equation} \label{eq:lpost}
 \begin{aligned}
 \ln \pi(\boldsymbol\alpha|X) &=(n+1) \ln \Gamma(\alpha_0 )-(n+1) \sum_{j=1}^{P} \ln \Gamma(\alpha_j ) \\
&\qquad+\sum_{j=1}^{P} \left\{\left(\alpha_j-1\right)\left[\left(\psi(\alpha_j)-\psi(\alpha_0)\right)+ \sum_{i=1}^n \ln x_{ij} \right]\right\} + c 
\end{aligned}
\end{equation}
where $c$ is an unknown constant (independent of the parameters).

\citet{Wicker2008} explain that the likelihood is globally concave. The addition of the prior does not appear to change this property, simplifying the process of finding the mode greatly, as the mode is simply the set of parameter values where the posterior reaches its peak. A simple gradient ascent algorithm can be employed, which only requires the first derivatives:
\begin{equation}
 \begin{aligned}
\frac{\partial \ln \pi}{\partial \alpha_j}&=(n+1)[\psi(\alpha_0 )-\psi(\alpha_j )]+\sum_{i=1}^n \ln  x_{ij} +(\psi(\alpha_j )-\psi(\alpha_0 ))\\
&\qquad +(\alpha_j-1) \psi' (\alpha_j )-\sum_{i=1}^{P}(\alpha_i-1) \psi' (\alpha_0 )\\
&=n[\psi(\alpha_0 )-\psi(\alpha_j )]+\sum_{i=1}^n \ln  x_{ij} +(\alpha_j-1) \psi' (\alpha_j ) - \psi' (\alpha_0 )(\alpha_0-P).
\end{aligned}
\end{equation}

One can simulate a sample from the posterior (\autoref{eq:lpost}) using the Metropolis-Hastings (MH) algorithm. See \citet[267--301]{Robert2004} for an in-depth discussion of this algorithm.

The MH algorithm requires the specification of a suitable jump distribution for choosing candidate values. The multivariate normal distribution is an obvious choice, \textit{i.e.\ }${\boldsymbol\alpha}^c \sim N_{P} ({\boldsymbol\alpha}^i,\Sigma)$ where $\Sigma$ is a diagonal matrix of constants that affect the acceptance rate and required burn-in period. As is standard practice, these constants are chosen by crudely matching the peak of the posterior to the peak of a multivariate normal distribution.

In this way the algorithm is simplified to:

Accept candidate ${\boldsymbol\alpha}^c$ as a new observation ${\boldsymbol\alpha}^{i+1}$ from the posterior if and only if

$\ln \pi({\boldsymbol\alpha}^c |\mathbf{X})-\ln \pi({\boldsymbol\alpha}^i |\mathbf{X})>\ln u$ where $u$ is drawn randomly from the standard uniform distribution.

It is worth noting that the above method works very well for any finite number of dimensions, including the beta case $(P=2)$.

\subsubsection{Jeffreys Prior} \label{sec:jeffreys}

An alternative to the MDI prior is the Jeffreys prior \citep{jeffreys1998} which is derived as follows:

Let $g=-\ln  f(\mathbf{x})=-\ln \Gamma(\alpha_0 )+\sum_{j=1}^{P} \ln \Gamma(\alpha_j )
-\sum_{j=1}^{P} [(\alpha_j-1)  \ln  x_j ] $. Then
\begin{equation}
 \begin{aligned}
&\frac{\partial g}{\partial \alpha_j }=-\psi(\alpha_0 )+\psi(\alpha_j )-\ln x_j ,\\
&\frac{\partial^2 g}{\partial \alpha_j^2 }=-\psi'(\alpha_0 )+\psi'(\alpha_j ),\\
&\frac{\partial^2 g}{\partial \alpha_j \partial \alpha_i }=-\psi'(\alpha_0 ),
\end{aligned}
\end{equation}
and thus
\begin{equation}
\begin{aligned}
&\text{Jeffreys Prior} \propto \\
&\left|\begin{matrix}
       \psi'(\alpha_1 )-\psi'(\alpha_0 ) & -\psi'(\alpha_0 ) & \cdots & -\psi'(\alpha_0 )         \\
       -\psi'(\alpha_0 ) & \psi'(\alpha_2 )-\psi'(\alpha_0 )  &   & \vdots \\
       \vdots    &  & \ddots & \vdots\\
	-\psi'(\alpha_0) & \cdots & \cdots & \psi'(\alpha_{P} )-\psi'(\alpha_0 )
     \end{matrix}\right|^{\frac{1}{2}}.
\end{aligned}
\end{equation}

In the beta case $(P=2)$:
\begin{equation}
\sqrt{\psi'(\alpha_1 ) \psi'(\alpha_2 )-\psi'(\alpha_0 )\left[\psi'(\alpha_1 )+\psi'(\alpha_2 )\right]}
\end{equation}

When $P=3$:
\begin{equation}
\left\{\psi'(\alpha_1 ) \psi'(\alpha_2 ) \psi'(\alpha_3 ) - \psi'(\alpha_0)\times\left[\psi'(\alpha_1 )\psi'(\alpha_2 ) + \psi'(\alpha_1 )\psi' (\alpha_3 )\psi'(\alpha_2 )\psi'(\alpha_3 )\right]\right\}^{\frac{1}{2}}
\end{equation}

In general the Jeffreys prior is proportional to the square root of
\begin{equation}
 \begin{aligned}
&\prod_{l=1}^{P} \psi'(\alpha_l ) -\psi'(\alpha_0 ) \sum_{l=1}^{P} \prod_{\substack{j=1\\j \ne l}}^{P}\psi' (\alpha_j )\\
&=\prod_{l=1}^{P}\psi' (\alpha_l )\left[1-\psi'(\alpha_0 ) \sum_{l=1}^{P} \frac{1}{\psi'(\alpha_l )}\right]
\end{aligned}
\end{equation}

Thus, the log Jeffreys prior (ignoring any constants for now) is given by
\begin{equation}
 \begin{aligned}
&0.5 \ln \prod_{l=1}^{P}\psi' (\alpha_l )\left[1-\psi'(\alpha_0 ) \sum_{l=1}^{P} \frac{1}{\psi'(\alpha_l )}\right]\\
&=0.5 \sum_{l=1}^{P} \ln \psi' (\alpha_l )+0.5 \ln\left[1-\psi'(\alpha_0 ) \sum_{l=1}^{P} \frac{1}{\psi'(\alpha_l )}\right]
\end{aligned}
\end{equation}

so that the log posterior is equal to
\begin{equation}
 \begin{aligned}
0.5 \sum_{j=1}^{P} \ln \psi' (\alpha_j ) &+0.5 \ln \left[ 1 - \psi'(\alpha_0 ) \sum_{j=1}^{P} \frac{1}{\psi'(\alpha_j )} \right] + n \ln \Gamma(\alpha_0 )\\
&-n\sum_{j=1}^{P} \ln \Gamma(\alpha_j ) + \sum_{j=1}^{P}\left[(\alpha_j-1) \sum_{i=1}^n \ln  x_{ij}\right] + c.
\end{aligned}
\end{equation}

The first derivatives (needed to find the mode) are found to be:
\begin{equation}
 \begin{aligned}
\frac{\partial \ln  \pi}{\partial \alpha_j } =&n[\psi(\alpha_0 )-\psi(\alpha_j )]+\sum_{i=1}^n \ln  x_{ij} + 0.5 \frac{\psi''(\alpha_j )}{\psi'(\alpha_j ) }\\
&-0.5\left[\psi'' (\alpha_0 ) \sum_{l=1}^{P} \frac{1}{\psi'(\alpha_l ) }-\psi' (\alpha_0 ) \psi''(\alpha_j ) (\psi'(\alpha_j ))^{-2}\right]\\
&\times\left[1-\psi' (\alpha_0 ) \sum_{l=1}^{P} \frac{1}{\psi'(\alpha_l ) }\right]^{-1}.
\end{aligned}
\end{equation}

The method of simulation is the same as for the MDI prior (using the Metropolis algorithm).

For an overview of the derivation of objective priors and simpler examples see \citet{Priors}.

\section{The Dirichlet Type II Distribution} \label{sec:type2}

In univariate analysis, the beta distribution is popular as it is a conjugate prior for the Bernoulli parameter $p$. Should the researcher wish to analyse the odds of success rather than the probability of success, then a conjugate prior is the beta prime (also called beta type II) distribution \parencite{broderick2014}. It is a transform of the beta distribution \parencite{johnson1970}, according to the following:
\begin{equation}
\text{If }A\sim Beta(a,b)\text{ then }B=\frac{A}{1-A}\sim Beta'(a,b) \text{ and }C=\frac{b}{a}B \sim F(2a,2b).
\end{equation}

A similar transformation can be done in the multivariate case. This transformation was first introduced by \textcite{tiao1965}. Consider again $\mathbf{X} \sim Dirichlet(\alpha_1,\ldots,\alpha_{P} )$. To be clear, $\mathbf{X}$ refers to an unobserved Dirichlet vector and $X_j$ to a specific element of this random vector. Transforming this distribution using
\begin{equation} \label{eq:type2transform}
 Y_j=\frac{X_j}{X_{P}}   ,j=1,\ldots,(P-1) \ \text{noting that}\ X_{P}=1-\sum_{l=1}^{P-1} X_l 
\end{equation}
yields $\mathbf{Y}$ (size $P-1$ by 1) $\sim Dirichlet\ Type\ II(\alpha_1,\ldots,\alpha_{P} )$.

% The above transformation (\autoref{eq:type2transform}}) is called the additive logratio transform (ALR).

In the above definition the last component is used as the reference component in the transformation as a matter of convenience. While a reference component is required for this transformation, it does not have to be the last component. In practical problems, careful consideration should be given to choosing the reference component. In the absence of natural reference component, a rule of thumb might be to choose the component with the largest minimum, resulting in the most stable transformation.

According to \textcite{tiao1965}, the Dirichlet Type II Distribution also arises as follows: \\
Starting with $P\ \chi^2$ random variables with $2\alpha_1,2\alpha_2,\ldots,2\alpha_P$ degrees of freedom respectively, say $Q_1,\ldots,Q_p$, then the ratios $Y_j=\frac{Q_j}{Q_{P}}   ,j=1,\ldots,(P-1)$ also form a $Dirichlet\ Type\ II(\alpha_1,\ldots,\alpha_{P} )$ distribution.

The density is now
\begin{equation}
 \begin{aligned}
f(\mathbf{y})=&\frac{\Gamma(\alpha_0 )}{\prod_{j=1}^{P}\Gamma(\alpha_j ) } (\prod_{j=1}^{P-1} y_j^{\alpha_j-1} ) (1+y_1+\cdots+y_{P-1} )^{-\alpha_0 }\\
&  \text{where}\ y_j>0 ,j=1,\ldots,(P-1).
\end{aligned}
\end{equation}
and the log density:
\begin{equation}
 \begin{aligned}
\ln f(\mathbf{y})=&\ln\Gamma(\alpha_0 )-\sum_{j=1}^{P}\ln\Gamma(\alpha_j ) - \alpha_0  \ln (1+\sum_{j=1}^{P-1} y_j )\\
&+\sum_{j=1}^{P-1} [(\alpha_j-1)  \ln y_j] 
\end{aligned}
\end{equation}

Thus, given a sample $(\mathbf{y}_1,\ldots,\mathbf{y}_n )$, the log likelihood is
\begin{equation}
 \begin{aligned}
\ell=&n \ln \Gamma(\alpha_0 )-n\sum_{j=1}^{P} \ln\Gamma(\alpha_j ) -\alpha_0 \sum_{i=1}^n \ln \left(1+\sum_{j=1}^{P-1} y_{ij} \right)\\
&+\sum_{i=1}^n \sum_{j=1}^{P-1} [(\alpha_j-1)  \ln  y_{ij} ] .
\end{aligned}
\end{equation}

Consider again the MDI prior. The log prior is:
\begin{equation}
 \begin{aligned}
E[\ln f(\mathbf{Y})]=&\ln \Gamma(\alpha_0 )-\sum_{j=1}^{P}\ln \Gamma(\alpha_j ) - \alpha_0 E\left[\ln \left(1+\sum_{j=1}^{P-1} Y_j \right)\right]\\
&+\sum_{j=1}^{P-1} [(\alpha_j-1)E(\ln Y_j)] ,
\end{aligned}
\end{equation}
which can be simplified as follows:

With $X_j \sim Beta(\alpha_j,\alpha_0-\alpha_j) \text{ and } X_P \sim Beta(\alpha_P,\alpha_0-\alpha_P)$, then
\begin{equation}
 \begin{aligned}
E[\ln  Y_j ]&=E[\ln X_j -\ln X_P ] \\
&=\psi(\alpha_j )-\psi(\alpha_P )
\end{aligned}
\end{equation}
and
\begin{equation}
 \begin{aligned}
E\left[\ln \left(1+\sum_{i=1}^{P-1} Y_i \right) \right] &= -E\left[\ln \left(1-\sum_{i=1}^{P-1} X_i \right) \right]\\
&=-E[\ln X_{P}] \\ 
&=\psi(\alpha_0 )-\psi(\alpha_{P} ).
\end{aligned}
\end{equation}

The log prior thus simplifies to 
\begin{equation}
 \begin{aligned}
\ln \Gamma(\alpha_0 )&-\sum_{j=1}^{P} \ln \Gamma(\alpha_j ) - \alpha_0 [\psi(\alpha_0 )-\psi(\alpha_{P} )]\\
&+\sum_{j=1}^{P-1} (\alpha_j-1)[\psi(\alpha_j )-\psi(\alpha_P )] 
\end{aligned}
\end{equation}
and the log posterior to
\begin{equation}
 \begin{aligned}
\ln \pi(\boldsymbol\alpha|\mathbf{Y})=&(n+1)\left[\ln \Gamma(\alpha_0 )-\sum_{j=1}^{P} \ln \Gamma(\alpha_j ) \right]\\
&-\alpha_0 \left\{\left[\psi(\alpha_0 )-\psi(\alpha_{P} )\right]+\sum_{i=1}^n \ln \left(1+\sum_{j=1}^{P-1} y_{ij} \right) \right\}\\
&+\sum_{j=1}^{P-1} (\alpha_j-1)[\psi(\alpha_j )-\psi(\alpha_P)] \\
&+\sum_{i=1}^n \sum_{j=1}^{P-1} [(\alpha_j-1)  \ln  y_{ij} ] 
\end{aligned}
\end{equation}

While it is possible to work directly with the above distribution, including simulating from the posterior as before, it is generally preferred to transform a sample from the Dirichlet Type II distribution to a sample from the Dirichlet Type I distribution and apply the methods discussed in \autoref{sec:fitting}. The transformation is the inverse of the previous transformation (\autoref{eq:type2transform}), \textit{i.e.\ } 

\begin{equation}
 X_j=\frac{Y_j}{1+\sum_{l=1}^{P-1} Y_l }  ,j=1,\ldots,(P-1).
\end{equation}

% The above transformation is called the inverse ALR. This approach is illustrated in \autoref{sec:examples}.

\section{Examples} \label{sec:examples}

As an example of composition data that should fit the assumptions of the Dirichlet Type I distribution, consider the hair colours of boys in different counties of Scotland. This dataset is given in \citet[p. 405]{Aitchison1986}. What identifies this data as Dirichlet is that no additional relationships are expected between hair colours beyond that which is induced by the constraint that the proportions sum to one. This expectation does not hold for the Dark and Black hair colours as these seem to have a strong positive relationship, hence consider Black hair as part of the Dark hair category prior to applying Dirichlet modelling. After this adjustment the expectation appears to hold.

This paper assesses the model fit in three ways. The first is to simulate a replicate sample and compare it to the data visually. This is done in \autoref{fig:hair}. The fit seems reasonable in all 4 dimensions.

\begin{figure}
\begin{center}
\includegraphics[width=\figwidth]{shadedballs.png}
\end{center}
\caption{Hair colour data with simulated replicates}

Dark balls are the observed values and the lighter balls are the simulated values. \label{fig:hair}
\end{figure}

The second is to compare the observed and simulated covariance structures (\autoref{tbl:covhair}). These seem remarkably similar, especially when one considers that the Dirichlet has no parameters to explicitly model dependence. The only difference is the apparent positive correspondence between the Fair and Red hair colours, which the Dirichlet model cannot capture.

\begin{table}
\caption{Correlation structure of observed and simulated data}
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
\textbf{Observed} & Fair  & Red   & Medium & Dark \\
\hline
Fair  & 1     & 0.22  & -0.61 & -0.41 \\
\hline
Red   & 0.22  & 1     & -0.11 & -0.33\\
\hline
Medium & -0.61 & -0.11 & 1     & -0.45 \\
\hline
Dark  & -0.41 & -0.33 & -0.45 & 1 \\
\hline
\textbf{Simulated} & Fair  & Red   & Medium & Dark \\
\hline
Fair  & 1     & -0.14 & -0.52 & -0.39 \\
\hline
Red   & -0.14 & 1     & -0.2  & -0.14 \\
\hline
Medium & -0.52 & -0.2  & 1     & -0.47 \\
\hline
Dark  & -0.39 & -0.14 & -0.47 & 1 \\
\hline
\end{tabular}
   \label{tbl:covhair}
\end{center}
\end{table}

Lastly, consider the QQ-plots of the marginal fits (\autoref{fig:qqmargins}). Again the fit seems reasonable except for the Red hair colour.

\begin{figure}
\begin{center}
\includegraphics[width=\figwidth]{QQplot.png}
\end{center}
\caption{QQ Plots for margins of hair colour data}

\label{fig:qqmargins}
\end{figure}


As an example of Dirichlet Type II fitting, consider Fisher’s Iris flower data. This is a well-known dataset with 4 variables and three populations of size 50 each. It contains observations on the four variables $Y_1$= sepal length, $Y_2$= sepal width, $Y_3$= petal length and $Y_4$= petal width of the species Setosa, Versicolor and Virginica. This dataset is available in the R package \citep{RCore}.

First assume that the distribution of a set of measurements
\begin{equation}
 \begin{aligned}
&(X_{i1} X_{i2},X_{i3},X_{i4} )=\\
&\left(\frac{Y_{i1}}{1+\sum_{j=1}^4 Y_{ij} },\frac{Y_{i2}}{1+\sum_{j=1}^4 Y_{ij} },\frac{Y_{i3}}{1+\sum_{j=1}^4 Y_{ij} },\frac{Y_{i4}}{1+\sum_{j=1}^4 Y_{ij} }\right)
\end{aligned}
\end{equation}
can be described through a Dirichlet distribution $D(\alpha_1,\ldots,\alpha_5)$ for each population, say $\mathbf{X}^{(p)}\sim\\ D\left(\alpha^{(p)}  =(\alpha_j^{(p)},j=1,\ldots,5)\right),\ p=1,2,3$.

Now proceed to fit Dirichlet distributions on the four variables for each species by transforming as above. Begin by applying the method of moments to obtain initial parameter estimates, and then apply numeric optimisation to the posterior distribution using the MDI prior (\autoref{eq:lpost}), and obtain at the following estimates of the Dirichlet parameters:
\begin{equation}
 \begin{aligned}
&\hat{\alpha}^{(SET)}  = {604.71,412.57,175.49,27.78,120.02}\\
&\hat{\alpha}^{(VER)}  = {378.31,176.39,270.87,84.17,64.15}\\
&\hat{\alpha}^{(VIR)}  = {427.61,193.26,360.36,131.62,65.37}.
\end{aligned}
\end{equation}

A basic test of goodness-of-fit is to simulate replicate samples and compare them to the original data visually \citep[see for example][159--171]{gelman2013}. Doing this with the Iris data produces graphs such as \autoref{fig:irisballs}. There appears to be sufficient correspondence between the observed data and the simulated data. The correspondence is in all four dimensions. More advanced goodness-of-fit procedures do not yet appear to exist for the Dirichlet distribution.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.6\figwidth]{setosaD2.png}\includegraphics[width=0.6\figwidth]{versicolorD2.png}

\includegraphics[width=0.6\figwidth]{virginicaD2.png}
\end{center}
\caption{Iris data with simulated replicate samples on Dirichlet type II scale}

Green indicates the original data and blue a simulated replicate sample. \label{fig:irisballs}
\end{figure}

All methods were implemented in the statistical package R \citep{RCore}. Packages used include the `parallel' package for increasing the speed of the simulation study in \autoref{sec:simstudy} \citep{RCore} and the `rgl' package for drawing the graphs in Figures \ref{fig:hair}, \ref{fig:irisballs} and \ref{fig:bihist} \citep{RGL}.

\section{Simulation Study} \label{sec:simstudy}

It is of interest to compare the accuracy of the different approaches with respect to estimating the parameters of a Dirichlet distribution. To this end, this paper simulates 5000 samples from this distribution for each of a variety of combinations of sample size, dimension and parameter values. For each sample, this paper estimates the parameters in six ways and compare the estimates to the true values. The discrepancies are summarised using the root mean square percentage error (RMSPE) in order to attempt to compare accuracy across varying parameter values.

Simulation of a value from the Dirichlet distribution can be done in at least three ways. The simplest method (used in this study) is to consider a Dirichlet observation as the relative sizes of gamma draws, \textit{i.e.\ }draw $y_j\sim Gamma(\alpha_j,1)$ for $j=1,\ldots,P$ and then set $x_j=\frac{y_j}{\sum_{j=1}^{P}y_j }$. This gives a draw $ \mathbf{x}\sim D(\alpha_1,\ldots,\alpha_{P} )$.

A second method relies on the fact that the univariate marginal and conditional distributions are beta and simulates values one dimension at a time conditional on the previous dimensions. A third method by \citet{Matsui2010} uses a path coupling approach. These methods are interesting in the theory of simulation, but not relevant to the current study.

The results of the simulation study are given in \autoref{tbl:simdir}.

\begin{table}
  \caption{Performance of parameter estimation methods as the sample size, parameter magnitude and number of dimensions are varied}
 \footnotesize
\begin{center}
\begin{tabular}{|l|l|l|l|l|l|l|l|}
    \hline
    \textbf{No.} & \textbf{RMSPE} & \textbf{Mode -} & \textbf{ML} & \textbf{Mean -} & \textbf{Mode -} & \textbf{Mean -} & \textbf{MOM} \\
     &  & \textbf{Jeffreys} &  & \textbf{Jeffreys} & \textbf{MDI} & \textbf{MDI} &  \\     \hline
    1     & \textbf{n=25} & 21.90\% & 25.20\% & 25.90\% & 27.20\% & 34.10\% & 38.60\% \\
        & $\alpha=(3,3,3)$ &  &  &  &  &  &  \\    \hline
    2     & \textbf{n=25} & 16.80\% & 19.00\% & 19.70\% & 20.90\% & 26.00\% & 43.00\% \\
        & $\alpha=(3,3,3,3,3)$ &  &  &  &  &  &  \\  \hline
    3     & \textbf{n=50} & 15.10\% & 16.30\% & 16.60\% & 17.00\% & 19.80\% & 23.90\% \\
        & $\alpha=(3,3,3)$ &  &  &  &  &  &  \\ \hline
    4     & \textbf{n=50} & 11.70\% & 12.40\% & 12.70\% & 13.00\% & 15.00\% & 25.70\% \\
        & $\alpha=(3,3,3,3,3)$ &  &  &  &  &  &  \\  \hline
    5     & \textbf{n=90} & 11.10\% & 11.60\% & 11.80\% & 11.90\% & 13.20\% & 17.10\% \\
        & $\alpha=(3,3,3)$ &  &  &  &  &  &  \\  \hline
    6     & \textbf{n=90} & 8.70\% & 9.10\% & 9.20\% & 9.30\% & 10.20\% & 18.30\% \\
        & $\alpha=(3,3,3,3,3)$ &  &  &  &  &  &  \\  \hline
    7     & \textbf{n=120} & 6.50\% & 6.70\% & 6.80\% & 6.90\% & 7.40\% & 16.30\% \\
        & $\alpha=(3,3,3,3,3,3,3,3)$ &  &  &  &  &  &  \\  \hline
    8     & \textbf{n=25} & 21.20\% & 24.50\% & 25.30\% & 26.70\% & 34.20\% & 38.20\% \\
        & $\alpha=(4,4,4)$ &  &  &  &  &  &  \\      \hline
    9     & \textbf{n=25} & 35.30\% & 31.50\% & 30.10\% & 36.00\% & 34.90\% & 38.90\% \\
        & $\alpha=(100,100,100)$ &  &  &  &  &  &  \\  \hline
    10    & \textbf{n=25} & 31.60\% & 27.40\% & 25.70\% & 33.60\% & 34.20\% & 41.60\%\\
        & $\alpha=(4,20,100)$ &  &  &  &  &  &  \\  \hline
    11    & \textbf{n=50} & 15.10\% & 16.30\% & 16.80\% & 17.10\% & 20.00\% & 24.00\% \\
        & $\alpha=(4,4,4)$ &  &  &  &  &  &  \\  \hline
    12    & \textbf{n=50} & 19.30\% & 19.20\% & 18.90\% & 19.90\% & 21.00\% & 23.60\% \\
        & $\alpha=(100,100,100)$ &  &  &  &  &  &  \\  \hline
    13    & \textbf{n=50} & 16.70\% & 17.20\% & 17.10\% & 18.10\% & 20.00\% & 26.10\% \\
        & $\alpha=(4,20,100)$ &  &  &  &  &  &  \\  \hline
    14    & \textbf{n=90} & 10.90\% & 11.40\% & 11.60\% & 11.70\% & 12.90\% & 16.40\% \\
        & $\alpha=(4,4,4)$ &  &  &  &  &  &  \\  \hline
    15    & \textbf{n=90} & 12.00\% & 13.20\% & 13.10\% & 12.40\% & 14.20\% & 16.40\% \\
        & $\alpha=(100,100,100)$ &  &  &  &  &  &  \\  \hline
    16    & \textbf{n=90} & 11.10\% & 12.00\% & 12.10\% & 11.80\% & 13.40\% & 18.40\%\\
        & $\alpha=(4,20,100)$ &  &  &  &  &  &  \\  \hline
    \end{tabular}
\end{center}

\normalsize
\vspace{2mm}
The effect of sample size can be seen by moving from row 1 to 7 or 8 to 16. The effect of parameter magnitude is evident in rows 8 vs 9 vs 10, 11 vs 12 vs 13 or 14 vs 15 vs 16. The effect of the number of dimensions can be seen by comparing rows 1 vs 2, 3 vs 4 or 5 vs 6 vs 7. Lower values indicate better estimation.
 \label{tbl:simdir}
\end{table}

There are many things worth noting in \autoref{tbl:simdir}. The first is that in almost all cases the errors increase moving from left to right in the table. This means that, in general, the posterior mode using the Jeffreys prior appears to perform the best, followed by the maximum likelihood estimator.

The second is that as some or all parameter values become large then the mean of the posterior simulations (where the mean is taken on each dimension individually) performs better. Thus, if another method is applied and the suggested parameter values are large then this method is recommended.  Special care must be taken when parameter values become small.

The third is that the accuracy of all estimates appears to be heavily dependent on sample size. Increased sample size can be expected to yield increased accuracy but the extent to which it does so is remarkable.

The fourth is that increasing the dimension appears to increase the accuracy of the estimates. Increasing the dimension while keeping the sample size constant increases the total information available from which to estimate parameters. The increase in the information is greater than the increase in parameters, in a relative sense.

The performance of these methods for the beta distribution are specifically investigated, varying the sample size and parameter values. The results are given in \autoref{tbl:simbeta}.

\begin{table}[hbt]
\caption{Performance of parameter estimation methods in the case of the beta distribution}
\begin{center}
\begin{small}
\begin{tabular}{|l|l|l|l|l|l|l|l|}
    \hline
    \textbf{No.} & \textbf{RMSPE} & \textbf{Mode -} & \textbf{ML} & \textbf{Mean -} & \textbf{Mode -} & \textbf{Mean -} & \textbf{MOM} \\
     &  & \textbf{Jeffreys} &  & \textbf{Jeffreys} & \textbf{MDI} & \textbf{MDI} &  \\      \hline
    1     & \textbf{n=60} & 18.10\% & 19.50\% & 20.40\% & 20.10\% & 23.10\% & 19.60\% \\
    & $\alpha=(6,6)$ &  &  &  &  &  &  \\      \hline
    2     & \textbf{n=60} & 21.10\% & 21.90\% & 24.10\% & 22.40\% & 26.70\% & 21.90\% \\
        & $\alpha=(60,60)$ &  &  &  &  &  &  \\      \hline
    3     & \textbf{n=60} & 20.70\% & 20.70\% & 22.30\% & 20.70\% & 23.10\% & 20.70\% \\
        & $\alpha=(600,600)$ &  &  &  &  &  &  \\      \hline
    4     & \textbf{n=60} & 19.80\% & 21.30\% & 22.80\% & 21.80\% & 26.00\% & 22.60\% \\
        & $\alpha=(6,60)$ &  &  &  &  &  &  \\      \hline
    5     & \textbf{n=25} & 30.50\% & 36.60\% & 38.30\% & 39.00\% & 50.00\% & 36.60\% \\
        & $\alpha=(3,3)$ &  &  &  &  &  &  \\    \hline
    6     & \textbf{n=50} & 20.80\% & 23.00\% & 23.30\% & 23.80\% & 28.40\% & 23.10\% \\
        & $\alpha=(3,3)$ &  &  &  &  &  &  \\    \hline
    7     & \textbf{n=75} & 16.10\% & 17.10\% & 17.50\% & 17.60\% & 20.10\% & 17.20\% \\
        & $\alpha=(3,3)$ &  &  &  &  &  &  \\    \hline
    8     & \textbf{n=100} & 14.30\% & 15.00\% & 15.20\% & 15.30\% & 17.00\% & 15.30\% \\
        & $\alpha=(3,3)$ &  &  &  &  &  &  \\    \hline
    \end{tabular}
\end{small}    
\end{center}

\vspace{1mm}
Parameter magnitude is varied over rows 1 to 4, while sample size is varied over rows 5 to 8.
Lower values indicate better estimation.
  \label{tbl:simbeta}
\end{table}

The patterns are the same as was observed in \autoref{tbl:simdir} --- in all rows it can be seen that the posterior mode using the Jeffreys prior produced the most accurate estimates. 

It should be noted that neither simulation analysis was extended to cover the case where the parameters are smaller than one. In that situation the shape of the posterior distribution may change and force a change in the simulation approach. Parameter values smaller than one correspond with extreme variance, with observations more spread out than one would expect from a uniform distribution. Extra care must be taken if such data are encountered in practice.

\section{Predictive posterior conditional distribution} \label{sec:ppcd}

Consider the situation where a number of components are jointly distributed Dirichlet and the researcher wishes to predict the values of a subset of components, given values for the remainder of the components. This can be done using the predictive posterior conditional distribution.

\subsection{Conditional distribution}
\begin{sloppypar}Initially assume that the parameter values are known, \textit{i.e.\ } $(X_1,\ldots,X_P )|\boldsymbol\alpha \sim D(\alpha_1,\ldots,\alpha_{P} )$ and there is interest in the distribution of $\left(X_1,\ldots,X_{P-1} \middle| X_{P} \right)$. Now condition on the last component for convenience only (the result holds for any component). The key result that the last component follows a $Beta(\alpha_P,\sum_{j=1}^{P-1}\alpha_j)$ marginal distribution, as proved in \textcite{Frigyik2010} is used.\end{sloppypar}

\begin{equation}
 \begin{aligned}
f(x_1,\ldots,x_{P-1}|X_P=x_P )&=\frac{f(x_1,\ldots,x_P )}{f(x_P)} \\
&=\frac{\frac{\Gamma\left(\sum_{j=1}^{P} \alpha_j \right)}{\prod_{j=1}^{P} \Gamma(\alpha_j ) } \left(\prod_{j=1}^P x_j^{\alpha_j-1} \right)}{\frac{\Gamma\left(\sum_{j=1}^{P} \alpha_j\right) }{\Gamma(\alpha_P ) \Gamma\left(\sum_{j=1}^{P-1} \alpha_j\right)} \left( x_P^{\alpha_P-1} (1-x_P)^{\left(\sum_{\alpha=1}^{P-1} \alpha_j\right) -1} \right)}\\
&=\frac{\Gamma\left(\sum_{j=1}^{P-1} \alpha_j\right)}{\prod_{j=1}^{P-1} \Gamma(\alpha_j ) } \left(\prod_{j=1}^{P-1} \left(\frac{x_j}{1-x_P}\right)^{\alpha_j-1} \right) (1-x_P)^{1-P}
\end{aligned}
\end{equation}

Set $W_j=\frac{X_j}{1-x_P}$, where $j=1,\ldots,(P-1)$ and $x_P$ is a constant in this situation, then $f(w_1,\ldots,w_{P-1}|X_{P}=x_{P} )$ becomes
\begin{equation}
\frac{\Gamma\left(\sum_{j=1}^{P-1} \alpha_j\right)}{\prod_{j=1}^{P-1} \Gamma(\alpha_j ) } \left(\prod_{j=1}^{P-1} \left(\frac{w_j}{1-x_P}\right)^{\alpha_j-1} \right)
\end{equation}
since the absolute Jacobian is $(1-x_P)^{P-1}$.

\begin{sloppypar}From this expression, $(W_1,\ldots,W_{P-1} |X_{P}=x_{P} ) \sim D(\alpha_1,\ldots,\alpha_{P-1} )$.\end{sloppypar}

The result above can be applied recursively, each time eliminating one component by dropping its parameter and rescaling the remaining components by their sum to ensure they again sum to one, as required by the Dirichlet distribution.

Thus, given a set of parameters and a set of known values for some components one can calculate quantities of interest. For example, one can simulate from the conditional distribution to calculate prediction intervals.

\subsection{Implementation of the predictive posterior conditional distribution}
In practise the parameters are unknown and must be estimated from data. This carries uncertainty across to the predictions, uncertainty that must be taken into account. A way to adjust for the uncertainty is using the predictive posterior conditional distribution.

Practically, the process for a single vector observation works as follows:
\begin{enumerate}
 \item Given a complete data matrix $X^{(c)}$ of size $n$ by $P$, simulate a large number of sets of parameters from the posterior distribution. Specifically, simulate $T$ sets and store them in a $T$ by $P$ parameter matrix $A$.
 \item Now consider a new matrix $X^{(n)}$ of size $m$ by $P$, where $Q$ components require prediction $(1\leq Q <P)$. For convenience, arrange the components of the new matrix, as well as the simulated parameter matrix, such that the $Q$ components are grouped together on the left.
 \item For each individual vector of simulated parameters $({\boldsymbol\alpha}_t)$, simulate a set of $w_1,\ldots,w_{Q}$ from a $D(A_{t1},\ldots,A_{tQ})$, where $t=1,\ldots,T$.
 \item Multiply each simulated set of $w_1,\ldots,w_{Q}$ by $\prod_{j=Q+1}^P (1-X^{(n)}_j)$.
 \item Combine the predictions over simulations  $t=1,\ldots,T$ to calculate values of interest, since these simulations constitute an approximation of the predictive posterior conditional distribution.
\end{enumerate}

\subsection{Conditional distribution in the Type II case}

In the case where the given values are on the Dirichlet Type II scale, a different approach is followed where the aggregation property is applied directly:

It is known that $ \left(X_1,\ldots,X_j,\left[1-\sum_{l=1}^j X_l \right]\right)  \sim D\left(\alpha_1,\ldots,\alpha_j,\sum_{l=j+1}^{P} \alpha_l \right)$ and that $X_j=\frac{Y_j}{1+\sum_{l=1}^{P-1} Y_l }$, where $\left(Y_1,\ldots,Y_{P-1}\right) \sim D2(\alpha_1,\ldots,\alpha_{P})$.

Thus, given $(Y_{j+1},\ldots,Y_{P-1} )$,
\begin{equation}
\frac{\frac{(Y_1,\ldots,Y_j )}{1+\sum_{l=1}^{P-1} Y_l }}{1-\sum_{i=1}^j \frac{Y_i}{1+\sum_{l=1}^{P-1} Y_l } } \sim D2\left(\alpha_1,\ldots,\alpha_j,\sum_{l=j+1}^{P} \alpha_l \right).	
\end{equation}

This expression simplifies to
\begin{equation}
\frac{(Y_1,\ldots,Y_j )}{1+\sum_{l=j+1}^{P-1} Y_l }  ,
\end{equation}
so, if one defines $\beta=\left(1+\sum_{l=j+1}^{P-1} Y_l \right)$ and $ V_i=\frac{Y_i}{\beta} $then
\begin{equation}
(V_1,\ldots,V_j |Y_{j+1}=y_{j+1},\ldots,Y_{P-1} = y_{P-1}) \sim D2\left(\alpha_1,\ldots,\alpha_j,\sum_{l=j+1}^{P} \alpha_l \right).
\end{equation}

Also define $U_i=\frac{V_i}{1+\sum_{l=1}^j V_l }$ and note that
\begin{equation}
\left(U_1,\ldots,U_j,1-\sum_{l=1}^j U_l \middle|Y_{j+1}=y_{j+1},\ldots,Y_{P-1}=y_{P-1}\right) \sim D\left(\alpha_1,\ldots,\alpha_j,\sum_{l=j+1}^{P} \alpha_l \right).
\end{equation}

Use either of these and appropriately modify steps 2 and 3 in the given process to calculate desired quantities. For example, if one simulates a set of $u_1,\ldots,u_j$ then note that $y_l=\frac{\beta u_l}{1-\sum_{h=1}^j u_h }$.

As an example, let us again consider the Iris data from \autoref{sec:examples}. Suppose a new flower is measured from the Versicolor family with a sepal length of 6 and a sepal width of 3, but it is not feasible to measure the petals. Apply the above process with $T=9999$ to find that the predictive posterior conditional distribution has the form depicted in \autoref{fig:bihist}. By taking means, this process predicts a petal length of 4.38 and a petal width of 1.36, which seems reasonable looking at \autoref{fig:irisballs}.

\begin{figure}[hbt]
\begin{center}
\includegraphics[width=\figwidth]{bihist.png}
\end{center}
\caption{Bivariate histogram to illustrate the predictive posterior conditional distribution}

Distribution of the petal size given that the sepal length and width are 6 and 3 respectively, for the Versicolor species. \label{fig:bihist}
\end{figure}

\section{Multiple Imputation of Dirichlet Type II Data} \label{sec:imputation}

As an additional example, an application of the methods discussed is given. Consider a set of monthly total rainfall figures collected at five stations in central South Africa. Only consider months where measurable amounts of rain fell. There are 37 months where all five stations produced measurements and 58 months where a station’s data are missing and needs to be imputed.

While many imputations methods exist, it is generally agreed that a flexible and efficient solution is multiple imputation. For an overview of multiple imputation see \citet{Schafer1998}. 

Assume that the rainfall measurements follow a Dirichlet Type II distribution with unknown parameters. Proceed through the following steps:

\begin{enumerate}
  \item Estimate the six parameters of the distribution using only the 37 complete months.
  \item Pick random values for all the missing measurements from the respective predictive posterior conditional distribution of each. This is done one at a time since each distribution is different.
 \item Estimate the six parameters of the distribution using all 95 months.
 \item Pick random values for all the missing measurements from the respective predictive posterior conditional distribution of each.
 \item Repeat steps 3 and 4 a number of times, preferably until the parameter estimates converge in distribution.
 \item Store the last set of completed data and last set of parameter estimates.
  \item Repeat steps 1 to 6 a number of times (min. 5) to obtain a collection of imputed datasets.
\end{enumerate}

The parameter estimates that are produced as a result of this approach are given in \autoref{tbl:imputes}. It is clear that parameters $\alpha_2, \alpha_4$ and $\alpha_5$ would be underestimated if only the complete data were used.

An extract of imputed data values is given in \autoref{tbl:imputespred} for illustration.

It is clear that this approach provides a valid alternative to the commonly used multivariate normal imputation \citep[see][for an explanation of why this is important]{vonHippel2013}. %Another, perhaps even more important use of this method is that it can address the problem of zeros in composition data. Zeros can be made missing and then multiply imputed.

\begin{table}[bh]
    \caption{Parameter estimates following multiple imputation} 
\begin{center}
    \begin{tabular}{|l|r|r|r|r|r|r|}
    \hline
    \textbf{Parameter estimates} & \multicolumn{1}{c|}{$\alpha_1$} & \multicolumn{1}{c|}{$\alpha_2$} & \multicolumn{1}{c|}{$\alpha_3$} & \multicolumn{1}{c|}{$\alpha_4$} & \multicolumn{1}{c|}{$\alpha_5$} & \multicolumn{1}{c|}{$\alpha_6$} \\
    \hline
    Complete data only & 3.57  & 3.49  & 2.72  & 2.8   & 3.62  & 0.42 \\    \hline
    Imputation 1 & 3.43  & 3.83  & 2.78  & 3.13  & 3.97  & 0.43\\    \hline
    Imputation 2 & 3.43  & 3.89  & 2.65  & 2.91  & 4     & 0.43 \\    \hline
    Imputation 3 & 3.5   & 3.82  & 2.77  & 3.03  & 4.03  & 0.43 \\    \hline
    Imputation 4 & 3.46  & 3.83  & 2.75  & 3.12  & 4.07  & 0.43 \\    \hline
    Imputation 5 & 3.55  & 3.94  & 2.8   & 3.17  & 4.05  & 0.43 \\    \hline
    Imputation 6 & 3.46  & 3.81  & 2.78  & 3.15  & 3.98  & 0.43 \\    \hline
    Imputation 7 & 3.63  & 4.04  & 2.84  & 3.33  & 4.17  & 0.44 \\    \hline
    Imputation 8 & 3.2   & 3.49  & 2.54  & 2.83  & 3.7   & 0.42 \\    \hline
    Imputation 9 & 3.28  & 3.65  & 2.65  & 2.93  & 3.66  & 0.42 \\    \hline
    Imputation 10 & 3.59  & 4.02  & 2.85  & 2.98  & 4.12  & 0.44 \\    \hline
    \end{tabular}
\label{tbl:imputes}
\end{center}
\end{table}

\begin{table}[thb]
  \caption{Extract from final round of Imputation 2 showing typical imputed values}
\begin{center}
    \begin{tabular}{|l|r|r|r|r|r|}
    \hline
    \textbf{Month} & \multicolumn{1}{l|}{\textbf{Station 1}} & \multicolumn{1}{l|}{\textbf{Station 2}} & \multicolumn{1}{l|}{\textbf{Station 3}} & \multicolumn{1}{l|}{\textbf{Station 4}} & \multicolumn{1}{l|}{\textbf{Station 5}} \\    \hline
    2003m4 & 5.8   & 6.5   & 8.5   & \textbf{6.4} & 9.4 \\    \hline
    2003m8 & 6.3   & 3.4   & 2.1   & \textbf{2.6} & 4.6 \\    \hline
    2003m9 & 30.4  & 18    & 20.2  & \textbf{28.3} & 20.3 \\    \hline
    2003m11 & 47.6  & 49.7  & 16.3  & \textbf{55.1} & 42.8 \\    \hline
    2003m12 & 24.3  & 52.4  & 0.7   & \textbf{7.9} & 45.2 \\    \hline
    2004m1 & 5.6   & 60.6  & 3.2   & \textbf{2.4} & 49.2 \\    \hline
    2004m3 & 73.8  & 85.3  & \textbf{75.7} & 46.7  & 151.5 \\    \hline
    \end{tabular}
   \label{tbl:imputespred}
\end{center}
\end{table}

\section{Conclusion} \label{sec:conclusion}
It is found that the Dirichlet distribution family is a useful distribution when modelling random variables that take on positive values and especially for composition data with negative correlations, particularly as the number of dimensions increases. It scales easily from the univariate (beta) case to the multivariate (Dirichlet) case with the number of parameters only increasing by one per dimension. The Bayesian approach can add accuracy and flexibility when working with this distribution family. Conditioning within this distribution family and how this can be extended to allow imputation is discussed.

\textit{Note:} This chapter has been edited from the stand-alone version, which may be viewed at this link: 
\url{https://arxiv.org/abs/1801.02962}

\printbibliography[title=List of References for Paper 1]
\end{refsection}




























% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{refsection}
\chapter{Paper 2 --- Multiple Discrimination Between Dirichlet Populations} \label{ch:p2}
\pagestyle{C3}
\begin{LARGE}
Daan de Waal and Sean van der Merwe
\end{LARGE}
\vspace{10mm}

\textbf{Abstract}

This paper shows applications of the discriminant function of the Dirichlet distribution in the multivariate and matrix cases. The first is on discriminating between unknown Dirichlet populations with an application on the Fisher's Iris data with four variables and three populations. This is a case where the Dirichlet parameters have to be estimated for each population. The second application is on constructing rejection regions for matrix Dirichlet variables with known parameters. The correspondence between the likelihood ratio test statistic for testing equality of covariance matrices and the discriminant function between matrix Dirichlet populations is discussed. The application is on volatility of the Southern Oscillation Index (SOI) for specific months among high, medium and low groupings.

\textit{Keywords}: Classification,\ Composition data,\ Dirichlet distribution,\ Discrimination.

\textit{Note:} This chapter has been edited from the stand-alone version, which may be viewed at this link: 
\url{https://journals.co.za/content/sasj/48/2/EJC159423}

\section{Introduction} \label{sec:Intro}
The Dirichlet distribution is a well-known multivariate distribution \citep{Gupta2001} and is defined as the joint density of $X_j=\frac{Y_j}{\sum_{j=1}^{P} Y_j}, j = 1,\ldots,P$ where $Y_j$ is distributed independently $Gamma(\alpha,1)$. It is clear that $0 < X_j <1$ and $\sum_{j=1}^{P} X_j=1$. It can therefore be seen as a $(P-1)$-variate distribution since $X_{P}=1-\sum_{j=1}^{P-1} X_j$. It is referred to as the compositional Dirichlet by Aitchison \parencite*{Aitchison1980,Aitchison1986} defined on the simplex $0 < x_j < 1$ and $\sum_{j=1}^{P} x_j = 1$ in $P$ dimensions. \textcite{Aitchison1986} discussed the advantages of the Dirichlet of which the number of parameters is one more than the number of dimensions. Composition data are usually transformed to the Euclidian space where standard statistical methods can be used \parencite[see][]{Filzmoser2012}. Due to the need for discriminant analysis in the context of composition data where the Dirichlet distribution is applicable, this paper concentrates on the Dirichlet distribution. The classification of a Dirichlet observation $\mathbf{x_0}=(x_{01},x_{02},\ldots,x_{0P})$ to one of $q > 1$ Dirichlet populations through the discriminant function is considered against the alternative that the observation is just a random observation. An application on Fisher's Iris data is given.

In \autoref{sec:matrixD} the assignment of a matrix Dirichlet observation to matrix Dirichlet populations is considered. The correspondence between the likelihood ratio test statistic for equivalence of covariance matrices and the discriminant function between two Dirichlet populations is shown. An application is shown on analysing the volatility of the Southern Oscillation Indexes (SOI) for several months jointly among three groups of indices.

\section{The Dirichlet Distribution} \label{sec:Dirichlet}
\subsection{Introduction} \label{sec:DirichletIntro}
Discussions on many of the properties of the Dirichlet distribution can be found in \textcite{Wilks1962}, \textcite{deGroot1970} and \textcite{Aitchison1986}. A few are mentioned.

If $X=(X_1,X_2,\ldots,X_P)$ is distributed $Dirichlet(\alpha_1,\ldots,\alpha_{P})$ or denoted by $D(\boldsymbol\alpha)$ on the positive simplex in $P$ dimensions, then the joint density is given by
\begin{equation} \label{eq:211jointd}f(x) = \frac{\prod_{j=1}^{P} \Gamma(\alpha_j)}{\Gamma(\alpha_0)}\left\{\prod_{j=1}^P x_j^{\alpha_j-1}\right\},\ 0<x_j<1,\ \sum_{j=1}^P x_j = 1,\ \alpha_0=\sum_{j=1}^{P} \alpha_j. \end{equation}

\begin{enumerate}
\item The means and covariances are
\begin{equation} \begin{aligned}\label{eq:212meancov}
&E(X_j)=\mu_j=\frac{\alpha_j}{\alpha_0},\ j = 1,\ldots,P\\
&Cov(X_j,X_j)=\left\{\begin{array}{l}
\sigma_{ij}=\frac{-\alpha_j \alpha_j}{\alpha_0^2 (1+\alpha_0)},\ j\ne j\\
\sigma_{ii}=\frac{\alpha_j (\alpha_0-\alpha_j)}{\alpha_0^2 (1+\alpha_0)},\ j=j
\end{array}\right.
\end{aligned} \end{equation}
$\mu=(\mu_1,\ldots,\mu_P)$ denotes the mean of the distribution and $\Sigma=(\sigma_{ij})\ j,j=1,\ldots,P$ the covariance matrix. Note the negative correlation between variables.

\item Let $Y_j\sim Gamma(\alpha_j,1),\ j = 1,\ldots,(P)$ independently, then $X_j = Y_j/\sum_{j=1}^{P} Y_j\ (j=1,\ldots,P) \sim D(\boldsymbol\alpha)$. The joint distribution of $X_j=Y_j/Y_{P}\ (j = 1,\ldots,(P-1))$ is referred to as the Dirichlet type II distribution on the positive (P-1)-dimensional space. \label{rem:2ofIntro}

\item The marginal distribution of $X_j$ is $Beta(\alpha_j,\alpha_0-\alpha_j)$.

\item Note that if $\boldsymbol\alpha$ is a multiple of $\boldsymbol\beta$ for two Dirichlet distributions $D(\boldsymbol\alpha)$ and $D(\boldsymbol\beta)$, the means are the same, but the covariance matrices differ.
\item The negative differential entropy $L = E\left[\ln f(X)\right]$ of the Dirichlet distribution (\autoref{eq:211jointd}) is given in \autoref{eq:mdi} as 
\begin{equation}
L=\ln \Gamma(\alpha_0)-\sum_{j=1}^{P}\ln\Gamma(\alpha_j)-\sum_{j=1}^{P}(\alpha_j-1)[\psi(\alpha_0)-\psi(\alpha_j)].
\label{eq:213entropy}
 \end{equation}
$\psi(.) \textnormal{ is the digamma function, }\frac{d}{da}\ln\Gamma(a)=\psi(a).$ 

\item \begin{sloppypar}Consider two Dirichlet populations, $D(\alpha_1,\ldots,\alpha_{P})$ (parent population) and $D(\beta_1,\ldots,\beta_{P})$, then the Kullback-Leibler measure of divergence between the two follows from \autoref{eq:213entropy} as\begin{equation} \label{eq:214kl}K\!L = E\left(\ln\frac{f(X)}{f_0(X)}\right)=\ln\frac{\Gamma(\alpha_0)}{\Gamma(\beta_0)}-\sum_{j=1}^{P}\ln\frac{\Gamma(\alpha_j)}{\Gamma(\beta_j)}-\sum_{j=1}^{P}(\alpha_j-\beta_j)[\psi(\alpha_0 )-\psi(\alpha_j)] \end{equation}  
The extent to which $K\!L\ge 0$ gives a measure of the deviance between $f_0$ and $f$.\end{sloppypar}

\item The exponent of $L$ given in \autoref{eq:213entropy} is the MDI prior of the Dirichlet distribution (\autoref{eq:211jointd}). See \textcite{zellner1997} for information on the MDI prior. \label{rem:7MDIP}

\item Complete composition independence exists between subsets of X \parencite[318--319]{Aitchison1986}.
\end{enumerate}

\subsection{The Uniform Distribution on P dimensional simplex} \label{sec:unidist}

If $\alpha_j=1,\ j=1,\ldots,P$, \autoref{eq:211jointd} becomes a uniform distribution. The joint density is given by
\begin{equation} \label{eq:221uniform}f(\mathbf{x})=\frac{1}{(P-1)!},\ 0<x_j<1,\ \sum_{j=1}^{P-1} x_j<1. \end{equation}
The distribution is defined on the P-dimensional simplex above.

The mean becomes $\mu=E(X)=\left[\begin{array}{c}
\frac{1}{P}\\ \vdots\\ \frac{1}{P}\end{array}\right]$
 and the covariance matrix
 
 $\Sigma = Cov(X,X') = \frac{1}{(P)^2(P+1)}\left[\begin{array}{ccc}
 (P-1)&\cdots&-1\\ \vdots&\ddots&\vdots\\ -1&\cdots&(P-1)\end{array}\right]$ with $\Sigma^{-1}=(P)(P+1)\left[\begin{array}{ccc}
 2&\cdots&1\\ \vdots&\ddots&\vdots\\ 1&\cdots&2\end{array}\right]$.
  
Let $P=3$, then given an arbitrarily chosen $x_2$, an $x_1$ can be simulated using $x_1=(1-x_2)u_1$ for a random $0<u_1<1$ using the property of the conditional distribution.
 
The observations are not attracted to the mean at all and the correlation $\rho_{ij}=-\frac{1}{P-1}$ between any two variables $X_i$ and $X_j$ becomes smaller as the dimension increases.

Visualisation of a dataset is shown in \autoref{fig:Uniform} up to 5 dimensions. The figure attempts to show that, while data is normally drawn in two dimensions, we must introduce depth, then width, then height to expand to more dimensions. Alternatively, one must use colour, size, or other tricks to add additional information on the same graph.

\begin{figure}[tbh]
\begin{center}
\includegraphics[width=\figwidth]{UniformDimensions.png}
\end{center}
\caption{A visualisation of the sampled data from the uniform distribution in 2 dimensions to 5 dimensions}

\label{fig:Uniform}
\end{figure}

The point to note is that the observations will be attracted to the mean, irrespective of the number of dimensions, unless the Dirichlet becomes uniform (the parameters all drop to one). While the uniform distribution is a special case of the Dirichlet distribution, it is also a contrast to the Dirichlet distribution in that it lacks the patterns associated with the Dirichlet distribution.

The uniform distribution is used going forward to represent the idea of randomness in observations. It is uniquely suited for this task as it lacks the patterns that a researcher might wish to model and is in the same model class as the models that will be applied.

\subsection{Discriminant function}\label{sec:discrimfunc}
Fisher's discriminant function in the case of multivariate normal populations is well known \parencite{fisher1936,anderson,Muirhead}. Bayesian discriminant rules resulting in linear or quadratic group separations based on mainly normal populations also exist \parencite{Filzmoser2012}. This paper considers only discrimination among Dirichlet populations.

In the case of discriminating between two populations with density functions $f_1$ and $f_2$ respectively, an observation $\mathbf{x_0}$ will be assigned to population 1 if $g(\mathbf{x_0})=\ln\frac{f_1(\mathbf{x_0})}{f_2(\mathbf{x_0})} > c$ for some constant $c$. In \autoref{sec:twopops} the assignment of such an observation is evaluated if it belongs to one of two Dirichlet populations $D_1$ or $D_2$. In \autoref{sec:multidiscrim} the analysis to more than two populations.
For $\mathbf{X}_0$ random, $E[g(\mathbf{X}_0)]$ becomes the Kullback-Leibler measure $K\!L$ defined in \autoref{eq:214kl} as the divergence measure between $f_1$ and $f_2$.

\subsubsection{Two populations}\label{sec:twopops}
Suppose the random variable $X$ is distributed either $D_1(\boldsymbol\alpha)$ or $D_2(\boldsymbol\beta)$. If $\mathbf{x_0}$ is observed, then according to the above ratio of the two densities, $\mathbf{x_0}$ will be assigned to $D_1$ if
\begin{equation} g(\mathbf{x_0})=(\boldsymbol\alpha-\boldsymbol\beta)\ln(\mathbf{x_0})>c. \end{equation}
This is a linear function in the elements of $\ln(\mathbf{x_0})$.

The constant $c$ is chosen such that the probabilities of misclassification
\begin{equation} \label{eq232defc} 
P[g(\mathbf{X})<c|D_1(\boldsymbol\alpha)]\textnormal{ and }P[g(\mathbf{X})>c |D_2(\boldsymbol\beta)] \end{equation}
are equal.

The distribution of $g(\mathbf{X})=(\boldsymbol\alpha-\boldsymbol\beta)\ln(\mathbf{X})$ is unknown, but can be simulated. The mean is
\begin{equation}
 \mu_g=E(g(\mathbf{X})|D_1)=(\boldsymbol\alpha - \boldsymbol\beta)E[\ln(\mathbf{X})|D_1]=\sum_{j=1}^{P}(\alpha_j-\beta_j)[\psi(\alpha_j)-\psi(\alpha_0)].
\end{equation}       

Also $E\ln(X_j)=E\ln(Y_j)-E\ln(\sum_{j=1}^{P} Y_j)$ according to the definition of $\mathbf{X}$ where $Y_j$ is distributed $Gamma(\alpha_j,1)$ (see \autoref{sec:Intro}). It is known that $E\ln(Y_j)=\psi(\alpha_j)$. The variance $Var(\ln(Y_j))=\psi'(\alpha_j)$. $\psi'$ is called the trigamma function.

The variance of $g(X)$ is
\begin{equation}
Var(g(X)|D_1)=(\boldsymbol\alpha-\boldsymbol\beta)\Phi(\boldsymbol\alpha-\boldsymbol\beta)',
\end{equation}
where $\Phi=Cov(\ln X,\ln X'|D_1)$.

Since $E[\ln(X_i)\ln(X_j)]=E\left[\ln(Y_i)-\ln\left(\sum_{i=1}^{P}Y_i\right)\right]\left[\ln(Y_j)-\ln\left(\sum_{i=1}^{P}Y_i\right)\right]$ and \\$Var(\ln(Y_i))=\psi'(\alpha_i)$ \parencite{Chan1993}, it follows that the $i,j^{\ th}$ element of $\Phi$ is
 \begin{equation} \begin{aligned}
\phi_{ij}&=Cov(\ln X_i,\ln X_j|D_1) \\
&\approx\left\{\ \begin{array}{l}
\psi'(\alpha_i)+{\psi(\alpha_i)}^2-2\psi(\alpha_i)\ln(\alpha_0)-[\psi(\alpha_i)-\psi(\alpha_0)]^2,\ i = j\\
\psi(\alpha_i)\psi(\alpha_j)-\ln(\alpha_0)\ln(\alpha_i)-\ln(\alpha_0)\ln(\alpha_j)\\
\ \ -[\psi(\alpha_i)-\psi(\alpha_0)][\psi(\alpha_j)-\psi(\alpha_0)],\ i\ne j,
\end{array}\right.
\end{aligned} \end{equation}
taking $E\left[\ln(Y_j)\ln\left(\sum_{j=1}^{P}Y_j\right)\right]\approx\ln(\alpha_j)\ln(\alpha_0)$.

\textbf{Remark:} Note that the joint distribution of $V_j=\frac{X_j}{X_{P}},\ j=1,\ldots,(P-1)$ is a Dirichlet type II and it follows easily that $Cov(\ln V_i,\ln V_j)=\psi'(\alpha_{P})+\psi(\alpha_{P})^2-\psi(\alpha_{P})[\psi(\alpha_i)+\psi(\alpha_j)]$. 

\paragraph{Example 2.3.1}
Consider as an example $n=30$ observations simulated from two independent Dirichlet populations $D_1(\mathbf{k}=(5,15,10))$ and $D_2(\mathbf{l}=(7,6,0.6))$. The observations are shown in \autoref{fig:232simdata2pops} with dots and their respective means with $\oplus$. According to \autoref{eq232defc}, $c=-20$. Suppose the researcher wishes to assign an observation $\mathbf{x_0}=(0.4,0.4,0.2)$ to one of the two populations. The distribution function of $g(X)$ given $D_1$ and survival function of $g(X)$ given $D_2$ defined in \autoref{eq232defc} as simulated are shown in \autoref{fig:231probsOfMis}. The observation is assigned to $D(\mathbf{l})$ since the observation falls on that side of the discriminant function in \autoref{fig:232simdata2pops}. 
 
\begin{figure}
\caption{Probabilities of misclassification given $D(\mathbf{k})$ and $D(\mathbf{l})$ respectively}
\begin{center}

\includegraphics[width=\figwidth]{probsOfMis.png}
\end{center}

Labels $\mathbf{x_1}$ and $\mathbf{x_2}$ refer to values of $g(X)$ and probabilities respectively.
\label{fig:231probsOfMis}
\end{figure}
 
\begin{figure}
\begin{center}
\includegraphics[width=\figwidth]{simdata2pops.png}
\end{center}
\caption{Simulated data from two Dirichlet populations $D_1$ and $D_2$}

The with the discriminant function $g(\mathbf{x_1},\mathbf{x_2})=-20$ and the observation $\mathbf{x_0}$.
\label{fig:232simdata2pops}
\end{figure}

Suppose the observation $\mathbf{x_0}=(0.4,0.1)$ is to be classified, then rather classify it as a random observation than to assign it to $D(\mathbf{l})$. This is the approach followed in the next section when one wants to assign an observation to one of several populations against the alternative that it is a random observation.

\subsubsection{Multiple discrimination}\label{sec:multidiscrim}
Consider classifying an observation $\mathbf{x_0}\in D(\boldsymbol\alpha^{(i)}),\ i=1\textnormal{ or }2,\ldots,\textnormal{ or }q$ if
\begin{equation} (\boldsymbol\alpha^{(i)}-\boldsymbol\alpha)\ln\mathbf{x_0})>c_i,\ i=1,\ldots,q, \end{equation}
else consider $\mathbf{x_0}$ as random, namely $\mathbf{x_0}\in D(\boldsymbol\alpha=(1,\ldots,1))$.

\emph{Criterion:} Select quantiles $c_i$ as critical values such that $P($Type I error$)\approx 0.05$, say,
\begin{equation} P(X_0\in D(\boldsymbol\alpha)|D(\boldsymbol\alpha^{(i)},\ i=1\textnormal{ or }\ldots\ q))\approx 0.05.\label{eq:237criterion} \end{equation}

The expectation of the discriminant function for a random $X\in D_i$, against randomness, can be written as
\begin{equation} \begin{aligned}
E[g(X)|D_i]&=(\boldsymbol\alpha^{(i)}-\boldsymbol\alpha)E\ln X\\
&=\sum_{j=1}^{P}\left(\alpha_j^{(i)}-1\right)\left[\psi\left(\alpha_j^{(i)}\right)-\psi\left(\alpha_0^{(i)}\right)\right].
\end{aligned} \end{equation}

\paragraph{Example 2.3.2}
Consider $q=4$ Dirichlet populations with $P=3$ and $\boldsymbol\alpha^{(1)}=(10,20,3);\ \boldsymbol\alpha^{(2)}=(1,0.2,13);\ \boldsymbol\alpha^{(3)} =(16,2,3);\ \boldsymbol\alpha^{(4)}=(160,20,30)$. Select the quantiles $c_i,\ i=1,\ldots,4$ such that the probability of the type I error
\begin{equation} P\left(X_0\in D(\boldsymbol\alpha)|D\left(\boldsymbol\alpha^{(i)},\ i=1\textnormal{ or }\ldots\ q\right)\right) \end{equation}
is small. A Type I error is to assign the observation $x_0$ as random given it belongs to one of the $q$ populations. Selecting $c=(c_i,\ i=1,\ldots,4)=(-30.4437, 1.8616, -160.2341, -14.8681)$, the probability is estimated through simulation, $P(I) = 0.0522$. The discriminating contours
\begin{equation} (\boldsymbol\alpha^{(i)}-\boldsymbol\alpha)\ln(x)=c_i \end{equation}
are obtained through simulating the x-values such that \autoref{eq:237criterion} is true (or almost true) for all $i=1,\ldots,4$. The simulated data are shown in \autoref{fig:233simdata4pops}, $n=1000$ for each population, indicated by blue (for included observations) and red (for excluded observations) dots. The data points on the right, indicated by blue, green and red dots, are two overlapping populations, $i=3,4$, with the same means but different covariance matrices. An observation $(0.6,0.1)$ will be assigned to population $i=3$ while another $(0.7,0.1)$ to population $i=4$. An observation $(0.2,0.1)$ will be considered as random since it cannot be assigned to any of the four populations.
 
\begin{figure}
\begin{center}
\includegraphics[width=\figwidth]{simdata4pops.png}
\end{center}
\caption{Simulated data from $q=4$ Dirichlet populations with the discriminant contours}

\label{fig:233simdata4pops}
\end{figure}

\paragraph{Example 2.3.3}
Fisher's Iris flower dataset is a well-known dataset with 4 variables and 3 populations. It contains observations on the four variables $Y_1=$ sepal length, $Y_2=$ sepal width, $Y_3=$ petal length and $Y_4=$ petal width of Iris-setosa, Iris-versicolor and Iris-virginica respectively measured in $cm$. Assume that $(Y_1,\ldots,Y_4)$ are jointly Dirichlet type II distributed (see Remark \ref{rem:2ofIntro} in \autoref{sec:DirichletIntro}). Then the distribution of $X=\left(\frac{Y_j}{1+\sum_{j=1}^4 Y_j},\ j=1,\ldots,4\right)$ can be described through a Dirichlet distribution $D(\alpha_1,\ldots,\alpha_5)$ for each population, say $X^{(i)}\sim D_i\left[\boldsymbol\alpha^{(i)}=\left(\alpha_j^{(i)},\ j=1,\ldots,5\right)\right],\ i=1,2,3$.
 
\paragraph{Parameter estimation}
There are moment, maximum likelihood and Bayes estimates of the parameters $\boldsymbol\alpha^{(i)}$ to consider, but the Bayes estimates are preferred since through many simulations they perform the best. The marginal posterior modes from the joint posterior of $\boldsymbol\alpha^{(i)}$ were used as the estimates. The maximal data information (MDI) prior of Zellner (see remark \ref{rem:7MDIP} of \autoref{sec:DirichletIntro}) is used as joint prior of the parameters.
 
The original data (available in the free package R, \parencite*{RCore}, or from Wikipedia at \url{http://en.wikipedia.org/wiki/Iris_flower_data_set}) is shown in \autoref{fig:234irisscatter}.
 
\begin{figure}
\begin{center}
\begin{minipage}[t]{\figwidth}
\vspace{0pt}
\includegraphics[width=\figwidth]{irisscatter.pdf}
\end{minipage}
\begin{minipage}[t]{0.25\figwidth}
\vspace{9mm}
\includegraphics[width=0.25\figwidth]{irisscatterLegend.pdf}
\end{minipage}
\end{center}
\caption{Fisher's Iris flower data}

\label{fig:234irisscatter}
\end{figure}

Dirichlet distributions are fitted on the four variables for each species, assuming each species follows a Dirichlet type II distribution and transforming the data to Dirichlet type I accordingly. As a measure of goodness-of-fit, all observations are assigned to its corresponding population under the multiple discriminant rules.
From the data, the estimates of the Dirichlet parameters are
\begin{equation} \begin{array}{l}
\hat{\alpha}^{(SET)} = (604.71, 412.57, 175.49, 27.78, 120.02)\\
\hat{\alpha}^{(VER)} = (378.31, 176.39, 270.87, 84.17, 64.15)\\
\hat{\alpha}^{(VIR)} = (427.61, 193.26, 360.36, 131.62, 65.37).
\end{array} \end{equation}

A basic test of goodness-of-fit is to simulate replicate samples and compare them to the original data visually \parencite[see for example][159--171]{gelman2013}. Doing this with the Iris data produces graphs such as \autoref{fig:235irisballs}. There appears to be sufficient correspondence between observed and simulated data.
 
\begin{figure}
\begin{center}
\includegraphics[width=0.8\figwidth]{setosaD2.png}\includegraphics[width=0.8\figwidth]{versicolorD2.png}

\includegraphics[width=0.8\figwidth]{virginicaD2.png}
\end{center}
\caption{Iris data with simulated replicate samples on Dirichlet type II scale}

Green indicates the original data and blue a simulated replicate sample.
\label{fig:235irisballs}
\end{figure}

\autoref{fig:235irisballs} shows a display of the 50 simulated observations with the 50 observed on each population in 4 dimensions. The fourth variable is indicated by the size of the balls. 

\textbf{Remark:} Instead of assuming a Dirichlet type II distribution on $Y_j,\ j=1,\ldots,4$, one can think of transforming the data to say $Z_j,\ j=1, \ldots,4$ such that it is reasonable to assume that $Z_j\sim Gamma(\alpha_j,1),\ j=1,\ldots,4$ independently. $X_j=\frac{Z_j}{\sum_{j=1}^4 Z_j},\ j=1,\ldots,3$ will then follow a Dirichlet distribution on a 3 dimensional simplex. The advantage is that the number of dimensions is reduced to three.

\section{The Matrix Dirichlet Distribution} \label{sec:matrixD}

The matrix Dirichlet distribution, denoted here as $MD()$, came about when \textcite{wilks1932} derived the likelihood ratio test statistic for testing the equality of the covariance matrices of multivariate normal populations \parencite[see also][291--293]{Muirhead}. An application of the discriminant function of this distribution is considered in this section.

\subsection{Introduction}\label{sec:matrixDintro}
Let $A_j\sim Wishart(n_j,I_p)$ independent of $B\sim Wishart(m,I_p)$.

The joint density of $L_j=\left(\sum_{j=1}^q{A_j+B}\right)^{-\frac{1}{2}}A_j\left(\sum_{j=1}^q{A_j+B}\right)^{-\frac{1}{2}},\ j=1,\ldots,q$ is given by
\begin{equation} \label{eq:311wishartjoint} 
\begin{aligned}
f(L_1,\ldots,L_q)=&\frac{\Gamma_p(\frac{1}{2}(n+m))}{\prod_{j=1}^q{\Gamma_p(\frac{1}{2} n_j)}\Gamma_p(\frac{1}{2}m)}\prod_{j=1}^q{|L_j|^{\frac{1}{2}(n_j-p-1)}|I-\sum_{j=1}^q L_j|^{\frac{1}{2}(m-p-1)}},\\
&\ 0<\sum_{j=1}^q L_j < I_p,\ \Gamma_p(a)=\pi^{p(p-1)/4}\prod_{i=1}^p{\Gamma(a+\frac{1}{2}(1-i))} ,
\end{aligned}
\end{equation}
and is called the matrix Dirichlet distribution.

The expected value of $L_j$ is given by $E(L_j)=\frac{n_j}{(n+m)}I_p$ where $n=\sum_{j=1}^q n_j$.

\subsubsection{Matrix Uniform}\label{sec:matrixunif}
If $n_j=p$ for all $j$ and $m=p$, then \autoref{eq:311wishartjoint} becomes the matrix uniform distribution in $q+1$ dimensional hyperspace with density
\begin{equation} f(L_1,\ldots,L_q)=\frac{\Gamma_p(\frac{1}{2}(q+1)(p))}{\left[\Gamma_p(\frac{1}{2}(p))\right]^{q+1}},\ 0<\sum_{j=1}^q L_j < I_p. \end{equation}
The expected value now becomes $E(L_j)=\frac{1}{q+1}I_p,\ j=1,\ldots,q$.

\subsubsection{The Matrix Dirichlet type II distribution}\label{sec:matrixDtype2}

\begin{sloppypar}Alternatively, define $V_j=B^{-\frac{1}{2}}A_jB^{-\frac{1}{2}}$ where $A_j\sim Wishart(n_j,I)$ independently of $B\sim Wishart(m,I)$, then the joint density of $V_1,\ldots,V_q$ becomes the matrix Dirichlet type II density with density function
\begin{equation} f(V_1,\ldots,V_q)=\frac{\Gamma_p(\frac{1}{2}(n+m))}{\prod_{j=1}^q \Gamma_p(\frac{1}{2} n_j)\Gamma_p(\frac{1}{2} m)} \prod_{j=1}^q |V_j|^{\frac{1}{2}(n_j-p-1)}|I+\sum_{j=1}^q V_j|^{-\frac{1}{2}(n+m)}, \end{equation}
$n=\sum_{j=1}^q n_j,\ V_j>0$. Denote it by $MDirich2(n_j,j=1,\ldots,q;m)$.\end{sloppypar}

The expected value of $V_j$ is $E(V_j)=\frac{n_j}{m-p-1}I$ which exists for $m>p$. If $p=1$, then the Dirichlet Type II case follows and $\frac{m}{2} > 1$ has been mentioned. This can be derived from the definition $V_j=B^{-\frac{1}{2}}A_jB^{-\frac{1}{2}}$ where $A_j\sim Wishart(n_j,I)$ independent of $B\sim Wishart(m,I)$ and using the property $E[E(V_j|B)]=E(V_j)$. Since $B^{-1}\sim InvWishart(m,I)$, the expected value is $E(B^{-1})=\frac{1}{m-p-1}I$ and $E(A_j)=n_j I$.
 
\subsection{The discriminant function}\label{sec:matrixdiscrimf}

\begin{sloppypar}
The log discrimination function for assigning a matrix Dirichlet observation $L_0=(L_{0,1},\ldots,L_{0,q+1})$ to one of two populations $MD(n^{(1)},m^{(1)})$ or $MD(n^{(2)},m^{(2)})$ is fairly easy to derive by taking the log of the ratio between the two densities at $L_0$ with parameters $\mathbf{a}=(n^{(1)},m^{(1)})$ and $\mathbf{b}=(n^{(2)},m^{(2)})$ as\end{sloppypar}
\begin{equation} \label{eq:321matDdiscrim} D = (\mathbf{a} - \mathbf{b})[\ln |L_{0,1}|,\ldots,\ln|L_{0,q+1}|]'. \end{equation}

The critical value $c$ can be chosen according to \autoref{eq232defc}.

If, however, one would like to assign an observation to a given Dirichlet population or classify it as random, the distribution function of $D$ in \autoref{eq:321matDdiscrim} can be simulated with $b=(p,\ldots,p)$ and $c$ can be calculated from the distribution function. The following example illustrates this.

\paragraph{Example 3.2.1}
Let $\mathbf{n}=(15,13)$ and $m=5$ be the parameters of a matrix Dirichlet where $p=3$ and $q=2$. It is of interest whether the observation $(|L_{0,1}|,|L_{0,2}|,|L_{0,3}|)=(0.1353, 0.0498, 0.0025)$ can be considered as an observation from the $MD$ or must be assigned as a random observation. The empirical DF of $D$ with $\mathbf{a}=(15,13,5)$ and $\mathbf{b}=(p,p,p)$ from simulated data is shown in \autoref{fig:321ECDF}. The critical 5\% value is $c=-71.7046$. The observed value is $d=-55$ which means the observation is seen as belonging to MDir$(a)$.
 
\begin{figure}
\begin{center}
\includegraphics[width=\figwidth]{ECDF.png}
\end{center}
\caption{Empirical DF of $D$}

The 5\% critical value is $c=-71.7046$.
\label{fig:321ECDF}
\end{figure}

\paragraph{Example 3.2.2}
This is an example there are $q=2$ matrix Dirichlet populations $MD_1(15,3,5)$ and $MD_2(7,7,7)$ with $p=3$ variables. An observation needs to be assigned to one of the two populations. The probabilities of misclassification are equal if $c=23$. \autoref{fig:322ECDFmd} shows the DF of the discriminant function $D$ defined in \autoref{eq:321matDdiscrim} if $MD_1$ is true and the survival function of $D$ if $MD_2$ is true from simulated data. An observation with $D=40$ will be assigned to $MD_1$ and with $D=0$ will be assigned to $MD_2$.
 
\begin{figure}
\begin{center}
\includegraphics[width=\figwidth]{ECDFmd.png}
\end{center}
\caption{Distribution function of $D$ given $MD_1$ and survival function of $D$ given $MD_2$}

The green dashed line indicates the survival function.
\label{fig:322ECDFmd}
\end{figure}

\subsubsection{Multiple matrix discrimination}
In the case of multiple matrix Dirichlet discrimination, \autoref{eq:321matDdiscrim} can be applied similarly to the procedure described in \autoref{sec:multidiscrim} by taking $\mathbf{b}=(p,p,\ldots,p)$ and classifying an observation as belonging to one of the populations or randomness.

\paragraph{Example 3.2.3}
\begin{sloppypar}
Consider $q=3$ matrix Dirichlet populations with $p=3$, namely $MD_1(15,8,5)$, $MD_2(7,7,7)$ and $MD_3(7,9,15)$. 1000 Matrices were simulated from each population and indicated with dots (blue, red and green respectively) in \autoref{fig:323logLscatter3D}. From the empirical DF's of $D$ defined in \autoref{eq:321matDdiscrim} with $\mathbf{a}$ given for each dataset from $MD_i(\mathbf{a})$ above, 53.5\% critical values have been calculated, namely $\mathbf{c}=(-46.2085, -35.5028, -62.8077)$ respectively. In each case $\mathbf{b}=(p, p, p)$. The observations that will be considered as random are encircled in \autoref{fig:323logLscatter3D}. The total number of misclassified observations (out of 3000) are 155, giving an error probability of misclassification of $\frac{155}{3000}=0.0517$. One can adjust the 53.5\% to get a desirable error probability.
\end{sloppypar}

\begin{figure}[hbt]
\begin{center}
\includegraphics[width=\figwidth]{logLscatter3D.png}
\end{center}
\caption{Scatterplot of simulated $\ln|L_1|,\ln|L_2|,\ln|L_3|$ values from each Dirichlet population}

Observations encircled are considered random.
\label{fig:323logLscatter3D}
\end{figure}

The examples considered in this section are for known parameters of the Dirichlet populations. If the parameters are unknown and the Dirichlet matrices observed, the parameters have to be estimated. This problem has not been considered.

The likelihood ratio test statistic for testing the equivalence between covariance matrices can be written in terms of matrix Dirichlet variables. It is of interest to know what the difference between this statistic and the matrix Dirichlet discriminant function is. This is considered in the next section.

\subsection{Testing equality of covariance matrices}\label{sec:equalcovmats}
Let $L_j=\left(\sum_{j=1}^q A_j\right)^{-\frac{1}{2}}A_j\left(\sum_{j=1}^q A_j\right)^{-\frac{1}{2}}$ where $A_j\sim W(\Sigma_j,n_j),\ j=1,\ldots,q$ independent of each other. Then the joint distribution of $L_j,\ j=1,\ldots,q$ is, according to \autoref{sec:matrixDintro}, $Dirichlet(\mathbf{n},\ j=1, \ldots,q)$ if $\Sigma_j=I_p\ \forall\ j$. The log likelihood ratio test statistic for testing the equality of $\Sigma_j=\Sigma,\ j=1,\ldots,q$, ignoring the constant, is given by \parencite{wilks1932,anderson,Muirhead}
\begin{equation} \Lambda=\sum_{j=1}^q{n_j\ln|L_j|}.\label{eq:331lambda} \end{equation}

This statistic is very similar to the discriminant function (\autoref{eq:321matDdiscrim}) if the covariance matrices are equal, and the discriminant function between the $Dirichlet(n_j,\ j=1,\ldots,q)$ and randomness, a $Dirichlet(p,\ldots,p)$, namely
\begin{equation} g(L_1,\ldots,L_q)=\sum_{j=1}^q{(n_j-p-1)\ln|L_j|}.\label{eq:332discrim} \end{equation}

The distribution of $\Lambda$ can be simulated fairly easily and a critical value $c$ be calculated for small values of $\Lambda$ below which the assumption of equal covariance matrices will be rejected. The two statistics \autoref{eq:331lambda} and \autoref{eq:332discrim} are very similar.

\paragraph{Example 3.3.1}
The southern oscillation index (SOI) is an index based on the air pressure difference between Darwin and Tahiti. This index is used as a leading indicator of rain/drought especially in the southern hemisphere. The monthly SOI values from 1876 are available on the Web at \url{http://www.bom.gov.au/climate/current/soihtm1.shtml}. The values are scaled to have mean zero and can be considered to be normally distributed. It has been found that high October SOI values (above 10) indicate good summer rains in the central part of South Africa while a low October SOI (below -10) gives an indication of draught. The summer rain season runs from October to April. It is of interest whether the covariance matrices of the SOI of the months November to April differ if their SOI values are grouped into $A$, $B$ and $C$ with $A\equiv October SOI > 10,\ B\equiv -10 < Oct SOI < 10$ and $C\equiv Oct SOI < -10$. This is considered a case where $p=6,\ q=3$ and
\begin{equation} 
M_A\sim Wishart(\Sigma_A,n_A),\ M_B\sim Wishart(\Sigma_B,n_B)\text{ and } M_C \sim Wishart(\Sigma_C,n_C) 
\end{equation}  
independent of each other. The degrees of freedom are $n_A=19$, $n_B=88$, $n_C=28$. It is of interest whether $\Sigma_A=\Sigma_B=\Sigma_C$. If this is true,  $(L_A,L_B,L_C)\sim MD(n_A,n_B,n_C)$ where $L_A=(M_A+M_B+M_C)^{-\frac{1}{2}}M_A(M_A+M_B+M_C)^{-\frac{1}{2}}$. Similarly for $L_B$ and $L_C$. The empirical distribution function of the simulated values of $\Lambda$ given in \autoref{eq:331lambda} is shown in \autoref{fig:ECDFlambda} with the 5\% critical value $c=-778.6722$. The observed value is $-779.2858$ which means one rejects the hypothesis of equal covariance matrices.

\begin{figure}
\begin{center}
\includegraphics[width=\figwidth]{ECDFlambda.png}
\end{center}
\caption{Empirical distribution function of $\Lambda$ with 5\% critical value}

\label{fig:ECDFlambda}
\end{figure}

\section{Conclusion} \label{sec:Conclusion}
Discriminating between Dirichlet populations is considered. The extension to matrix Dirichlet populations is shown and the relationship with the test statistic for equality of covariance matrices with the discriminant function illustrated.

\textit{Note: The statistical packages that were used to produce the results illustrated in the figures were MATLAB\textregistered\ \parencite{MATLAB} and R \parencite{RCore}. \autoref{fig:235irisballs} was produced using the RGL package \parencite{RGL}.}

\printbibliography[title=List of References for Paper 2]
\end{refsection}
















% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{refsection}
\chapter{Paper 3 --- Classification of Multiple Dirichlet Observations under a Multinomial Model} \label{ch:p3}
\pagestyle{C4}
\begin{LARGE}
Daan de Waal, Roelof Coetzer, and Sean van der Merwe
\end{LARGE}
\vspace{10mm}

\textbf{Abstract}

The amount of gas produced from a coal gasification facility depends crucially on the properties and the size distribution of the coal being used in the process. The particle size distribution and the composition of the coal are measured as composition data. This paper applies the Dirichlet distribution for the inputs and presents a new classification scheme for yielding low, medium and high gas production. The approach presented is a linear partitioning of the Dirichlet simplexes that can also be extended to high dimensional cases. An alternative clustering approach based on a distance measure is also presented.

{\it Keywords:}  Coal gasification, Composition data, Dirichlet distribution, Multinomial Model

\textit{Note:} This chapter has been edited from the stand-alone version, which may be viewed at this link: 
\url{https://www.sciencedirect.com/science/article/pii/S0169743915002634}

\section{Introduction} \label{sectionIntro}

Sasol's coal-to-liquids facility delivers nearly 29\% of the fuel requirements in South Africa, and the continuous improvement and optimisation of the coal gasification plant is of critical importance to the company to ensure stable supply of high quality synthetic gas to the downstream units. The amount of gas produced from the coal gasification facility depends crucially on the properties and the size distribution of the coal being used in the process. Therefore, in order to optimise gas production, the relationship between the properties of the coal and the gas produced must be understood and quantified \citep{Coetzer2003}. This paper develops and discusses a new classification scheme that discriminates between coal properties and size distributions for yielding high gas production. 

The particle size distribution (PSD) and the composition of the coal are measured as composition data. The Dirichlet distribution has been widely accepted in literature for modelling composition data, subject to the constraint that all the correlations between variables are negative \citep{Minka}. The coal PSD and composition conforms to the negative correlation constraint, and therefore, the Dirichlet distribution is fitted to these input variables.  A wider class of distributions which allows for positive correlations defined on the same sample space is the logistic-normal (LN) \citep{Aitchison1980}. However, the LN distribution has many parameters to estimate due to the unknown covariance parameter matrix. In contrast, the Dirichlet distribution has only $p$ unknown parameters to estimate for $p$ composition variables.  

The current study is motivated by a coal gasification industrial problem where the amount of gas produced is a function of the quality of the coal being fed to the reactors, and it is crucially important to the business to know which coal sizes and compositions will yield optimal or sub-optimal production. The coal-to-liquids production facility gasifies no less than 40 million tons of bituminous coal per year which is delivered to the factory from 5 to 7 coal sources. Therefore, the coal sizes and compositions can be very different from source to source and from week to week. Furthermore, the coal blending schedule is updated weekly subject to feed availability. Since the coal gasification process is a continuous process and the performance of the reactors is monitored in real time, it is also desired to employ a real-time monitoring model for the coal feed to the reactors using data obtained from appropriate analytical equipment. Therefore, the classification scheme presented in this study was developed for this purpose. This classification scheme can, for example, be implemented for real-time monitoring of the coal properties composition and the size distribution of the coal being fed to the reactors for early detection of coal blends which may yield low gas production \citep{Coetzer2014}. Real-time detection of sub-optimal feed allows the factory to be pro-active in changing the coal blends to minimize the losses in gas production and to sustain profits.

In this study, three size fractions (expressed as proportions) $x_1,\ x_2$ and $x_3$ and three compositions (also expressed as proportions) $c_1,\ c_2$ and $c_3$ were considered. The dependent variable $(Y)$ depending on covariates $C = (C_1,C_2,C_3)$ and $X = (X_1,X_2,X_3)$ is divided into three categories namely high (H), standard (S) and low (L), representing 20\%, 68\% and 12\% of production respectively. These boundaries were chosen to illustrate the methodology only, and have no relevance to current operational performance.  It is assumed that $Y$ conditional on $C=\mathbf{c}$ or $X=\mathbf{x}$ is distributed $Multinomial(n,\ p_\mathbf{x})$ where $p_\mathbf{x}=(p_H,p_S,p_L)$ represents the multinomial probabilities. The aim of this paper is to classify a coal sample as H, S or L under the desired constraint $p_0=(0.20,0.68,0.12)$, given its composition $(\mathbf{c})$ and sizes $(\mathbf{x})$. Here $P_X$ denotes the random variable and $p_\mathbf{x}$ the observed proportion given $\mathbf{x}$ (or $\mathbf{c}$).  
 
The classification scheme based on a single covariate $(X)$ can be summarised as follows:
\begin{enumerate}
\item $Y|P_X \sim f(\mathbf{y}|P_X) \equiv Multinomial(n,P_X)$
\item $P_X|\mathbf{x} \sim k(p_\mathbf{x} |g(X))$ where $k$ is some selected distribution of $P_X$ and $g(X)$ refers to the information on $X$ used to specify $k$.
\item $X|\boldsymbol\alpha \sim h(\mathbf{x}|\boldsymbol\alpha) \equiv Dirichlet(\boldsymbol\alpha)$
\item The classification scheme: Partition the Dirichlet sample space into subspaces H, S and L such that $p_x = (p_H,p_S,p_L)$ coincides with $p_0$. 
\end{enumerate}

This paper demonstrates and discusses the application of the above classification scheme using simulated data from the specified probability distributions. The results are generated through simulation since the integrals cannot be solved analytically. This paper illustrates and discusses the performance of two alternative classification schemes:
\begin{enumerate}
\item a linear partitioning of the Dirichlet sample spaces, and
\item a clustering approach based on minimising a chosen metric.
\end{enumerate}  

The paper is outlined as follows: In \autoref{sectionDirichlet}, a brief introduction on applicable Dirichlet properties is given. In \autoref{secClas} the linear method is described and applied on a simulated dataset. In \autoref{secHighDim} the proposed method is extended to higher dimensions. In \autoref{secClus} the clustering approach is discussed. Conclusions and future work are discussed in \autoref{sectionConclusion}.

\section{The Dirichlet distribution} \label{sectionDirichlet}

\citet{Wilks1962} and \citet{deGroot1970} provided detailed discussions on many of the properties of the Dirichlet distribution. A few are mentioned. If $X$ of order $P$ by 1 is distributed $Dirichlet(\alpha_1,\ldots,\alpha_{P})$, denoted by $D(\boldsymbol\alpha)$, then the joint density is given by 

\begin{equation} \label{D1}
\begin{aligned}
f(\mathbf{x})= &\frac{\prod_{j=1}^{P}\Gamma(\alpha_j)}{\Gamma(\alpha_0)}\left\{ \prod_{j=1}^P x_j^{\alpha_j-1}\right\} ,\\
&\ 0 < x_j < 1,\ \sum_{j=1}^{P-1} x_j < 1,\ \alpha_0 = \sum_{j=1}^{P} \alpha_j .
\end{aligned}
\end{equation}

\citet{Aitchison1986} calls this distribution the compositional Dirichlet defined on the specified simplex.
\begin{enumerate}
\item The means and covariances are
\begin{equation} \label{cov}
\begin{aligned}
E(X_j)&= \mu_j=\frac{\alpha_j}{\alpha_0}, j=1,\ldots,p \\
\sigma_{ij}&=\left\{\begin{array}{l}
\dfrac{-\alpha_i\alpha_j}{\alpha_0^2(1+ \alpha_0)}, i\neq j\vspace{1ex} \\
\dfrac{\alpha_i(\alpha_0 - \alpha_i)}{\alpha_0^2(1+ \alpha_0)}, i=j
\end{array}\right.
\end{aligned}
\end{equation}
 
$\mu = (\mu_1,\ldots,\mu_P)$ denotes the mean of the distribution and $\Sigma=(\sigma_{ij})\ i,j=1,\ldots,P$ the covariance matrix.

\item The marginal distribution of $X_i$ is $Beta(\alpha_i,\alpha_0-\alpha_i)$. The marginal distribution of any subset $J$ of $X$, with the remaining components aggregated, is again $Dirichlet(\alpha_J, \alpha_0 - \sum_{i\in J} \alpha_i)$.
\item Note in \autoref{cov} that if $\boldsymbol\alpha$ is a multiple of $\boldsymbol\beta$ for two Dirichlet distributions $D(\boldsymbol\alpha)$ and $D(\boldsymbol\beta)$, the means are the same, but the covariance matrices differ.
\item The negative differential entropy $L = E(\ln{f(X)})$ of the Dirichlet distribution (\autoref{D1}) is given by 
 
\begin{equation} \label{entropy}
\begin{aligned}
 L=&\ln{\Gamma(\alpha_0)} - \sum_{j=1}^{P}\ln \Gamma(\alpha_j )\\
&\ - \sum_{j=1}^{P}(\alpha_j-1) \left\{\psi(\alpha_0 )-\psi(\alpha_j)\right\}
\end{aligned}
\end{equation}

(\autoref{eq:mdi}). $\psi(.)$ refers to the digamma function.

\item The log Jeffreys prior \citep{jeffreys1998} for the $Dirichlet(\boldsymbol\alpha)$ is

\begin{equation}
0.5\sum_{l=1}^{P}\ln \psi'(k_l) + 0.5\ln \left[1 - \psi'(k_0)\sum_{l=1}^{P}\frac{1}{\psi'(k_l)}\right].
\end{equation}  

$\psi'(.)$ refers to the trigamma function.

This prior is used later to obtain a posterior distribution of parameters. The derivation of the above is given in \autoref{sec:post}.
\item \begin{sloppypar}Suppose there are two Dirichlet populations, $D(\alpha_1,\ldots,\alpha_{P})$ (parent population) and $D(\beta_1,\ldots,\beta_{P})$, then the Kullback-Leibler measure of divergence between the two follows from \autoref{entropy} as\end{sloppypar}

\begin{equation}
\begin{aligned}
K\!L=&E\ln\frac{f(X)}{f_0(X)} = \ln\frac{\Gamma(\alpha_0)}{\Gamma(\beta_0)}- \sum_{j=1}^{P}\ln\frac{\Gamma(\alpha_j)}{\Gamma(\beta_j)} \\
&\ - \sum_{j=1}^{P}(\alpha_j-\beta_j)\left\{\psi(\alpha_0)- \psi(\alpha_j)\right\} 
\end{aligned}
\end{equation}

The size of $K\!L \geq 0$ gives a measure of divergence between $f_0$ and $f$.
\end{enumerate}

\section{Classification of Dirichlet sample spaces} \label{secClas}

\begin{sloppypar}A dataset is constructed on compositions $(C)$ and sizes $(X)$ by simulating 200 $Dirichlet(\boldsymbol\alpha)$ and $Dirichlet(\boldsymbol\beta)$ observations independently. To each random vector of observations $(C_1,C_2,X_1,X_2)$, the production level H, S or L is assigned according to the linear functions partitioning the two sample spaces (see \autoref{compsizeslines}) such that 20\% of the observations are classified as H, 68\% as S and 12\% as L. Note these boundaries were chosen to illustrate the methodology only, and have no relevance to current operational performance.\end{sloppypar}

\begin{figure}
\begin{center}
\includegraphics[width=\figwidth]{FIG1_JPEG.jpg}
\end{center}
\caption{Simulated Dirichlet observations on C and X with linear lines classifying each observation}

Classifications are H (green), S (blue) or L (red). \label{compsizeslines}
\end{figure}

The classification rules are specified as follows:

Classify a sample as 
\begin{itemize}
\item[] H if $C_2-1+c_1 C_1<0$ and $X_2-1+\frac{1}{x_2}  X_1>0$,
\item[] L if $C_2-1+c_2 C_1>0$ and $X_2-1+\frac{1}{x_1}  X_1<0$,
\item[] S otherwise.
\end{itemize}
 
Note that $c_i$ and $x_i , i = 1,2$ denote chosen end points of the lines with origins $(0,1)$ and $(1,0)$ respectively such that the production criteria $(20,68,12)$ are met. \autoref{dataextract1} shows part of such a dataset.

\begin{table}
\caption{Extract of the dataset}
\begin{center}

\begin{tabular}{|c|c|c|c|c|}
\hline
  &$C_1$&$C_2$& $X_1$&$X_2$\\ \hline
L	&0.222	&0.475	&0.075	&0.546\\ \hline
S	&0.086	&0.740	&0.492	&0.508\\ \hline
S	&0.100	&0.395	&0.415	&0.466\\ \hline
S	&0.038	&0.553	&0.196	&0.509\\ \hline
H	&0.203	&0.061	&0.413	&0.579\\ \hline
L	&0.501	&0.086	&0.199	&0.726\\ \hline
S	&0.014	&0.326	&0.313	&0.625\\ \hline
S	&0.160	&0.582	&0.416	&0.401\\ \hline
H	&0.034	&0.066	&0.307	&0.679\\ \hline
H	&0.063	&0.233	&0.141	&0.855\\ \hline
S	&0.690	&0.028	&0.396	&0.569\\ \hline
\end{tabular}
\label{dataextract1}
\end{center}

\end{table}

In practice, only such data is available and the aim is to set up the classification scheme to classify future observations of $C$ and $X$ as H, S or L .

\subsection{Partitions with linear functions} \label{secPartLinear}

The partitioning linear functions shown in \autoref{compsizeslines} partition the sample spaces into subspaces H, S and L with probabilities  $p_H,\ p_S$ and $p_L$ such that they equal $p_0=(0.2,0.68,0.12)$ approximately. Various options exist to draw these lines, but due to the nature of the three regions in both $C$ and $X$ Dirichlet spaces, it was decided to let the three lines originate at the point $(0,1)$ on compositions and $(1,0)$ on sizes. The regions may be such that the lines start from a centre point $(a,b)$ which will require careful thinking. The above choices make the analysis fairly simple, since the only variables needed to fully specify the lines are their slopes.

Let the slopes be $b_{c1}$ and $b_{c2}$ on the $C$ space and $b_{x1}$ and $b_{x2}$ on the $X$ space, and the corresponding abscissa points be $0 < c_1 < c_2 < 1$ on $C$ and $0 < x_1 < x_2 < 1$ on $X$. The only values to be selected are $c_1,\ c_2,\ x_1$ and $x_2$. These points are selected independently for $C$ and $X$.

The distribution function of the random variable $P_X$ can be written in terms of the marginal $C_1\sim Beta(a,b)$ of $C$. Explicitly
\begin{equation}\begin{aligned}
P(P_X<p_x|\theta)&= P(P_{X_1}<p_{x_1},P_{X_2}<p_{x_2}) \\
&=\int\!\!\int_Sg(x_1,x_2)\,dx_1dx_2\ ,
\end{aligned}
\end{equation}  

where S is defined by the two linear functions such that $p_0$ holds approximately for determining $c_1$ and $c_2$. Similarly one can write the distribution of $P_X$ on the $X$ space.

The classification rules are specified as follows:

\begin{itemize}
\item Composition: Classify a sample as H if $C_2-1+c_1 C_1<0,$ L if $C_2-1+c_2 C_1>0$ and S otherwise

\item Sizes: Classify a sample as H if $X_2-1+\frac{1}{x_2} X_1>0$, L if $X_2-1+\frac{1}{x_1} X_1<0$ and S otherwise.

\item Combined rules: These rules can be combined by numerically selecting $c_i$ and $x_i$ simultaneously or separately as presented in \autoref{clasrules}. For example, classify a sample as H if Composition is H and Size is S, as L if Composition is S and Size is L, \textit{etc.}
\end{itemize}

\begin{table}
\caption{Classification rules to match the goal of $(20\%,68\%,12\%)$ of production}
\begin{center}

\begin{tabular}{|c|c|c|c|}
\hline
Sizes\textbackslash\textbf{Compositions}	&\textbf{H}	&\textbf{S}	&\textbf{L} \\ \hline
H	&H	&H	&H\\ \hline
S	&S	&S	&L\\ \hline
L	&L	&L	&L\\ \hline
\end{tabular}

\label{clasrules}
\end{center}

\end{table}

The combined rules in \autoref{clasrules} are obtained from the simulations where the outcomes match the target probabilities $p_0$ closely.

\subsection{Validation of the linear partitioning functions} \label{secValLinear}

To validate the rules in \autoref{clasrules}, a sample of $m=200$ observations was drawn on $C$ and $X$ and classified according to the rules in \autoref{clasrules}. The resulting misclassification matrix, expressed in percentages, is given in \autoref{valtable1}. The success rate of correct classification was 82\%. In addition the success rate of correct classification for the low gas production was 92\%, which is significant for the commercial plant since it highlights the ability to detect sub-optimal coal feeds.

\begin{table}
\caption{Misclassification matrix according to rules in \autoref{clasrules}}
\begin{center}

\begin{tabular}{|c|c|c|c|}
\hline
Classified\textbackslash\textbf{True}	&\textbf{H}	&\textbf{S}	&\textbf{L} \\ \hline
H	&4	&1	&0\\ \hline
S	&14	&67	&1\\ \hline
L	&0	&2	&11\\ \hline
\end{tabular}

\label{valtable1}
\end{center}

\end{table} 

The results in \autoref{valtable1} were tested against the assumption that the generating process is uniform (this is equivalent to say that the observations come from a $Multinomial(n=200,p_o)$ distribution) and the outcome in \autoref{valtable2mn} was obtained. The success rate is now only 55\%. This is the success rate one would expect if the explanatory variables add no information on the production output.

\begin{table}[hbt]
\caption{Misclassification matrix according to uniform Dirichlet inputs}
\begin{center}

\begin{tabular}{|c|c|c|c|}
\hline
Classified\textbackslash\textbf{True}	&\textbf{H}	&\textbf{S}	&\textbf{L} \\ \hline
H	&4	&22	&5\\ \hline
S	&16	&50	&1\\ \hline
L	&0	&1	&1\\ \hline
\end{tabular}

\label{valtable2mn}
\end{center}

\end{table}

This process was repeated a number of times and similar results were obtained. This shows that, through linear partitioning of the sample spaces with chosen anchor points, classification is improved by almost 30\% using the Dirichlet assumptions.

However, ignoring the distribution of sizes and only concentrating on the distribution of compositions in the above simulated data, an average success rate of about 90\% is obtained. This is an improvement of 35\% above a random classification from the multinomial distribution. %due to the large extent of overlapping of the high (green), standard (blue) and low (red) subpopulations in the sizes sample. Unless we can find a better anchor point for the linear partitioning of the simplex, it is better to omit it.

In the commercial application it is most important to identify those coal types which will yield S to L gas production. In this case, the results are promising and illustrate that the classification scheme can be implemented as a real-time process monitoring tool for alarming the factory of the sub-optimal coal feed.

\section{Higher dimensions} \label{secHighDim}

The method described can be applied for classifying high dimensional Dirichlet observations in the multinomial setup. The number of multinomial cells can also be larger than three. Let $X \sim Dirichlet(\boldsymbol\alpha)$ be defined on the $k$ dimensional simplex describing the independent variables and $Y|\mathbf{x}\sim Multinomial(n, p_\mathbf{x})$. In a regression model $E(Y|\mathbf{x})=np_\mathbf{x}$. To estimate $p_\mathbf{x}$, partition the $k$ dimensional simplex with linear functions through the anchor point $(0,0,0,\ldots,0,1)$. The anchor point can be any other point in the simplex depending on the data. The criteria for classifying $\mathbf{x}$ as L, S or H is given by the following linear functions:

Classify $\mathbf{x}$ as L if $x_k - c_{k,L} + b_{1,L} x_1 + \dots + b_{k-1,L} x_{k-1} < 0,$
as H if  $x_k - c_{k,H} + b_{1,H} x_1 + \dots + b_{k-1,H} x_{k-1} > 0$,
otherwise classify $\mathbf{x}$ as S. 

$b_{i,j}$ is the slope of $x_i$ with respect to $x_k$ of the linear function $j$ and $c_{k,j}$ the intercept of the function $b_{(j)} \mathbf{x}$ with the axis $x_k$. The intercepts of $b_{(j)} \mathbf{x}$ with axis $x_{i\neq k}$ is denoted by $c_{(j)}$ where $b_{(j)} = (b_{1,j} ,\ldots, b_{k-1,j}) = -(\frac{c_{k,j}}{c_{1,j}} ,\ldots,\frac{c_{k,j}}{c_{k-1,j}})$. Estimates of $c_{i,j}$ can be obtained through simulating a large number of observations from the marginal distribution of $(X_1,\ldots,X_{k-1}) \sim Dirichlet(\alpha_1,\ldots,\alpha_{k-1},\alpha_t)$ where $\alpha_t = \alpha_{k+1} + \alpha_k$ minimizes the sum of squares $SSE = \sum_{i=1}^{3}(\hat{p}_i- p_{0i})^2$. $p_0$ is the-pre assigned target for the yield. The estimates can be validated similarly to the three dimensional case if a dataset of $2n$ observations has been simulated. The first $n$ observations on $\mathbf{y}$ and $\mathbf{x}$ are treated as the data and then the next $n$ observations are classified and counted.  

\subsection{Illustration on 4-dimensional case} \label{sec4Dillus}

Let $Y \sim Multinomial(n,P_X)$ with density function
\begin{equation} \label{eqMultiDensity}
\begin{aligned}
f(\mathbf{y}|P_X)=&\frac{n!}{y_1!y_2!y_3!} P_{X_1}^{y_1} P_{X_2}^{y_2} P_{X_3}^{y_3},\\
&\ y_i > 0,\ y_1+y_2+y_3 =n.
\end{aligned}
\end{equation}

\begin{sloppypar}The parameter vector $P_X$ is considered random conditional on the composition $X\sim Dirichlet(\boldsymbol\alpha)$ defined on the simplex $0<X_j<1,\  j=1,2,3$ and $X_1+X_2+X_3<1$ \citep{Aitchison1986}. The density function of $X$ is given by\end{sloppypar}

\begin{equation}
g(\mathbf{x})= \frac{\Gamma(\alpha_0)}{\prod_{j=1}^4 \Gamma(\alpha_j)} \prod_{j=1}^4x_j^{\alpha_j-1},\ 0<x_j<1,\ \sum_{j=1}^4x_j=1,\ \alpha_0 = \sum_{j=1}^4 \alpha_j.
\end{equation}

Partitioning the simplex $0<x_j<1,\ j=1,2,3$ and $x_1+x_2+x_3<1$ with linear functions $a\mathbf{x}=c$ and $b\mathbf{x}=d$ such that $(0,0,1)$ is considered an anchor point in the simplex for both functions, the distribution function of $P_X$ is given by

\begin{equation}
\begin{aligned}
P(P_X<p_\mathbf{x}|\theta)&= P(P_{X_1}<p_{x_1},P_{X_2}<p_{x_2},P_{X_3}<p_{x_3})\\
&= \int\!\!\int_Sg(x_1,x_2)\,dx_1dx_2
\end{aligned}
\end{equation}
 
where S is defined by the two linear functions such that $p_0$ holds. This implies that the distribution of $P_X$ becomes the joint distribution of the four end points of the two linear functions. A discrete approach is used here and estimate the distribution of $P_X$ from the joint density of $X_1$ and $X_2$, namely Dirichlet$(\alpha_1,\alpha_2,\alpha_0\!-\!\alpha_1\!-\!\alpha_2)$, through simulation given $p_\mathbf{x}$. It is a property of the Dirichlet that any subset is again Dirichlet \citep{Aitchison1986,Frigyik2010}, as explained in \autoref{sec:intro}.
 
\subsection{Application} \label{secAppli}

Let $X\sim $Dirichlet$(\boldsymbol\alpha)$ where $\boldsymbol\alpha = (9,11,57,80)$. A dataset of 1000 observations was simulated from the Dirichlet and a yield vector $y$ constructed with two linear functions through an anchor point $(0,0,1)$ such that 13\% of observations are classified as H, 60\% as S and 27\% as L. The dataset is shown in \autoref{3dzoomed}.

\begin{figure}[hbt]
\begin{center}
\includegraphics[width=\figwidth]{FIG2_JPEG.jpg}
\end{center}
\caption{Dirichlet data of different yields}

Red dots show the L yield data, green the H yield data and blue the S yield data.
\label{3dzoomed}

\end{figure}
 
The constants are chosen for constructing the data and will be unknown in practice. Since one only has to consider the marginal distribution of $X_1$ and $X_2$ as shown in \autoref{eqMultiDensity} to estimate the distribution of $P_X$, the abscissa on the $X_3$ axis is irrelevant.
The classification scheme is:

Classify $(x_1,x_2,x_3)$ as L if $x_2-c_2+\frac{c_2x_1}{c_1} < 0$, \\H if $x_2-d_2+\frac{d_2x_1}{d_1} > 0$ and S otherwise.

The slopes of the two linear functions from lower to upper are $\frac{-c_2}{c_1} $ and $\frac{-d_2}{d_1}$  respectively, $c_1$ and $d_1$ being the abscissae. There is a need to estimate the parameters $c_1, c_2, d_1$ and $d_2$ from the data. The first $n=500\ y$ values form the yield data while the first $500\ \mathbf{x}$ values form the independent Dirichlet observations. \autoref{dataextract2} shows a part of this dataset.
 
\begin{table}[hbt]
\caption{Extract of higher dimensional dataset}
\begin{center}

\begin{tabular}{|c|c|c|c|c|} \hline
$y$	&$x_1$	&$x_2$	&$x_3$	&$x_4$\\ \hline
L	&0.046	&0.054	&0.398	&0.502\\ \hline
H	&0.045	&0.119	&0.397	&0.438\\ \hline
S	&0.070	&0.079	&0.325	&0.527\\ \hline
S	&0.057	&0.102	&0.371	&0.471\\ \hline
S	&0.071	&0.071	&0.381	&0.478\\ \hline
S	&0.058	&0.041	&0.418	&0.483\\ \hline
S	&0.042	&0.091	&0.361	&0.507\\ \hline
S	&0.032	&0.072	&0.420	&0.477\\ \hline
S	&0.054	&0.061	&0.406	&0.479\\ \hline
L	&0.065	&0.024	&0.387	&0.525\\ \hline
\end{tabular}

\label{dataextract2}

\end{center}
\end{table}

A scatterplot of $(x_1,x_2)$ is shown in \autoref{scatterplot}. To be able to classify a new observation $x_0$ as H, S or L, the linear functions need to be estimated.
 
\begin{figure}[hbt]
\begin{center}
\includegraphics[width=\figwidth]{FIG3_JPEG.jpg}
\end{center}
\caption{Scatterplot of $n=500$ observations on $x_1$ and $x_2$ with two linear functions}

The linear functions pass through $(0,c_2)$ and $(c_1,0)$ for the lower function, and $(0,d_2)$ and $(d_1,0)$ for the upper function.
\label{scatterplot}
\end{figure}

To estimate the parameters $c=(c_1,c_2)$ and $d=(d_1,d_2)$ the joint density of $X_1$ and $X_2$ needs to be integrated, which has to be done numerically. This paper uses simulation for this. The joint marginal of $X_1$ and $X_2$ is $Dirichlet(\alpha_1,\alpha_2,\alpha_0\!-\!\alpha_1\!-\!\alpha_2)$ with estimates of the parameters  $(8.5269,10.5715,132.1404)$ respectively. Simulating a large number of observations from this distribution and selecting $c_i, d_i,i=1,2$ such that our goal of a yield of $(27,60,13)$\% is achieved, this method obtains $\hat{c} = (0.1132,0.1040)$ and $\hat{d} = (0.1307,0.1889)$. These estimates provide the information for drawing the linear functions. 

\subsection{Validation} \label{secVali}

To validate the classification criteria, classify each observation from the 500 simulated observations left out from the original sample by substituting the $(x_1,x_2)$ values in the two equations. 76\% correct classifications was obtained. The summary of the results is depicted in \autoref{clasoutcomes}.

\begin{table}
\caption{Classification outcomes}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
Classified\textbackslash\textbf{True}	&\textbf{H}	&\textbf{S}	&\textbf{L} \\ \hline
H	&15	&54	 &0\\ \hline
S	&11	&268 &29\\ \hline
L	&0	&29	 &98\\ \hline
\end{tabular}
\end{center}
\label{clasoutcomes}
\end{table}

From \autoref{clasoutcomes} it is noted that 15 out of 26 observations were classified as H, 268 out of 351 were classified as S and 98 out of 127 were classified as L. The rate of correct classification for the low gas production was 77\%, which is slightly lower than for the initial example of three dimensions. 

\section{Clustering approach} \label{secClus}

The aim of this approach is to classify an observation as H, S or L based on which of these three groups it is closest to, given some chosen metric. For the purpose of illustration this paper uses the Mahalanobis distance. The Mahalanobis distance is based on the mean vector and covariance matrix of each group. In order to visualise this concept, consider every point in the simplex for the composition data and one by one calculate the distance from the point to the group centres. Classify each point based on which of the three distances is the smallest. Using simulated data as in \autoref{secClas}, the results are illustrated in \autoref{malfig1}.

\begin{figure}[hbt]
\begin{center}
\includegraphics[width=\figwidth]{FIG4_PDF.pdf}
\end{center}
\caption{Mahalanobis classification of each point to a group}

\label{malfig1}
\end{figure}

There is a choice of calculating the distances empirically, or by fitting Dirichlet distributions to each group and using the theoretical mean and covariance matrices (\autoref{cov}), which was done here. 

The advantage of this approach is that every point in the domain is classified and the calculations needed in order to classify any new observation are straightforward. For example, one can validate the approach similarly to \autoref{secPartLinear} by overlaying the points on the graph, as seen in \autoref{malfig2}. Surprisingly few observations are misclassified, as can be seen in \autoref{malclas1}. Specifically, note the excellent classification of the low production cases (92\%) which is very important for the commercial implementation of the classification scheme.

\begin{figure}[hbt]
\begin{center}
\includegraphics[width=\figwidth]{FIG5_PDF.pdf}
\end{center}
\caption{Composition sample overlaid on Mahalanobis classification}

\label{malfig2}
\end{figure}

\begin{table}[hbt]
\caption{Misclassification matrix for Mahalanobis classification of compositions}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
Classified\textbackslash\textbf{True}	&\textbf{H}	&\textbf{S}	&\textbf{L} \\ \hline
H	&29	&10	 &0\\ \hline
S	&11	&112 &2\\ \hline
L	&0	&14	 &22\\ \hline
\end{tabular}
\end{center}
\label{malclas1}
\end{table}

Repeating the above procedure for the size data yields the results in \autoref{malfig3} and \autoref{malclas2}. Here the rate of correct classification on the low gas production is only 50\%. It is clear that classification of the sizes is a more difficult enterprise, as was indicated earlier with the linear classification. Since the discrimination power of the two types of data differ so greatly under this analysis this paper did not investigate combining the analyses using decision rules. It is worth noting that the option to do so exists and may have value in more balanced scenarios.

\begin{figure}[hbt]
\begin{center}
\includegraphics[width=\figwidth]{FIG6_PDF.pdf}
\end{center}
\caption{Sizes sample overlaid on Mahalanobis classification}
\label{malfig3}
\end{figure}

\begin{table}[hbt]
\caption{Misclassification matrix for Mahalanobis classification of sizes}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
Classified\textbackslash\textbf{True}	&\textbf{H}	&\textbf{S}	&\textbf{L} \\ \hline
H	&23	&42	 &1\\ \hline
S	&13	&59  &11\\ \hline
L	&4	&35	 &12\\ \hline
\end{tabular}
\end{center}

\label{malclas2}
\end{table}

In the case of the Mahalanobis (or similar metric) classification scheme, it is straightforward to extend the approach to higher dimensions because in this approach the dimension is irrelevant. In \autoref{malfig3d} classification of the 4-component data shown in \autoref{3dzoomed} and \autoref{dataextract2} is illustrated to aid the reader in visualising this approach as the dimension increases.

\begin{figure}[hbt]
\begin{center}
\includegraphics[width=\figwidth]{FIG7_JPEG.jpg}
\end{center}
\caption{Mahalanobis classification of higher dimensional Dirichlet data}

The data are illustrated in \autoref{3dzoomed} and \autoref{dataextract2}.
\label{malfig3d}
\end{figure}

\section{Violations of assumptions} \label{secRobustness}

Most of the assumptions are intuitive but the assumption of the Dirichlet distribution is worth scrutinising. This is done by simulating data which vaguely resembles our data, but is not Dirichlet. Note that the clustering method in \autoref{secClus} yields similar results as it did for the Dirichlet data, suggesting that the researcher need not be overly concerned with satisfying the Dirichlet assumption in this case.

As an alternative to the Dirichlet distribution one can simulate new data from the logistic-normal (LN) \citep{Aitchison1980}. To understand to what extent the new dataset differs, generate QQ plots, as given in \autoref{fig:QQplots}. Plot first the best Dirichlet fit on Dirichlet data and then the best Dirichlet fit on the LN data. This illustration uses only the Low data (our primary focus) and only the first dimension, which yields plots that are typical of this type of data. One may conclude from \autoref{fig:QQplots} that the Dirichlet distribution fits both sets of data adequately.

\begin{figure}[hbt]
\begin{center}
\includegraphics[width=0.7\figwidth]{FIG9_PDF.pdf}
\includegraphics[width=0.7\figwidth]{FIG10_PDF.pdf}
\end{center}
\caption{Best Dirichlet fit on Dirichlet data and LN data}

Best Dirichlet fit (line) on Dirichlet data (circles) in the first figure, and best Dirichlet fit (line) on the LN data (circles) in the second figure.
\label{fig:QQplots}
\end{figure}

Considering the Mahalanobis classification in \autoref{fig:MalLN} a similar pattern to \autoref{malfig2} is seen. The resulting classification matrix (\autoref{malclasLN}) is also very similar, yielding results that are almost, but not quite as good as for the Dirichlet data, as is to be expected.

\begin{figure}[hbt]
\begin{center}
\includegraphics[width=\figwidth]{FIG11_PDF.pdf}
\end{center}
\caption{Logistic-normal sample overlaid on Mahalanobis classification}

\label{fig:MalLN}
\end{figure}

\begin{table}[hbt]
\caption{Misclassification matrix for Mahalanobis classification of logistic-normal sample}
\begin{center}

\begin{tabular}{|c|c|c|c|}
\hline
Classified\textbackslash\textbf{True}	&\textbf{H}	&\textbf{S}	&\textbf{L} \\ \hline
H	&30	&14	 &0\\ \hline
S	&10	&92  &3\\ \hline
L	&0	&30	 &21\\ \hline
\end{tabular}
\end{center}

\label{malclasLN}
\end{table}

The correct classification for the Low data is now 88\% instead of 92\%. This difference is surprisingly small considering the difference in the data visible in the QQ plots (\autoref{fig:QQplots}) and there is no reason to expect anything different for the other scenarios. The authors are thus confident that the clustering approach is robust to violations of its assumptions.

\section{Comparison to normal approach} \label{secNormal}

It is worth comparing the performance of these methods to that of some common approaches. Often researchers will ignore the restricted nature of the observations and apply well known statistical techniques as is, and it is argued that there are advantages to using the Dirichlet assumption instead.

A key point to bear in mind is that the Dirichlet is fitted with fewer parameters, so for comparable within-sample performance one would favour the Dirichlet out-of-sample, as it should be less prone to over-fitting and thus generalise better. This is particularly important in industrial applications such as ours as the out-of-sample sample size is often many orders of magnitude larger than what is used to fit the model.

Even if out-of-sample performance is comparable, one would still favour the Dirichlet model as it takes the nature of the data into account making it more likely to correctly classify fringe cases (near the corners of the domain). It also cannot produce predictions and intervals outside the domain of the data.

For the situation described in the previous sections a large number of simulated experiments were run where the following steps are followed for each experiment: simulate a small training sample, fit the Dirichlet Mahalanobis approach, the normal Mahalanobis approach, and linear discriminant analysis (LDA), and then apply the three fitted models to a large testing sample. The rates of correct classification were averaged and it is noted that the methods produced rates in the region of 75\%, 76\% and  81\% respectively, with some variation.

However, when focusing on the more extreme observations (near the corners of the domain) the Dirichlet approach performs the best, with classification rates now averaging 76\%, 74\% and 73\% respectively.


\section{Conclusions} \label{sectionConclusion}

This paper developed and illustrated a classification scheme using the Dirichlet distribution for composition inputs, and applied it to discriminate between the coal properties and size distributions for yielding low, medium and high gas production for a coal-to-liquids facility. The paper illustrates and discusses the performance of two alternative classification schemes:  
\begin{enumerate}
\item a linear partitioning of the Dirichlet sample spaces, and
\item a clustering approach based on minimising a chosen metric.
\end{enumerate}  

The classification scheme based on a single covariate $(X)$ can be summarised as follows:
\begin{enumerate}
\item $Y|P_X  \sim f(\mathbf{y}|P_X) \equiv Multinomial(n,P_X)$
\item $P_X|\mathbf{x} \sim k(p_\mathbf{x} |g(X))$ where $k$ is some selected distribution of $P_X$ and $g(X)$ refers to the information on $X$ used to specify $k$.
\item $X|\boldsymbol\alpha \sim h(\mathbf{x}|\boldsymbol\alpha) \equiv Dirichlet(\boldsymbol\alpha)$
\item The classification scheme: Partition the Dirichlet sample space into subspaces H, S and L such that $p_\mathbf{x} = (p_H, p_S, p_L)$ coincides with $p_0$. 
\end{enumerate}

The linear partitioning of the Dirichlet simplex has to be treated with care in high dimensional cases. The choice of the anchor point plays an important role in the classification scheme and if not one of the corner points, it becomes complicated for higher than 3 dimensional cases. The aim of the clustering  approach is to classify an observation as H, S or L based on which of these three groups it is closest to, given some chosen metric. For the purpose of illustration the Mahalanobis distance was used. The Mahalanobis distance is based on the mean vector and covariance matrix of each group. The advantage of this approach is that every point in the domain is classified and the calculations needed in order to classify any new observation are straightforward. The validation results were very promising for both methods, and the methods described here are ideal when it is clear that the data arises from distinct groups or clusters.

Furthermore, for the current coal gasification problem under study, the results showed that the classification scheme developed can be implemented as a real-time monitoring model for the coal feed to the reactors. The scheme can be implemented by using either the linear partitioning of the Dirichlet sample spaces, or the clustering approach based on minimising a chosen metric. Therefore, deviations in the composition and the size distribution of the coal being fed to the reactors can be detected early on by classifying the coal for yielding high, standard or low gas production. This predictive process monitoring approach may yield significant prevention of losses in gas production for the business.

\textit{Note: The statistical packages that were used to produce the results illustrated in the figures and tables were MATLAB\textregistered\ \citep{MATLAB} and R \citep{RCore}. \autoref{malfig3d} was produced using the RGL package \citep{RGL}.}

\printbibliography[title=List of References for Paper 3]
\end{refsection}

















% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{refsection}
\chapter{Paper 4 --- A Multivariate Gamma Distribution Applied To Composition Data Analysis} \label{ch:p4}
\pagestyle{C5}
\begin{LARGE}
Daan de Waal, Roelof Coetzer, and Sean van der Merwe
\end{LARGE}
\vspace{10mm}

\textbf{Abstract}

Parametric composition data analysis in a high dimensional simplex can be performed by employing the Dirichlet distribution, or alternatively, through the logistic normal distribution if the Dirichlet is not appropriate. In this paper, a multivariate gamma (MGAM) distribution is proposed as an alternative distribution for composition data. In addition, the MGAM distribution is extended to a multivariate extreme value (MEV) distribution and goodness-of-fit statistics are calculated for comparison against the logistic normal distribution. An application is considered where the amount of gas produced from a coal gasication facility depends crucially on the size distribution of the coal, which is measured as composition data and characterized by six variables. The observed sample space is divided into three regions of high (H), standard (S) and low (L) gas production by choosing appropriate thresholds, and new observations are classified among the regions.

\textit{Keywords}: Composition data,\ Multivariate Gamma distribution,\ Multivariate extreme value distribution,\ Classification.

\section{Introduction}

The amount of gas produced from a coal gasication facility depends crucially on the properties and the size distribution of the coal being used in the process. Therefore, in order to optimise gas production, the relationship between the type of coal and the gas produced must be understood and quantified. In the current application, data on six coal sizes were observed and it is expressed as compositions in a five dimensional simplex. The data are not Dirichlet since positive correlations have been found amongst the variables. Therefore, the logistic normal (LN) distribution as an alternative parametric model to apply \parencite{Aitchison1986}. This paper proposes a multivariate gamma distribution (MGAM) as an alternative to the LN. There exists a number of multivariate gamma distributions in the literature, see for example \textcite{prekopa1978,patil1984dictionary}.  However, this distribution presented here is new, and has the advantage that it is fairly simple to fit to data. It also has the advantage that if the data are heavy tailed, the gamma is a more appropriate distribution than the logistic Normal. It has fewer number of parameters than the LN and estimates can easily be obtained. A goodness-of-fit test is considered using the largest trace of a difference square matrix between the observed matrix and the simulated data from the specified distribution. In the case of the LN, the posterior predictive distribution was used \parencite{geisser1983prediction} instead of the plug-in estimates as in the case of the MGAM. The posterior predictive distribution of the MGAM is not considered since the distribution is not available explicitly.

The sample of size $n = 125$ is divided into three groups, namely $n_H = 23$ yielding high (H) gas production, $n_S = 84$ yielding standard (S) production and $n_L = 18$ yielding low (L) production.  The objective is to classify a new observation as H, S or L. The sample space of the MGAM is partitioned into the three regions by partitioning the sample space of the trace of the distance measure between an observation and its simulated equivalent from the model with thresholds derived from the observed outcomes such that 18\% of the observations belong to H, 68\% to S and 14\% to L. 

Zeros do occur, but  were considered to be outliers. Therefore, observations containing zeros were deleted.

\section{Multivariate Gamma distribution} \label{sec:MGAM}
There exist many definitions of multivariate gamma distributions. The version proposed here is a generalisation of the one defined by \textcite{deWaal2004}. The definition is as follows:

\begin{sloppypar}\textbf{Definition:} The random variable $V$ of size $p$ by 1 is distributed multivariate gamma (MGAM) with shape parameter $\alpha$ ($p$ by 1) and correlation structure $(\rho_{ij}), i,j=1, \ldots , p$ if the elements of $Z = \exp(W)$ are distributed independently $Gamma(\alpha_i,1), i= 1 , \ldots , p$ where $W  = H^{1/2} D^{-1} \{ logV - \psi(\alpha)\} + \psi(\alpha)$. $D$ is defined as a matrix square root of $\Lambda = \left( \rho_{ij} \sqrt{\psi'(\alpha_i) \psi'(\alpha_j)}\right), i,j = 1, \ldots , p$. Explicitly, $\Lambda = D'D$ and $H = diag(\psi'(\alpha))$. $\psi(.)$ refers to the standard  digamma function applied independently to each element of its argument. $\psi'(.)$ refers to its first derivative, the trigamma function, and $diag(.)$ refers to the diagonal matrix.\end{sloppypar}

\textbf{Remarks:}
\begin{enumerate} %[(i)]
\item $E(W) = \psi(\alpha)$ and Cov$(W,W') = H$.
\item This paper refers to the distribution of $\ln V$ as a multivariate extreme value (MEV) distribution with $E(logV) = \psi(\alpha)$ and Cov$(\ln V) = \Lambda$. Cov(.) refers to the covariance matrix. The special case $p = 1$ and $\alpha = 1$ reduces the MGAM of $V$ to the standard exponential distribution with the distribution of $\ln V$ known as the Extreme Value Distribution \parencite{Kotz2000}.
\item It is known that the $\ln \Gamma(\alpha , \beta)$ distribution belongs to the Pareto class \parencite{beirlant2004statistics} with extreme value index $1/\alpha$. Therefore, it may be specified that $W = \exp(V)$ belongs to a class of multivariate Pareto distributions (the Fréchet domain).
\item Let us consider the marginal distributions of $V_i,\ i=1,\ldots,p$ :
\begin{enumerate}
\item If $p=1$ then $V_1\sim Gamma(\alpha_1,1)$ as the matrix inverse becomes a simple reciprocal and $\rho_{11}=1$.
\item In the case where $\rho_{ij} = 0, i \ne j$, the $V_i$ are distributed independent $Gamma(\alpha_i,1)$, as all matrices become diagonal and the transformation then cancels out trivially.
\item In the general case where $\rho_{ij} \ne 0$ for some $i \ne j$, the marginal distributions of $V_i$ and $V_j$ do not appear to follow any standard distribution. See Appendix A for an explanation.
\end{enumerate}
\item Let $\ln V = \Sigma^{1/2}Y + \mu$, then $Y$ is distributed as a generalized multivariate gamma distribution with parameters $\mu , \Sigma , \alpha , \Lambda$ denoted by $GMGAM(\mu , \Sigma , \alpha , \Lambda)$, with  $E(Y) = \mu$, Cov$(Y,Y') = \Sigma$.
\item If $\mu = 0$ and $\Sigma =I_p$, refer to $Y \sim MGAM(\alpha,\rho)$ as the multivariate gamma distribution.
\item If $\rho_{ij} = 0, i \ne j$, then $Y \sim GMGAM(\alpha , \mu , \Sigma)$ defined in \textcite{deWaal2004}
\item To fit a MGAM distribution to composition data on $X$ in the $p$ dimensional simplex, let $Y = -\ln X$ and assume $Y \sim MGAM(\alpha, \rho)$.
\item This paper refer to the distribution of $\ln(Y)$ as a multivariate extreme value (MEV) distribution.
\item To simulate an $x$ value on $X$, one needs to go backwards in the definition above. The steps are:
\begin{enumerate} %[(a)]
\item Generate $w_i \in$ $Gamma(\alpha_i, 1)$ independently for each $i = 1, \ldots , p$.
\item Let $\mu_i = \ln(w_i)$ and standardise by subtracting the mean $\psi(\alpha_i)$ and dividing by $\sqrt{\psi'(\alpha_i)}$, namely $\mu_s = diag^{-1} (\sqrt{\psi'(\alpha )}) \{\mu - \psi(\alpha) \}$
\item Transform $\mu_s$ to have covariance $\Lambda$ and mean $\psi(\alpha)$, namely $v = \Lambda^{-1/2} \mu_s + \psi(\alpha)$.
\item Let $x = \exp(-v)$
\item Rescale $x$ such that $\sum_{i=1}^p x_i < 1$. One way to do this is to divide each element of $x$ by $p + \sum_{i=1}^p x_i . \;  x_1, \ldots , x_p$ becomes the required generated composition observation in the $p$ dimensional simplex.

\end{enumerate} 
\end{enumerate}

\section{Two dimensional MGAM}

As an illustration data from a MGAM distribution in two dimensions is simulated. \autoref{fig:MEVscatter1} shows scatter plots of two sets of data simulated from a MEV distribution with different $\alpha$'s and $\rho$'s. The data are heavy tailed and the marginals are skew distributed. The figures highlight the ability of this distribution to capture some non-traditional relationships.

\begin{figure}[hbt]
\begin{center}

   \includegraphics[width=\figwidth]{1}
\end{center}
\caption{Simulations from two dimensional MGAM distributions}

Two sets of size $n = 200$ simulated from two dimensional MGAM distributions with $\alpha$ = [2 4], and $\rho_{12} = 0.7$ and $-0.7$ respectively.
\label{fig:MEVscatter1}
\end{figure}

\section{Application of the MGAM}

Consider a $MGAM(\alpha,\rho)$ fit to the composition data on sizes of coal in a five dimensional simplex. Let $y = -\ln(x)$ where $x$ is the $(n = 130, p = 5)$ data matrix and it is desired to fit a $MGAM(\alpha,\rho)$ to $y$. It is however advisable to rather fit a $MEV(\alpha, \Lambda)$ distribution to $z = \ln(y)$. The moment estimates are simple to obtain. Since $E\ln(Y)  = \psi(\alpha)$, it can be estimated by the mean of $z = \ln(y)$, say $\tilde{z}$. From the inverse function $\hat{\alpha} =  \psi^{-1} (\tilde{y})$, estimates of the shape parameters can be obtained by numerically calculating the inverse values. Also, since the covariance of  $Z = \ln(Y)$ is $\Lambda$, it is easy to estimate. An estimate of $\rho$ can be obtained instead. The estimates are:

\[
\begin{split}
\hat{\alpha} &= \begin{bmatrix} 1.4700  &  2.0700  &  1.8800 &   4.2100 &   2.7700 \end{bmatrix} \\
\hat{\rho} &= \begin{bmatrix}
1.0000 &  -0.7716 &  -0.4760 &  -0.0475 &   -0.3735 \\
   -0.7716 &   1.0000 &  -0.0954 &   0.0808 &   0.2949 \\
   -0.4760  & -0.0954  &  1.0000 &  -0.1252 &  -0.1743 \\
   -0.0475  &  0.0808  & -0.1252  &  1.0000 &  -0.1160 \\
   -0.3735  &  0.2949 &  -0.1743 &  -0.1160 &   1.0000
\end{bmatrix}
\end{split}
\]

The estimate of $\Lambda$ becomes $\hat{\Lambda} = \left(\hat{\rho}_{ij} \sqrt{\tilde{z}_i \tilde{z}_j}\right),  i,j = 1, \ldots ,p$.

According to the definition of the MGAM, ones needs to transform the data $y$ to have marginals $Gamma(\alpha,1)$. Therefore, transform $Y_i \sim Gamma(a,b)$ to $Y_i \sim Gamma(\alpha_i,1)$ independently for all $i = 1, \ldots , p$, and note that $Z = \ln Y \sim MEV(\alpha,\rho)$. Thereafter, proceed with the analysis on $z$.

This allows the researcher to generate $y$ from $MGAM(\alpha,\rho)$ or $z$ from $MEV(\alpha,\rho)$. The Mahalanobis distance measure 

\begin{equation}
D(z) = \frac{1}{n} trace[(z_0 - z)' (z_0 - z)].
\end{equation}

between an observed matrix $z$ and a generated $z_0$ from a given model can be used to select $\hat{\alpha}$. The closer $z_0$ gets to $z$, that is, the smaller $D(z)$, the better the $\hat{\alpha}$.

Since $D(z)$ varies from simulation to simulation, 500 $z_0$ values were generated and the distribution of $D(z)$ considered. \autoref{fig:EDFsDs} shows the empirical cdf of such $D(z)$ values for the above selected $\hat{\alpha}$. It also serves as a goodness-of-fit check in the sense that it provides the $\hat{\alpha}$ where the generated $z_0$ is close to the observed $z$.

\begin{figure}
\begin{center}

   \includegraphics[width=\figwidth]{2}
\end{center}
\caption{Empirical cdf's of repetitions on distances between $z$ and a simulated dataset $z_0$ with  simulations from MEV}

Empirical cdf's of 30 repetitions on distances $D$ between $z$ and a simulated dataset $z_0$ with 500 simulations from $MEV(\alpha,\rho$). Dotted cdf on $D$ is obtained for the observed $z$.
 \label{fig:EDFsDs}
\end{figure}

\subsection{Classification of new observations under the MGAM model.}

Consider now the classification of a new observation $z$ under the MEV model between three groups H, S and L on coal sizes. A measure to classify a new observation as H, S or L is needed. As a measure, consider the linear function $w = \beta_0 + \beta_1 z_1 + \ldots +\beta_p z_p$ with $z_i$ the MEV variables and $w$ the percentiles of the cdf of the observed gas production. The sample space of $w$ will then be divided into three regions representing H, S and L.

Consider the classification of a new observation $z$ on coal sizes obtained through the transformations discussed in \autoref{sec:MGAM} under the assumption that $Z \sim MEV(\alpha, \rho)$. This study wishes to classify the observation as H, S or L.

Estimates of $\beta$'s defined in the linear function $w = \beta [1 \, z_1 \ldots z_p]^T$ are obtained through least squares minimization. From the dataset of $n = 112$  observations, one obtains
$$\hat{\beta} = \begin{bmatrix} 1.5467  &  0.2394 &  -0.2877 &   0.4000  &  0.1341 &   0.2767 \end{bmatrix}$$
The $R^2$ is equal to 0.48 for the linear model. This is increased to 0.56 if interaction terms are included in the linear function. The specific value of $R^2$ for this example is not meaningful on its own; it is used merely as a starting point for model comparison. Fitting more complex models in an attempt to raise this value is possible but outside the scope of this article.

The empirical cdf of the calculated $w$ values is shown in \autoref{fig:EDFW}. Selecting thresholds on the sample space of $w$, one can specify classification rules: say if $w > 2.3$ classify the observation as L, $w < 1.6$ classify the observation as H and otherwise classify it as S. The thresholds can be chosen optimally by striving for the outcome on the classification accuracy to be as close as possible to the real outcomes H = 20, S = 76 and L = 16. The classification results are H = 24, S = 67 and L = 21 of which 13 H were classified correctly, 55 S were correct and 11 L observations classified correctly. Adjusting the thresholds one may get improved classifications.  The success rate is 71\% . Although the multiple linear regression model does not fit that well on the MEV variables, it does contain information on the relationship between the coal sizes and the gas production.

\begin{figure}[hbt]
\begin{center}

   \includegraphics[width=\figwidth]{3}
\end{center}
\caption{Empirical cdf of $W$ with thresholds dividing the sample space}

\label{fig:EDFW}
\end{figure}

\section{Logistic normal}

Under the logistic normal model, assume that $X \sim Logistic\ Normal(\mu,\Sigma)$ or $Y \sim N(\mu,\Sigma)$ where $Y_i = \ln\left( \frac{X_i}{1-\sum_{i=1}^{P-1} X_i} \right).$ It is known that the predictive distribution of a future $Y_0$ is multivariate $t$ with $n-p$ degrees of freedom, mean $\bar{y}$ and precision matrix $T = \frac{n(n-p)}{n^2 -1} S^{-1}$ \parencite{geisser1983prediction}. $S = \frac{1}{n-1} \sum_{j=1}^n (y_j - \bar{y})'(y_j - \bar{y})$ is the estimated covariance matrix of $Y$.

Considering the dataset $y$ of size 124 by 5 on coal sizes, the estimates of $\mu$ and $\Sigma$ are

$$\bar{y} = \begin{bmatrix} 2.7617 &   2.1730 &  2.3919  &  0.0584  &  1.5043  \end{bmatrix}$$
and
$$ S =  \begin{bmatrix} 
0.6416 &   0.1496 &   0.3571 &   0.3134 &   0.4040 \\
  0.1496  &  0.6737 &   0.3751 &   0.2859 &   0.5625 \\
  0.3571 &   0.3751 &   0.7714 &   0.3239 &   0.4942 \\
 0.3134  &  0.2859 &   0.3239  &  0.8963 &   0.3551 \\
  0.4040  &  0.5625 &   0.4942 &   0.3551  &  1.1189
\end{bmatrix} $$

A goodness-of-fit test is employed with measure $D(y) = \frac{1}{n} trace[(y_0 - y)' (y_0 - y)]$. 

If $y$ denotes the data and $y_0$ simulated from a multivariate normal with parameters  above and repeated 500 times, the cdf of $D(y)$ is shown in \autoref{fig:EDFDsim}. 

\begin{figure}
\begin{center}

   \includegraphics[width=\figwidth]{4}
\end{center}
\caption{Empirical cdf of $D(y)$ shown by dots based on 500 simulations}

Lines show a few repetitions if normal data are used instead of real data.
\label{fig:EDFDsim}
\end{figure}

The logistic normal fit is also acceptable, as shown in \autoref{fig:EDFDsim} from a few repetitions of simulated data instead of the data.

Alternatively, similar to the distance measure above, one can specify
$$U(Y_0) = \frac{1}{p} (Y_0-\bar{y})' T(Y_0-\bar{y}) \sim F(p,n-p).$$

The statistic $U$ can also be considered as a distance measure between $Y_0$ and the mean $\bar{y}$. From the distribution of $U$, 95\% critical values can be obtained to be able to make a decision on judging the validity of the model.

One cannot conclude whether the MGAM fits better than the LN. The objective of the paper is to show that the MGAM is an alternative distribution to consider and especially for heavier tailed data it may be more appropriate. For the classification outcomes, the logistic normal model also gives good results in the following application.

\subsection{Classifying a new observation under the LN model}

This paper considers the classification of a new observations under the logistic normal model, using the linear function $w = \beta_0 + \beta_1 y_1 + \ldots + \beta_p y_p$ for classification. \autoref{fig:EDFwThresholds} shows the empirical cdf of $w$ for the observed data with thresholds chosen to classify an observation between H, S and L. As a classification rule, if $w < 1.7$, the observation is classified as H, if $w > 2.26$ as L and as S otherwise. These thresholds can be optimised by minimizing the difference between the outcomes of the classification on all the data and the observed, namely  20 H, 76 S and 16 L observations.

The two cdf's of $w$, the one for the MGAM (\autoref{fig:EDFW}) and the other for the logistic normal (\autoref{fig:EDFwThresholds}), are very similar except that the cdf corresponding to the MGAM model has a heavier tail for H observations.

\begin{figure}[hbt]
\begin{center}
   \includegraphics[width=\figwidth]{5}
\end{center}
\caption{Empirical cdf of $w$ with lines showing the thresholds}

\label{fig:EDFwThresholds}
\end{figure}

The classification matrix in \autoref{tab:clasA} shows the outcomes.

\begin{table}[hbt]

 \begin{center}
 \caption{Classification outcome matrix A}
 \label{tab:clasA}
   \begin{tabular}{| l | c  | c | c |}
   \hline
   Observed\ \textbackslash\ Assigned & H & S & L \\ \hline
 H &  13 & 7 & 0 \\ \hline
  S & 13 & 49 & 4 \\ \hline
  L & 0 & 4 & 11 \\
   \hline
   \end{tabular}
 \end{center}

\end{table}

The percentage correct classification is 65\%. This is lower than the outcome under the MGAM model, but by changing the thresholds, this can be higher.

\section{Conclusion}

In conclusion, the MGAM and logistic normal give very similar classification results on our dataset as a whole. Both models seem to fit well. The largest eigenvalues on the square distances between observed and simulated values as a classification measure is giving satisfactorily results, but other measures can be used. The choice of the thresholds is also debatable and other choices than squared determinant differences may be appropriate. The assignment of new observations to H, S or L can  be done under both models with at least 55\%  confidence on properties and sizes.

\subsection*{Notes}

All analyses in this chapter were programmed in MATLAB\textregistered\ \parencite{MATLAB}.

This chapter has been edited from the stand-alone version, which may be viewed at this link: 
\url{https://journals.co.za/content/sasj/50/2/EJC194760}

\section*{Appendix A: Marginal and conditional distributions}

\subsection*{Marginal distributions in the general case}
To understand the case where there is non-zero correlation, consider $p=2$ and proceed through the following derivation:

First, let R be the matrix of correlation parameters $[\rho_{ij}]=\begin{bmatrix}1&\rho\\ \rho & 1\end{bmatrix}$, then $\Lambda=H^{1/2} R H^{1/2}$.
$$\therefore D= R^{1/2} H^{1/2}\text{, where }R^{1/2}=\frac{1}{\sqrt{2+2\sqrt{1-\rho^2}}}\begin{bmatrix}
1 + \sqrt{1-\rho^2} & \rho \\
\rho & 1 + \sqrt{1-\rho^2}
\end{bmatrix} $$
And so, defining $r=1 + \sqrt{1-\rho^2}$,
$$\begin{aligned}
H^{1/2} D^{-1}&=(2r)^{-0.5}\begin{bmatrix}
r & -\rho \\
-\rho & r
\end{bmatrix}\\
\end{aligned} $$
$$\therefore W = (2r)^{-0.5}\begin{bmatrix}
r(\ln V_1-\psi(\alpha_1)) -\rho(\ln V_2-\psi(\alpha_2)) \\
-\rho(\ln V_1-\psi(\alpha_1)) + r(\ln V_2-\psi(\alpha_2))
\end{bmatrix} + \begin{bmatrix} \psi(\alpha_1) \\ \psi(\alpha_2) \end{bmatrix}
$$
$$\begin{aligned}
\text{and } Z &= \begin{bmatrix}
V_1^{(2r)^{-0.5}r} V_2^{-(2r)^{-0.5}\rho} e^{(2r)^{-0.5}[-r\psi(\alpha_1)+\psi(\alpha_2)\rho ]+\psi(\alpha_1)} \\
V_1^{-(2r)^{-0.5}\rho} V_2^{(2r)^{-0.5}r} e^{(2r)^{-0.5}[-r\psi(\alpha_2)+\psi(\alpha_1)\rho ]+\psi(\alpha_2)}
\end{bmatrix} \\
& = \begin{bmatrix}
V_1^{(2r)^{-0.5}r} V_2^{-(2r)^{-0.5}\rho} c_1 \\
V_1^{-(2r)^{-0.5}\rho} V_2^{(2r)^{-0.5}r} c_2
\end{bmatrix}
\end{aligned} $$
The joint density of $Z$ is $f_Z(\mathbf{z})=[\Gamma(\alpha_1)\Gamma(\alpha_2)]^{-1}z_1^{\alpha_1-1}z_2^{\alpha_2-1}\exp\{-z_1-z_2\}$, so

$$\begin{aligned}
f_V(\mathbf{v})&= c_3 v_1^{(2r)^{-0.5}[r(\alpha_1-1) - (\alpha_2-1)\rho]} v_2^{(2r)^{-0.5}[r(\alpha_2-1) - (\alpha_1-1)\rho]} \\
&\quad \times \exp\{-v_1^{(2r)^{-0.5}r} v_2^{-(2r)^{-0.5}\rho} c_1 - v_1^{-(2r)^{-0.5}\rho} v_2^{(2r)^{-0.5}r} c_2\} \\
&\quad \times  c_1 c_2 v_1^{(2r)^{-0.5}r - 2 -(2r)^{-0.5}\rho} v_2^{(2r)^{-0.5}r - 2 -(2r)^{-0.5}\rho} \times |c_4| \\
&= c_5 v_1^{(2r)^{-0.5}[r\alpha_1 - \alpha_2\rho]-1} v_2^{(2r)^{-0.5}[r\alpha_2 - \alpha_1\rho]-1} \\
&\quad \times \exp\{-c_1v_1^{(2r)^{-0.5}r} v_2^{-(2r)^{-0.5}\rho} - c_2v_1^{-(2r)^{-0.5}\rho} v_2^{(2r)^{-0.5}r}\}
\end{aligned}$$
$\therefore f_{v_1}(v_1)=\int_0^\infty av_2^b e^{-cv_2^d-gv_2^h}\ dv_2$ where $a,b,c,g,h$ are functions of $v_1$ and the parameters; while $d,h$ are just functions of $\rho$. This integral does not have a closed form solution for $\rho\ne 0,1$.

Note that if $\rho=0$ then $r=2$ and one again arrives at the product of independent gamma densities.

\subsection*{Conditional distributions}

Since the marginals do not have a closed form, neither do the conditional distributions. However, given specific values of one or more components, it is possible to evaluate the conditional densities empirically, up to an unknown constant. It is thus possible to simulate from the conditional distributions and then calculate desired quantities from the simulations.

\printbibliography[title=List of References for Paper 4]
\end{refsection}

















% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{refsection}
\chapter{Paper 5 --- A Method for Bayesian regression modelling of composition data} \label{ch:p5}
\pagestyle{C6}
\begin{LARGE}
Sean van~der~Merwe
\end{LARGE}
\vspace{10mm}

\textbf{Abstract}

Many scientific and industrial processes produce data that is best analysed as vectors of relative values, often called compositions or proportions. The Dirichlet distribution is a natural distribution to use for composition or proportion data. It has the advantage of a low number of parameters, making it the parsimonious choice in many cases. This paper considers the case where the outcome of a process is Dirichlet, dependent on one or more explanatory variables in a regression setting. The paper explores some existing approaches to this problem, and then introduces a new simulation approach to fitting such models, based on the Bayesian framework. The paper illustrates the advantages of the new approach through simulated examples and an application in sport science. These advantages include: increased accuracy of fit, increased power for inference, and the ability to introduce random effects without additional complexity in the analysis.

\noindent%
{\it Keywords:}  Composition data, proportions, Dirichlet distribution, regression, Bayes, simulation

\section{Introduction} \label{secIntro}

When vectors are measured in whole numbers, say based on a classification process, analysis is often based on the multinomial distribution, with the total count being seen as a nuisance parameter. In cases where the total count is not relevant it seems more natural to work directly with the observed proportions. Sometimes the proportions themselves are observed directly, rather than counts and totals. Any situation where the quantity surveyed or analysed does not affect the expected vector (but may affect the precision) falls in this class of problem.
 
The Dirichlet distribution has been widely accepted in literature for modelling composition data, subject to the constraint that all the correlations between variables are negative \citep{Maier2014}. A wider class of distributions which allows for positive correlations defined on the same sample space is the logistic-normal (LN) \citep{Aitchison1980}. However, the LN distribution has many parameters to estimate due to the unknown covariance parameter matrix. The MGAM model, discussed in \autoref{ch:p4}, suffers the same drawbacks. In contrast, the Dirichlet distribution has only $P$ unknown parameters to estimate for $P$ composition components. Further, when the model is expanded to the regression framework then the requirement of negative correlations no longer applies \citep{Maier2014}. 

The only remaining restriction of relevance is that observed values should be strictly positive. The occasional zero will not affect results, but the method cannot handle binary data in the dependent variable vectors. No restrictions apply to the explanatory variables --- they may have any structure.

As examples of applications, consider the work of \citet{boukal2015analyses}, who look at the development of insects, or the work of \citet{EspinGarcia2014}, who look at genetic analysis problems, or the work of \citet{smithson2006better}, who discuss applications in psychology. For an industrial application see \citet{DEWAAL2016} where the composition of coal is analysed. Another area in which Dirichlet regression may be useful is in politics, where only the proportion of voters supporting a set of candidates is of interest, possibly dependent on district or demographics. In ecology, preferences of animals for specific types of prey is of interest. To generalise, any situation where people choose between a set of options, and the researcher is interested in the choice and not the number of people making a choice, then the Dirichlet regression model may be of use.

The paper is outlined as follows: In \autoref{secDirichlet}, a brief introduction on applicable Dirichlet properties is given. After that \autoref{secDirReg} considers existing approaches to the problem of regression with dependent variables that follow the Dirichlet distribution (conditional on explanatory data). The new methodology is then introduced, in \autoref{newmethod}, and its usefulness is illustrated via simulated examples (\autoref{sims}). As an example based on observed data, movement data arising from a school netball tournament is analysed (\autoref{secRealdata}). Conclusions and future work are discussed in \autoref{secConclusion}.

\section{The Dirichlet distribution} \label{secDirichlet}

If $\mathbf{Y}$ is distributed Dirichlet$(\alpha_1,\ldots,\alpha_{P})$, denoted by $D(\boldsymbol\alpha)$, then the joint density is given by 
\begin{equation} \label{eq:D1}
f(\mathbf{y})= \frac{\prod_{i=1}^{P}\Gamma(\alpha_i)}{\Gamma(\alpha_0)}\left\{ \prod_{i=1}^P y_i^{\alpha_i-1}\right\}\ ,\ 0 < y_i < 1,\ \sum_{i=1}^P y_i = 1,\ \alpha_0 = \sum_{i=1}^{P} \alpha_i 
\end{equation}
\citet{Aitchison1986} calls this distribution the compositional Dirichlet defined on the specified simplex.

\citet{Wilks1962} and \citet{deGroot1970} provide detailed discussions on many of the properties of the Dirichlet distribution. Some relevant properties follow.

\begin{enumerate}
\item The means and covariances are
\begin{equation} \label{eq:cov}
\begin{aligned}
E[Y_i]&= \mu_i=\frac{\alpha_i}{\alpha_0}, i=1,\ldots,P \\
\sigma_{ij}&=\left\{\begin{array}{l}
\dfrac{-\alpha_i\alpha_j}{\alpha_0^2(1+ \alpha_0)}, i\neq j\vspace{1ex} \\
\dfrac{\alpha_i(\alpha_0 - \alpha_i)}{\alpha_0^2(1+ \alpha_0)}, i=j
\end{array}\right.
\end{aligned}
\end{equation}
 
$\boldsymbol\mu = (\mu_1,\ldots,\mu_P)$ %changed small p to P
 denotes the mean of the distribution and $\Sigma=(\sigma_{ij})\ i,j=1,\ldots,P$ the covariance matrix.

\item The marginal distribution of any subset $\mathbf{W}$ of $\mathbf{Y}$ is again Dirichlet$(\alpha_{\mathbf{W}}, \alpha_0 - \sum_{i\in {\mathbf{W}}} \alpha_i)$ \citep{Aitchison1986}.

\item In \autoref{eq:cov}, if $\boldsymbol\alpha$ is a multiple of $\boldsymbol\beta$ for two Dirichlet distributions $D(\boldsymbol\alpha)$ and $D(\boldsymbol\beta)$, the means are the same, but the covariance matrices differ.

\item Fitting the Dirichlet distribution to a sample of reasonable size can be done by following the method of moments or method of maximum likelihood (ML) as described in \citet{Minka}, or by using the Bayesian approach \autoref{sec:post}. Note that the ML approach is iterative, while the Bayesian approach uses simulation. Both methods are fast in most cases due to the neat convex nature of the likelihood.

\end{enumerate}

\section{Dirichlet Regression} \label{secDirReg}

In \citet{Carmargo2012} they define the problem as follows:

Let $Y=(\by_{1\cdot };\by_{2\cdot};\dots ;\by_{n\cdot})$ be a sample of vector observations. They use the $i\cdot$ notation to indicate that the vectors are arranged in rows of the matrix $Y$, for practical convenience. Let $X=(\bx_{1\cdot};\bx_{2\cdot};\dots ;\bx_{n\cdot})$ be $Q$ explanatory variables arranged the same way (each column of $X$ is a component and each row of $X$ corresponding with the same row of $Y$). Recall that $\sum_{j=1}^P y_{ij}=1,\ y_{ij}>0$; while the values of $X$ could be any real numbers. Their notation is used going forward.

Based on the work of \citet{campbell1987multivariate}, \citet{Hijazi2009} and \citet{Carmargo2012} each parameter is modelled as a function of the explanatory variables. In terms of a single observation $\by_{i\cdot}=(y_{i1},\dots,y_{iP})\sim D(\alpha_{i1},\ldots,\alpha_{iP})$,
\begin{equation} \label{eq:basicDirReg}
\alpha_{ij}=x_{i1}\beta_{1j}+\dots+x_{iQ}\beta_{Qj}=\bx_{i\cdot}\boldsymbol\beta_{\cdot j}
\end{equation}

Thus, the parameters to be estimated are all the $\beta_{kj},\ k=1,\ldots,Q,\ j=1,\ldots,P$, subject to the constraint $\alpha_{ij}>0\ \forall\ i=1,\ldots,n,\ j=1,\ldots,P$. They describe a custom optimisation procedure to estimate these parameters under these constraints. Finally, \citet{Carmargo2012} explain an approach to testing $\beta_{kj}=0$, which is useful in many problems.

\citet{gueorguieva2008} propose using a log-link in each dimension, thus eliminating constraints in the optimisation procedure.

However, each $\beta_{kj}$ does not have a clean interpretation in the above models as $E[Y_{ij}]$ is a function of all $\beta_{kj}$. The difficulty of interpretation is seen as a major drawback by many researchers, and led to the investigation of alternatives.

\citet{Maier2014} applies a multivariate transformation to the parameters of the Dirichlet distribution, arriving at an alternative formulation that has the advantage of modelling the expected value of an observation separately from its precision. He begins by defining a parameter $\phi=\alpha_0$ to denote the precision of an observation. Looking at \autoref{eq:cov} it can be seen that $\phi$ is not exactly the precision but acts like the precision in the sense that, for large values, if the value of $\phi$ doubles while the mean vector is unchanged then the variance halves. Note the relationship $\boldsymbol\alpha=\boldsymbol\mu \phi$. \citet{Maier2014} applies a log-link function to model $\phi$, \textit{i.e.}\ $\ln \phi_i = \mathbf{w}_{i\cdot}\boldsymbol\beta_{\cdot\phi}$, where $W$ is a matrix of explanatory variables for the precision.

For the purpose of modelling the mean \citet{Maier2014} uses a multivariate logit link. This involves choosing a base category and setting all coefficients to zero for this category, then using a log-link to model the other categories and rescaling the results so that the means sum to one. In the notation defined previously, and using dimension 1 as the base, for every observation $i=1,\ldots,n$,
\begin{equation} \label{eq:multilogit}
\begin{aligned}
\mu_{ij}&=\frac{\exp (\bx_{i\cdot}\boldsymbol\beta_{\cdot j})}{\sum_{k=1}^P \exp (\bx_{i\cdot}\boldsymbol\beta_{\cdot k})} \ ,\ j=2,\ldots, P \\
\mu_{i1}&=\frac{\exp \bx_{i\cdot}\mathbf{0}}{\sum_{k=1}^P \exp \bx_{i\cdot}\boldsymbol\beta_{\cdot k}} = \frac{1}{1+\sum_{k=2}^P \exp \bx_{i\cdot}\boldsymbol\beta_{\cdot k}}
\end{aligned}
\end{equation}

\citet{Maier2014} explains that using the transformation above results in coefficients that are interpretable as odds ratios if exponentiated.

However, each $\mu_{ij}$ is still a function of all $\beta_{kj}$, and there are no coefficients for the base dimension. These limitations inhibit interpretation. The next section introduces a new methodology that incorporates the best aspects of the approaches described above.

\section{New methodology} \label{newmethod}

The first change is the use of a univariate logit transformation for each mean parameter individually, thus allowing all $\beta_{kj}$ coefficients to be unrestricted real numbers. The second change that ties in with the first is that this new methodology abandons the idea of a reference category or dimension.

In theory, one dimension is redundant since it is a linear combination of the others, but in practice in is of interest to know the relationship between the explanatory variables and the outcomes in all dimensions. Often the dimensions are equal in the view of the researcher and it is not sensible to relegate one to reference status. It is for this reason that researchers such as \citet{chen2016} resort to modelling each dimension individually as beta distributed random variables, but that in turn ignores the multivariate nature of the data.

It is desirable to have each $\beta_{kj}$ relate directly to a single dimension in a way that can be directly interpreted. By modelling all dimensions the researcher moves closer to this ideal.

The third change is a move to the Bayesian framework. This involves introducing vague normal priors on all $\beta$ parameters. All other parameters are defined in terms of these, so no further priors are necessary at this stage.

Specifically, this method uses the Bayesian simulation framework, which holds many advantages. It allows us to directly quantify uncertainty in both the coefficients and the means. Also, when moving to a predictive framework, construction of predictive densities is relatively straightforward.

The method is still bound by the conditions 
\begin{equation} \label{eq:restrict}
\sum_{j=1}^P \mu_{ij} = 1,\ \forall\ i=1,\ldots ,n
\end{equation}
which impede the standard simulation approach greatly. In order to have the simulation process run smoothly, the researcher must introduce a source of flexibility into the model. This paper adds flexibility by replacing the restriction (\autoref{eq:restrict}) by a penalty on the likelihood:
\begin{equation} \label{eq:penalty}
L^{*}\propto L\times \exp\left\{-\frac{1}{\xi}\sum_{i=1}^n \left[\left(\sum_{j=1}^P \mu_{ij}\right) -1\right]^2\right\}
\end{equation}

The hyperparameter $\xi$ must be chosen large enough to allow the simulation procedure to run smoothly, but small enough to have minimal impact on the simulation results. Here minimal impact implies mean deviations that can easily be corrected. The valid region for $\xi$ to meet these criteria seems surprisingly large. One may consider $\xi$ as a hyperparameter and choose its value manually, or, more conveniently, apply a prior distribution to $\xi$ and have the value vary as part of the simulation process.

Then, for further simulation flexibility a second penalty parameter $(\xi^{*})$ is introduced in the relationship between $\boldsymbol\alpha$, $\boldsymbol\mu$ and $\phi$. For both penalty parameters a simple exponential prior works well in all scenarios tested.

Assuming explanatory data captured in matrices $X$ and $W$ (which may overlap), the model is defined in a hierarchical fashion:
\begin{equation} \label{eq:hierarchy}
\begin{aligned}
\by_{i\cdot}|\bx_{i\cdot},\mathbf{w}_{i\cdot} &\sim D(\boldsymbol\alpha_{i\cdot}) \\ % added idot to alpha
\ln \alpha_{ij} &\sim N\left(\ln \mu_{ij} + \ln \phi_i\ ,\ \frac{1}{\xi^{*}}\right) \\
\ln \phi_i &= \mathbf{w}_{i\cdot}\boldsymbol\beta_{\cdot\phi} \\
\text{logit}(\mu_{ij}) &= \bx_{i\cdot}\boldsymbol\beta_{\cdot j} \\
\sum_{j=1}^P \mu_{ij} &\sim N\left(1,\frac{1}{\xi} \right) \\
\beta_{ij}\ ,\ \beta_{i\phi} &\sim N(0,10000) \\
\xi^{*} &\sim Exp(\mu=100/P) \\
\xi &\sim Exp(\mu=1000/P) \\
\end{aligned}
\end{equation}
Note that in the expressions above, the models for the mean $(\text{logit}(\mu_{ij}) = \bx_{i\cdot}\boldsymbol\beta_{\cdot j})$ and precision $(\ln \phi_i = \mathbf{w}_{i\cdot}\boldsymbol\beta_{\cdot\phi})$ are in linear form for easy of understanding only --- these models can be extended as needed by the researcher.

Looking only at the likelihood it might be said that the model is over-parameterised or lacking identifiability; however, the model is identifiable in the Bayesian sense \citep{rannala2002}. All parameters are given proper priors, ensuring a proper posterior, and use of informative priors on the penalty parameters ensures that the posterior distributions of the parameters of interest reflect the information in the data. In a practical sense the identifiability of the model is improved over previous models as there is now exactly one coefficient joining each component of the explanatory data to each component of the observation vectors. 

This paper implements the model using Gibbs sampling \citep{Gelfand1990} via the OpenBUGS program (\url{http://www.openbugs.net/}). Implementation is done indirectly through the R2OpenBUGS package \citep{R2OpenBUGS} for R \citep{RCore}. All pre- and post-calculations are done in R. The MASS and parallel packages supplied with R were also used to facilitate calculations.

Post-simulation, the following corrections are applied to each simulated parameter set $k$ individually to ensure that fitted expected values sum to one for each observation:
\begin{equation} \label{eq:corrections}
\begin{aligned}
\mu^{adj}_{ijk} &= \frac{\mu^{sim}_{ijk}}{\sum_{j=1}^P \mu^{sim}_{ijk}} \\
\alpha^{adj}_{ijk} &= \mu^{adj}_{ijk}*\phi^{sim}_{ik}
\end{aligned}
\end{equation}


\section{Simulation Study} \label{sims}

Since the different methods discussed in the previous sections use different transformations (identity, log, multivariate logit, and univariate logit), the methods are compared on a single scale. In general the researcher fitting Dirichlet regression models is interested in three things: the significance of the coefficients, the direction of any significant relationships, and the accuracy of the model fit on the observed data. The first point of focus is model accuracy.

As a measurement of error, consider the average \textit{sum of composition errors} (SCE), explained by \citet{Hijazi2009}. It is the sum of the Aitchison distances \citep{Aitchison1986} between estimated compositions and the target values.

Further, it is sensible to calculate the intervals for each expected value individually and then report the average coverage, along with the average width (standardised by dividing by the expected values). Lastly, the posterior predictive distribution allows calculation of coverage.

Datasets are constructed from a model exactly in line with the `alternative' specification of \citet{Maier2014}. Models are then fitted using his DirichletReg package, as well as the new methodology. Models are correctly specified in all cases --- model misspecification is beyond the scope of this work.

\subsection{Simulation Scenario A} \label{sec:scenA}

Scenario A is a simple analogy to the MANOVA problem. Consider a single explanatory variable that is a factor with three levels.  A researcher might be interested in whether the mean vectors differ between the three groups created by the factor levels, under the assumption of constant variance.

Let the observed vectors have three dimensions. Set the coefficients for the first dimension to zero, and then use coefficients of (-0.9,0.6,1.2) for the second % changed from first to second
 dimension and (0.8,-1,0.5) for the third. % changed from second to third
  The inverse multivariate logit transformation is then applied to obtain the `true' expected values for every observation, and expanded to 20 observations per factor level ($n=60$ in total). The first step is creating a matrix of binary variables $(X)$ from an expansion of the explanatory factor. The transformation (explained in detail in \cite{Maier2014}) involves multiplying $X$ by each set of coefficients, and then exponentiating to obtain raw expected values, which are then standardised to sum to one for each observation.

The next step is to multiply by a chosen value for $\phi$. This paper uses $\phi=1$ for illustration. The effect of changing this value will be discussed after the results. Multiplying the expected values by $\phi$ yields a matrix of $\alpha_{ij}$ values, which is used to simulate hundreds of Dirichlet samples in the standard way.

After every sample is modelled and the results summarised, various statistics are produced. See \autoref{tb:ScenA} for the most important values. It is clear that the new methodology produces better fit.

\begin{table}[htb]\begin{center}
\caption{Fit statistics for Scenario A}
 \label{tb:ScenA}
\begin{tabular}{|c|c|c|c|} \hline
Scenario A & Target & Maier approach & New approach  \\ \hline
Error  & 0 & 19.59 & 18.38 \\ \hline
Coverage & 0.95 & 0.87 & 0.94 \\ \hline
Std. Width & 0 & 0.70 & 0.75 \\ \hline
\end{tabular}\end{center}
\end{table}

The next question of interest would be the effect of varying the chosen parameter values. It appears that the important parameter is the underlying precision, which is closely related to the $\phi$ parameter. When the data are measured with high precision (say $\phi\geq 5$) then there is little difference between the fits obtained by the two methods (both methods fit very well). As the precision drops and the underlying relationships become more obfuscated then the accuracy of the previous methodology falls away rapidly, while the new methodology loses accuracy slowly, resulting a the difference observed in \autoref{tb:ScenA}.

Another question of interest would be the effect of increasing the dimensionality of the problem. If the dimension is increased to 8 and the coefficients are chosen as random $U(-1.5,1.5)$ values, then the relative differences between the methods become even more exaggerated. The new methodology adapts easily to having many categories in the dependent variable.

\subsection{Simulation Scenario B} \label{sec:scenB}

Scenario B is a more complex scenario where a linear term is introduced in every mean vector as well as the precision model, in addition to the factor explained in \autoref{sec:scenA}. The goal is to determine whether the model can identify both relationships simultaneously in all categories.

The explanatory factor $(X_1)$ is given two levels with coefficients (0,0), (-0.9,0.6), (1.8,-1) in three dimensions respectively. The explanatory variable with a linear relationship $(X_2)$ is given real values between 4.5 and 7.5 with 40 values per level of the explanatory factor. The linear relationship is created by adding $0.75X_2$ to the third dimension and then correcting the means to add to one, thus creating a positive relationship in the third dimension and implicitly creating negative relationships in the first two dimensions. As $X_2$ increases $y_1$ and $y_2$ will tend to decrease, while $y_3$ will tend to increase. The expression for $\ln \phi$ used to generate the data is $-1+0.5X_2$.

The results are averaged from hundreds of samples and summarised. As indicated in \autoref{tb:ScenB}, the results from the methods are similar, but the new methodology is more accurate. 

\begin{table}[htb]\begin{center}
\caption{Fit statistics for Scenario B}
 \label{tb:ScenB}
\begin{tabular}{|c|c|c|c|} \hline
Scenario B & Target & Maier approach & New approach  \\ \hline
Error  & 0 & 19.19 & 18.81 \\ \hline
Coverage & 0.95 & 0.85 & 0.86 \\ \hline
Std. Width & 0 & 0.52 & 0.54 \\ \hline
\end{tabular}\end{center}

\end{table}

% There were two paragraphs here, they have been combined into one for clarity.
Concerning inference, the new method shows a major improvement over the existing method. In the precision model the median p-value over simulated samples for the existing method is 0.1\%, while for the new method it is approximately 0\%. In the second and third dimensions the median p-values for the existing method are 50\% and 24\%, while the new method yields 1\% and 0.1\%. The new method also reports a p-value for the first dimension $(\text{median }0.4\%)$, while the existing method does not. Thus, for almost all simulations done, the new method correctly identifies the direction of the linear relationships and marks all of them as significant. The existing method only identifies the linear relationship in the precision model.

\section{Observed data example from sport science} \label{secRealdata}

During a school netball tournament, scholars were tracked accurately as they move across the field. One of the resulting measurements was the proportion of time spent standing/walking/running during the course of a match. The goal is to investigate the relationship between these measurements and the playing position. 

A major complicating factor is the fitness and behaviour variation between players. Some players were observed for only one match, while others were observed for up to nine matches. This suggests an unbalanced mixed effects model, with position as fixed effect and player as random effect.

The model is defined explicitly as follows:
\begin{equation} \label{eq:netballmodel}
\begin{aligned}
\by_{i\cdot}|x_{i},w_{i} &\sim D(\alpha_{i1},\alpha_{i2},\alpha_{i3}) \\
x_i&=\text{position for observation }i \\
w_i&=\text{player for observation }i \\
\ln \alpha_{ij} &\sim N\left(\ln \mu_{ij} + \ln\phi_{x_i}\ ,\ \frac{1}{\xi^{*}}\right) \\
\text{logit}(\mu_{ij}) &= \beta_{x_i j}+\beta_{w_i j} \\
\sum_{j=1}^3 \mu_{ij} &\sim N\left(1,\frac{1}{\xi} \right) \\
\end{aligned}
\end{equation}
With priors:
\begin{equation} \label{eq:netballmodelp2}
\begin{aligned}
\ln\phi_{x_i} &\sim N(0,10000) \\
\beta_{x_i j} &\sim N(0,10000) \\
\beta_{w_i j} &\sim N(0,\tau^{-1}_{j}) \\
\tau_j &\sim Gamma(0.001,0.001) \\
\xi^{*} &\sim Exp(\mu=100/3) \\
\xi &\sim Exp(\mu=1000/3) \\
\end{aligned}
\end{equation}

For additional clarity, where there is a second subscript or the letter $j$ these refer to the category: standing=1, walking=2, running=3. The model as implemented has exactly one precision parameter for each playing position, but three intercepts for each playing position corresponding to the three categories. The model has three random effects for each player, corresponding to the three categories. These random effects have a single common variance within each category (three in total).

The results are summarised in terms of posterior means and 95\% credibility intervals for each position, by category. \autoref{fig:netball} shows this result with the positions sorted from highest to lowest expected running proportion. This sorting makes it easier to identify differences between positions. The information in this figure may be of use to netball coaches going forward.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.9\textwidth]{Netball11.pdf}
\end{center}
\caption{Differences in movement between playing positions in a netball tournament as suggested by Bayesian unbalanced mixed effects model.}
\label{fig:netball}
\end{figure}

Code with detailed comments for this example is available at \url{http://seanvdm.co.za/files/ProportionsExample20180704.pdf}, along with its accompanying dataset (anonymised): \url{http://seanvdm.co.za/files/NetballData.csv}.

%***** Most of example changed

%We begin by implementing simple approaches that might be used to analyse this data. This is to illustrate the practical differences between approaches. First, we fit linear models on the logit scale to each dimension separately. This is essentially the same as using descriptive statistics to analyse the data. It fails to account for the multivariate nature of the observations and for the player effect. Second, we fit mixed effects models in each dimension (on the logit scale) to account for the player effect. The mixed effects model fitted did not account for the multivariate nature of the data.
%
%Third, we implement the method of \citet{Maier2014}. This accounts for the multivariate structure, but not the player effect; and then finally we implement the new approach, which takes into account both the structure of the data and the player effect at once.
%
%In \autoref{fig:netball} we see the different methods discussed and the differences between them. We note that the univariate mixed effects model appears to distort the expected values slightly, while the models without the random effects may distort the intervals. The results produced by the new methodology appear sensible in both the expected values and intervals. It appears to combine the best aspects of each of the simpler models, with none of the weaknesses, and allows for unbalance intervals which are often shorter for the same coverage.
%
%\begin{figure}
%\caption{Differences in movement between playing positions in a netball tournament as suggested by different modelling approaches}
%\begin{center}
%
%\includegraphics[width=0.9\textwidth]{Netball20180611.pdf}
% %Graph labels updated

%\end{center}
%
% \label{fig:netball}
%\end{figure}

\section{Conclusions} \label{secConclusion}

In this paper the goal was a generic solution to the problem where the outcome of a process is Dirichlet, dependent on one or more explanatory variables in a regression setting. Existing approaches were discussed and a new methodology introduced. The new methodology was directly compared to the latest of the existing approaches and found to perform well. At worst the performance is in line with existing tools, but in many cases the improvement is remarkable, especially when the data has high variance. Advantages of the new methodology were discussed, including ease of interpretation and prediction, with accurate intervals, as well as the ability to introduce random effects.

\section*{Acknowledgements}

The netball data are from a study by Michael Shaw and Derik Coetzee of the Department of Exercise and Sports Science, University of the Free State, South Africa. The author wishes to thank Robert Schall and Michael von Maltitz for useful advice during the development of this work.

\textit{Note:} This chapter has been edited from the stand-alone version, which may be viewed at this link: 
\url{https://arxiv.org/abs/1801.02954}

\printbibliography[title=List of References for Paper 5]
\end{refsection}














% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{refsection}
\chapter{Conclusion} \label{ch:conclusion}
\pagestyle{C7}

\section{General conclusion}

In the preceding chapters, advances were made in the analysis of composition data. Advances were made in direct analysis, transformed analysis, classification analysis and regression, of data vectors summing to one.

Composition data is prevalent in many applications, but the theory is far from fully developed. Data and related problems were presented by researchers in industry and academia. In attempting to address these problems, the inadequacies in the theory of analysing such data were realised and new theory and methods were developed. The hope is that future researchers working with composition data will be able to refer to this work as a starting point in their efforts in solving their own problems. 

\section{Conclusions from each chapter}

\autoref{ch:intro} introduces the subject area of focus. It is clear from the sources cited that the applications of the work in the later chapters are numerous and widely spread.

\autoref{ch:p1} expands the theory related to the use of the Dirichlet distribution. Contributions made by this work include the full derivation of objective prior distributions for both the standard Dirichlet distribution as well as the first adaptation (the Dirichlet Type II distribution). The conditional distributions are derived and their uses elaborated on. A technique of imputing missing values using these distributions is also developed and explained. A simulation study found that the Bayesian estimation using the Jeffreys prior was superior.

In \autoref{ch:p2}, discriminating between Dirichlet populations is considered. The extension to matrix Dirichlet populations is shown. A relationship between the discriminant function and the test statistic for equality of covariance matrices is illustrated.

\autoref{ch:p3} develops two classification schemes using the Dirichlet distribution. The first is a linear partitioning of the Dirichlet sample space. The second is a clustering approach based on minimising a chosen metric. The ease of use and flexibility of the schemes are shown.

\autoref{ch:p4} continues to expand the theory of classification of composition data. Multivariate gamma and multivariate extreme value distributions are introduced and developed. The application of the new theory to industry problems is explained.

\autoref{ch:p5} develops and explains a new approach to the regression problem involving composition data. The new approach yields increased prediction accuracy due to the inclusion of a small number of additional parameters. The additional parameters allow additional flexibility and ease of interpretation, without weaker generalisation. This new approach does not suffer any of the disadvantages of the previous approaches, and yet is more straightforward to implement and use.

\section{Future research}

Future research in the fields discussed in these articles may include expanding on generalisations of the Dirichlet distribution, such as that of \textcite{wong1998}, or mixture models such as that of \textcite{bouguila2003}; but there are many unexplored topics in these areas that will no doubt gain interest over time as the analysis tools improve. One such topic is that of appropriate incorporation of subjective prior information into the models discussed in the previous chapters, be it from previous studies or expert knowledge.

Another avenue of future research is goodness-of-fit assessments for the Dirichlet distribution, both frequentist and Bayesian. Currently, composition models are chosen for their implementation properties, not for optimal fit, and there is room for new developments in these model choices.

\printbibliography[title=List of References for the concluding chapter]
\end{refsection}

\paragraph{Final note:} 
some code for this work is available at \url{http://seanvdm.co.za/project/thesisstuff/}. This includes the most important code for the objective Bayes methods and the regression example.

\end{document}
