Bayesian density regression for count data
Despite the increasing popularity of quantile regression models for continuous responses, models for count data have so far received little attention. The main quantile regression technique for count data involves adding uniform random noise or ``jittering'', thus overcoming the problem that the conditional quantile function is not a continuous function of the parameters of interest. Although jittering allows estimating the conditional quantiles, it has the drawback that, for small values of the response variable $Y,$ the added noise can have a large influence on the estimated quantiles. In addition, quantile regression can lead to ``crossing'' quantiles. We propose a Bayesian Dirichlet process (DP)-based approach to quantile regression for count data. The approach is based on an adaptive DP mixture (DPM) of COM-Poisson regression models and determines the quantiles by estimating the density of the data, thus eliminating all the aforementioned problems. Taking advantage of the exchange algorithm, the proposed MCMC algorithm can be applied to distributions on which the likelihood can only be computed up to a normalising constant.
compstat2014.tex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[11pt]{article}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\usepackage{fancyhdr}
%
\usepackage{authblk}

%\pagestyle{fancy}

\usepackage{enumitem}
\usepackage{graphicx,epsfig,amssymb,amstext,amsmath,amsthm}
% Packages graphicx,epsfig, [english]{babel}, amssymb,amstext,amsmath, amsthm are included by default
\makeatletter
\newenvironment{tablehere}
{\def\@captype{table}}
{}
\newenvironment{figurehere}
{\def\@captype{figure}}
{}
\makeatother

%\usepackage[]{natbib}
\usepackage{tikz}
\usepackage{graphics}
\usepackage{amsmath}
\usepackage{multicol}
\usepackage{amssymb}
\usepackage{float}
\usepackage{booktabs}
\newcommand{\bftab}{\fontseries{b}\selectfont}
\usepackage{multirow,tabularx}
\newcolumntype{Y}{>{\centering\arraybackslash}X}
\renewcommand{\arraystretch}{2}
%\renewcommand\refname{}
\newcommand{\bo}{\boldsymbol}
\newcommand{\argmin}[1]{\underset{#1}{\operatorname{arg}\,\operatorname{min}}\;}
\usepackage{mathtools}
\DeclarePairedDelimiter{\ceil}{\lceil}{\rceil}
\DeclarePairedDelimiter{\floor}{\lfloor}{\rfloor}
\begin{document}

%\mainmatter
%\pagestyle{proceedings}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Here comes your COMPSTAT 2014 manuscript %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% ----------------------------------------------------------------------------------------------------------------------
% Declaration of COMPSTAT Manuscript for HEADERS and MAIN TITLE (please select the appropriate example and modify it)
% \COMPSTATManuscript{Authors}{Long title (at most 15 words)}{Short title (for header)}

% Examples: (please modify the corresponding data as you wish)
% Single Author
%\COMPSTATManuscript{John Smith}{This is the long manuscript title: at most 15 words}{Short manuscript title (for header)}

% Multiple authors
%\COMPSTATManuscript{John Smith \textit{et al.}}{This is the long manuscript title: at most 15 words}{Short manuscript title (for header)}
% or (if enough room)
%\COMPSTATManuscript{Charalampos Chanialidis, Ludger Evers and Tereza Neocleous}{Bayesian density regression for count data}{Bayesian density regression for count data}
%-----------------------------------------------------------------------------------------------------------------------


%-----------------------------------------------------------------------------------------------------------------------
% Declaration of the different authors (please modify the corresponding data as you wish)
%\COMPSTATAuthor{Author Name}{Author Surname}{Institution}{Email}

% Example: (use a line for each author)
%\COMPSTATAuthor{Charalampos}{Chanialidis}{University of Glasgow}{c.chanialidis1@research.gla.ac.uk} %Add a line for each author
%\COMPSTATAuthor{Ludger}{Evers}{University of Glasgow}{ludger.evers@glasgow.ac.uk} %Add a line for each author
%\COMPSTATAuthor{Tereza}{Neocleous}{University of Glasgow}{tereza.neocleous@glasgow.ac.uk} %Add a line for each author
%-----------------------------------------------------------------------------------------------------------------------


%-----------------------------------------------------------------------------------------------------------------------
% Declaration of the abstract
%\COMPSTATAbstract{Abstract paragraph}
\title{Bayesian density regression for count data}
\author[1]{Charalampos Chanialidis}
\author[2]{Ludger Evers}
\author[3]{Tereza Neocleous}
\affil[1]{University of Glasgow, c.chanialidis1@research.gla.ac.uk}
\affil[2]{University of Glasgow, ludger.evers@glasgow.ac.uk}
\affil[3]{University of Glasgow, tereza.neocleous@glasgow.ac.uk}

%
\maketitle
% Example: (please modify the corresponding data as you wish)
\begin{abstract}
Despite the increasing popularity of quantile regression models for continuous responses, models for count data have so far received little attention. The main quantile regression technique for count data involves adding uniform random noise or ``jittering'', thus overcoming the problem that the conditional quantile function is not a continuous function of the parameters of interest. Although jittering allows estimating the conditional quantiles, it has the drawback that, for small values of the response variable $Y,$ the added noise can have a large influence on the estimated quantiles. In addition, quantile regression can lead to ``crossing'' quantiles. We propose a Bayesian Dirichlet process (DP)-based approach to quantile regression for count data. The approach is based on an adaptive DP mixture (DPM) of COM-Poisson regression models and determines the quantiles by estimating the density of the data, thus eliminating all the aforementioned problems. Taking advantage of the exchange algorithm, the proposed MCMC algorithm can be applied to distributions on which the likelihood can only be computed up to a normalising constant.
\end{abstract}
%-----------------------------------------------------------------------------------------------------------------------


%-----------------------------------------------------------------------------------------------------------------------
% Declaration of the keywords
%\COMPSTATKeywords{Keyword1, Keyword2, etc.}

% Example: (please modify the corresponding data as you wish)
%\COMPSTATKeywords{Quantile regression, Bayesian density regression, Bayesian nonparametrics, Dirichlet processes, COM-Poisson distribution}
%-----------------------------------------------------------------------------------------------------------------------



%-----------------------------------------------------------------------------------------------------------------------
%\section{Introduction}
% It is important to provide always a unique label (name of author:name of label) to avoid crashes with other authors labels.
%Say something
% \label{john.smith:basic}

% Please avoid to insert two consecutive headers, instead let every header to be followed by at least a short introductory paragraph. Note that subsections and subsubsections are not numbered, thus it is advisable to use only sections and subsections at most.

% Use the template \verb"author_COMPSTAT2014.tex" together with the COMPSTAT2014 proceedings class \verb"proceedings_COMPSTAT2014.cls" to style the various elements of your manuscript content to be incorporated in COMPSTAT 2014 proceedings.

% Please note that the first line of text after each heading is not indented, whereas the first lines of all subsequent paragraphs are indented. Please, do not attempt to modify that.

\section{Quantile regression}
Quantile regression was introduced as a nonparametric method for modelling a variable of interest as a function of covariates \cite{KB1978}. By estimating the conditional quantiles rather than the mean, it gives a more complete description of the conditional distribution of the response variable than least squares regression, and is especially relevant in certain types of applications. 
% Quantile regression has been applied to a lot of diverse areas such as economics, growth curve  estimation, modelling of ecological processes, and observing the effects of reducing class size on students; see \cite{cade2003,C1994,levin2001,wei2006}.

Consider a random variable $Y$ with cumulative distribution function $F(y)$. The $p$th quantile function of $Y$ is defined as \begin{equation}
Q(p)=\inf\{y \in \mathbb{R}:p\le F(y)\}\end{equation}
 and can be obtained by minimising the expected loss $E[\rho_p(Y-u)]$ with respect to $u$, where $\rho _{p}(y)=|y( p-I(y<0))|$. The $p$th sample quantile is obtained in a similar way by minimising $\sum_{i=1}^n \rho_p(y_i-u)$. 

Suppose that the $p$th conditional quantile function, $Q_Y(p|X=\boldsymbol{x})$, is a linear function of the predictors so that $ Q_Y(p|X=\bo{x})=X'\bo{\beta}_p.$
The parameter estimates $\bo{\hat{\beta}}_p$ are then obtained as 
\begin{equation}\bo{\hat{\beta}}_p=\argmin{\bo{\beta}_p \in \mathbb{R}^k}\sum_{i=1}^{n}\rho _{p}(Y-X'\bo{\beta}_p).\end{equation} 
A closed-form solution for this minimisation problem does not exist since the objective function is not differentiable at the origin, and it is solved using linear programming techniques \cite{B1998}. %Linear quantile regression can be generalised to non-linear quantile regression using the same approaches as used for classical least-squares regression.


%\paragraph{Paragraph.} Additionally you can add a header to a simple paragraph like this if you wish. Note that a paragraph is not indented but

%\subparagraph{A subparagraph} always is indented.

\subsection{Quantile regression for count data}
%Please use the \LaTeX\ automatism for all your cross-references. As an example of cross-reference see Section \ref{john.smith:basic}.
 The problem with applying quantile regression to count data is that the cumulative distribution function of the response variable is not continuous, resulting in quantiles that are not continuous, and which thus can not be expressed as a continuous function of the covariates. One way to overcome this problem is by adding uniform random noise (``jittering'') to the counts \cite{MS2002}. The general idea is to construct a continuous variable whose conditional quantiles have a one-to-one relationship with the conditional quantiles of the counts. Defining the new continuous variable $Z= Y + U$ where $Y$ is the count variable and $U$ is a uniform random variable in the interval $[0, 1)$, the conditional quantiles $Q_Z(p|X=\bo{x})= p + \exp(X'\bo{\beta}_p).$

The variable $Z$ is transformed in such a way that the new quantile function is linear in the parameters, i.e.$Q_{T(Z;p)}(p|X=\bo{x})=X'\bo{\beta}_p$ where

\begin{equation}
T(Z;p) = \left\{
\begin{array}{rl}\label{quantilestransform}
\log(Z-p) & \text{for } Z > p,\\
\log(\varsigma) & \text{for } Z \le p,
\end{array} \right.
\end{equation}
with $\varsigma$ being a small positive number. The parameters $\bo{\beta}_p$ are estimated by running a linear quantile regression of $T(Z;p)$ on $x$. % This is possible, since quantiles are invariant to monotonic transformations and to censoring from below up to the quantile of interest.
Finally, the conditional quantiles of interest, $Q_Y(p|X=\bo{x})$ can be obtained from  the previous quantiles as
 \begin{equation}Q_Y(p|X=\bo{x})=\ceil*{Q_Z(p|X=\bo{x})-1}\end{equation} where $\ceil*{p}$ denotes the ceiling function which returns the smallest integer greater than, or equal to, $p$.  %For more information in this method, see .

While the jittering approach eliminates the problem of a discrete response distribution, for small values of the response variable $Y$, the mean and the variance in the transformed variable $Z$ will be mainly due to the added noise, resulting in poor estimates of the conditional quantiles  $Q_Y(p|X=\bo{x})$. As an example, when $Y=0$ the term $\log(Z-p)=\log(U-p)$ could go from $-\infty$ to $0$, simply due to the added noise.  In addition, quantile regression can suffer from the problem of crossing quantile curves, which is usually seen in sparse regions of the covariate space. This happens due to the fact that the conditional quantile curve for a given $X=\bo{x}$ will not be a monotonically increasing function of $p$. % Specifically, we may get that the $50\%$ quantile of the conditional distribution given the covariates is smaller than the $25\%$ quantile, which is, of course, implausible.

 Another approach would be to view the counts as ordinal variables with fixed thresholds and then model the new latent variable by an infinite mixture of normal densities \cite{KW2012}. Instead of using the aforementioned methods, we propose an adaptive Dirichlet process mixture approach which estimates the conditional density of the data. The approach is based on an adaptive Dirichlet Process mixture (DPM) of COM-Poisson regression models.
\section{COM-Poisson distribution}
%Our alternative approach is based on mixtures of COM-Poisson distributions. The section gives a brief summary of key properties of this distribution.
 The COM-Poisson distribution \cite{CM1962,SMKBB2005} is a two-parameter generalisation of the Poisson distribution that allows for different levels of dispersion. The probability mass function of the COM-Poisson($\lambda$, $\nu$) distribution is
 \begin{equation}P(Y=y|\lambda, \nu)=\frac{\lambda^y}{(y!)^\nu}\frac{1}{Z(\lambda, \nu)} \ \  y=0, 1, 2, \ldots\end{equation}
where $Z(\lambda, \nu)=\displaystyle \sum_{j=0}^{\infty}\frac{\lambda^j}{(j!)^\nu}$ and $\lambda > 0$ and $\nu \ge 0$, where the normalisation constant does not have a closed form and has to be approximated numerically. The extra parameter $\nu$ allows the distribution to model under- ($\nu>1$) or over-dispersed ($\nu<1$) data, having the Poisson distribution as a special case ($\nu=1$). %The normalisation constant $Z(\lambda, \nu)=\displaystyle \sum_{j=0}^{\infty}\frac{\lambda^j}{(j!)^\nu}$ does not have a closed form and has to be approximated. % An upper bound on $Z(\lambda, \nu)$ and an asymptotic approximation that is accurate for $\lambda>10^{\nu}$ can be used, see \cite{MSKBB2003}. Arbitrarily precise bounds can also be found based on piecewise geometric approximations to the series. Methods for estimating the parameters of the COM-Poisson and showing its flexibility in fitting count data, compared to other distributions, can be seen in \cite{SMKBB2005}.
% The mean and variance of the COM-Poisson are approximated by  
% \begin{equation}
% \mathbb{E}[Y]\approx \lambda^{\frac{1}{\nu}}, \ \ \  \mathbb{V}[Y]\approx \frac{\lambda^{\frac{1}{\nu}}}{\nu}. \label{formulation}
% \end{equation}

The above formulation of the COM-Poisson does not have a clear centering parameter since the parameter $\lambda$ is close to the mean only when $\nu$ takes values close to $1$, which makes it difficult to interpret for under- or over-dispersed data. Substituting the parameter $\lambda$ with $\mu=\lambda^{\frac{1}{\nu}}$, where $\floor{\mu}$ is the mode of the distribution
\begin{equation}
\mathbb{E}[Y]\approx \mu, \ \ \  \mathbb{V}[Y]\approx \frac{\mu}{\nu} \label{reform}
\end{equation}
and the new probability mass function is 
  \begin{equation}P(Y=y|\mu, \nu)=\left(\frac{\mu^y}{y!}\right)^\nu\frac{1}{Z(\mu, \nu)} \  \ y=0, 1, 2, \ldots \label{reformpmf}\end{equation}
where $Z(\mu, \nu)=\displaystyle \sum_{j=0}^{\infty}\left(\frac{\mu^j}{j!}\right)^\nu$.
%with  $Z(\mu, \nu)=\displaystyle \sum_{j=0}^{\infty}\left(\frac{\mu^j}{j!}\right)^\nu$. %For more information on the COM-Poisson distribution see .
%\subsection{Approximating over/under- dispersed distributions}
\subsection{Mixtures of COM-Poisson distributions}
% Unlike the Poisson distribution, the COM-Poisson is flexible enough to approximate distributions with any kind of dispersion. In contrast to this, a mixture of Poisson distributions is sufficient only in the case of
% %  \begin{figure}
% % \centering
% %  \includegraphics[height=11cm]{approximation.pdf}
% %  \caption{Approximating a Binomial}\label{fig:approxbinomforhighp}
% %  \end{figure}
%   over-dispersion. A Poisson (or a mixture of Poisson distributions) approximation is not able to approximate an under-dispersed distribution.
  % This can be seen in Figure \ref{fig:approxbinomforhighp} where we approximate a binomial (with sample size $n=60$ and probability parameter $p=0.99$) with a Poisson, a COM-Poisson and a mixture of COM-Poisson distributions. The COM-Poisson and Poisson approximations plotted are the ones that minimize the Kullback-Leibler divergence. The COM-Poisson distribution is much close to the target distribution than the Poisson distribution.
The COM-Poisson is flexible enough to approximate distributions with any kind of dispersion in contrast to a Poisson or a mixture of Poisson distributions which can only deal with overdispersion.

The two parameters of the COM-Poisson distribution allow it to have arbitrary (positive) mean and variance; one can obtain a point mass by letting the variance parameter $\nu$ tend to infinity. Thus one can show that mixtures of COM-Poisson distributions can provide an arbitrarily precise approximation to any discrete distribution with support $\mathbb{N}_0$, which is why COM-Poisson distributions are used by our method. All other generalisations of the Poisson distribution we are aware of do not have this property. 

\subsection{COM-Poisson regression}
A regression model can be defined based on \eqref{reformpmf}, in which both the mean and the variance parameter are modelled as a function of covariates:
 \begin{align}
% P(Y_i=y_i|\mu_i, \nu_i)=&\left(\frac{\mu_i^{y_i}}{y_i!}\right)^{\nu_i}\frac{1}{Z(\mu_i, \nu_i)}\\
% Z(\mu_i, \nu_i)&=\sum_{j=1}^{\infty}\left(\frac{\mu_i^j}{j!}\right)^{\nu_i}\\
 \log{\bo{\mu_i}}&= \bo{x_i}^\intercal\bo{\beta}\\
 \log{\bo{\nu_i}}&= \bo{x_i}^\intercal\bo{c}
  \end{align}
% \begin{align}
% P(Y_i=y_i|\mu_i, \nu_i)=&\left(\frac{\mu_i^{y_i}}{y_i!}\right)^{\nu_i}\frac{1}{Z(\mu_i, \nu_i)}, &\log{\bo{\mu}}&= X\bo{\beta}, \\
% Z(\mu_i, \nu_i)=&\sum_{j=1}^{\infty}\left(\frac{\mu_i^j}{j!}\right)^{\nu_i},&\log{\bo{\nu}}&= X\bo{c}
% %\log{\bo{\nu}}&= X\bo{c}
%  \end{align}
where $Y$ is the response variable being modelled, and $\bo{\beta}, \bo{c}$ are the regression coefficients for the centering link function and the shape link function respectively. The parameters in this formulation have a direct link to either the mean or the variance, providing insight into the behaviour of the response variable. Notably,
% \begin{align}
% \mathbb{E}[Y_i]&\approx \exp(\bo{x}_i'\bo{\beta}),\\
% \mathbb{V}[Y]&\approx \frac{\exp(\bo{x}_i'\bo{\beta})}{\exp(\bo{x}_i'\bo{c})} =\exp(\bo{x}_i'(\bo{\beta}-\bo{c})).
% \end{align}
\begin{equation}
\mathbb{E}[Y_i]\approx \exp(\bo{x}_i'\bo{\beta}), \ \ \ \mathbb{V}[Y]\approx \frac{\exp(\bo{x}_i'\bo{\beta})}{\exp(\bo{x}_i'\bo{c})} =\exp(\bo{x}_i'(\bo{\beta}-\bo{c})).
\end{equation}
The calculation of the normalisation constant of the COM-Poisson distribution is the computationally most expensive part of the proposed regression model. It can be seen, in the next subsection, that this calculation is redundant.

\subsection{Exchange algorithm}
Any probability density function $p(y|\theta)$ can be written as
\begin{equation}p(y|\theta)=\frac{q_\theta(y)}{Z(\theta)}\end{equation} where $q_\theta(y)$ is the unnormalised density and the normalising constant $Z(\theta)=\int p(y,\theta) \, \mathrm{d}y $  is unknown. In this case the acceptance ratio of the Metropolis-Hastings algorithm is
\begin{align}
\alpha&=\min \left(1,\frac{q_{\theta^*}(y) \pi(\theta^*) Z(\theta)h(\theta|\theta^*)}{q_\theta(y) \pi(\theta) Z(\theta^*)h(\theta ^*|\theta)}\right)\label{defaultmh}
\end{align}
where $\pi (\theta)$ is the prior distribution of $\theta$. The acceptance ratio in \eqref{defaultmh} involves computing unknown normalising constants. Introducing auxiliary variables $\theta^*, y^*$ and sampling from an augmented distribution
\begin{equation}\pi (\theta ^*,y^*,\theta|y) \propto p(y|\theta)\pi (\theta) p(y^*|\theta ^*) h(\theta ^*|\theta) \end{equation}
results in
\begin{align}
\alpha&=\min \left(1,\frac{p(y|\theta^*)\pi (\theta^*) p(y^*|\theta) h(\theta|\theta^*)}{p(y|\theta)\pi (\theta) p(y^*|\theta ^*) h(\theta ^*|\theta)}\right)\\
&=\min \left(1,\frac{q_\theta(y^*)\pi (\theta^*) h(\theta|\theta^*)q_{\theta^*}(y) Z(\theta) Z(\theta^*) }{q_\theta(y) \pi (\theta) h(\theta ^*|\theta) q_{\theta^*}(y^*) Z(\theta^*)  Z(\theta) }\right)\\
&=\min \left(1,\frac{q_\theta(y^*)\pi (\theta^*) q_{\theta^*}(y) }{q_\theta(y) \pi (\theta) q_{\theta^*}(y^*) }\right) \label{exchangemh}
\end{align}
where the normalising constants cancel out and $h()$ is a symmetric distribution \cite{MGM06,MPRB2006}. In order to be able to use this algorithm one has to be able to sample from from the unnormalised density which in the case of the COM-Poisson distribution can be done efficiently using rejection sampling. 

Updating the parameter $\mu$ of the COM-Poisson we have $\theta=(\mu,\nu)$ and $\theta^*=(\mu^*,\nu)$
where $\mu^*$ follows a Normal distribution centered at $\mu$ and 
\begin{align}
q_\theta(y^*)&=\left(\frac{\mu_i^{y_i^*}}{y_i^*!}\right)^{\nu_i} & q_{\theta^*}(y)&=\left(\frac{(\mu_i^*)^{y_i}}{y_i!}\right)^{\nu_i}\\
q_\theta(y)&=\left(\frac{\mu_i^{y_i}}{y_i!}\right)^{\nu_i}&q_{\theta^*}(y^*)&=\left(\frac{(\mu_i^*)^{y_i^*}}{y_i^*!}\right)^{\nu_i}
\end{align}
Likewise for updating the parameter $\nu$.
%The idea behind the exchange algorithm is that you propose a parameter ($\theta^*$) and at the same time sample data ($y^*$) from the COM-Poisson distribution having this parameter. If the proposed parameter ($\theta^*$) ``explains'' the data ($y$) and the proposed data ($y^*$) better than the current parameter ($\theta$); then you accept the proposed parameter ($\theta^*$).

\section{Bayesian density regression}
Density regression is similar to quantile regression in that it allows flexible modelling of the response variable $Y$ given the covariates $\boldsymbol{x}=(x_1, \ldots, x_p)'$. Features (mean, quantiles, spread) of the conditional distribution of the response variable, vary with  $\boldsymbol{x}$, so, depending on the predictor values, features of the conditional distribution can change in a different way than the population mean. The difference between density regression and quantile regression is that density regression models the probability density function or probability mass function rather than directly modelling the quantiles. %Bayesian methods for density regression are considered in \cite{DPP2007} where the conditional distribution of the response variable is expressed as a mixture of regression models where the mixing weights vary with covariates

% Bayesian methods for density regression are considered in \cite{DPP2007} where the conditional distribution of the response variable is expressed as a mixture of regression models where the mixing weights vary with covariates. The prior for the uncountable collection of the mixture distributions is a weighted mixture of Dirichlet process priors where the weights are dependent on the distance between covariate values. Covariates that are ``close'' to each other tend to have higher weight. This procedure results in a generalisation of the  P\'{o}lya urn scheme that does not rely on the exchangeability of the subjects.
\subsection{Bayesian density regression for count data}
This paper focuses on the following mixture of regression models:
\begin{equation}f(y_i|\boldsymbol{x}_i)=\int f(y_i|\boldsymbol{x}_i,\phi_i)G_{\boldsymbol{x}_i}(\mathrm{d}\phi_i)\end{equation}
where \begin{equation}f(y_i|\boldsymbol{x}_i,\phi_i)=\text{COM-P}(y_i;\exp(\boldsymbol{x}_i'\boldsymbol{b}_i),\exp(\boldsymbol{x}_i'\boldsymbol{c}_i))\end{equation}
%where
%\begin{equation}f(y_i|\boldsymbol{x}_i,\phi_i)=\text{COM-Poisson}(y_i;\exp(\boldsymbol{x}_i'\boldsymbol{b}_i),\exp(\boldsymbol{x}_i'\boldsymbol{c}_i)),\end{equation}
the conditional density of the response variable given the covariates is expressed as a mixture of COM-Poisson regression models  with $\phi_i=(\boldsymbol{b}_i,\boldsymbol{c}_i)$ and $G_{\boldsymbol{x}_i}$ is an unknown mixture distribution that changes according to the location of $\boldsymbol{x}_i$. 


\subsection{MCMC algorithm}
Let $\boldsymbol{\theta}=(\theta_1, \ldots, \theta_k)'$ denote the $k\le n$ distinct values of $\phi$ and let $\boldsymbol{S}=(S_1, \ldots, S_n)'$ be a vector of indicators denoting the global configuration of subjects to distinct values $\boldsymbol{\theta}$, with $S_i=h$ indexing the location of the $i$th subject within the $\boldsymbol{\theta}$. In addition, let $\boldsymbol{C}=(C_1, \ldots, C_k)'$  with $C_h=j$ denoting that $\theta_h$ is an atom from the basis distribution, $G_{\boldsymbol{x}_j}^*$. Hence $C_{S_i}=Z_i=j$ denotes that subject $i$ is drawn from the $j$th basis distribution.

Excluding the $i$th subject, $\boldsymbol{\theta}^{(i)}=\boldsymbol{\theta}\backslash \{\phi_i\}$ denotes the $k^{(i)}$ distinct values of $\boldsymbol{\phi}^{(i)}=\boldsymbol{\phi}\backslash \{\phi_i\}$, $\boldsymbol{S}^{(i)}$ denotes the configuration of subjects $\{1, \ldots, n\}\backslash \{i\}$ to these values and  % $\boldsymbol{C}^{(i)}=(C_1^{(i)}, \ldots, C_{k^{(i)}}^{(i)})'$ 
$\boldsymbol{C}^{(i)}$ indexes the DP component numbers for the elements of  $\boldsymbol{\theta}^{(i)}$.

% Conditioning on $\boldsymbol{Z}^{(i)}$ but marginalizing over $Z_i$, the following conditional prior for $\phi_i$ is obtained:
% \begin{equation}
% (\phi_i|\boldsymbol{Z}^{(i)}, \phi_i,\boldsymbol{X},a) \sim \sum_{i=1}^n \frac{ab_{ij}}{a + \sum_{l\ne i}\boldsymbol{1}(Z_l=j)}G_0  \\+ \sum_{m\ne i}\Big\{ \sum_{j=1}^n \frac{b_{ij}\boldsymbol{1}(Z_m=j)}{a+ \sum_{l\ne i}\boldsymbol{1}(Z_l=j)} \Big\} \delta_{\phi_m}. \end{equation}
% Grouping the subjects in the same cluster, we obtain the expression
% \begin{equation}(\phi_i|\boldsymbol{S}^{(i)},\boldsymbol{C}^{(i)},\boldsymbol{\theta}^{(i)} ,\boldsymbol{X},a) \sim w_{i0}G_0 + \sum_{h=1}^{k^{(i)}} w_{ih} \delta_{\theta_h^{(i)}},\end{equation}
% where the probability weights on the various components are defined as 
% % \begin{align}
% % w_{i0}&=\sum_{j=1}^n \frac{ab_{ij}}{a + \sum_{l\ne i}\boldsymbol{1}(C_{S_l^{(i)}}^{(i)}=j)} \\
% % w_{ih}&=\frac{b_{i,C_h^{(i)}}\sum_{m\ne i}\boldsymbol{1}(S_m^{(i)}=h)}{a+ \sum_{l\ne i}\boldsymbol{1}((C_{S_l^{(i)}}^{(i)}=C_h)},\  h=1, \ldots , k^{(i)}.
% % \end{align}
% \begin{equation}
% w_{i0}=\sum_{j=1}^n \frac{ab_{ij}}{a + \sum_{l\ne i}\boldsymbol{1}(C_{S_l^{(i)}}^{(i)}=j)} \ \ \,
% w_{ih}=\frac{b_{i,C_h^{(i)}}\sum_{m\ne i}\boldsymbol{1}(S_m^{(i)}=h)}{a+ \sum_{l\ne i}\boldsymbol{1}(C_{S_l^{(i)}}^{(i)}=C_h)},\  h=1, \ldots , k^{(i)}.
% \end{equation}

% Updating the prior with the likelihood for the data $\boldsymbol{y}$, we obtain the conditional posterior:
% \begin{equation}(\phi_i|\boldsymbol{S}^{(i)},\boldsymbol{C}^{(i)},\boldsymbol{\theta}^{(i)} ,\boldsymbol{X},a) \sim q_{i0}G_{i,0} + \sum_{h=1}^{k^{(i)}} q_{ih}\delta_{\theta_h^{(i)}},\end{equation}
% where $G_{i,0}(\phi)$ is the posterior obtained by updating the prior $G_0(\phi)$ and the likelihood $f(y_i|\boldsymbol{x}_i,\phi)$:
Grouping the subjects in the same cluster and updating the prior with the likelihood for the data  $\boldsymbol{y}$, we obtain the conditional posterior
\begin{equation}(\phi_i|\boldsymbol{S}^{(i)},\boldsymbol{C}^{(i)},\boldsymbol{\theta}^{(i)} ,\boldsymbol{X},a) \sim q_{i0}G_{i,0} + \sum_{h=1}^{k^{(i)}} q_{ih}\delta_{\theta_h^{(i)}},\end{equation}
where $G_{i,0}(\phi)$ is the posterior obtained by updating the prior $G_0(\phi)$ and the likelihood $f(y_i|\boldsymbol{x}_i,\phi)$:
\begin{align}
  G_{i,0}(\phi)&= \frac{G_0(\phi)f(y_i|\boldsymbol{x}_i,\phi)}{h_i(y_i|\boldsymbol{x}_i)}, & &\\
  q_{i0}&= cw_{i0}h_i(y_i|\boldsymbol{x}_i),& q_{ih}&=cw_{ih}f(y_i|\boldsymbol{x}_i,\theta_h),\\
  w_{i0}&=\sum_{j=1}^n \frac{ab_{ij}}{a + \sum_{l\ne i}\boldsymbol{1}(C_{S_l^{(i)}}^{(i)}=j)}, &w_{ih}&=\frac{b_{i,C_h^{(i)}}\sum_{m\ne i}\boldsymbol{1}(S_m^{(i)}=h)}{a+ \sum_{l\ne i}\boldsymbol{1}(C_{S_l^{(i)}}^{(i)}=C_h)},
 % &[x \mapsto s](t_1 \; t_2)&&= ([x \mapsto s]t_1)\;([x \mapsto s]t_2) &&
\end{align}

% \begin{alignnat}{2}
% &G_{i,0}(\phi) &=  \frac{G_0(\phi)f(y_i|\boldsymbol{x}_i,\phi)}{h_i(y_i|\boldsymbol{x}_i)}&&,
% &q_{i0}&=cw_{i0}h_i(y_i|\boldsymbol{x}_i), q_{ih}&=cw_{ih}f(y_i|\boldsymbol{x}_i,\theta_h)&,
% &w_{i0}&=\sum_{j=1}^n \frac{ab_{ij}}{a + \sum_{l\ne i}\boldsymbol{1}(C_{S_l^{(i)}}^{(i)}=j)}, w_{ih}&=\frac{b_{i,C_h^{(i)}}\sum_{m\ne i}\boldsymbol{1}(S_m^{(i)}=h)}{a+ \sum_{l\ne i}\boldsymbol{1}(C_{S_l^{(i)}}^{(i)}=C_h)},\  h=1, \ldots , k^{(i)}&.
% \end{alignnat}
% \begin{equation}
% G_{i,0}(\phi) =  \frac{G_0(\phi)f(y_i|\boldsymbol{x}_i,\phi)}{h_i(y_i|\boldsymbol{x}_i)},
% \end{equation}
% \begin{equation}\hspace{-7cm}q_{i0}=cw_{i0}h_i(y_i|\boldsymbol{x}_i), q_{ih}=cw_{ih}f(y_i|\boldsymbol{x}_i,\theta_h),\end{equation}
% \begin{equation}
% w_{i0}=\sum_{j=1}^n \frac{ab_{ij}}{a + \sum_{l\ne i}\boldsymbol{1}(C_{S_l^{(i)}}^{(i)}=j)} \ \ \,
% w_{ih}=\frac{b_{i,C_h^{(i)}}\sum_{m\ne i}\boldsymbol{1}(S_m^{(i)}=h)}{a+ \sum_{l\ne i}\boldsymbol{1}(C_{S_l^{(i)}}^{(i)}=C_h)},\  h=1, \ldots , k^{(i)}.
% \end{equation}
%\begin{equation}(\phi_i|\boldsymbol{S}^{(i)},\boldsymbol{C}^{(i)},\boldsymbol{\theta}^{(i)} ,\boldsymbol{X},a) \sim w_{i0}G_0 + \sum_{h=1}^{k^{(i)}} w_{ih} \delta_{\theta_h^{(i)}},\end{equation}
%where the probability weights on the various components are defined as 
% \begin{align}
% w_{i0}&=\sum_{j=1}^n \frac{ab_{ij}}{a + \sum_{l\ne i}\boldsymbol{1}(C_{S_l^{(i)}}^{(i)}=j)} \\
% w_{ih}&=\frac{b_{i,C_h^{(i)}}\sum_{m\ne i}\boldsymbol{1}(S_m^{(i)}=h)}{a+ \sum_{l\ne i}\boldsymbol{1}((C_{S_l^{(i)}}^{(i)}=C_h)},\  h=1, \ldots , k^{(i)}.
% \end{align}
%Updating the prior with the likelihood for the data $\boldsymbol{y}$, we obtain the conditional posterior:
where $b_{ij}$ are weights that depend on the distance between subjects' predictor values, $c$ is a normalising constant and $h=1, \ldots , k^{(i)}$. Since there is no closed form expression for the posterior distribution, approximation of the probability $q_{i0}=cw_{i0}h_i(y_i|\boldsymbol{x}_i)$ is difficult. 

We overcome this problem by bridging: i) an MCMC algorithm for sampling from the posterior distribution of a Dirichlet process model, with a non-conjugate prior, found in \cite{N2000}; ii) the MCMC algorithm found in \cite{DPP2007}; and iii) a variation of the MCMC exchange algorithm. %We use a variation of algorithm $5$ from \cite{N2000}.




% \begin{align}
% G_{i,0}(\phi) &= \frac{G_0(\phi)f(y_i|\boldsymbol{x}_i,\phi)}{ \int f(y_i|\boldsymbol{x}_i,\phi)\mathrm{d} G_0(\phi)}\\
% &=  \frac{G_0(\phi)f(y_i|\boldsymbol{x}_i,\phi)}{h_i(y_i|\boldsymbol{x}_i)},
% \end{align}



% \cite{N2000} introduces auxilliary parameters that exist only temporarily and define a valid MCMC sampler. The idea is that if  $$h_i(y_i|\boldsymbol{x}_i)= \int f(y_i|\boldsymbol{x}_i,\phi)\mathrm{d} G_0(\phi)$$ is approximated by an average of $m$ values sampled from $G_0(\phi)$, we could also approximate a draw from the posterior distribution by sampling from the $m$ values with probabilities proportional to their likelihoods. The probability of allocating an observation to a new cluster, in this case, is $$q_{i0j}=cw_{i0}\frac{f(y_i|\boldsymbol{x}_i,\boldsymbol{\theta_j})}{m} \ j=1,2, \ldots, m$$ while the probability of allocating an observation to an already existing cluster \begin{equation}q_{ih}=cw_{ih}f(y_i|\boldsymbol{x}_i,\boldsymbol{\theta_h}) \ h=m+1,m+2, \ldots,m+k\end{equation}  where $\boldsymbol{\theta_h}=(b_h,c_h)$ represent the regression coefficients for the $h$ cluster.




 The MCMC algorithm alternates between the following steps:
\begin{description}
\item[Step 1:] Update $S_i$ for $i=1, \ldots, n,$ by proposing, from the conditional prior,  a move to a new cluster or an already existing cluster with probabilities proportional to $w_{i0}$ and $w_{ih}$ for $h=1, \ldots,k^{(i)}.$ 
\begin{enumerate}
\item[a)] If the proposed move is to go to a new cluster we draw parameters ($\mu_0,\nu_o$) for that cluster from $G_0$ and at the same time sample an observation $y^*$ from the COM-Poisson($\mu_0,\nu_0$). The acceptance ratio of the Metropolis-Hastings algorithm is 
\begin{equation}\min \left(1,\frac{q_\theta(y^*)q_{\theta^*}(y) }{q_\theta(y) q_{\theta^*}(y^*) }\right) \label{exchangemh2}
\end{equation}
%It is slightly different to \ref{exchangemh}, because in this case we choose the conditional prior of $S_i$ as the proposal distribution. % As a result the factors $\pi (\theta^*), \pi (\theta)$ and $h(\theta ^*|\theta), h(\theta|\theta^*)$ cancel out.
If the proposal is accepted,  $C_{S_i} \sim$ multinomial $(\{1, \ldots, n\},\boldsymbol{b}_i)$.
\item[b)] If the proposed move is to an already existing cluster $h$, we sample an observation $y^*$ from the COM-Poisson($\mu_h,\nu_h$) and accept with the same probability as in \eqref{exchangemh2}. If the proposal is accepted $C_{S_i}=C_h$.
  % sampling from the multinomial conditional posterior $\textit{Pr}(S_i=h)=q_{ih}, h=1, \ldots,k^{(i)}.$ When $S_i=0$, this means that the observation will go to a new cluster ($S_i=k+1$) with the parameters that we have already drawn ( $\boldsymbol{\theta}_{S_{i}}$) for that cluster.
 \end{enumerate}


\item[Step 2:] Update the parameters $\theta_h$, for $h=1, \ldots, k$ by sampling from the conditional posterior distribution
\begin{equation}(\theta_h|\boldsymbol{S},\boldsymbol{C},\boldsymbol{\theta}^{(h)} ,k,\boldsymbol{y},\boldsymbol{X}) \sim \prod_{i:S_i=h} f(y_i|\boldsymbol{x}_i,\theta_h) \} G_0(\theta_h),\end{equation}
using the  Metropolis-Hasting algorithm with acceptance probability as in \eqref{exchangemh}.
\item[Step 3:] Update $C_h$, for $h=1, \ldots, k$, by sampling from the multinomial conditional with 
\begin{equation}(C_h|\boldsymbol{S},\boldsymbol{C}^{(h)},\boldsymbol{\theta} ,k,\boldsymbol{y},\boldsymbol{X}) \sim \frac{\prod_{i:S_i=h} b_{ij}}{\sum_{l=1}^n \prod_{i:S_i=h} b_{il}}, \ j=1, \ldots, n\end{equation} and location weights $\gamma_j$ for $j=1, 2, \ldots, n$, using an approach used in \cite{DS2005}. %For more information on choosing the proposed weight function $\boldsymbol{b}(\boldsymbol{x})$ and the location weights $\gamma$ see \cite{DPP2007}.
%$ \boldsymbol{b}(\boldsymbol{x})=\gamma_j\textit{K}(\boldsymbol{x},\boldsymbol{x}_j)/\sum_{l=1}^n \gamma_l\textit{K}(\boldsymbol{x},\boldsymbol{x}_j) \ \ j=1, \ldots, n, \ \ \forall \boldsymbol{x}$
%and $\gamma=(\gamma_1, \ldots,\gamma_n)'$ represent location weights and K is a kernel function, such as
%$\textit{K}(\boldsymbol{x},\boldsymbol{x}')=\exp\{-\psi \|\boldsymbol{x}-\boldsymbol{x}'\|^2 \}.$
%\vspace{-1cm}
%\item[Step 4:] Update the location weights $\gamma_j$ for $j=1, 2, \ldots, n$, by using a data augmentation approach used in \cite{DS2005}.
\end{description}



%Say something
% \subsection{Simulation}
% In this section we consider two simulated data sets to compare the proposed discrete Bayesian density regression method to the ``jittering'' method in \cite{MS2002}. The first data set is simulated from 
% \begin{equation}P(Y_i=y_i|x_i)= \text{Binomial}(10,0.3x_i)\end{equation} where $x_i \sim \text{Unif}(0,1)$.

% The second simulated data set is more complicated, as it is a mixture of three distributions with different levels of dispersion,
% \begin{equation}
% P(Y_i=y_i|x_i)= 0.4\text{Pois}(\exp(1+x_i)) +\\ 0.2\text{Binomial}(10,1-x_i) +0.4 \text{Geom}(0.2),\end{equation}
%  where $x \sim \text{Unif}(0,1)$. We compare our method to the ``jittering'' method (both using a linear model and a non-linear spline model). For our method we use the above MCMC algorithm based on $5000$ iterations with the first $3000$ used as a burn-in period.

% Table~\ref{difference} shows the absolute mean errors obtained using both methods. If $q_p$ is the true conditional quantile when $x=p$ and $\hat{q_p}$ is the estimated conditional quantile, the mean absolute error is defined as
% \begin{equation}\mathbb{E}[|q_p - \hat{q_p}|]\end{equation}

% The discrete Bayesian density regression (BDR) estimates outperform the ``jittering'' method  being used. % Figure~\ref{fig:quantiles} shows the true and estimated conditional quantiles for the second simulation (for $q=0.1, 0.5, 0.95$) across the covariate space, where the

% % \begin{figure}
% % \centering
% % \includegraphics[height=11cm]{quantiles.pdf}
% % \caption{$10\%$, $50\%$ and $95\%$ quantiles obtained using our method (BDR) and the ``jittering'' approach for the simulated mixture data.}\label{fig:quantiles}
% % \end{figure}
% %  estimated conditional quantiles are drawn with a dashed line.

%  \begin{table}[h]
% \centering
% %\begin{center}
% %\tabcolsep=0.04cm
% \scalebox{0.9}{
% \begin{tabular}{llccc}
% %\hline
% %\hline\noalign{\smallskip}
% Data&Method& \multicolumn{3}{c}{\centering Number of Observations} \\ \hline\noalign{\smallskip}
%  & &\multicolumn{1}{c}{20} & \multicolumn{1}{c}{100} & \multicolumn{1}{c}{500}\\ 
% %&\cline{2-4} & &\multicolumn{1}{c}{20} & \multicolumn{1}{c}{100} & \multicolumn{1}{c}{500}\\ 
% \hline\noalign{\smallskip}
% &Dens. Regression & 0.5576 &   \bftab{0.2820} & \bftab{0.2421}\\
% Binomial&Jittering (linear) & \bftab{0.5256}  & 0.8461  & 0.4765\\
% &Jittering (splines) &0.7820 &0.5128 & 0.3020\\
% \hline\noalign{\smallskip}
% &Dens. Regression & \bftab{0.7435} &   \bftab{0.5833} & \bftab{0.3589}\\
% Mixture&Jittering (linear) & 1.1923 & 0.6666  & 0.4294\\
% &Jittering (splines) &1.9487 & 0.8269 & 0.3910\\
% \noalign{\smallskip}\hline\noalign{\smallskip}
% \end{tabular}}
% \caption{Mean absolute error obtained using the different density/quantile regression methods for the two simulated datasets.} \label{difference}
% %\end{center}
% \end{table}

% In almost all cases the ``jittering'' method leads to crossing quantiles, except when the number of observations was $n=500$.



\section{Simulations and application}
We consider two simulated data sets to compare the proposed discrete Bayesian density regression method to the ``jittering'' method. These are
\begin{align}
Y_i|X_i=x_i&\sim \text{Binomial}(10,0.3x_i)\\
Y_i|X_i=x_i&\sim 0.4\text{Pois}(\exp(1+x_i)) + 0.2\text{Binomial}(10,1-x_i) +0.4 \text{Geom}(0.2)
\end{align} 
 where $x_i \sim \text{Unif}(0,1)$. Table~\eqref{difference} shows the absolute mean errors obtained using both methods. If $q_p$ is the true conditional quantile when $x=p$ and $\hat{q_p}$ is the estimated conditional quantile, the mean absolute error is defined as $\mathbb{E}[|q_p - \hat{q_p}|]$. The discrete Bayesian density regression (BDR) estimates outperform the ``jittering'' method and in almost all cases the ``jittering'' method leads to crossing quantiles (except when $n=500$).


% The second simulated data set is more complicated, as it is a mixture of three distributions with different levels of dispersion,
% \begin{multline*}
% P(Y_i=y_i|x_i)= 0.4\text{Pois}(\exp(1+x_i)) +\\ 0.2\text{Binomial}(10,1-x_i) +0.4 \text{Geom}(0.2),\end{multline*}
%  where $x \sim \text{Unif}(0,1)$
% We compare our method to the ``jittering'' method (both using a linear model and a non-linear spline model). For our method we use the above MCMC algorithm based on $5000$ iterations with the first $3000$ used as a burn-in period.
% \begin{tablehere}
% \centering
% \tabcolsep=0.05cm
% \scalebox{0.8}{
% \begin{tabular}{llccc}
% %\hline
% %Data&Method& \multicolumn{3}{c}{\centering Number of Observations} \\ \hline\noalign{\smallskip}
% %&\cline{2-4} & &\multicolumn{1}{c}{20} & \multicolumn{1}{c}{100} & \multicolumn{1}{c}{500}\\ \hline\noalign{\smallskip}
% Data&Method& \multicolumn{3}{c}{\centering Number of Observations} \\ \hline\noalign{\smallskip}
%  & &\multicolumn{1}{c}{20} & \multicolumn{1}{c}{100} & \multicolumn{1}{c}{500}\\ \hline\noalign{\smallskip}  
% &Dens. Regression & 0.5576 &   \bftab{0.2820} & \bftab{0.2421}\\
% Binomial&Jittering (linear) & \bftab{0.5256}  & 0.8461  & 0.4765\\
% &Jittering (splines) &0.7820 &0.5128 & 0.3020\\
% \hline\noalign{\smallskip}
% &Dens. Regression & \bftab{0.7435} &   \bftab{0.5833} & \bftab{0.3589}\\
% Mixture&Jittering (linear) & 1.1923 & 0.6666  & 0.4294\\
% &Jittering (splines) &1.9487 & 0.8269 & 0.3910\\
% \noalign{\smallskip}\hline\noalign{\smallskip}
% %\hline
% \end{tabular}}
% \caption{Mean absolute error obtained using the different density/quantile regression methods. % for the two simulated datasets.
% } \label{difference}
% \end{tablehere}
\begin{tablehere}
\centering
\tabcolsep=0.05cm
\scalebox{0.9}{
\begin{tabular}{lcccccc}
%\hline
%Data&Method& \multicolumn{3}{c}{\centering Number of Observations} \\ \hline\noalign{\smallskip}
%&\cline{2-4} & &\multicolumn{1}{c}{20} & \multicolumn{1}{c}{100} & \multicolumn{1}{c}{500}\\ \hline\noalign{\smallskip}
Method& \multicolumn{6}{c}{\centering Number of Observations} \\ \hline\noalign{\smallskip}
&\multicolumn{3}{c}{\centering Binomial}&\multicolumn{3}{c}{\centering Mixture}  \\
  &\multicolumn{1}{c}{20} & \multicolumn{1}{c}{100} & \multicolumn{1}{c}{500}&\multicolumn{1}{c}{20} & \multicolumn{1}{c}{100} & \multicolumn{1}{c}{500}\\ \hline\noalign{\smallskip}  
Density Regression & 0.5576 &   \bftab{0.2820} & \bftab{0.2421} & \bftab{0.7435} &   \bftab{0.5833} & \bftab{0.3589}\\
Jittering (linear) & \bftab{0.5256}  & 0.8461  & 0.4765 & 1.1923 & 0.6666  & 0.4294\\
Jittering (splines) &0.7820 &0.5128 & 0.3020 &1.9487 & 0.8269 & 0.3910\\
%\hline\noalign{\smallskip}
%&Dens. Regression & \bftab{0.7435} &   \bftab{0.5833} & \bftab{0.3589}\\
%Mixture&Jittering (linear) & 1.1923 & 0.6666  & 0.4294\\
%&Jittering (splines) &1.9487 & 0.8269 & 0.3910\\
\noalign{\smallskip}\hline\noalign{\smallskip}
%\hline
\end{tabular}}
\caption{Mean absolute error obtained using the different density/quantile regression methods. % for the two simulated datasets.
} \label{difference}
\end{tablehere}



We apply the discrete density regression technique to data on housebreakings in Greater Glasgow (Scotland). The data consist of the number of housebreakings in each of the 127 intermediate geographies in Greater Glasgow in 2010.  We aim to relate the number of housebreakings to the deprivation score of the intermediate geography area, as measured by the Scottish Index of Multiple Deprivation (SIMD). The deprivation score is standardised by considering the difference of each intermediate geography's deprivation from the average deprivation in Greater Glasgow e.g. low values relate to affluent areas, large values to deprived areas. % The estimated quantiles  are shown as dashed lines .  The number of housebreakings is higher in affluent areas and lower in deprived areas.
The solid and dashed lines in figure~\ref{fig:glm2} show the quantiles (for $p=0.1, 0.5, 0.95$) obtained for the standard Poisson regression model and the COM-Poisson model respectively. The first model is not able to capture the overdispersion of the data, nor the skewness of the distribution. %This highlights the benefits of using a density-based or quantile-based regression  approach, compared to the classical (generalised) mean regression.

\begin{figure}
\centering
\includegraphics[height=7.2cm]{glm2.pdf}
%\vspace{-1.cm}
\caption{Estimated quantiles for housebreaking data, using discrete Bayesian density regression (dashed lines) and derived from a Poisson model.}\label{fig:glm2}
\end{figure}


 
 



% Table~\ref{difference} shows the absolute mean errors obtained using both methods. If $q_p$ is the true conditional quantile when $x=p$ and $\hat{q_p}$ is the estimated conditional quantile, the mean absolute error is defined as $\mathbb{E}[|q_p - \hat{q_p}|]$. The discrete Bayesian density regression (BDR) estimates outperform the ``jittering'' method  being used and in almost all cases the ``jittering'' method leads to crossing quantiles (except when $n=500$).


% In addition, we apply our discrete density regression technique to data on housebreakings in Greater Glasgow (Scotland). The data consist of the number of housebreakings in each of the 127 intermediate geographies in Greater Glasgow in 2010.  We aim to relate the number of housebreakings to the deprivation score of the intermediate geography area, as measured by the Scottish Index of Multiple Deprivation (SIMD). Before using it in the model, the deprivation score is standardised by considering the difference of each intermediate geography's deprivation from the average deprivation in Greater Glasgow. Low values of the normalised deprivation score relate to affluent aeras, large values to deprived areas.


%  The estimated quantiles (for $p=0.1, 0.5, 0.95$) are shown as dashed lines in figure~\ref{fig:glm2}. The number of housebreakings is higher in affluent areas and lower in deprived areas. The solid lines in Figure~\ref{fig:glm2} show the quantiles obtained when using a standard Poisson regression model. The standard Poisson regression model is neither able to capture the overdispersion of the data, nor the skewness of the distribution.
%  This highlights the benefits of using a density-based or quantile-based regression  approach, compared to the classical (generalised) mean regression.

\section{Conclusions and further research}
In this manuscript we have proposed a novel Bayesian density regression technique for discrete data which is based on mixing COM-Poisson distributions. The new method takes advantage of the exchange algorithm and updates the cluster allocations by drawing a new allocation for an auxiliary observation and then accepting or rejecting it. As a result the MCMC samples from the target distribution without the need to estimate the normalisation constant of the likelihood. The method overcomes the two main drawbacks of the ``jittering'' method for discrete quantile regression, namely that it does not require the addition of artificial additional noise and that it does not suffer from the problem of crossing quantiles. We have illustrated the method in a real world application as well as simulated examples in which our method compared favourably to the ``jittering'' method. Further research efforts will be devoted in improving the computational speed and efficiency of the MCMC algorithm to make it an even more attractive alternative to ``jittering''.

%We would like to thank the anonymous referee for his suggestions and comments.%  \cite{KW2012}, which models the mixture weights found in \eqref{Karabatsos} by an ordered-category probit regression for the infinite component labels, and }
 % {\color{red} and exploring different procedures for applying Bayesian density regression to count data. }
% (1) Why not perform a log transformation of the count dependent variable (Y) data,
% such as Y* = log(Y + 1), and then perform a usual continuous Bayesian density regression
% analysis of Y* ?  (which may already be computationally faster, anyways)

% (2) Another Bayesian density regression approach for a count dependent variable Y,
% which was suggested by Karabatsos & Walker (2012, Elec J Stat, Appendix),
% involves viewing Y as an ordinal variable with fixed thresholds -Inf,0,1,2,...
% and then modeling the latent variable Y* by a infinite mixture
% of normal densities, with covariate-dependent mixture weights.
% Then, Y* is a continuous function of the covariates, interpretable on a counts scale
% (by virtue of the chosen thresholds).
% Specifically, in each stage of the MCMC algorithm for this model, for each observation i=1,...,n,
% the full conditional posterior density of Y*_i
% is a normal density truncated to the threshold interval (Y_i - 1,Y_i) when Y_i>0,
% and otherwise is truncated into the interval (-Inf, 0) when Y_i = 0. 







\begin{thebibliography}{99}

% \bibitem{smithetal:2006} Doe, J., Smith, A. and Smith, J. (2006) \textit{Hand book on first useful conclusions on the well-known interesting research topic.} Interesting topics editions.

% \bibitem{smith:1999} Smith, J. (1999) \textit{A first introductory research on some interesting topics.} Local journal of interesting topics research, \textbf{1}, 1--34.

% \bibitem{smith:2005} Smith, J. (2005) \textit{A second and not so introductory research on the same interesting topics.} Local journal of interesting topics research, \textbf{6}, 2345--2358.

% \bibitem{smithdoe:1999} Smith, A. and Doe, J. (1999) \textit{Discussion on: A first introductory research on some interesting topics.} Local journal of interesting topics research, \textbf{1}, 35--41.

\bibitem{B1998} Buchinsky, Moshe (1998) \textit{Recent advances in quantile regression models: A practical guideline for empirical research.} The journal of human resources, \textbf{33}, 88--126.

%\bibitem{cade2003} Cade, Brian S. and Noon, Barry R. (2003) \textit{Frontiers in Ecology and the Environment} A gentle introduction to quantile regression for ecologists, \textbf{1}, 412--420.

%\bibitem{C1994} Chamberlain, Gary (1994) \textit{Advances in Econometrics} Quantile regression, censoring and the structure of wages.


\bibitem{CM1962} Conway, Richard W. and Maxwell, William L. (1962) \textit{A queuing model with state dependent service rate.} Journal of industrial engineering, \textbf{12}, 132--136.


\bibitem{DPP2007} Dunson, David B. and Pillai, Natesh and Park, Ju-Hyun. (2007) \textit{Bayesian density regression.} Journal of the royal statistical society: Series B, \textbf{69}, 163--183.


\bibitem{DS2005} Dunson, David B. and Stanford, Joseph B. (2005) \textit{Bayesian inferences on predictors of conception probabilities.} Biometrics, \textbf{61}, 126--133.


%\bibitem{GCJG2008} Guikema, Seth D. and Coffelt, Jeremy P. (2008) \textit{Risk analysis : an official publication of the Society for Risk Analysis} A flexible count data regression model for risk analysis, \textbf{28}, 213--223.

\bibitem{KW2012} Karabatsos, George and Walker, Stephen G. (2012) \textit{Adaptive-modal Bayesian nonparametric regression.} Electronic journal of statistics, 6, 2038--2068.


\bibitem{KB1978} Koenker, Roger and Bassett, Gilbert (1978) \textit{Regression quantiles.} Econometrica, \textbf{46}, 33--50.


%\bibitem{levin2001} Levin, Jesse (2006) \textit{For whom the reductions count: A quantile regression analysis of class size and peer effects on scholastic achievement} Empirical Economics, \textbf{1}, 221--246.


\bibitem{MS2002} Machado, Jos\`e Ant\'onio Ferreira and Santos Silva, Jo\~ao M.C.. (2005) \textit{Quantiles for counts.} Journal of the american statistical association, \textbf{100}, 1226--1237.

%\bibitem{MSKBB2003} Minka, Thomas P. and Shmueli, Galit and Kadane,  Joseph B. and Borle, Sharad and Boatwright, Peter (2003) \textit{Technical report} Computing with the {COM}-{P}oisson distribution.

\bibitem{MPRB2006} M{\o}ller, J. and Pettitt, A. N. and Reeves, R. and Berthelsen, K. K. (2006) \textit{An efficient Markov chain Monte Carlo method for distributions with intractable normalising constants.} Biometrika, 93\textbf{2}, 451--458.

\bibitem{MGM06} Murray, Ian and Ghahramani, Zoubin and MacKay, David J. C. (2006) \textit{MCMC for doubly-intractable distributions.} Proceedings of the 22nd Annual UAI Conference, 359--366.


\bibitem{N2000} Neal, Radford M. (2000) \textit{Markov chain sampling methods for {D}irichlet process mixture models.} Journal of computational and graphical statistics, \textbf{9}, 249--265.


%\bibitem{S2004} Schulze, Niels (2004) \textit{PhD thesis} Applied Quantile Regression: Microeconometric, Financial, and Environmental Analyses.


\bibitem{SMKBB2005} Shmueli, Galit and Minka, Thomas P. and Kadane, Joseph B. and Borle, Sharad and Boatwright, Peter (2008) \textit{A useful distribution for fitting discrete data: revival of the {C}onway-{M}axwell-{P}oisson distribution.} Journal of the royal statistical society: Series C, \textbf{54}, 127--142.


%\bibitem{wei2006} Wei, Ying  and He, Xuming (2006) \textit{Annals of Statistics} Conditional growth charts, \textbf{1}, 2069--2097.


%\bibitem{W2006} Winkelmann, Rainer (2006) \textit{Journal of Industrial Engineering} Reforming health care: Evidence from quantile regressions for counts, \textbf{25}, 131--145.





\end{thebibliography}


% \begingroup
% \let\clearpage\relax
% \bibliographystyle{abbrv}
%  \bibliography{compstat}{}
% \endgroup

% \section{Citations and lists}
% %Please use the \LaTeX\ automatism for your citations. For instance, see \cite{smith:1999} or see \cite{smith:2005} and \cite{smithdoe:1999} or see \cite{smithetal:2006}, \cite{smith:1999} and \cite{smithdoe:1999}. You can also use joint citations as in \cite{smith:1999,smith:2005} if you prefer.

% \section{Lists}

% Note that \verb|enumeration| environment uses numbers in first level and letters in the nested one.
% \begin{enumerate}
% \item External level, first item.
% \begin{enumerate}
% \item First item of the nested level.
% \item Second item of the nested level.
% \end{enumerate}
% \item External level, second item.
% \end{enumerate}

% For unnumbered list we recommend to use the \verb|itemize| environment. Note that by default external level uses dots while first nested level uses -.

% \begin{itemize}
% \item External level, first item.
% \begin{itemize}
% \item First item of the nested level.
% \item Second item of the nested level.
% \item[$\diamond$] Third item with specified mark
% \end{itemize}
% \item External level, second item.
% \end{itemize}

% \section{Equations}
% Use the standard \verb"equation" environment to typeset your equations,
% \begin{equation}
% a \times b = c\;,
% \end{equation}

% For multiline equations it is advisable to use the \verb|eqnarray| environment
% \begin{eqnarray}
% a \times b = c \nonumber\\
% \sum_{i=1}^k \prod_{j=1}^n A_i \times B_j \\
% \int_{[0,1]} x^p (1-x)^q \lambda(dx)
% \end{eqnarray}

% Please avoid inline equations (within text) whenever it is possible.

% See Section \ref{john.smith:theorems} for the way of handling Theorems, Definitions, Lemmas and so on.


% \section{Tables}
% Please insert tables by using the standard \verb"table" environment as in Table \ref{john.smith:table}. Add a caption to each table to explain it clearly. As in the case of figures (see Section \ref{john.smith:figures}), please add a \verb"center" environment inside the \verb"table" one to center all the elements. Please check the example for a way of adding extra vertical space between horizontal lines and elements in the table.

% \begin{table}[h]
% 	\centering
% 		\begin{tabular}{clcr}
% 		Column1 & Column2 & Column3 & Column4  \\ \hline\noalign{\smallskip}
% 		A & Left & 123 & Right justified\\
% 		A & Left & 222 & Right justified\\
% 		B & Left & 12--21 & Right justified\\
% 		B & Left & 1000 & Right justified\\
% 		\noalign{\smallskip}\hline\noalign{\smallskip}
% 		\end{tabular}
% 		\caption{Please write your table caption here.}
% 		\label{john.smith:table}       % Give a unique label
% \end{table}
% %


% \section{Figures}
% \label{john.smith:figures}
% Please, insert figures by using the standard \verb"figure" environment as in Figure \ref{john.smith:figure}. You can use both \verb"\includegraphics" or \verb"\epsfig" for inserting the graphics in the figure environment. Eps figures are preferable. Add a caption to each figure to explain it clearly.

% % \begin{figure}[h]
% % 	\centering
% % 		\includegraphics[scale=.45]{figure.eps}
% % 		\caption{Inside each figure environment add a center environment}
% % 		\label{john.smith:figure}      % Give a unique label as usual
% % \end{figure}

% \section{Theorems, Definitions, etc}
% \label{john.smith:theorems}

% You can insert Theorems, Lemmas, Definitions, Proofs and so on in the following way.

% \begin{theorem}
% \label{john.smith:theorem1}
% This is the first theorem.
% \end{theorem}

% \begin{proof}
% This is the proof of Theorem \ref{john.smith:theorem1}.
% \end{proof}

% \begin{corollary}
% \label{john.smith:corollary1}
% This is a Corollary to Theorem \ref{john.smith:theorem1}.
% \end{corollary}

% \begin{definition}
% This is a first definition. Please do not try to modify the counters for Theorems, Definitions and the rest of elements, which always are consecutive and related to each section.
% \end{definition}

% \begin{lemma}
% \label{john.smith:lemma1}
% Lemma text goes here.
% \end{lemma}

% \begin{proof}
% This is the proof of Lemma \ref{john.smith:lemma1}.
% \end{proof}

% \begin{algorithm}
% \label{john.smith:algorithm1}
% \indent Introduction to the Algorithm if any.
% \begin{description}
%  \item[Step 1] This is the first step
%  \item[Step 2] This is the second step
%  \item[Step 3] This is the third step
%  \item[Step 4] This is the fourth step
% \end{description}
% \end{algorithm}

% Please, note the usage of a \verb|description| environment for the different steps in an algorithm instead of an \verb|itemize| environment.


% \section*{Acknowledgement}
% You can include acknowledgements at the end of your manuscript by using the \verb|acknowledgement| environment with no section numbering as in this example.


% \section*{Appendix}
% When placed at the end of your manuscript (as opposed to at the end of the whole book), the numbering of tables, figures, and equations in the appendix section continues on from that in the main text. You can use the \verb|appendix| environment in your manuscript if you wish again with no section numbering as in this example.

% Please, make sure that all references from the list are cited in the text. Those not cited should be moved to a separate
% \textit{Further Reading} section. The reference list should be \textit{sorted} in alphabetical order by author and coauthors. If there are several works by the same author (or same author and coauthors) they should be \textit{ordered chronologically by year of publication}.
%-----------------------------------------------------------------------------------------------------------------------


%-----------------------------------------------------------------------------------------------------------------------

%-----------------------------------------------------------------------------------------------------------------------


\end{document} 
1406.1882.extracted_bib
{smithetal:2006 Doe, J., Smith, A. and Smith, J. (2006) Hand book on first useful conclusions on the well-known interesting research topic. Interesting topics editions. %
1406.1882.extracted_bib
{smith:1999 Smith, J. (1999) A first introductory research on some interesting topics. Local journal of interesting topics research, 1, 1--34. %
1406.1882.extracted_bib
{smith:2005 Smith, J. (2005) A second and not so introductory research on the same interesting topics. Local journal of interesting topics research, 6, 2345--2358. %
1406.1882.extracted_bib
{smithdoe:1999 Smith, A. and Doe, J. (1999) Discussion on: A first introductory research on some interesting topics. Local journal of interesting topics research, 1, 35--41.
1406.1882.extracted_bib
{B1998 Buchinsky, Moshe (1998) Recent advances in quantile regression models: A practical guideline for empirical research. The journal of human resources, 33, 88--126. %
1406.1882.extracted_bib
{cade2003 Cade, Brian S. and Noon, Barry R. (2003) Frontiers in Ecology and the Environment A gentle introduction to quantile regression for ecologists, 1, 412--420. %
1406.1882.extracted_bib
{C1994 Chamberlain, Gary (1994) Advances in Econometrics Quantile regression, censoring and the structure of wages.
1406.1882.extracted_bib
{CM1962 Conway, Richard W. and Maxwell, William L. (1962) A queuing model with state dependent service rate. Journal of industrial engineering, 12, 132--136.
1406.1882.extracted_bib
{DPP2007 Dunson, David B. and Pillai, Natesh and Park, Ju-Hyun. (2007) Bayesian density regression. Journal of the royal statistical society: Series B, 69, 163--183.
1406.1882.extracted_bib
{DS2005 Dunson, David B. and Stanford, Joseph B. (2005) Bayesian inferences on predictors of conception probabilities. Biometrics, 61, 126--133. %
1406.1882.extracted_bib
{GCJG2008 Guikema, Seth D. and Coffelt, Jeremy P. (2008) Risk analysis : an official publication of the Society for Risk Analysis A flexible count data regression model for risk analysis, 28, 213--223.
1406.1882.extracted_bib
{KW2012 Karabatsos, George and Walker, Stephen G. (2012) Adaptive-modal Bayesian nonparametric regression. Electronic journal of statistics, 6, 2038--2068.
1406.1882.extracted_bib
{KB1978 Koenker, Roger and Bassett, Gilbert (1978) Regression quantiles. Econometrica, 46, 33--50. %
1406.1882.extracted_bib
{levin2001 Levin, Jesse (2006) For whom the reductions count: A quantile regression analysis of class size and peer effects on scholastic achievement Empirical Economics, 1, 221--246.
1406.1882.extracted_bib
{MS2002 Machado, Jos\`e Ant\'onio Ferreira and Santos Silva, Jo\~ao M.C.. (2005) Quantiles for counts. Journal of the american statistical association, 100, 1226--1237. %
1406.1882.extracted_bib
{MSKBB2003 Minka, Thomas P. and Shmueli, Galit and Kadane, Joseph B. and Borle, Sharad and Boatwright, Peter (2003) Technical report Computing with the {COM-{Poisson distribution.
1406.1882.extracted_bib
{MPRB2006 M{\oller, J. and Pettitt, A. N. and Reeves, R. and Berthelsen, K. K. (2006) An efficient Markov chain Monte Carlo method for distributions with intractable normalising constants. Biometrika, 932, 451--458.
1406.1882.extracted_bib
{MGM06 Murray, Ian and Ghahramani, Zoubin and MacKay, David J. C. (2006) MCMC for doubly-intractable distributions. Proceedings of the 22nd Annual UAI Conference, 359--366.
1406.1882.extracted_bib
{N2000 Neal, Radford M. (2000) Markov chain sampling methods for {Dirichlet process mixture models. Journal of computational and graphical statistics, 9, 249--265. %
1406.1882.extracted_bib
{S2004 Schulze, Niels (2004) PhD thesis Applied Quantile Regression: Microeconometric, Financial, and Environmental Analyses.
1406.1882.extracted_bib
{SMKBB2005 Shmueli, Galit and Minka, Thomas P. and Kadane, Joseph B. and Borle, Sharad and Boatwright, Peter (2008) A useful distribution for fitting discrete data: revival of the {Conway-{Maxwell-{Poisson distribution. Journal of the royal statistical society: Series C, 54, 127--142. %
1406.1882.extracted_bib
{wei2006 Wei, Ying and He, Xuming (2006) Annals of Statistics Conditional growth charts, 1, 2069--2097. %
1406.1882.extracted_bib
{W2006 Winkelmann, Rainer (2006) Journal of Industrial Engineering Reforming health care: Evidence from quantile regressions for counts, 25, 131--145.