# Multivariate and regression models for directional data based on projected
Pólya trees

Luis E. Nieto-Barajas
Department of Statistics, ITAM, Mexico

###### Abstract

Projected distributions have proved to be useful in the study of circular and
directional data. Although any multivariate distribution can be used to
produce a projected model, these distributions are typically parametric. In
this article we consider a multivariate Pólya tree on $\mathbb{R}^{k}$ and
project it to the unit hypersphere $\mathbb{S}^{k}$ to define a new Bayesian
nonparametric model for directional data. We study the properties of the
proposed model and in particular, concentrate on the implied conditional
distributions of some directions given the others to define a directional-
directional regression model. We also define a multivariate linear regression
model with Pólya tree error and project it to define a linear-directional
regression model. We obtain the posterior characterisation of all models and
show their performance with simulated and real datasets.

Keywords: Bayesian nonparametrics, circular data, directional data, projected
normal.

##  1 Introduction

Directional data arise from the observation of unit vectors in a
$k$-dimensional space, $\mathbb{S}^{k}$, which is also called $k$-dimensional
unit hypersphere. These data can also be represented as $k-1$ angles in polar
coordinates, also known as direction cosines. Particular cases arise when
$k=2$ producing circular data, and when $k=3$ producing spherical data. If
data are only observed in a portion of a hypersphere, they are known as axial
data. In general, directional data occur in different disciplines like
biology, geophysics, meteorology, ecology and environmental sciences. Specific
applications include the study of wind directions, orientation data in
biology, direction of birds migration, directions of fissures propagation in
concrete and other materials, orientation of geological deposits, and the
analysis of mammalian activity patterns in ecological reserves, among others.
For a survey on the area, the reader is referred to Mardia,  (1972), Fisher,
(1987) and Mardia and Jupp,  (1999).

The study of directional data started in the early years of the XX century.
For instance, the von Mises distribution (Von Mises, , 1918), wrapped
distributions (Lévy, , 1939) and projected, offset or displaced distributions
(Klotz, , 1964). In the XXI century, thanks to the advances of the statistical
science, new more elaborated models have been proposed, specially in the
context of bayesian nonparametrics. For instance, Núñez-Antonio et al.,
(2015) proposed a Dirichlet process mixture of wrapped and projected normals,
and Nieto-Barajas and Núñez-Antonio,  (2021) introduced a projected Pólya tree
for the analysis of circular data.

Following Jammalamadaka and SenGupta,  (2001), regression models for circular
response variables can be classified according to whether the explanatory
variables are linear (linear-circular) or circular (circular-circular). For
the linear-circular case Gould,  (1969) considered a von Mises distribution
for the response with mean parameter expressed in terms of a linear predictor
of explanatory variables, later Fisher and Lee,  (1992) generalised the
previous model by including a link function of the explanatory variables and
they further proposed a concentration parameter regression. Presnell et al.,
(1998) considered a multivariate normal regression model with linear responses
and linear predictors and projected the responses to the unit hypersphere. For
the circular-circular case, Jammalamadaka and Sarma,  (1993) defined the
conditional first trigonometric moment of the response in terms of functions
of the explanatory variable and Downs and Mardia,  (2002) considered a von
Mises distribution for the response and proposed a particular link function of
a circular explanatory variable.

The idea of this article is to introduce a multivariate model for directional
data of any dimension $k$ based on a projected Pólya tree. By considering the
implied conditional distributions we also define circular-circular, or better
said directional-directional regression models. If before projection we
consider a multivariate median regression model with Pólya tree errors and
with linear explanatory variables, the projected model for the responses is a
linear-directional regression model.

The rest of the paper is organized as follows. In Section 2 we set the
notation and present basic ideas about multivariate Pólya trees. In Section 3
we introduce the directional model as a projected Pólya tree prior and define
the two types of regression models. In Section 4 we describe how to perform
posterior inference via a data augmentation technique and consider hyper
priors for some of the parameters. We illustrate the performance of our
proposal in Section 5, for the spherical case ($k=3$), via simulations and the
analysis of real datasets with and without covariates. We conclude with some
remarks in Section 6.

Before proceeding we introduce notation: $\mbox{Ga}(\alpha,\beta)$ denotes a
gamma density with mean $\alpha/\beta$; $\mbox{N}(\mu,\tau)$ denotes a normal
density with mean $\mu$ and precision $\tau$;
$\mbox{N}_{k}(\boldsymbol{\mu},{\bf T})$ denotes a $k$-variate normal density
with mean vector $\boldsymbol{\mu}$ and precision matrix ${\bf T}$; and
$\mbox{Dir}(\boldsymbol{\alpha})$ denotes a dirichlet density with parameter
vector $\boldsymbol{\alpha}$.

##  2 Multivariate Pólya tree

In this section we recall the definition of a multivariate Pólya tree and set
notation. Let $(\mathbb{R}^{k},\mathcal{B}^{k})$ be a measurable space for
$k\in\mathbb{N}$. There are several ways of defining and denoting the nested
partition $\Pi$ (Hanson, , 2006; Jara et al., , 2009; Filippi and Holmes, ,
2017). For simplicity, we define the partition as the cross product of
univariate partitions and use the notation of Nieto-Barajas and Müller,
(2012). In other words, $\Pi=\\{B_{m,\boldsymbol{j}_{m}}\\}$ such that
$B_{m,\boldsymbol{j}_{m}}=B_{m,j_{m,1}}\times\cdots\times B_{m,j_{m,k}}$, for
$j_{m,l}=1,\ldots,2^{m}$, $l=1,\ldots,k$ and $m=1,2,\ldots$. The index $m$
denotes the level of the tree and the vector
$\boldsymbol{j}_{m}=(j_{m,1},\ldots,j_{m,k})$ locates the partitioning subset
within the level. In general, the set $B_{m,\boldsymbol{j}_{m}}$ splits into
$2^{k}$ disjoint subsets $B_{m+1,\boldsymbol{j}_{m+1}}$ such that
$\boldsymbol{j}_{m+1}\in\\{\times_{l=1}^{k}(2j_{m,l}-1,2j_{m,l})\\}$. At each
level $m$ we will have a partition of size $2^{km}$. We associate random
branching probabilities $Y_{m,\boldsymbol{j}_{m}}$ with every set
$B_{m,\boldsymbol{j}_{m}}$ such that, for example,
$Y_{m+1,2j_{1}-1,\ldots,2j_{k}-1}=F(B_{m+1,2j_{1}-1,\ldots,2j_{k}-1}\mid
B_{m,\boldsymbol{j}_{m}})$, where $F$ denotes a cumulative distribution
function (cdf) or a probability measure, indistinctively.

###### Definition 1

Let $\mathcal{A}=\\{\alpha_{m,\boldsymbol{j}_{m}}\\}$,
$\boldsymbol{j}_{m}=(j_{m,1},\ldots,j_{m,k})$, $j_{m,l}=1,\ldots,2^{m}$,
$l=1,\ldots,k$, $m=1,2,\ldots$ be a set of nonnegative real numbers. A random
probability measure $F$ on $(\mathbb{R}^{k},\mathcal{B}^{k})$ is said to have
a multivariate Pólya tree prior with parameters $(\Pi,\mathcal{A})$ if there
exists random vectors ${\bf
Y}_{m,\boldsymbol{j}_{m}}=\\{Y_{m+1,\boldsymbol{j}_{m+1}},\,\boldsymbol{j}_{m+1}\in\\{\times_{l=1}^{k}(2j_{m,l}-1,2j_{m,l})\\}\\}$
of dimension $2^{k}$ such that the following hold:

  1. (a)

All random vectors ${\bf Y}_{m,\boldsymbol{j}_{m}}$, for
$\boldsymbol{j}_{m}=(j_{m,1},\ldots,j_{m,k})$, $j_{m,l}=1,\ldots,2^{m}$,
$l=1,\ldots,k$, $m=0,1,\ldots$ are independent

  2. (b)

For every $m=0,1,\ldots$ and every $\boldsymbol{j}_{m}$, ${\bf
Y}_{m,\boldsymbol{j}_{m}}\sim\mbox{Dir}(\boldsymbol{\alpha}_{m,\boldsymbol{j}_{m}})$,
where
$\boldsymbol{\alpha}_{m,\boldsymbol{j}_{m}}=\\{\alpha_{m+1,\boldsymbol{j}_{m+1}},\,\boldsymbol{j}_{m+1}\in\\{\times_{l=1}^{k}(2j_{m,l}-1,2j_{m,l})\\}\\}$
such that $\dim(\boldsymbol{\alpha}_{m,\boldsymbol{j}_{m}})=2^{k}$

  3. (c)

For every $m=1,2,\ldots$ and every $\boldsymbol{j}_{m}$,

$F(B_{m,\boldsymbol{j}_{m}})=\prod_{h=1}^{m}Y_{m-h+1,j_{m-h+1,1},\ldots,j_{m-h+1,k}},$

where $j_{h-1,l}=\left\lceil\frac{j_{h,l}}{2}\right\rceil$ for $l=1,\ldots,k$
are recursive decreasing formulae that locate the set
$B_{m,\boldsymbol{j}_{m}}$ with its ancestors upwards in the tree.

Note that in the previous definition ${\bf
Y}_{0,1,\ldots,1}=\\{Y_{1,\boldsymbol{j}_{1}},\boldsymbol{j}_{1}\in\\{\times_{l=1}^{k}(1,2)\\}\\}$
and
$\boldsymbol{\alpha}_{0,1,\ldots,1}=\\{Y_{1,\boldsymbol{j}_{1}},\boldsymbol{j}_{1}\in\\{\times_{l=1}^{k}(1,2)\\}\\}$
are the vectors associated to the partition elements at level $m=1$.

It is desired to center the multivariate Pólya tree around a parametric
probability measure $F_{0}$. For simplicity, let us assume that
$F_{0}(x_{1},\ldots,x_{k})=\prod_{l=1}^{k}F_{0,l}(x_{l})$. Non-independence
$F_{0}$ could also be considered but a suitable transformation of the
partition sets $B_{m,\boldsymbol{j}_{m}}$ would be required (e.g. Jara et al.,
, 2009). Therefore, we proceed by matching the marginal partitions
$\left\\{B_{m,j_{m,l}}\right\\}$ with the dyadic quantiles of the marginals
$F_{0,l}$, i.e.,

$B_{m,j_{m,l}}=\left(F_{0,l}^{-1}\left(\frac{j_{m,l}-1}{2^{m}}\right),F_{0,l}^{-1}\left(\frac{j_{m,l}}{2^{m}}\right)\right],$
(1)

for $j_{m,l}=1,\ldots,2^{m}$ and $l=1,\ldots,k$. We further define
$\boldsymbol{\alpha}_{m,\boldsymbol{j}_{m}}=(\alpha\varphi(m+1),\ldots,\alpha\varphi(m+1))$
where $\alpha>0$ can be interpreted as the precision parameter (Walker and
Mallick, , 1997) and the function $\varphi$ controls the speed at which the
variance of the branching probabilities moves down in the tree. As suggested
by Watson et al.,  (2017) we take $\varphi(m)=m^{\delta}$ with $\delta>1$ to
define an absolutely continuous multivariate Pólya tree (Kraft, , 1964). It is
not difficult to prove that a multivariate Pólya tree, defined in this way,
satisfies
$\mbox{E}\\{F(B_{m,\boldsymbol{j}_{m}})\\}=F_{0}(B_{m,\boldsymbol{j}_{m}})=(1/2)^{km}$.

In practice we need to stop partitioning the space at a finite level $M$ to
define a finite tree process. At the deepest level $M$, we can spread the
probability within each set $B_{M,\boldsymbol{j}_{M}}$ according to $f_{0}$,
the density associated to $F_{0}$. In this case the random probability measure
defined will have a multivariate density of the form

$f({\bf x})=\left\\{\prod_{m=1}^{M}Y_{m,\boldsymbol{j}_{m}^{({\bf
x})}}\right\\}2^{kM}f_{0}({\bf x}),$ (2)

where ${\bf x}^{\prime}=(x_{1},\ldots,x_{k})\in\mathbb{R}^{k}$, and with
$\boldsymbol{j}_{m}^{({\bf x})}=(j_{m,1}^{(x_{1})},\ldots,j_{m,k}^{(x_{k})})$
identifying the set at level $m$ that contains ${\bf x}$. This maintains the
condition $\mbox{E}(f)=f_{0}$. By taking $M\to\infty$ we recover the
(infinite) multivariate Pólya tree of Definition 1.

It is well known (e.g. Lavine, , 1992) that Pólya tree densities are
discontinuous at the boundaries of the partitions. To overcome this feature,
an extra mixture with respect to the parameters of the centring measure is
imposed, that is, $f_{0}$ is replaced by $f_{0}(\cdot\mid\boldsymbol{\eta})$
and a prior $f(\boldsymbol{\eta})$ is placed to induce smoothness.

##  3 Projected tree models

###  3.1 Directional and circular models

We are now in a position to construct the projected Pólya tree. Let us assume
a multivariate random vector ${\bf X}^{\prime}=(X_{1},\ldots,X_{k})$ such that
${\bf X}\mid f\sim f$ and $f$ is given in (2). We project the random vector
${\bf X}$ to the unit circle by defining ${\bf U}={\bf X}/R$, where $R=||{\bf
X}||$ is the resultant length of the vector. Alternatively, we can work with
the multivariate polar coordinates transformation ${\bf
X}\leftrightarrow(\boldsymbol{\Theta},R)$, where
$\boldsymbol{\Theta}=(\Theta_{1},\ldots,\Theta_{k-1})$ are the $k-1$ angles
such that $\boldsymbol{\theta}\in{\bf H}$ with ${\bf
H}=\\{\theta_{l}\in[0,\pi],l=1,\ldots,k-2,\theta_{k-1}\in[0,2\pi],\\}$ and
$R\in\mathbb{R}^{+}$. The inverse transformation is $X_{1}=R\cos\Theta_{1}$,
$X_{2}=R\sin\Theta_{1}\cos\Theta_{2}$, $\ldots$,
$X_{k-1}=R(\prod_{l=1}^{k-2}\sin\Theta_{l})\cos\Theta_{k-1}$,
$X_{k}=R\prod_{l=1}^{k-1}\sin\Theta_{l}$. Thus, after some algebra, we obtain
that the corresponding Jacobian becomes
$J=(-r)^{k-1}\prod_{l=1}^{k-2}(\sin\theta_{l})^{k-l-1}$. Then, the induced
marginal density for the $k-1$ angles $\boldsymbol{\Theta}$ has the form

$f(\boldsymbol{\theta})=\int_{0}^{\infty}\left\\{\prod_{m=1}^{M}Y_{m,\boldsymbol{j}_{m}^{({\bf
x}(\boldsymbol{\theta},r))}}\right\\}2^{kM}f_{0}({\bf
x}(\boldsymbol{\theta},r))\,|J|\,\mbox{d}r,$ (3)

where we have made explicitly that ${\bf x}$ has to be written in terms of the
multivariate polar coordinates $(\boldsymbol{\theta},r)$. We will refer to the
density $f(\boldsymbol{\theta})$, given in (3), as the $(k-1)$-directional
model based on a projected Pólya tree, and will be denoted by
$\mbox{PPT}_{k-1}(\alpha,\varphi,f_{0},M)$.

In contrast to Pólya tree densities, the projected Pólya tree (3) is not
discontinuous at the boundaries of the partitions. This is a consequence of
the marginalisation when passing from the joint density
$f(\boldsymbol{\theta},r)$ to the marginal $f(\boldsymbol{\theta})$, which can
also be seen as a mixture of the form $f(\boldsymbol{\theta})=\int
f(\boldsymbol{\theta}\mid r)f(r)\mbox{d}r$. A specific angle
$\boldsymbol{\theta}_{0}$ might come from many points in $\mathbb{R}^{k}$
defined by different resultants in polar coordinates say
$(\boldsymbol{\theta}_{0},r_{h})$, $h=1,2,\ldots$. Each of these points might
belong to different partition sets, which are added (integrated) in the
marginalisation.

In particular, we can center our projected Pólya tree on the projected normal
distribution, considered by Wang and Gelfand,  (2013), by taking $f_{0}({\bf
x})=\mbox{N}_{k}({\bf x}\mid\boldsymbol{\mu},{\bf I})$, that is, a
multivariate normal density with mean vector
$\boldsymbol{\mu}^{\prime}=(\mu_{1},\ldots,\mu_{k})$ and precision the
$k\times k$ identity matrix ${\bf I}$. In this case, the projected Pólya tree
becomes

$\displaystyle f(\boldsymbol{\theta})=$
$\displaystyle\,2^{kM}(2\pi)^{-k/2}e^{-\frac{1}{2}\boldsymbol{\mu}^{\prime}\boldsymbol{\mu}}\,\left\\{\prod_{l=1}^{k-2}(\sin\theta_{l})^{k-l-1}\right\\}I_{{\bf
H}}(\boldsymbol{\theta})$
$\displaystyle\times\int_{0}^{\infty}\left\\{\prod_{m=1}^{M}Y_{m,\boldsymbol{j}_{m}^{{\bf
x}(\boldsymbol{\theta},r)}}\right\\}r^{k-1}\exp\left[-\frac{1}{2}\left\\{r^{2}-2r\boldsymbol{\mu}^{\prime}{\bf
x}(\boldsymbol{\theta},1)\right\\}\right]\mbox{d}r.$ (4)

Marginal densities for each $\Theta_{l}$, $l=1,\ldots,k-1$ can be obtained by
standard marginalisation procedures. If we define
$\boldsymbol{\Theta}_{-l}=\\{\Theta_{h}:h\neq l\in\\{1,\ldots,k-1\\}\\}$ the
vector of all angles except $\Theta_{l}$, and ${\bf H}_{-l}$ the corresponding
support, then marginal densities have the expression

$f(\theta_{l})=\int_{{\bf
H}_{-l}}f(\boldsymbol{\theta})\mbox{d}\boldsymbol{\theta}_{-l}.$ (5)

Expression (5) becomes a circular model (Nieto-Barajas and Núñez-Antonio, ,
2021) if $l=k-1$, and becomes an axial model for $l=1,\ldots,k-2$. In either
case we will denote this model as $\mbox{PPT}_{1}(\alpha,\varphi,f_{0},M)$.

Densities for circular variables on $[0,2\pi]$ are periodic, that is
$f(\theta_{l}+2\pi)=f(\theta_{l})$, therefore moments are not defined in the
usual way (e.g. Jammalamadaka and Umbach, , 2010). Instead, the $q^{th}$
trigonometric moment of a random variable $\Theta_{l}$, $l=1,\ldots,k-1$, is a
complex number of the form
$\varphi_{l,q}=\mbox{E}(e^{iq\Theta_{l}})=a_{l,q}+ib_{l,q}$, for any integer
$q$, where $a_{l,q}=\mbox{E}\\{\cos(q\Theta_{l})\\}$ and
$b_{l,q}=\mbox{E}\\{\sin(q\Theta_{l})\\}$. The mean direction of $\Theta_{l}$,
denoted by $\nu_{l}$, and a concentration measure around the mean, denoted by
$\varrho_{l}$, are defined in terms of the moment of order $q=1$ as

$\nu_{l}=\arctan(b_{l,1}/a_{l,1})\quad\mbox{and}\quad\varrho_{l}=\sqrt{a_{l,1}^{2}+b_{l,1}^{2}},$
(6)

where $\varrho_{l}\in[0,1]$ for $l=1,\ldots,k-1$. A value of $\varrho_{l}$
close to one means that $\Theta_{l}$ is highly concentrated around its mean
$\nu_{l}$, whereas a value of $\varrho_{l}$ close to zero means that
$\Theta_{l}$ is highly disperse.

Association between two random directional variables has been described via
generalisations of the Pearson correlation measure. One proposal is that of
Jammalamadaka and Sarma,  (1988) given by

$\rho(\Theta_{l},\Theta_{h})=\frac{\mbox{E}\\{\sin(\Theta_{l}-\nu_{l})\sin(\Theta_{h}-\nu_{h})\\}}{\sqrt{\mbox{Var}\\{\sin(\Theta_{l}-\nu_{l})\\}\mbox{Var}\\{\sin(\Theta_{h}-\nu_{h})\\}}},$
(7)

which maintains the properties $|\rho(\Theta_{l},\Theta_{h})|\leq 1$ and
$\rho(\Theta_{l},\Theta_{h})=0$ if $\Theta_{l}$ and $\Theta_{h}$ are
independent.

###  3.2 Regression models

The simple linear regression model for any two (non circular) random
variables, say $Y$ and $Z$, arises from a joint normal model for the bivariate
vector $f(y,z)$ and defining the conditional distribution of the random
variable of interest, say $Y$, given $Z$ of the form $f(y\mid z)$. We extend
this idea to circular variables to define a regression model for the variables
of interest, say, without loss of generality, let
$\boldsymbol{\Theta}_{I}=\\{\Theta_{1},\ldots,\Theta_{l}\\}$ be a subset of
the $(k-1)$-dimensional vector $\boldsymbol{\Theta}$, and
$\boldsymbol{\Theta}_{II}=\\{\Theta_{l+1},\ldots,\Theta_{k-1}\\}$ the
complement. If $\boldsymbol{\Theta}_{I}$ is the set of responses (variables of
interest) and $\boldsymbol{\Theta}_{II}$ the set of explanatory variables,
then the regression model
$f(\boldsymbol{\theta}_{I}\mid\boldsymbol{\theta}_{II})$ is given by

$f(\boldsymbol{\theta}_{I}\mid\boldsymbol{\theta}_{II})=\frac{f(\boldsymbol{\theta})}{\int_{{\bf
H}_{I}}f(\boldsymbol{\theta})\mbox{d}\boldsymbol{\theta}_{I}}$ (8)

where ${\bf H}_{I}$ is the cross product of the corresponding marginal
supports $H_{l}=[0,\pi]$ if $l=1,\ldots,k-2$ or $H_{l}=[0,2\pi]$ if $l=k-1$;
and $f(\boldsymbol{\theta})$ is either (3) for the general case, or (4) for
the projected normal case. Expression (8) could be seen as a directional-
directional regression model emphasising the fact that the response as well as
the explanatory variables are directional.

However, linear-directional regression models can also be defined with
projected trees. In particular, we consider a median regression model for the
unprojected variables ${\bf X}$ given a $p$-set of (linear) explanatory
variables ${\bf Z}^{\prime}=(Z_{1},\ldots,Z_{p})$, where one of them could be
an intercept, of the form ${\bf X}=\Gamma\,{\bf Z}+\boldsymbol{\epsilon}$,
with $\Gamma$ a $(k\times p)$-matrix of coefficients and
$\boldsymbol{\epsilon}^{\prime}=(\epsilon_{1},\ldots,\epsilon_{k})$ random
errors such that $\boldsymbol{\epsilon}\mid f\sim f$ where $f$ is a finite
multivariate Pólya tree process, as in (2), centred around $f_{0}({\bf
x})=N_{k}({\bf x}\mid{\bf 0},{\bf I})$. To make this model identifiable, we
need the additional constraints in the first level of the tree (Hanson and
Johnson, , 2002). In (1), for $m=1$ and $l=1,\ldots,k$, we require
$B_{1,1}=(-\infty,0]$ and $B_{1,2}=(0,\infty)$, together with fixed branching
probabilities ${\bf Y}_{0,1,\ldots,1}=(1/2^{k},\ldots,1/2^{k})$. These
constraints ensure that the marginal medians of ${\bf X}$ are $\Gamma\,{\bf
Z}$. If we finally consider the polar coordinates transformation
$\boldsymbol{\epsilon}\leftrightarrow(\Theta,R)$ and marginalise with respect
to $R$, the regression model $f(\boldsymbol{\theta}\mid{\bf z})$ is of the
form (4) with ${\bf x}(\boldsymbol{\theta},r)$ replaced by
$\boldsymbol{\epsilon}(\boldsymbol{\theta},r)={\bf
x}(\boldsymbol{\theta},r)-\Gamma\,{\bf z}$ and with $\boldsymbol{\mu}={\bf
0}$.

Depending on the dimension $k-1$ of the directional response vector
$\boldsymbol{\Theta}$, if $k-1\geq 2$, we can partition
$\boldsymbol{\Theta}=(\boldsymbol{\Theta}_{I},\boldsymbol{\Theta}_{II})$ and
define a (directional+linear)-directional regression model
$f(\boldsymbol{\theta}_{I}\mid\boldsymbol{\theta}_{II},{\bf z})$, as in (8)
with $f(\boldsymbol{\theta})$ replaced by $f(\boldsymbol{\theta}\mid{\bf z})$.
Here $\boldsymbol{\Theta}_{I}$ would be the response variable,
$\boldsymbol{\Theta}_{II}$ the directional and ${\bf Z}$ the linear
explanatory variables.

##  4 Posterior inference

Let
$\boldsymbol{\Theta}_{1},\boldsymbol{\Theta}_{2},\ldots,\boldsymbol{\Theta}_{n}$
be a sample of size $n$ of dimension $(k-1)$ vectors such that
$\boldsymbol{\Theta}_{i}\mid f\sim f$, independently, and
$f\sim\mbox{PPT}_{k-1}(\alpha,\varphi,f_{0},M)$, as in (3). We consider a data
augmentation approach (Tanner, , 1991) by defining latent resultant lengths
$R_{1},R_{2},\ldots,R_{n}$ such that $(\boldsymbol{\Theta}_{i},R_{i})$ define
the polar coordinate transformation of the multivariate linear vectors ${\bf
X}_{i}=(X_{1,i},\ldots,X_{k,i})$, for $i=1,\ldots,n$.

Then, the likelihood for $\\{{\bf Y}_{m,j_{m}}\\}$,
$\boldsymbol{j}_{m}=(j_{m,1},\ldots,j_{m,k})$, $j_{m,l}=1,\ldots,2^{m}$,
$l=1,\ldots,k$ and $m=0,1,\ldots,M-1$, given the extended data, is

$\mbox{lik}({\bf
Y}\mid\mbox{data})=\prod_{i=1}^{n}\prod_{m=1}^{M}Y_{m,\boldsymbol{j}_{m}^{({\bf
x}_{i}(\boldsymbol{\theta}_{i},r_{i}))}}=\prod_{m=1}^{M}\prod_{j_{m,1}=1}^{2^{m}}\cdots\prod_{j_{m,k}=1}^{2^{m}}Y_{m,j_{m,1},\ldots,j_{m,k}}^{N_{m,j_{m,1},\ldots,j_{m,k}}},$

where $N_{m,\boldsymbol{j}_{m}}=\sum_{i=1}^{n}I\left({\bf
x}_{i}(\boldsymbol{\theta}_{i},r_{i})\in B_{m,\boldsymbol{j}_{m}}\right)$.

Recalling from Definition 1 that the prior distribution of the vectors ${\bf
Y}_{m,\boldsymbol{j}_{m}}$ is Dirichlet with parameter
$\boldsymbol{\alpha}_{m,\boldsymbol{j}_{m}}$, and noting that the prior is
conjugate with respect to this likelihood, then the posterior distribution for
the branching probability vectors is

${\bf Y}_{m,\boldsymbol{j}_{m}}\mid\boldsymbol{\theta},{\bf
r}\sim\mbox{Dir}(\boldsymbol{\alpha}_{m,\boldsymbol{j}_{m}}+{\bf
N}_{m,\boldsymbol{j}_{m}}).$ (9)

where ${\bf
N}_{m,\boldsymbol{j}_{m}}=\\{N_{m+1,\boldsymbol{j}_{m+1}},\,\boldsymbol{j}_{m+1}\in\\{\times_{l=1}^{k}(2j_{m,l}-1,2j_{m,l})\\}\\}$.

Remember that this posterior depends on an extended version of the data. The
latent resultant lengths $R_{i}$’s have to be sampled from their corresponding
posterior predictive distribution which is simply

$f(r_{i}\mid{\bf
Y},\boldsymbol{\theta}_{i})\propto\left\\{\prod_{m=1}^{M}Y_{m,\boldsymbol{j}_{m}^{({\bf
x}_{i}(\boldsymbol{\theta}_{i},r_{i}))}}\right\\}f_{0}({\bf
x}_{i}(\boldsymbol{\theta}_{i},r_{i}))\,r_{i}^{k-1},$ (10)

for $i=1,\ldots,n$.

There are two parameters that have to be defined a-priori, one is the
precision parameter $\alpha$ and the other is the centring measure $f_{0}$.
Instead of fixing them, we suggest to assign hyper prior distributions. Say
$\alpha\sim\mbox{Ga}(a_{\alpha},b_{\alpha})$ and for $f_{0}$, in the case that
this is a multivariate normal with mean $\boldsymbol{\mu}$ and precision
matrix ${\bf I}$, we further take $\mu_{l}\sim\mbox{N}(\mu_{0},\tau_{\mu})$,
independently for $l=1,\ldots,k$. Conditional posterior distributions for
these parameters are:

$f(\alpha\mid{\bf Y})\propto\mbox{Ga}(\alpha\mid
a_{\alpha},b_{\alpha})\prod_{m=0}^{M-1}\prod_{\boldsymbol{j}_{m}}\mbox{Dir}({\bf
y}_{m,\boldsymbol{j}_{m}}\mid\boldsymbol{\alpha}_{m,\boldsymbol{j}_{m}})$ (11)

and

$f(\mu_{l}\mid\boldsymbol{\theta},{\bf
r})=\mbox{N}\left(\mu_{l}\left|\,\frac{n\bar{x}_{l}(\boldsymbol{\theta},{\bf
r})+\tau_{\mu}\mu_{0}}{n+\tau_{\mu}},n+\tau_{\mu}\right.\right),$ (12)

where $\bar{x}_{l}(\boldsymbol{\theta},{\bf
r})=\frac{1}{n}\sum_{i=1}^{n}x_{l,i}(\boldsymbol{\theta}_{i},r_{i})$, for
$l=1,\ldots,k$.

In the linear-directional regression setting, together with each
$\boldsymbol{\Theta}_{i}$ we observe a set of linear covariates ${\bf Z}_{i}$,
such that $\boldsymbol{\Theta}_{i}\mid{\bf Z}_{i},f\sim
f(\boldsymbol{\theta}_{i}\mid{\bf z}_{i})$ for $i=1,\ldots,n$, as in Section
3.2. Augmenting the data, as before, with the latent resultant lengths, we
define the one to one transformation of polar coordinates
$(\boldsymbol{\theta}_{i},R_{i})$ to the multivariate linear vectors
$\boldsymbol{\epsilon}_{i}$, where $\boldsymbol{\epsilon}_{i}={\bf
X}_{i}-\Gamma\,{\bf Z}_{i}$ for $i=1,\ldots,n$. Then, conditional posterior
distribution of the branching probabilities ${\bf Y}_{m,\boldsymbol{j}_{m}}$
is as in (9) with
$N_{m,\boldsymbol{j}_{m}}=\sum_{i=1}^{n}I\left(\boldsymbol{\epsilon}_{i}(\boldsymbol{\theta}_{i},r_{i})\in
B_{m,\boldsymbol{j}_{m}}\right)$, where the notation emphasises the dependence
on $\boldsymbol{\theta}_{i}$ and $r_{i}$.

To complete the full bayesian regression model, we need to specify the prior
for the regression coefficients $\Gamma=(\gamma_{l,h})$ for $l=1,\ldots,k$ and
$h=1,\ldots,p$. For simplicity we take
$\gamma_{l,h}\sim\mbox{N}(0,\tau_{\gamma})$ independently $\forall(l,h)$.
Therefore, the conditional posterior distribution for each $\gamma_{l,h}$,
given the observed extended data, has the form

$f(\gamma_{l,h}\mid{\bf Y},\boldsymbol{\theta},{\bf
r})\propto\mbox{N}(\gamma_{l,h}\mid
0,\tau_{\gamma})\prod_{i=1}^{n}\left\\{\prod_{m=1}^{M}Y_{m,\boldsymbol{j}_{m}^{(\boldsymbol{\epsilon}_{i}(\boldsymbol{\theta}_{i},r_{i},\Gamma))}}\right\\}f_{0}(\boldsymbol{\epsilon}_{i}(\boldsymbol{\theta}_{i},r_{i},\Gamma)),$
(13)

where now the notation also emphasises the dependence on $\Gamma$.

With equations (9) – (12), plus (13) (in the case of linear-directional
regression), we can implement a Gibbs sampler (Smith and Roberts, , 1993).
Sampling from (9) and (12) is straightforward, however, to sample from (10),
(11) and (13) we will require Metropolis-Hastings (MH) steps (Tierney, ,
1994). In those cases we suggest random walk proposal distributions with the
following specifications. Dropping the indexes, at iteration $(t+1)$: sample
$r^{*}$ from $\mbox{Ga}(\kappa_{r},\kappa_{r}/r^{(t)})$ and accept it with
probability

$\pi(r^{*},r^{(t)})=\frac{f(r^{*}\mid{\bf
Y},\boldsymbol{\theta})\mbox{Ga}(r^{(t)}\mid\kappa_{r},\kappa_{r}/r^{*})}{f(r^{(t)}\mid{\bf
Y},\boldsymbol{\theta})\mbox{Ga}(r^{*}\mid\kappa_{r},\kappa_{r}/r^{(t)})};$

sample $\alpha^{*}$ from
$\mbox{Ga}(\kappa_{\alpha},\kappa_{\alpha}/\alpha^{(t)})$ and accept it with
probability

$\pi(\alpha^{*},\alpha^{(t)})=\frac{f(\alpha^{*}\mid{\bf
Y})\mbox{Ga}(\alpha^{(t)}\mid\kappa_{\alpha},\kappa_{\alpha}/\alpha^{*})}{f(\alpha^{(t)}\mid{\bf
Y})\mbox{Ga}(\alpha^{*}\mid\kappa_{\alpha},\kappa_{\alpha}/\alpha^{(t)})};$

and sample $\gamma^{*}$ from $\mbox{N}(\gamma^{(t)},\kappa_{\gamma})$ and
accept it with probability

$\pi(\gamma^{*},\gamma^{(t)})=\frac{f(\gamma^{*}\mid{\bf
Y},\boldsymbol{\theta},{\bf
r})\mbox{N}(\gamma^{(t)}\mid\gamma^{*},\kappa_{\gamma})}{f(\gamma^{(t)}\mid{\bf
Y},\boldsymbol{\theta},{\bf
r})\mbox{N}(\gamma^{*}\mid\gamma^{(t)},\kappa_{\gamma})}.$

In any case the acceptance probability $\pi$ is truncated to the interval
$[0,1]$. The parameters $\kappa_{r}$, $\kappa_{\alpha}$ and $\kappa_{\gamma}$
are nonnegative tuning parameters, chosen appropriately to produce good
acceptance probabilities. For the examples considered here we used
$\kappa_{r}=0.5$ and $\kappa_{\alpha}=10$ for the models without covariates,
and $\kappa_{r}=3$, $\kappa_{\alpha}=15$ and $\kappa_{\gamma}=50$ for the
models with covariates, to obtain acceptance rates between $0.2$ and $0.4$,
which according to Robert and Casella,  (2010) are optimal.

As a measure of goodness of fit, for each prior scenario we will compute the
logarithm of the pseudo marginal likelihood (LPML), originally suggested by
Geisser and Eddy,  (1979). In particular, we will choose the best value for
the hyper-parameters by comparing these measures.

In general, the projected Pólya tree directional model
$f(\boldsymbol{\theta})$, as in (3) or (4), does not have a closed analytic
expression. It depends on an integral with respect to the latent $R$. We
propose to approximate the density via numerical quadrature. Say, if
$0=r^{(0)}<r^{(1)}<\cdots<r^{(T)}<\infty$ is a partition of the positive real
line, then (3) is approximated by

$f(\boldsymbol{\theta})\approx\sum_{t=1}^{T}\left\\{\prod_{m=1}^{M}Y_{m,\boldsymbol{j}_{m}^{({\bf
x}(\boldsymbol{\theta},r^{(t)}))}}\right\\}2^{kM}f_{0}({\bf
x}(\boldsymbol{\theta},r^{(t)}))(r^{(t)})^{k-1}\left\\{\prod_{l=1}^{k-2}(\sin\theta_{l})^{k-l-1}\right\\}(r^{(t)}-r^{(t-1)}).$

Specifically we take $T=100$, $r^{(T)}=||\boldsymbol{\mu}||+4$ and the rest
$r^{(t)}$ equally spread.

##  5 Numerical studies

###  5.1 Simulated prior paths

Let us consider a scenario for spherical data, $k=3$. For this, let us
envision a sphere and imagine an arrow piercing the sphere in the center from
south to north, and a table cutting the sphere right in the equator. Consider
an arbitrary point on the surface of the sphere. The angle $\theta_{1}$,
called _colatitud_ , is the angle from the arrow to the point, measured from
north to south. The angle $\theta_{2}$, called _longitude_ , is the angle
measured anti clock-wise of the projection of the point to the table. The pair
$(\theta_{1},\theta_{2})$ characterises any point on the sphere surface.

We first illustrate some paths of the projected Pólya tree centred around the
projected normal, as in (4). We took three levels of the partition $M=3$, a
precision parameter $\alpha=1$, a function $\varphi(m)=m^{\delta}$ with
$\delta=1.1$, and different values of $\boldsymbol{\mu}$. For each setting we
sampled twenty paths (densities) from the model. The joint and marginal
densities of $\Theta_{l}$, $l=1,2$, were approximated numerically with $T=100$
points.

Figure 1 contains nine panels which correspond to
$\boldsymbol{\mu}^{\prime}=(-1,-1,-1)$ (first column),
$\boldsymbol{\mu}^{\prime}=(0,0,0)$ (middle column) and
$\boldsymbol{\mu}^{\prime}=(1,1,1)$ (third column). Marginal densities for
$\Theta_{1}$ and $\Theta_{2}$ are shown in the first and second rows,
respectively. For $\Theta_{1}$ we note that the support is bounded to
$[0,\pi]$ and that the densities are not periodic, $f(0)\neq f(\pi)$, whereas
for $\Theta_{2}$, the support is $[0,2\pi]$ and the densities are periodic,
$f(0)=f(2\pi)$. If for $\Theta_{1}$ we completed the plot up to $2\pi$, the
density in $[\pi,2\pi]$ would be a reflection of the density in $[0,\pi]$.
Third row shows some population moments as boxplots, the mean marginal
directions, $\nu_{1}$ and $\nu_{2}$, marginal concentrations, $\varrho_{1}$
and $\varrho_{2}$, given in (6), and the correlation measure $\rho$ given in
(7).

When we take $\boldsymbol{\mu}^{\prime}=(-1,-1,-1)$, paths for $\Theta_{1}$
and $\Theta_{2}$ seem to have a predominant mode around $3\pi/4$ and $2\pi/3$,
respectively, with the paths for $\Theta_{2}$ slightly more concentrated than
those for $\Theta_{1}$. For $\boldsymbol{\mu}^{\prime}=(0,0,0)$ our projected
Pólya tree is centred around the uniform model, so the paths for both
$\Theta_{1}$ and $\Theta_{2}$ are very chaotic around uniform lines. The mean
marginal direction $\nu_{2}$ shows a huge dispersion in the boxplot. When we
take $\boldsymbol{\mu}^{\prime}=(1,1,1)$ paths for $\Theta_{1}$ and
$\Theta_{2}$ have predominant modes around the same value, $\pi/3$, but there
is less concentration for $\Theta_{1}$ than for $\Theta_{2}$. In the paths for
$\Theta_{2}$ there seem to be a bimodality, close to the extremes of the
interval $[0,2\pi]$, however this effect is due to the choice of plot, if we
had plotted the paths around a circle, there would only be a single
predominant mode.

Although the correlation between the two circular variables $\Theta_{1}$ and
$\Theta_{2}$ in Figure 1 is not zero, it is highly concentrated around zero
(shown as narrow boxes), with few paths having a positive/negative
correlation. In order to find combinations of $\boldsymbol{\mu}$ that induce
stronger correlations, in Figure 2 we show simulations with
$\boldsymbol{\mu}^{\prime}=(0,-1,1)$ (first column),
$\boldsymbol{\mu}^{\prime}=(-1,0,1)$ (middle column) and
$\boldsymbol{\mu}^{\prime}=(-1,1,0)$ (third column). In these three cases,
correlation boxplots (third row, fifth plot) show wider boxes, specially that
for $\boldsymbol{\mu}^{\prime}=(-1,0,1)$ (middle column). Simulated paths when
$\boldsymbol{\mu}^{\prime}=(-1,1,0)$ for $\Theta_{2}$ show a u-shape, due
again to the choice of plot but have a single predominant mode. What is
interesting in this case is that the boxplot for the mean direction $\nu_{2}$
has a huge box, which can be interpreted as very diverse locations along the
whole interval $[0,2\pi]$. For $\boldsymbol{\mu}^{\prime}=(0,-1,1)$ (first
column) however, the boxplot for the concentration $\varrho_{1}$ is very
narrow, implying almost null variability.

To illustrate the shapes of the joint distributions, Figure 3 contains single
simulations for each of the six different values of $\boldsymbol{\mu}$
considered in Figures 1 and 2. Contours of the joint densities are shown as
heat maps with darker colour indicating higher density. Joint densities in the
first row correspond to cases where the correlation is low (close to zero). In
particular, in the joint density for $\boldsymbol{\mu}^{\prime}=(0,0,0)$ we
are able to see eight rectangles, which correspond to the partitions
$\\{B_{m,j_{m},l}\\}$ needed to define the Pólya tree, the projection to polar
coordinates does not seem to smooth away the discontinuities in this
particular case. In the second row of Figure 3 we show joint distributions
with higher correlations, the first two with negative dependence and the last
one with positive dependence. As for a circular variable, the best way of
plotting a density is on a circle, for a bivariate directional vector, the
best plot would be a cylinder, rolling back top and bottom edges of the plots.
By doing so, contour for $\boldsymbol{\mu}^{\prime}=(-1,1,0)$ that could be
seen having two separate clouds, would turn into a single cloud. Additionally
in the six panels, we also include the conditional means
$\mbox{E}(\Theta_{1}\mid\theta_{2})$ and $\mbox{E}(\Theta_{2}\mid\theta_{1})$,
which are functions of the conditioning variable $\theta_{2}$ and
$\theta_{1}$, respectively. This conditional means can be interpreted as the
regression curves. Broadly speaking, conditional means for the low correlation
densities (top row) are parallel to the axis, whereas those for higher
correlation (bottom row) show more diagonal lines, specially in the region
with high probability.

###  5.2 Real data analysis

The book Fisher,  (1987) presents a series of spherical datasets. We have
selected some of them to illustrate the performance of our model. Specifically
we have chosen the sets B15 and B19. The first dataset deals with the
orientation of the dendritic field at various sites in the retina of 6 cats,
in response to different visual stimuli. The spherical data are defined by
latitude ($\vartheta_{1}$) and longitude ($\vartheta_{2}$), which take values
in $[90,180]$ and $[0,360]$ degrees, respectively. The second dataset deal
with orientation of beddings of F1 folds in Ordovician turbidites. In plain
words, F1 is a category used by geologists to define the first folding,
Ordovician is the second period of the Paleozoic era, and a tubidite is the
sediment left by an underwater current. The spherical data in this context is
defined by two measures, dip ($\vartheta_{1}$) and dip direction
($\vartheta_{2}$), which take values in $[0,90]$ and $[0,360]$ degrees,
respectively. In both cases, we transform the data into radians by defining
$\theta_{1}=\pi\eta_{1}/180$ and $\theta_{2}=2\pi\eta_{2}/360$.

In either case we have a sample $(\Theta_{1,i},\Theta_{2,i})$ for
$i=1,\ldots,n$, where the sample sizes are 94 and 64, respectively for each
datasets. To analyse the data we use the direction (spherical) projected Pólya
tree model of Section 3 centred in a projected normal with the following prior
specifications: The depth of the tree is $M=3$, the function $\varphi$ is
defined by taking $\delta=1.1$. For the precision parameter we take a prior
$\alpha\sim\mbox{Ga}(1,2)$ and for the location parameter of the centring
measure we take $\mu_{l}\sim\mbox{N}(0,\tau_{\mu})$ for $l=1,2,3$, with
varying $\tau_{\mu}\in\\{1,100\\}$ to compare. Posterior inference was based
on an MCMC with 5,500 iterations with a burn-in of 500 and keeping one of
every 5th iteration after burn-in. Using an Intel Xenon at 2.20 GHz, 2 cores
and 12 GB of RAM, it took a little more than 9 hours each run.

Table 1 shows the LPML measures obtained with the two prior specifications for
the two datasets. For dataset B15 there is little difference between the
fittings for the two different values of $\tau_{\mu}$, whereas for B19 there
is a huge gain when taking $\tau_{\mu}=100$ as compared to $\tau_{\mu}=1$. So
we take $\tau_{\mu}=1$, for B15, and $\tau_{\mu}=100$, for B19, to report
further posterior inferences also included in Table 1. Posterior 95% credible
intervals (CI) for the precision parameter $\alpha$ are around $1$ for B15 and
around $0.5$ for B19. Recall that the parameter vector $\boldsymbol{\mu}$ also
plays the role of precision in the sense that if the magnitude
$||\boldsymbol{\mu}||$ is large, then the projected Pólya tree prior is highly
concentrated around its mean. On the other hand, the prior on
$\boldsymbol{\mu}$ induces a smoothness effect in the posterior density
estimates. For the two datasets, posterior estimates of $\mu_{2}$ and
$\mu_{3}$ are around zero (highest dispersion), but for $\mu_{1}$ posterior
estimates are negative for B15 and positive for B19.

Marginal, joint density estimates and moments are reported in Figures 4 and 5,
for each dataset, respectively. The observed values of the first coordinate
$\Theta_{1}$ lie in the subinterval $[\pi/2,\pi]$ for B15, and in $[0,\pi/2]$
for B19. The second (circular) coordinate $\Theta_{2}$, in both datasets,
takes values in $[0,2\pi]$, however we have decided to rotate it to the
interval $[-\pi,\pi]$ to improve visualisations. For both datasets, marginal
density estimates (top row) follow closely the shape of the histogram, but in
a smooth way.

For B15 (Figure 4), marginal density estimate for $\Theta_{1}$ is unimodal
with the mode around $2\pi/3$, whereas for $\Theta_{2}$ marginal density is
bimodal with modes around $-2\pi/3$ and $0$. Joint density estimate, shown as
a heat map (bottom left panel), indicates two predominant modes (darker
shadows) and perhaps one smaller mode around coordinates $(2\pi/3,\pi)$.
Regression curves, $\mbox{E}(\Theta_{1}\mid\theta_{2})$ and
$\mbox{E}(\Theta_{2}\mid\theta_{1})$, are also included as dashed and dotted
lines, respectively, over the heat map. Since these two lines are almost
parallel to the axes, we do not see any explanation power in any of the two
explanatory variables. Regarding the moments, posterior distributions are
included as box plots (bottom right panel). Posterior distributions of the
mean directions have totally different behaviour, $\nu_{1}$ shows almost null
dispersion, whereas $\nu_{2}$ has a lot more dispersion. Regarding
concentrations, $\varrho_{1}>\varrho_{2}$. Linear dependence measure $\rho$
shows little dispersion and is highly concentrated around zero. In fact 95% CI
for $\rho$ is $(-0.16,0.19)$, confirming that there is almost no dependence
between these two coordinates (latitude and longitude).

For B19 (Figure 5), observed data for $\Theta_{1}$ has a histogram with two
clear modes, however, posterior estimates of the marginal density for
$\Theta_{1}$ is almost uniform with little resemblance of bimodality. On the
other hand, histogram and marginal density estimates for $\Theta_{2}$ both
show a bimodality, with modes around $-\pi/3$ and $2\pi/3$. Joint density
estimate shows two predominant modes around coordinates $(\pi/2,-\pi/3)$ and
$(\pi/3,2\pi/3)$. Gross speaking, the heat map shows a letter “C” shape, which
suggests an interesting dependence between the two angles. Regression curves,
$\mbox{E}(\Theta_{1}\mid\theta_{2})$ and $\mbox{E}(\Theta_{2}\mid\theta_{1})$,
wiggle around the two modes. For instance, if $\theta_{1}<1.3$, mean values
for $\Theta_{2}$ are positive (north-east direction), whereas if
$1.3<\theta_{1}<1.8$, mean values for $\Theta_{2}$ are negative (south-east
direction), and if $\theta_{1}>1.8$ mean values for $\Theta_{2}$ are around
zero (east-direction). Complementing this, the box plot of dependence measure
$\rho$ is mainly located in negative values, with posterior mean $-0.25$ and
posterior probability of $\rho$ being negative $0.967$. This confirms a
negative dependence between dip and dip direction in the folds. Moreover,
marginal mean directions $\nu_{1}$ and $\nu_{2}$ are both positive around
$\pi/3$ and $\pi/6$, respectively, and marginal concentrations behave
similarly to dataset B15 with $\varrho_{1}>\varrho_{2}$.

###  5.3 Real data analysis with covariates

Considering again Fisher,  (1987), dataset B23 consists of measurements of
natural remanent magnetisation in Old Red Sandstones rocks from two sites, A
and B, in Pembrokeshire Wales. The spherical observed data are inclination
($\eta_{1}$) and declination ($\eta_{2}$), with values in $[-25,52]$ and
$[167,235]$, respectively. We transform the data into radians by defining
$\theta_{1}=\pi(\eta_{1}+90)/180$ and $\theta_{2}=2\pi\eta_{2}/360$. We define
indicator variables $Z_{1}=1$ if observation comes from site A and $Z_{2}=1$
if observation comes from site B.

Available data therefore consists of
$(\Theta_{1,i},\Theta_{2,i},Z_{1,i},Z_{2,i})$ for $i=1,\ldots,n$. Data size is
$n=48$, where $35$ of them come from site A and $13$ from site B. For this
dataset we fit our linear-directional regression model as the one described in
Section 3, with explanatory variables ${\bf Z}$ of dimension $p=2$ and no
intercept. The matrix of coefficients $\Gamma$ is of dimension $3\times 2$.
Prior specifications for the Pólya tree are as in the previous example, $M=3$,
$\delta=1.1$, $\alpha\sim\mbox{Ga}(1,2)$. For the regression coefficients we
took independent prior distributions of the form
$\gamma_{l,h}\sim\mbox{N}(0,0.1)$. Posterior inference was based on an MCMC
with 20,000 iterations with a burn-in of 10,000 and keeping one of every 10th
iteration after burn-in. Computational time was around 36 hours in a machine
with the characteristics defined above.

Figure 6 displays the joint posterior predictive for sites A (left panel) and
B (right panel). Note that in both graphs we did a close up around the region
of high probability because the data are highly concentrated. The contours
seem to be mainly parallel to the axis and the correlation between
$\Theta_{1}$ and $\Theta_{2}$, in both sites, is almost null with point
estimates $0.06$ and $0.10$, respectively. Moreover,
$\mbox{P}(\rho>0\mid\mbox{data})$ is $0.7$ and $0.8$ in sites A and B,
respectively. This low association between angles is also shown by the
conditional means (dashed and dotted lines), which are mainly parallel to the
axis with slight curves around the areas of high probability. In particular,
in site B, $\mbox{E}(\Theta_{1}\mid\theta_{2})$ (dashed line) shows a diagonal
shape in the center of the joint distribution.

One apparent difference between the joint distributions in the two sites is a
shift in the location of the variable $\Theta_{1}$ which is located around 1.5
radians for site A and around 2 radians for site B. In order to understand how
this shift is captured by our model, we show in Figure 7 the joint posterior
distributions of the regression coefficients $\gamma_{l,h}$, for $l=1,2,3$ and
$h=1,2$, as dispersion diagrams. We first note that the coefficients are all
away from zero with a marginal probability larger than $0.95$. Parameters
associated to site A are shown as light colours in the figure and have
posterior estimates $\hat{\gamma}_{1,1}=0.20$, $\hat{\gamma}_{2,1}=-3.43$ and
$\hat{\gamma}_{3,1}=-1.06$. Parameters associated to site B are shown as dark
colours in the figure and have posterior estimates $\hat{\gamma}_{1,2}=-1.92$,
$\hat{\gamma}_{2,2}=-4.30$ and $\hat{\gamma}_{3,2}=-1.52$. Therefore, it seems
that the shift in the location is controlled by $\gamma_{1,1}$ and
$\gamma_{1,2}$. This makes sense since in the projection there is a one to one
relationship between $\Theta_{1}$ and $X_{1}$ via $X_{1}=R\cos(\Theta_{1})$.

##  6 Concluding remarks

We have proposed a Bayesian nonparametric model for directional data of any
dimension. Our proposal is based on the projection of a multivariate Pólya
tree to the unit hypersphere. Random densities obtained from the model turned
out to be smoother than typical densities arising from Pólya trees, which are
discontinuous at the boundaries of the partitions. Further smoothing can be
obtained by taking hyper-prior distributions on the location parameter
$\boldsymbol{\mu}$ of the centring distribution $F_{0}$. As a consequence of
the joint multivariate modelling, we can define conditional distributions of
some of the directions given the others and define directional-directional
regression curves via conditional expected values. We emphasize that
regression models defined in this way are not linear, since they are non
parametric, they can capture any linear and non linear relationship. Inclusion
of linear covariates is also possible by considering median regression models
for the unprojected directions.

Posterior inference, for models with and without linear covariates, is done in
a similar way. We just augment the data with (unobserved) latent resultants
and update the multivariate Pólya tree. To simplify the posterior dependence
on the prior choice of $\alpha$, we suggest to place a hyper-prior on this
parameter, with minimal extra effort in sampling from its conditional
posterior distribution.

Although posterior inference is theoretically simple, it is computationally
very demanding, specially for dimensions $k\geq 3$. Multivariate Pólya trees
are defined via nested partitions of $\mathbb{R}^{k}$ into $2^{km}$ sets for
each level $m=1,\ldots,M$. For $k=2$, values of $M\leq 4$ produce a manageable
number of partitioning sets, $2^{kM}\leq 256$. For $k=3$, we can only go up to
$M=3$ to produce posterior inferences in a reasonable amount of time because
$2^{kM}\leq 512$. For larger values of $k$ and or $M$, inference turns almost
computationally impossible, due to the enormous amount of time needed. One
possible solution is to define asymmetric partition depths, that is, for each
coordinate $l$, define a different maximum depth, say $M_{l}$ for
$l=1,\ldots,k$, such that the final partition size is
$2^{\sum_{l=1}^{k}M_{l}}$. However it might not be easy to determine an
optimal maximum depth for each coordinate because we do not actually observe
data in $\mathbb{R}^{k}$, we observe data in the hypersphere $\mathbb{S}^{k}$,
which is characterised by $k-1$ directions and augment them with an extra
latent resultant to produce data in the whole space $\mathbb{R}^{k}$, so we
can not easily identify coordinates that would require finer partitions.

Proposing solutions to solve the computational burden, due to the curse of
dimensionality, is open and left for further study in a future work.

## Acknowledgements

The author acknowledges support from Asociación Mexicana de Cultura, A.C.

## References

  * Downs and Mardia, (2002) Downs, T.D. and Mardia, K.V. (2002). Circular regression. Biometrika 89, 683–697.
  * Filippi and Holmes, (2017) Filippi, S., Holmes, C.C. (2017). A Bayesian nonparametric approach to testing for dependence between random variables. Bayesian Analysis 12, 919–938.
  * Fisher, (1987) Fisher, N.I., Lewis, T. and Embleton, B.J.J. (1987). Statistics analysis of spherical data. Cambridge University Press, Cambridge.
  * Fisher and Lee, (1992) Fisher, N.I. and Lee, A.J. (1992). Regression models for an angular response. Biometrics 48, 665-677.
  * Geisser and Eddy, (1979) Geisser, S. and Eddy, W.F. (1979). A predictive approach to model selection. Journal of the American Statistical Association 74, 153–160.
  * Gould, (1969) Gould, A.L. (1969). A regression technique for angular data. Biometrics 25, 683–700.
  * Hanson, (2006) Hanson, T. (2006). Inference for mixtures of finite Pólya tree models. Journal of the American Statistical Association 101, 1548–1565.
  * Hanson and Johnson, (2002) Hanson, T. and Johnson, W. (2002). Modeling regression error with a mixture of Pólya trees. Journal of the American Statistical Association 97, 1020–1033.
  * Jammalamadaka and Sarma, (1988) Jammalamadaka, S.R. and Sarma, Y.R. (1988). A correlation coefficient for angular variables, In Statistical theory and data analysis II, K. Matusita (editor), pp. 349–364. North Holland, Amsterdam.
  * Jammalamadaka and Sarma, (1993) Jammalamadaka, S.R. and Sarma, Y.R. (1993). Circular regression. In Statistical Sciences and Data Analysis, K. Matusita (editor), pp. 109–128. VSP, Utrecht.
  * Jammalamadaka and SenGupta, (2001) Jammalamadaka, S.R. and SenGupta, A. (2001). Topics in circular statistics. World Scientific, Singapore.
  * Jammalamadaka and Umbach, (2010) Jammalamadaka, S.R. and Umbach, D. (2010). Some moment properties of skew-symmetric circular distributions. Metron 68, 265–273.
  * Jara et al., (2009) Jara, A., Hanson, T. and Lesaffre, E. (2009). Robustifying generalized linear mixed models using a new class of mixtures of multivariate Pólya trees. Journal of Computational and Graphical Statistics 18, 838–860.
  * Klotz, (1964) Klotz, J. (1964). Small sample power of the bivariate sign tests of Blumen and Hodges. Annals of Mathematical Statistics 35, 1576–82.
  * Kraft, (1964) Kraft, C. (1964). A class of distribution function processes which have derivatives. Journal of Applied Probability 1, 385–388.
  * Lavine, (1992) Lavine, M. (1992). Some aspects of Pólya tree distributions for statistical modelling. Annals of Statistics 20, 1222–1235.
  * Lévy, (1939) Levy, P. (1939). L’addition des variable aléatoires définies sur une circonférence. Bulletin de la Société Mathématique de France 67, 1–41.
  * Mardia, (1972) Mardia, K.V. (1972). Statistics of directional data. Academic Press, London.
  * Mardia and Jupp, (1999) Mardia, K.V. and Jupp, P.E. (1999). Directional Statistics. Wiley, Chichester.
  * Nieto-Barajas and Müller, (2012) Nieto-Barajas, L.E. and Müller, P. (2012). Rubbery Pólya tree. Scandinavian Journal of Statistics 39, 166–184.
  * Nieto-Barajas and Núñez-Antonio, (2021) Nieto-Barajas, L.E. and Núñez-Antonio, G. (2021). Projected Pólya tree. Journal of Computational and Graphical Statistics, DOI: 10.1080/10618600.2021.1923515
  * Núñez-Antonio et al., (2015) Núñez-Antonio, G., Ausín, M.C. and Wiper, M.P. (2015). Bayesian nonparametric models of circular variables based on Dirichlet process mixtures of normal distributions. Journal of Agricultural, Biological, and Environmental Statistics 20, 47–64.
  * Presnell et al., (1998) Presnell, B., Morrison, S.P. and Littell, R.C. (1998). Projected multivariate linear models for directional data. Journal of the American Statistical Association 93, 1068–1077.
  * Robert and Casella, (2010) Robert, C.P. and Casella, G. (2010). Introducing Monte Carlo methods with R. Springer, New York.
  * Smith and Roberts, (1993) Smith, A. and Roberts, G. (1993). Bayesian computations via the Gibbs sampler and related Markov chain Monte Carlo methods. Journal of the Royal Statistical Society, Series B 55, 3–23.
  * Tanner, (1991) Tanner, M.A. (1991). Tools for statistical inference: Observed data and data augmentation methods. Springer, New York.
  * Tierney, (1994) Tierney, L. (1994). Markov chains for exploring posterior distributions. Annals of Statistics 22, 1701–1722.
  * Von Mises, (1918) Von Mises, R. (1918). Über die ganzzahligkeit der atomgewicht und verwandte Fragen. Physikalische Zeitschrift 19, 490–500.
  * Wang and Gelfand, (2013) Wang, F. and Gelfand, A.E. (2013). Directional data analysis under the general projected normal distribution. Statistical Methodology 10, 113–127.
  * Watson et al., (2017) Watson, J. Nieto-Barajas, L. and Holmes, C. (2017). Characterising variation of nonparametric random probability measures using the Kullback-Leibler divergence. Statistics 51, 558–571.
  * Walker and Mallick, (1997) Walker, S. and Mallick, B. (1997). Hierarchical generalized linear models and frailty models with Bayesian nonparametric mixing. Journal of the Royal Statistical Society, Series B 59, 845–860.

Table 1: Data analysis of the three real datasets. LPML goodness of fit measures for different values of $\tau_{\mu}$ and 95% credible intervals (CI) for some model parameters for $\tau_{\mu}=100$. | B15 | B19
---|---|---
$\tau_{\mu}$ | LPML
$1$ | $-232.7$ | $-209.6$
$100$ | $-233.3$ | $-184.6$
Par. | 95% CI
$\alpha$ | $(0.42,2.01)$ | $(0.17,0.78)$
$\mu_{1}$ | $(-0.92,-0.46)$ | $(0.26,0.63)$
$\mu_{2}$ | $(-0.12,0.35)$ | $(-0.16,0.17)$
$\mu_{3}$ | $(-0.44,0.01)$ | $(-0.13,0.18)$

Figure 1: Simulated densities of the prior $\mbox{PPT}_{2}$ with $M=3$,
$\alpha=1$, $\delta=1.1$, for varying $\boldsymbol{\mu}^{\prime}$. Marginal
densities of $\Theta_{1}$ and $\Theta_{2}$ (first two rows) and moments (third
row).

Figure 2: Simulated densities of the prior $\mbox{PPT}_{2}$ with $M=3$,
$\alpha=1$, $\delta=1.1$, for varying $\boldsymbol{\mu}^{\prime}$. Marginal
densities of $\Theta_{1}$ and $\Theta_{2}$ (first two rows) and moments (third
row).

Figure 3: Simulated densities of the prior $\mbox{PPT}_{2}$ with $M=3$,
$\alpha=1$, $\delta=1.1$, for varying $\boldsymbol{\mu}^{\prime}$. Joint
density (contour plot) and conditional means
$\mbox{E}(\Theta_{1}\mid\theta_{2})$ (dashed line) and
$\mbox{E}(\Theta_{2}\mid\theta_{1})$ (dotted line).

Figure 4: Posterior inference for dataset B15. Marginal densities for
$\Theta_{1}$ (top-left) and $\Theta_{2}$ (top-right), posterior mean (solid
line) and 95% CI (dotted line). Contour plot of joint density (bottom-left)
and conditional means $\mbox{E}(\Theta_{1}\mid\theta_{2})$ (dashed line) and
$\mbox{E}(\Theta_{2}\mid\theta_{1})$ (dotted line). Boxplot of moments
(bottom-right).

Figure 5: Posterior inference for dataset B19. Marginal densities for
$\Theta_{1}$ (top-left) and $\Theta_{2}$ (top-right), posterior mean (solid
line) and 95% CI (dotted line). Contour plot of joint density (bottom-left)
and conditional means $\mbox{E}(\Theta_{1}\mid\theta_{2})$ (dashed line) and
$\mbox{E}(\Theta_{2}\mid\theta_{1})$ (dotted line). Boxplot of moments
(bottom-right).

Figure 6: Posterior inference for dataset B23. Contour plot of joint densities
for site A (left) and site B (right).

Figure 7: Posterior inference for dataset B23. Dispersion diagrams of
regression parameters: $\gamma_{1,1}$ vs $\gamma_{2,1}$ (light dots),
$\gamma_{1,1}$ vs $\gamma_{3,1}$ (light squares), $\gamma_{1,2}$ vs
$\gamma_{2,2}$ (dark dots), $\gamma_{1,2}$ vs $\gamma_{3,2}$ (dark squares).