---
abstract: 'One way of defining probability distributions for circular variables (directions in two dimensions) is to radially project probability distributions, originally defined on $\mathbb{R}^2$, to the unit circle. Projected distributions have proved to be useful in the study of circular and directional data. Although any bivariate distribution can be used to produce a projected circular model, these distribution are typically parametric. In this article we consider a bivariate Pólya tree on $\mathbb{R}^2$ and project it to the unit circle to define a new Bayesian nonparametric model for circular data. We study the properties of the proposed model, obtain its posterior characterisation and show its performance with simulated and real datasets.'
author:
- |
    [Luis Nieto-Barajas$^1$]{} and [Gabriel Nuñez-Antonio$^2$]{}\
    [*[$^1$Department of Statistics, ITAM, Mexico]{}*]{}\
    [*[$^2$Department of Mathematics, UAM-I, Mexico]{}*]{}
title: '[**Projected Pólya Tree**]{}'
---

[*Keywords*]{}: Bayesian nonparametrics, circular data, directional data, projected normal.

Introduction {#sec:intro}
============

Directional data arise from the observation of unit vectors in the $k$-dimensional space and, consequently, they can be represented through $k-1$ angles. Thus, the sample space associated with this type of data is the $k$-dimensional unit sphere, $\mathbb{S}^k$. The most common case is when $k=2$ producing the so called circular data. This type of data is especially common in biology, geophysics, meteorology, ecology and environmental sciences. Specific applications are 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, analysis of mammalian activity patterns in ecological reserves, among others.

For a survey on the area, the reader is referred to classic literature, e.g. [@mardia:72], [@fisher:95], and . For a more recent overview of applications of circular data analysis in ecological and environmental sciences see and [@lee:10].

In recent years the development of statistical methods to analyze directional data has had a new interest. considered the case of projected linear models, studied longitudinal circular data, and introduced an elliptically symmetric angular Gaussian distribution for the study of directional data on $\mathbb{S}^k$.

While there are several ways to define probability distributions for directional random vectors, one of the simplest ways to generate distributions on $\mathbb{S}^k$ is to radially project probability distributions originally defined on $\mathbb{R}^k$. A directional distribution which has received a lot of attention is the special case where the distribution to project is a $k$-variate Normal distribution, in this case it is said the corresponding directional variable has a *projected Normal distribution* . Within a Bayesian context, this model has been studied by and for the circular case, and for the $k$-dimensional case.

Althought any bivariate distribution can be used to produce a projected circular model, these distribution are typically parametric. However, in real situations it may be preferable to consider semiparametric or nonparametric models as an alternative to properly describe the behaviour of this kind of data. In a classical context, nonparametric modelling for circular data has been typically carried out using a circular kernel density such as the von Mises distribution [e.g. @fisher:89]. Within a Bayesian nonparametric approach, Dirichlet processes mixtures (DPM) of von Mises distributions and DPM of projected normal distributions have been proposed. Other semiparametric approaches are mixtures of triangular distributions and log-spline distributions . In this work we consider a bivariate Pólya tree and project it to the unit circle to produce a projected Pólya tree model. This new Bayesian nonparametric model for circular data will be shown to be competitive with respect to other Bayesian nonparametric proposals with the advantage of the simplicity in carrying out posterior inference.

The rest of the paper is organized as follows. In Section \[sec:pt\] we set the notation and present basic ideas about univariate Pólya trees. In Section \[sec:model\] we introduce the projected Pólya tree prior and study its properties. In Section \[sec:post\] we describe how to perform posterior inference via a data augmentation technique. We illustrate the performance of our proposal in Section \[sec:numerical\] via a simulation study and the analysis of a real data set. We conclude with some remarks in Section \[sec:concl\].

Pólya Tree {#sec:pt}
==========

In this section we recall the definition of a univariate Pólya tree and set notation. Consider $(\mathbb{R},{\mathcal{B}})$ the measurable space with $\mathbb{R}$ the real line and ${\mathcal{B}}$ the Borel sigma algebra of subsets of $\mathbb{R}$. We require a binary partition tree, which using notation from , is denoted by $\Pi=\{B_{mj}: m\in\mathbb{N}, j=1,\dots,2^m\}$, where the index $m$ specifies the level of the tree and $j$ the location of the partitioning subset within the level. In general, at level $m$, the set $B_{mj}$ splits into two disjoint sets $(B_{m+1,2j-1},B_{m+1,2j})$. For every set $B_{mj}$ there is associated a branching probability $Y_{mj}$ such that $Y_{m+1,2j-1} = F(B_{m+1,2j-1} \mid B_{mj})$, and $Y_{m+1,2j}=1-Y_{m+1,2j-1} = F(B_{m+1,2j} \mid B_{mj})$, where $F$ will be used to denote a cumulative distribution function (cdf) or a probability measure indistinctively.

\[def:uniPT\] [@lavine:92]. Let ${\mathcal{A}}_m=\{\alpha_{mj},\, j=1,\ldots,2^m\}$ be non-negative real numbers, $m=1,2,\ldots,$ and let ${\mathcal{A}}=\bigcup {\mathcal{A}}_m$. A random probability measure $F$ on $(\mathbb{R},{\mathcal{B}})$ is said to have a Pólya tree prior with parameters $(\Pi,{\mathcal{A}})$, if for $m=1,2,\ldots$ there exist random variables ${\mathcal{Y}}_m=\{Y_{m,2j-1}\}$ for $j=1,\ldots,2^{m-1}$, such that the following hold:

1.  All the random variables in ${\mathcal{Y}}=\cup_m\{{\mathcal{Y}}_m\}$ are independent.

2.  For every $m=1,2,\ldots$ and every $j=1,\ldots,2^{m-1}$, $Y_{m,2j-1}\sim\be(\alpha_{m,2j-1},\alpha_{m,2j})$.

3.  For every $m=1,2,\ldots$ and every $j=1,\ldots,2^m$ $$F(B_{mj})=\prod_{k=1}^m Y_{m-k+1,j_{m-k+1}^{(m,j)}},$$ where $j_{k-1}^{(m,j)}=\lceil j_k^{(m,j)}/2 \rceil$ is a recursive decreasing formula, whose initial value is $j_{m}^{(m,j)}=j$, that locates the set $B_{mj}$ with its ancestors upwards in the tree. $\lceil\cdot\rceil$ denotes the ceiling function, and $Y_{m,2j}=1-Y_{m,2j-1}$ for $j=1,\ldots,2^{m-1}$.

A Pólya tree prior can be centred around a parametric probability measure $F_0$. The simplest way consists of matching the partition with the dyadic quantiles of the desired centring measure and keeping $\alpha_{mj}$ constant within each level $m$. More explicitly, at each level $m$ we take $$\label{bmj}
B_{mj}=\left(F_0^{-1}\left(\frac{j-1}{2^m}\right),F_0^{-1}\left(\frac{j}{2^m}\right)\right],$$ for $j=1,\ldots,2^m$, with $F_0^{-1}(0)=-\infty$ and $F_0^{-1}(1)=\infty$. If we further take $\alpha_{mj}=\alpha_m$ for $j=1,\ldots,2^m$ we get $\E\{F(B_{mj})\}=F_0(B_{mj})$.

In particular, we take $\alpha_{mj}=\alpha\rho(m)$, so that the parameter $\alpha$ can be interpreted as a precision parameter of the Pólya tree , and the function $\rho$ controls the speed at which the variance of the branching probabilities moves down in the tree. As suggested by we take $\rho(m)=m^\delta$ with $\delta>1$ to ensure the process $F$ is absolutely continuous [@kraft:64].

Main model {#sec:model}
==========

Bivariate Pólya tree
--------------------

In this section we generalize the univariate Pólya tree to a bivariate one. Let $(\mathbb{R}^2,{\mathcal{B}}^2)$ be our measurable space. There are several ways of defining and denoting the nested partition $\Pi$ . For simplicity, we define the partition as the cross product of univariate partitions and use the notation of presented in Section \[sec:intro\]. In other words, $\Pi=\{B_{m,j,k}\}$ such that $B_{m,j,k}=B_{m,j}\times B_{m,k}$, for $j,k=1,\ldots,2^{m}$ and $m=1,2,\ldots$. The index $m$ denotes the level of the tree and the pair $(j,k)$ locates the partitioning subset within the level. In general, the set $B_{m,j,k}$ splits into four disjoint subsets $(B_{m+1,2j-1,2k-1},B_{m+1,2j-1,2k},B_{m+1,2j,2k-1},B_{m+1,2j,2k})$. At each level $m$ we will have a partition of size $4^m$. We associate random branching probabilities $Y_{m,j,k}$ with every set $B_{m,j,k}$ such that, for example, $Y_{m+1,2j-1,2k-1}=F(B_{m+1,2j-1,2k-1}\mid B_{m,j,k})$, where again $F$ denotes a cdf or a probability measure, indistinctively.

\[def:biPT\] Let ${\mathcal{A}}=\{\alpha_{m,j,k}\}$, $j,k=1,\ldots,2^m$, $m=1,2,\ldots$ be a set of nonnegative real numbers. A random probability measure $F$ on $(\mathbb{R}^2,{\mathcal{B}}^2)$ is said to have a bivariate Pólya tree prior with parameters $(\Pi,{\mathcal{A}})$ if there exists random vectors $\bY_{m,j,k}=(Y_{m+1,2j-1,2k-1},Y_{m+1,2j-1,2k},Y_{m+1,2j,2k-1},Y_{m+1,2j,2k})$ such that the following hold:

1.  All random vectors $\bY_{m,j,k}$, $j,k=1,\ldots,2^m$ and $m=0,1,2,\ldots$ are independent

2.  For every $m=0,1,\ldots$ and every $j,k=1,\ldots,2^m$, $\bY_{m,j,k}\sim\dir({\boldsymbol{\alpha}}_{m,j,k})$, where ${\boldsymbol{\alpha}}_{m,j,k}=(\alpha_{m+1,2j-1,2k-1},\alpha_{m+1,2j-1,2k},\alpha_{m+1,2j,2k-1},\alpha_{m+1,2j,2k})$

3.  For every $m=1,2,\ldots$ and every $j,k=1,\ldots,2^m$, $$F(B_{m,j,k})=\prod_{l=1}^m Y_{m-l+1,j_{m-l+1}^{m,j,k},k_{m-l+1}^{m,j,k}},$$ where $j_{l-1}^{(m,j,k)}=\left\lceil\frac{j_l^{(m,j,k)}}{2}\right\rceil$ and $k_{l-1}^{(m,j,k)}=\left\lceil\frac{k_l^{(m,j,k)}}{2}\right\rceil$ are recursive decreasing formulae, whose initial values are $j_m^{(m,j,k)}=j$ and $k_m^{(m,j,k)}=k$, that locate the set $B_{m,j,k}$ with its ancestors upwards in the tree.

Note that in the previous definition $\bY_{0,1,1}=(Y_{1,1,1},Y_{1,1,2},Y_{1,2,1},Y_{1,2,2})$ and ${\boldsymbol{\alpha}}_{0,1,1}=(\alpha_{1,1,1},$ $\alpha_{1,1,2},\alpha_{1,2,1},\alpha_{1,2,2})$ are the vectors associated to the partition elements at level $m=1$.

It is desired to center the bivariate Pólya tree around a parametric probability measure $F_0$. For simplicity, let us assume that $F_0(x_1,x_2)=F_{1_0}(x_1)F_{2_0}(x_2)$. Non-independence $F_0$ could also be considered but a suitable transformation of the partition sets $B_{m,j,k}$ would be required . Therefore, we proceed as in the univariate case by matching the partition $B_{m,j,k}=B_{m,j}\times B_{m,k}$ with the dyadic quantiles of the marginals $F_{1_0}$ and $F_{2_0}$, i.e., $$\label{bmj2}
B_{m,j}=\left(F_{1_0}^{-1}\left(\frac{j-1}{2^m}\right),F_{1_0}^{-1}\left(\frac{j}{2^m}\right)\right]\quad\mbox{and}\quad B_{m,k}=\left(F_{2_0}^{-1}\left(\frac{k-1}{2^m}\right),F_{2_0}^{-1}\left(\frac{k}{2^m}\right)\right],$$ for $j,k=1,\ldots,2^m$. We further define ${\boldsymbol{\alpha}}_{m,j,k}=(\alpha\rho(m+1),\ldots,\alpha\rho(m+1))$ where $\alpha>0$ is the precision parameter and $\rho(m)=m^\delta$ with $\delta>1$ to define an absolutely continuous bivariate Pólya tree. It is not difficult to prove that a bivariate Pólya tree, defined in this way, satisfies $\E\{F(B_{m,j,k})\}=F_0(B_{m,j,k})=1/4^m$.

In practice we need to stop partitioning the space at a finite level $M$ to define a finite tree process. At the lowest level $M$, we can spread the probability within each set $B_{M,j,k}$ according to $f_0$, the density associated to $F_0$. In this case the random probability measure defined will have a bivariate density of the form $$\label{eq:dPT}
f(\bx)=\left\{\prod_{m=1}^M Y_{m,j_m^{(x_1)},k_m^{(x_2)}}\right\}4^M f_0(\bx),$$ where $\bx=(x_1,x_2)\in\mathbb{R}^2$, and with $(j_m^{(x_1)},k_m^{(x_2)})$ identifying the set at level $m$ that contains $\bx$. This maintains the condition $\E(f)=f_0$. We denote a finite bivariate Pólya tree process as ${\mbox{PT}}_M(\alpha,\rho,F_0)$. By taking $M\to\infty$ we recover the (infinite) bivariate Pólya tree of Definition \[def:biPT\].

It is well known [e.g. @lavine:92] that univariate and multivariate 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 centering measure is imposed, that is, $f_0$ is replaced by $f_0(\cdot\mid\eta)$ and a prior $f(\eta)$ is placed to induce smoothness.

Projected tree
--------------

We are now in a position to construct the projected Pólya tree. Let us assume a bivariate random vector $ \bX=(X_1,X_2)$ such that $\bX\mid f \sim f$ and $f$ is given in . We project the random vector $\bX$ to the unit circle by defining $\bU=\bX/||\bX||$. Alternatively, we can work with the polar coordinates transformation $(X_1,X_2)\rightarrow (\Theta,R)$, where $\theta$ is the angle and $R=||\bX||$ is the resultant length of the vector in the plane. The inverse transformation becomes $X_1=R\cos\Theta$ and $X_2=R\sin\Theta$. Thus, the corresponding Jacobian is $J=R$. Then, the induced marginal distribution for the angle $\theta$ has the form $$\label{eq:PPT}
f(\theta)=\int_0^\infty\left\{\prod_{m=1}^M Y_{m,j_m^{(r\cos\theta)},k_m^{(r\sin\theta)}}\right\}4^M f_0(r\cos\theta,r\sin\theta)\,|J|\,\d r.$$ We will refer to the distribution $f(\theta)$, given in , as the *projected Pólya tree* and will be denoted by ${\mbox{PPT}}_M(\alpha,\rho,f_0)$.

In contrast to Pólya tree densities, the projected Pólya tree is not discontinuous at the boundaries of the partitions. The reason for the smoothing effect relies on the marginalisation when passing from the joint distribution $f(\theta,r)$ to the marginal $f(\theta)$, which can also be seen as a mixture of the form $f(\theta)=\int f(\theta\mid r)f(r)\d r$. A specific angle $\theta_0$ might come from many points in ${{\rm I}\!{\rm R}}^2$ defined by different resultants in polar coordinates say $(\theta_0,r_l)$, $l=1,2,\ldots$. Each of these points might belong to different partition sets, which are added (integrated) in the marginalisation. Therefore, no extra mixing is required to produce smooth densities.

In particular, we can center our projected Pólya tree on the projected normal distribution, considered by , by taking $f_0(\bx)=\no_2(\bx\mid{\boldsymbol{\mu}},\bI)$, that is, a bivariate normal density with mean vector ${\boldsymbol{\mu}}=(\mu_1,\mu_2)$ and variance-covariance the identity matrix $\bI$. In this case, the projected Pólya tree becomes $$\begin{aligned}
\nonumber
f(\theta)=\int_0^\infty & \left\{\prod_{m=1}^M Y_{m,j_m^{(r\cos\theta)},k_m^{(r\sin\theta)}}\right\}4^M (2\pi)^{-1}e^{-\frac{1}{2}{\boldsymbol{\mu}}'{\boldsymbol{\mu}}}\ r \\
\label{eq:PPTn}
&\times\exp\left[-\frac{1}{2}\left\{r^2-2r\left(\mu_1\cos\theta+\mu_2\sin\theta\right)\right\}\right]I_{(0,2\pi]}(\theta)\d r.\end{aligned}$$

In general, the marginal density $f(\theta)$ does not have an analytic expression. However, it can be computed numerically via quadrature. Say, if $0=r^{(0)}<r^{(1)}<\cdots<r^{(L)}<\infty$ is a partition of the positive real line, then $$f(\theta)\approx\sum_{l=1}^L f(r^{(l)}\cos\theta,r^{(l)}\sin\theta)\,|J|\left(r^{(l)}-r^{(l-1)}\right),$$ where $f(\cdot,\cdot)$ is given in . Alternatively, $f(\theta)$ can also be approximated via Monte Carlo.

Densities for circular variables are periodic, that is $f(\theta+2\pi)=f(\theta)$, therefore moments are not defined in the usual way . Instead, the $p^{th}$ trigonometric moment of a random variable $\Theta$ is a complex number of the form $\varphi_p=\E(e^{ip\Theta})=a_p+ib_p$, for any integer $p$, where $a_p=\E(\cos p\Theta)$ and $b_p=\E(\sin p\Theta)$. The mean direction of $\Theta$, $\nu_\theta$, and the concentration around the mean, $\varrho_\theta$, are defined as $$\label{eq:moments}
\nu_\theta=\arctan(b_1/a_1),\quad\varrho_\theta=\sqrt{a_1^2+b_1^2},$$ where $\varrho_\theta\in[0,1]$. A value of $\varrho_\theta$ close to one means that $\Theta$ is highly concentrated around its mean $\nu_\theta$, whereas a value of $\varrho_\theta$ close to zero means that $\Theta$ is highly disperse.

To illustrate how the paths of the projected Pólya tree look like, we consider the model centred around the projected normal, as in , with four levels of the partition ($M=4$), a precision parameter $\alpha=1$, a function $\rho(m)=m^{\delta}$ with $\delta=1.1$, and different values of ${\boldsymbol{\mu}}$. For each setting we sampled ten paths (densities) from the model. The marginal density of $\theta$ is approximated numerically with $L=100$ points.

Figure \[fig:priorppt1\] contains four panels which correspond to ${\boldsymbol{\mu}}=(0,1)$ (top left), ${\boldsymbol{\mu}}=(1,0)$ (top right), ${\boldsymbol{\mu}}=(0,-1)$ (bottom left) and ${\boldsymbol{\mu}}=(-1,0)$ (bottom right). These values of ${\boldsymbol{\mu}}$ represent specific values of the bivariate normal mean in locations around the unit circle at $\pi/2$, $2\pi$, $3\pi/2$, and $\pi$, respectively. We first note that the densities are connected in the sense that the value at $\theta=0$ coincides with the value at $\theta=2\pi$, as they should be. Within each panel we see the diversity of the paths, most of them present a multimodal behaviour. However, the predominant modes are located around the directions of the ${\boldsymbol{\mu}}$’s for each of the graphs in the four panels.

In a different scenario, we move the bivariate normal mean away from the origin to see the impact in the projected tree. This is presented in Figure \[fig:priorppt2\] that contains four panels which correspond to ${\boldsymbol{\mu}}=(0,0)$ (top left), ${\boldsymbol{\mu}}=(1,1)$ (top right), ${\boldsymbol{\mu}}=(2,2)$ (bottom left) and ${\boldsymbol{\mu}}=(5,5)$ (bottom right). The first panel corresponds to the projected tree centred around the uniform density, obtained when ${\boldsymbol{\mu}}=(0,0)$, however the simulated paths show a high variability around the centring density. As we move away from the origin (second to fourth panels) two things happen, there starts to appear a dominant mode around $\pi/4$, and the variability of the paths highly decreases. This is an interesting finding because in Pólya trees the variability is entirely controlled by the parameters $\alpha$ and $\rho(\cdot)$ [e.g. @hanson:06] and not by the centring measure. What we are seeing in this Figure \[fig:priorppt2\] is that the variability of the paths in this projected Pólya tree is also controlled by the centring measure and specifically by its location. In other words, the location parameter of the centring measure not only controls the shape of the densities but also the variability of the paths.

For the eight values of ${\boldsymbol{\mu}}$ studied, in Figure \[fig:priormoments\] we also show the prior distribution of the mean $\nu_\theta$ and the prior distribution of the concentration $\varrho_\theta$, given in . We based our prior distributions on 500 simulated paths of the corresponding projected Pólya tree. For varying ${\boldsymbol{\mu}}$ around the unit circle, we see that $\nu_\theta$ (top left panel) has a symmetric distribution with low variability and locations that move at $\pi/2$, $2\pi$, $3\pi/2$ (or $-\pi/2$), and $\pi$, respectively. However, the concentration parameter $\varrho_\theta$ (top right panel) has practically the same symmetric distribution for the different values of ${\boldsymbol{\mu}}$ around the unit circle. On the other hand, the distribution of $\nu_\theta$ for ${\boldsymbol{\mu}}=(0,0)$ is uniform, whereas this distribution is more and more less dispersed around $\pi/4$, when ${\boldsymbol{\mu}}$ moves from $(2,2)$ to $(5,5)$. Finally, the concentration parameter $\varrho_\theta$ (bottom right panel) has a dispersed distribution around 0.4, for ${\boldsymbol{\mu}}=(0,0)$, and moves to distributions less dispersed and locations that increase its values closer to one, when ${\boldsymbol{\mu}}$ moves from $(2,2)$ to $(5,5)$.

A typical concern in Bayesian nonparemtric priors is posterior consistency of the model. That is, we want to be sure that the posterior distribution concentrates around (weak) neighbours of a particular density, say $f^*(\theta)$, when the sample size $n$ goes to infinity. [@barron:98] proved that this property if satisfied as long as the prior $f$ puts positive mass around a Kullback-Leibler neighbour of $f^*$. That is, we want $\P\{\mbox{KL}(f^*,f)<\epsilon\}>0$ where $\mbox{KL}(f^*,f)=\int\log\left\{f^*(x)/f(x)\right\}f^*(x)\d x$. The following result states conditions for this to happen.

\[prop:consistency\] Let $f\sim{\mbox{PPT}}(\alpha,\rho,f_0)$ as in . Let $f^*(\theta)$ be an arbitrary density such that $\mbox{KL}(f^*,f_0)<\infty$. Then, if $\sum_{m=1}^\infty\rho(m)^{-1/2}<\infty$, as $n\to\infty$ $f$ achieves weak posterior consistency.

The idea of the proof is to prove posterior consistency for a bivariate density $f^*(\theta,r)=f^*(\theta)f^*(r)$, where $f^*(r)$ is an arbitrary density for a latent resultant $r$. Following proof of Theorem 3.1 in , by the martingale convergence theorem, there exists a collection of numbers $\{y_{m,j,k}\}$ in $[0,1]$ such that with probability one $f^*(\theta,r)=\lim_{M\to\infty}r\prod_{m=1}^M\left\{4y_{m,j_m^{(r\cos\theta)},k_m^{(r\sin\theta)}}\right\}$. Now, by and for $M\to\infty$ we have that $f(\theta,r)=\lim_{M\to\infty}r\prod_{m=1}^M\left\{4Y_{m,j_m^{(r\cos\theta)},k_m^{(r\sin\theta)}}\right\}$. The proof continues analogous to Ghosal’s. However, they show that $\eta(k)=\E\{|\log(2U_k)|\}=O(k^{-1/2})$ where $U_k\sim\be(k,k)$. In our case we need to prove that $\E\{|\log(4V_k)|\}=O(k^{-1/2})$ where $V_k\sim\be(k,3k)$. We note that $\E(2U_k)=1$ and $\V(2U_k)=1/(2k+1)$ and that $\E(4V_k)=1$ and $\V(4V_k)=3/(4k+1)$. Since both variances have the same rate of decay, our requirement is also true, proving the result.

In other words, what Proposition \[prop:consistency\] states is that, if $\rho(m)=m^\delta$, we need $\delta>2$ to satisfy the posterior consistency property. On the other hand, suggest $\delta$ close to one, say $\delta=1.1$, to maximise the dispersion of a finite tree and make the prior less informative. Since we will mostly be using finite trees we will follow this latter suggestion.

Posterior inference {#sec:post}
===================

Let $\theta_1,\theta_2,\ldots,\theta_n$ be a sample of size $n$ such that $\theta_i\mid f\sim f$, independently, and $f\sim{\mbox{PPT}}_M(\alpha,\rho,f_0)$, as in . We consider a data augmentation approach [@tanner:91] by defining latent resultant lengths $R_1,R_2,\ldots,R_n$ such that $(\Theta_i,R_i)$ define the polar coordinate transformation of the bivariate $(X_{1i},X_{2i})$ on the plane, for $i=1,\ldots,n$.

Then, the likelihood for $\{\bY_{m,j,k}\}$, $j,k=1,\ldots,2^m$ and $m=0,1,\ldots,M-1$, given the extended data, is $$\mbox{lik}(\bY\mid\data)=\prod_{i=1}^n\prod_{m=1}^M Y_{m,j_m^{(r_i\cos\theta_i)},k_m^{(r_i\sin\theta_i)}}=
\prod_{m=1}^M\prod_{j=1}^{2^m}\prod_{k=1}^{2^m}Y_{m,j,k}^{N_{m,j,k}},$$ where $N_{m,j,k}=\sum_{i=1}^n I(r_i\cos\theta_i\in B_{m,j})I(r_i\sin\theta_i\in B_{m,k})$.

Recalling from Definition \[def:biPT\] that the prior distribution of the vectors $\bY_{m,j,k}$ is Dirichlet with parameter ${\boldsymbol{\alpha}}_{m,j,k}$, and noting that the likelihood is conjugate with respect to this prior, then the posterior distribution for the branching probability vectors is $$\label{eq:posty}
\bY_{m,j,k}\mid\data\sim\dir({\boldsymbol{\alpha}}_{m,j,k}+\bN_{m,j,k}).$$ where $\bN_{m,j,k}=(N_{m+1,2j-1,2k-1},N_{m+1,2j-1,2k},N_{m+1,2j,2k-1},N_{m+1,2j,2k})$.

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 $$\label{eq:postr}
f(r_i\mid\bY,\theta_i)\propto\left\{\prod_{m=1}^M Y_{m,j_m^{(r_i\cos\theta_i)},k_m^{(r_i\sin\theta_i)}}\right\}f_0(r_i\cos\theta_i,r_i\sin\theta_i)\,r_i,$$ for $i=1,\ldots,n$.

With equations and we can implement a Gibbs sampler . Sampling from is straightforward and to sample from we will require a Metropolis-Hastings (MH) step [@tierney:94]. For this we propose a random walk proposal distribution such that, dropping the index $i$, at iteration $(t+1)$ we sample $r^*$ from $\ga(\kappa,\kappa/r^{(t)})$ and accept it with probability $$\pi(r^*,r^{(t)})=\frac{f(r^*\mid\bY,\theta)\ga(r^{(t)}\mid\kappa,\kappa/r^*)}{f(r^{(t)}\mid\bY,\theta)\ga(r^*\mid\kappa,\kappa/r^{(t)})}.$$ This latter is truncated to the interval $[0,1]$. The parameter $\kappa$ is a tuning parameter, chosen appropriately to produce good acceptance probabilities. For the examples considered here we used $\kappa=0.5$ and obtained acceptance rates between $0.2$ and $0.4$, which according to are optimal.

Numerical studies {#sec:numerical}
=================

Simulation study
----------------

We consider a model that is based on the projection of a bivariate normal mixture with four components. Specifically we define $f(\bx)=\sum_{j=1}^4\pi_j\no_2(\bx\mid{\boldsymbol{\eta}}_j,\bI)$, with ${\boldsymbol{\pi}}=(0.1,0.2, 0.4,0.3)$ and ${\boldsymbol{\eta}}_1=(1.5,1.5)$, ${\boldsymbol{\eta}}_2=(-1,1)$, ${\boldsymbol{\eta}}_3=(-1,-2)$, ${\boldsymbol{\eta}}_4=(1.5,-1.5)$, and project it to the unit circle. From this model we took two samples of sizes, $n=50$ and $n=500$. For each dataset we fitted our projected Pólya tree model ${\mbox{PPT}}_M(\alpha,\rho,f_0)$. To define the prior we took $f_0=\no_2({\boldsymbol{\mu}},\bI)$ such that the model is centred on the projected normal distribution, as in . We varied the value of the location of the bivariate normal to see the effect in the posterior estimation. In particular we took ${\boldsymbol{\mu}}\in\{(0,0),(1,1),(2,2)\}$ and a function $\rho(m)=m^\delta$, with $\delta=1.1$. We also played with different values of the precision parameter $\alpha\in\{0.5,1,2\}$. The depth of the tree was taken as $M=4$.

We ran our MCMC for $10,000$ iterations with a burn-in of $1,000$ and keeping one of every 5$^{th}$ iteration after burn-in to produce posterior inference. As a measure of goodness of fit, for each prior scenario we computed the logarithm of the pseudo marginal likelihood (LPML), originally suggested by . These numbers are summarised in Table \[tab:sim1\]. Additionally, in Figure \[fig:postsim1\] we present posterior estimates for most prior scenarios and for $n=500$. The solid line corresponds to the point estimate and the dotted lines form a 95% credible interval (CI). We accompany all graphs with a probability histogram of the data in the background.

From Table \[tab:sim1\] we can see that the fitting becomes worse (smaller LPML values) when the mean of the bivariate normal goes away from the origin. This behaviour was foreseen since the concentration (dispersion) of the projected Pólya tree prior highly increases (reduces) for larger $||{\boldsymbol{\mu}}||$ (see Figure \[fig:priormoments\]), being harder for the model to adjust to the data. Depending on the value of ${\boldsymbol{\mu}}$, some values of $\alpha$ provide better fitting than others. This latter parameter is usually interpreted as a precision parameter in Pólya trees [@hanson:06]. For smaller values of $\alpha$ the model becomes more nonparametric, and more parametric for larger values. Moreover, $\alpha$ also plays the role of a smoothing parameter. This smoothing effect can be appreciated in the top row in Figure \[fig:postsim1\], but not so much in the lower rows. What is interesting is that when ${\boldsymbol{\mu}}=(2,2)$ the posterior estimate is highly dependent on the prior and barely moves with the data, despite the large dataset of size $n=500$.

The best fitting, according to the LPML, is obtained when ${\boldsymbol{\mu}}=(0,0)$ and $\alpha=2$. This is regardless of the data size $n$. The fitting for $n=500$ is depicted in the top-right graph in Figure \[fig:postsim1\]. The posterior estimate follows smoothly the path of the data. In our experience we do not advise to go beyond $\alpha=2$ unless the data size is very large. All scenarios in Table \[tab:sim1\] were re-ran with a deeper Pólya tree with $M=6$, no real advantage was observed in terms of the LPML statistic, but the running time was a lot larger. Perhaps for the case when ${\boldsymbol{\mu}}=(2,2)$, posterior estimates are somehow better than with $M=4$, but still a lot worse than with ${\boldsymbol{\mu}}=(0,0)$.

Real data analysis
------------------

In this section, we apply our methodology to the analysis of a real dataset. A study of the interaction among species was carried out as part of a larger research project at El Triunfo biosphere reserve in Mexico in 2015. The use of camera-trapping strategies allowed ecologists to generate temporal activity information (time of the day) for three animal species, peccary, tapir and deer. The data sizes were 16, 35 and 115, respectively, and are reported in Table \[tab:data\].

This data sets has been previously analysed by , using DPM of projected normals to estimate overlapping coefficient among species. Here we fitted our projected Pólya tree model for each of the directions of the three animals. We centred our prior on a spherical bivariate normal with ${\boldsymbol{\mu}}=(0,0)$, which produces a very dispersed prior projected tree. The concentration function was $\rho(m)=m^{1.1}$ and the depth of the tree was $M=4$. We tried different values of the precision/smoothing parameter $\alpha\in\{0.5,1,2\}$ to compare. The MCMC specifications were the same as for the simulated data, and for each value of $\alpha$ we computed the statistic LPML.

The goodness of fit statistics are reported in Table \[tab:real\]. Interestingly, for the tapir and deer datasets the best fitting is achieved with $\alpha=2$, whereas for the peccary dataset the largest LPML value is obtained with $\alpha=0.5$. This is explained by the small data size of peccary which only has $n=16$ points. Alternatively, instead of selecting the best value for $\alpha$ from a range of values, we could place a hyper-prior distribution, say $\alpha\sim\ga(c_\alpha,d_\alpha)$, and update it with its corresponding conditional posterior distribution, which has the form $$f(\alpha\mid\bY)\propto\left\{\prod_{m=0}^{M-1}\prod_{j=1}^{2^m}\prod_{k=1}^{2^m}
\dir(\by_{m,j,k}\mid{\boldsymbol{\alpha}}_{m,j,k})\right\}\ga(\alpha\mid c_\alpha,d_\alpha),$$ with ${\boldsymbol{\alpha}}_{m,j,k}=(\alpha\rho(m+1),\ldots,\alpha\rho(m+1))$. Obtaining draws from this distribution requires a MH step. When taking $c_\alpha=1$ and $d_\alpha=2$ we obtain a-priori that $\P(\alpha<2)=0.98$, so we mainly concentrate on values smaller than 2. This produces posterior 95% CI for $\alpha$, for the three datasets: $(0.17,1.49)$ for peccary, $(0.40,3.00)$ for tapir, and $(0.45,2.61)$ for deer, which are consistent with the selected best values of $\alpha$. The corresponding LPML values are also reported in Table \[tab:real\].

Posterior density estimates with the best fitting settings are shown in Figure \[fig:postreal\]. Point estimates correspond to the solid lines and 95% CIs to the dotted lines. In all cases, density estimates are multimodal, perhaps for the tapir data the first mode is not so clear. For the peccaries, the directions where they move have a bounded support, mainly from 1.5 to 4.5 radians with a somehow uniform pattern. This range corresponds to the time of the day from 6:00 to 18:00 in a 24-hours clock. On the other hand, tapirs and deer appear everywhere. The mode direction where both tapirs and deer are seen is around 18:00 hrs. ($3\pi/2$ radians).

We use the mean $\nu_\theta$, as in , to summarise the preferred direction. We have assumed that the density $f$ of the directions $\theta$ is nonparametric, therefore the mean of $\theta$ is not a singe value, but a set of values whose probability distribution can be obtained. Posterior distribution for the mean direction of the three animals are presented in Figure \[fig:postmeans\] as boxplots. On a 24-hours clock, the mean direction for peccaries goes from 10:02 to 14:48 hours (2.6 to 3.9 radians) with 95% probability, the mean direction for tapirs goes from 18:10 to 23:09 hours (-1.5 to -0.2 radians) with 95% probability, and finally, the mean direction for deer goes from 16:37 to 21:35 hours (-1.9 to -0.6 radians) with 95% probability. Broadly speaking we can say that peccaries have a preferred activity-time around midday, which is totally different to the other two animals whose preferred times are around 20:30 hours, for tapirs, and around 19:00 hours, for deer.

Finally, we compare our results with alternative models, specifically we consider a parametric projected normal and a nonparametric DPM of projected normals . The corresponding LPML statistics are reported in the last two rows in Table \[tab:real\]. Now, comparing the best fit of our PPT model with the two competitors we have very interesting findings. For the peccary dataset our proposal is by far the best model with the DPM in second and the parametric model in third place. For the tapir dataset the three models practically achieve the same fit. In an attempt to explain why this happens, we recall that circular data have no beginning and end points, so histograms are better seen in a circle. The first block of points, close to zero, can be seen as a continuation of the larger points, close to $2\pi$, and therefore we could appreciate a single predominant mode characterising the data, thus a parametric (unimodal) model, like the projected normal, could do a good job describing this dataset. Lastly, for the deer data, the best fit is obtained by the DPM, followed closely by our PPT and the parametric model in a far third place.

Concluding remarks {#sec:concl}
==================

We have proposed a Bayesian nonparametric model for circular data. Our proposal is based on the projection of a bivariate Pólya tree to the unit circle. Random densities obtained from the model turned out to be smooth. This is in contrast to the bivariate densities obtained from a bivariate Pólya tree which are discontinuous at the boundaries of the partitions.

Posterior inference is simply done by augmenting the data with (unobserved) latent resultants and updating the bivariate 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. Comparing the performance of our model to other alternatives, for the three datasets studied here, we showed that our proposal is a good competitor, with the advantage of the simplicity in the posterior inference.

Generalising our model to directional data with more than two dimensions would require to project a multivariate Pólya tree on $\mathbb{R}^k$ to the unit sphere $\mathbb{S}^k$. This can be done straightforwardly by generalising the nested partition to include sets $\Pi=\{B_{m,j_1,\ldots,j_k}\}$, where at each level $m$ the partition size would be $2^{km}$. Studying the properties of this generalisation remains open and is left for future work. The inclusion of covariates in the projected (bivariate) Pólya tree also deserves study.

Acknowledgements {#acknowledgements .unnumbered}
================

This work was partial supported by the National System of Researchers, Mexico. The first author acknowledges support from *Asociación Mexicana de Cultura, A.C.* Finally, the authors are deeply thankful to professor Eduardo Mendoza from *Universidad Michoacana de San Nicolás de Hidalgo*, Mexico for the corresponding permission to use the data from project at El Triunfo biosphere reserve in Mexico.

[99]{}

Arnold, B. C. and SenGupta, A. (2006). Recent advances in the analyses of directional data in ecological and environmental sciences. [*Environmental and Ecological Statistics*]{} [**13(3)**]{}, 253–256.

Barron, A.R. (1998). The exponential convergence of posterior probabilities with implications for Bayes estimators of density functions. [*Technical Report*]{} [**7**]{}, Dept. Statistics, Univ. Illinois, Champaign.

D’Elia, A., Borgioli, C. and Scapini, F. (2001). Orientation of sandhoppers under natural conditions in repeated trials: an analysis using longitudinal directional data. [*Estuarine, Coastal and Shelf Science*]{} [**53**]{}, 839–847.

Ferreira, J.T.A.S., Juárez, M.A. and Steel, M.F.J. (2008). Directional log-spline distributions. [*Bayesian Analysis*]{} [**3**]{}, 297–316.

Filippi, S., Holmes, C.C. (2017). A Bayesian nonparametric approach to testing for dependence between random variables. [*Bayesian Analysis*]{} [**12**]{}, 919–938.

Fisher, N.I. (1989). Smoothing a sample of circular data. [*Journal of Structural Geology*]{}. [**11**]{}, 775–778.

Fisher, N.I. (1995). [*Statistical analysis of circular data.*]{} Cambridge, University Press.

Geisser, S. and Eddy, W.F. (1979). A predictive approach to model selection. [*Journal of the American Statistical Association*]{} [**74**]{}, 153–160.

Ghosal, S., Ghosh, W.F. and Ramamoorthi, R.V. (1999). Consistent semiparametric Bayesian inference about a location parameter. [*Journal of Statistical Planning and Inference*]{} [**77**]{}, 181–193.

Ghosh, K., Jammalamadaka, R. and Tiwari, R. (2003). Semiparametric Bayesian techniques for problems in circular data. *Journal of Applied Statistics* **30**, 145–161.

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.

Hanson, T. (2006). Inference for mixtures of finite Pólya tree models. [*Journal of the American Statistical Association*]{} [**101**]{}, 1548–1565.

Hernandez-Stumpfhauser, D., Breidt, F.J. and van der Woerd, M.J. (2017). The general projected normal distribution of arbitrary dimension: Modeling and Bayesian inference. [*Bayesian Analysis*]{} [**12**]{}, 113–133.

Jammalamadaka, S.R. and SenGupta, A. (2001). [*Topics in circular statistics.*]{} Singapore, World Scientific.

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.

Kraft, C. (1964). A class of distribution function processes which have derivatives. [*Journal of Applied Probability*]{} [**1**]{}, 385–388.

Lavine, M. (1992). Some aspects of Pólya tree distributions for statistical modelling. [*Annals of Statistics*]{} [**20**]{}, 1222–1235.

Lee, A. (2010). Circular data. [*Wiley Interdisciplinary Reviews: Computational Statistics*]{}, [**2(4)**]{}, 477–486.

Mardia, K.V. (1972). [*Statistics of Directional Data.*]{} London, Academic press.

Mardia, K.V. and Jupp, P.E. (2000). [*Directional Statistics.*]{} Chichester, Wiley.

McVinish, R. and Mengersen, K. (2008). Semiparametric Bayesian circular statistics. *Computational Statistics and Data Analysis* **52**, 4722–4730.

Nieto-Barajas, L.E. and Müller, P. (2012). Rubbery Pólya tree. [*Scandinavian Journal of Statistics*]{} [**39**]{}, 166–184.

Nuñez-Antonio, G. and Gutiérrez-Peña, E. (2005). A Bayesian analysis of directional data using the projected normal distribution. [*Journal of Applied Statistics*]{} [**32**]{}, 995–1001.

Nuñ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.

Nuñez-Antonio, G., Mendoza, M., Contreras-Cristán, A., Gutiérrez-Peña, E. and Mendoza, E. (2018). Bayesian nonparametric inference for the overlap of daily animal activity patterns. [*Environmental and Ecological Statistics*]{} [**25**]{}, 471–494.

Paddock, S.M. (2002). Bayesian nonparametric multiple imputation of partially observed data with ignorable nonresponse. [*Biometrika*]{} [**89**]{}, 529–538.

Paine, P.J., Preston, S.P., Tsagris, M. and Wood, A.T. (2018). An elliptically symmetric angular Gaussian distribution. [*Statistics and Computing*]{} [**28**]{}, 689–697.

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.

Rao Jammalamadaka, S. and Umbach, D. (2010). Some moment properties of skew-symmetric circular distributions. [*Metron*]{} [**68**]{}, 265–273.

Robert, C.P. and Casella, G. (2010). [*Introducing Monte Carlo methods with R*]{}. Springer, New York.

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, M.A. (1991). [*Tools for statistical inference: Observed data and data augmentation methods.*]{} Springer, New York.

Tierney, L. (1994). Markov chains for exploring posterior distributions. [*Annals of Statistics*]{} [**22**]{}, 1701–1722.

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.

Wang, F. and Gelfand, A.E. (2013). Directional data analysis under the general projected normal distribution. [*Statistical Methodology*]{} [**10**]{}, 113–127.

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.

  ---------------------- ---------- -------------- ---------------
                                                   
   ${\boldsymbol{\mu}}$   $\alpha$      $n=50$         $n=500$
         $(0,0)$           $0.5$       $-88.10$       $-849.73$
         $(0,0)$           $1.0$       $-87.07$       $-848.25$
         $(0,0)$           $2.0$     $\bf -86.98$   $\bf -847.70$
         $(1,1)$           $0.5$       $-89.13$       $-848.04$
         $(1,1)$           $1.0$       $-91.41$       $-848.07$
         $(1,1)$           $2.0$       $-95.82$       $-851.15$
         $(2,2)$           $0.5$      $-118.06$      $-1050.26$
         $(2,2)$           $1.0$      $-129.12$      $-1065.09$
         $(2,2)$           $2.0$      $-147.12$      $-1093.84$
  ---------------------- ---------- -------------- ---------------

  : LPML goodness of fit measures for simulated data.[]{data-label="tab:sim1"}

[llllllll]{}\
3.0757 & 2.7422 & 3.2214 & 0.8017 & 2.3065 & 2.6849 & 4.5517 & 4.3300\
2.3421 & 4.6541 & 2.2754 & 2.4580 & 3.3150 & 4.0887 & 4.4092 & 4.2632\
\
3.3352 & 4.6813 & 4.7835 & 5.4591 & 5.4929 & 3.6559 & 4.9567 & 4.5505\
3.7114 & 4.6214 & 5.5011 & 0.7815 & 0.4264 & 5.6929 & 4.6098 & 0.0712\
4.7340 & 4.7583 & 0.8511 & 4.5465 & 4.0871 & 1.3747 & 4.8558 & 0.9962\
4.9629 & 2.7328 & 5.9844 & 0.6099 & 5.9213 & 1.9393 & 6.2521 & 4.7322\
4.8155 & 5.1034 & 0.5203\
\
4.5338 & 4.9636 & 2.3963 & 0.1049 & 0.6435 & 1.6665 & 2.7504 & 0.5619\
5.2474 & 4.5670 & 4.4406 & 5.3001 & 4.6440 & 0.8320 & 1.5593 & 2.6858\
5.3614 & 1.5104 & 2.1596 & 4.5811 & 4.9057 & 6.1155 & 1.9216 & 3.6685\
4.7676 & 4.1158 & 3.3225 & 1.0981 & 4.7476 & 2.0472 & 4.0766 & 4.4075\
4.4901 & 5.6538 & 5.4914 & 2.0064 & 5.8532 & 0.0833 & 2.3170 & 0.6101\
5.3250 & 0.7459 & 3.4606 & 4.8188 & 4.4032 & 4.2024 & 1.5408 & 5.3556\
5.2969 & 5.9074 & 5.1198 & 4.7095 & 4.9927 & 1.5943 & 4.8544 & 0.9802\
4.7600 & 4.8139 & 4.9786 & 2.3377 & 5.0841 & 4.1202 & 6.2377 & 2.7648\
4.7023 & 4.3310 & 2.5126 & 6.0751 & 2.2459 & 1.2403 & 2.7941 & 5.0400\
5.3202 & 1.4342 & 3.2619 & 1.9663 & 4.7633 & 5.7232 & 2.1505 & 3.9069\
0.8642 & 3.5219 & 4.9393 & 2.3317 & 4.0359 & 2.0050 & 5.4570 & 4.6069\
6.0874 & 0.1445 & 0.9540 & 3.4935 & 1.6002 & 5.2741 & 0.5729 & 6.1006\
1.0324 & 4.8253 & 5.9624 & 3.5083 & 4.3276 & 4.6632 & 0.6040 & 0.7223\
3.4750 & 5.1140 & 4.9180 & 4.2155 & 4.5710 & 0.5368 & 5.1135 & 3.1823\
3.1831 & 4.4513 & 5.5457\

      $\alpha$         Peccary         Tapir           Deer
  ----------------- -------------- -------------- ---------------
        $0.5$        $\bf -23.05$     $-61.02$       $-208.31$
         $1$           $-23.22$       $-60.20$       $-206.92$
         $2$           $-24.10$     $\bf -59.57$   $\bf -205.68$
     $\ga(1,2)$        $-23.40$       $-60.15$       $-206.77$
     Proj.Normal       $-26.52$       $-59.43$       $-207.54$
   DPM Proj.Normal     $-24.64$       $-59.56$       $-204.31$

  : LPML goodness of fit measures for El Triunfo Reserve data. First three rows correspond to our PPT model.[]{data-label="tab:real"}

![[Ten simulated densities of the prior projected Pólya tree with $M=4$, $\alpha=1$, $\delta=1.1$, for varying ${\boldsymbol{\mu}}$. $(0,1)$ (top left), $(1,0)$ (top right), $(0,-1)$ (bottom left) and $(-1,0)$ bottom right.]{}[]{data-label="fig:priorppt1"}](SimDens_0__1_M4all "fig:") ![[Ten simulated densities of the prior projected Pólya tree with $M=4$, $\alpha=1$, $\delta=1.1$, for varying ${\boldsymbol{\mu}}$. $(0,1)$ (top left), $(1,0)$ (top right), $(0,-1)$ (bottom left) and $(-1,0)$ bottom right.]{}[]{data-label="fig:priorppt1"}](SimDens_1__0_M4all "fig:")

![[Ten simulated densities of the prior projected Pólya tree with $M=4$, $\alpha=1$, $\delta=1.1$, for varying ${\boldsymbol{\mu}}$. $(0,1)$ (top left), $(1,0)$ (top right), $(0,-1)$ (bottom left) and $(-1,0)$ bottom right.]{}[]{data-label="fig:priorppt1"}](SimDens_0__-1_M4all "fig:") ![[Ten simulated densities of the prior projected Pólya tree with $M=4$, $\alpha=1$, $\delta=1.1$, for varying ${\boldsymbol{\mu}}$. $(0,1)$ (top left), $(1,0)$ (top right), $(0,-1)$ (bottom left) and $(-1,0)$ bottom right.]{}[]{data-label="fig:priorppt1"}](SimDens_-1__0_M4all "fig:")

![[Ten simulated densities of the prior projected Pólya tree with $M=4$, $\alpha=1$, $\delta=1.1$, for varying ${\boldsymbol{\mu}}$. $(0,0)$ (top left), $(1,1)$ (top right), $(2,2)$ (bottom left) and $(5,5)$ bottom right.]{}[]{data-label="fig:priorppt2"}](SimDens_0__0_M4all "fig:") ![[Ten simulated densities of the prior projected Pólya tree with $M=4$, $\alpha=1$, $\delta=1.1$, for varying ${\boldsymbol{\mu}}$. $(0,0)$ (top left), $(1,1)$ (top right), $(2,2)$ (bottom left) and $(5,5)$ bottom right.]{}[]{data-label="fig:priorppt2"}](SimDens_1__1_M4all "fig:")

![[Ten simulated densities of the prior projected Pólya tree with $M=4$, $\alpha=1$, $\delta=1.1$, for varying ${\boldsymbol{\mu}}$. $(0,0)$ (top left), $(1,1)$ (top right), $(2,2)$ (bottom left) and $(5,5)$ bottom right.]{}[]{data-label="fig:priorppt2"}](SimDens_2__2_M4all "fig:") ![[Ten simulated densities of the prior projected Pólya tree with $M=4$, $\alpha=1$, $\delta=1.1$, for varying ${\boldsymbol{\mu}}$. $(0,0)$ (top left), $(1,1)$ (top right), $(2,2)$ (bottom left) and $(5,5)$ bottom right.]{}[]{data-label="fig:priorppt2"}](SimDens_5__5_M4all "fig:")

![[Prior Distribution (boxplot) of moments for the projected Pólya tree with $M=4$, $\alpha=1$, $\delta=1.1$, for varying ${\boldsymbol{\mu}}$. Mean $\nu_\theta$ (first column), and concentration $\varrho_\theta$ (second column). In each of the panels of top row $(0,1)$, $(1,0)$, $(0,-1)$ and $(-1,0)$. In each of the panels of bottom row $(0,0)$, $(1,1)$, $(2,2)$ and $(5,5)$.]{}[]{data-label="fig:priormoments"}](mt2 "fig:") ![[Prior Distribution (boxplot) of moments for the projected Pólya tree with $M=4$, $\alpha=1$, $\delta=1.1$, for varying ${\boldsymbol{\mu}}$. Mean $\nu_\theta$ (first column), and concentration $\varrho_\theta$ (second column). In each of the panels of top row $(0,1)$, $(1,0)$, $(0,-1)$ and $(-1,0)$. In each of the panels of bottom row $(0,0)$, $(1,1)$, $(2,2)$ and $(5,5)$.]{}[]{data-label="fig:priormoments"}](vt2 "fig:")

![[Prior Distribution (boxplot) of moments for the projected Pólya tree with $M=4$, $\alpha=1$, $\delta=1.1$, for varying ${\boldsymbol{\mu}}$. Mean $\nu_\theta$ (first column), and concentration $\varrho_\theta$ (second column). In each of the panels of top row $(0,1)$, $(1,0)$, $(0,-1)$ and $(-1,0)$. In each of the panels of bottom row $(0,0)$, $(1,1)$, $(2,2)$ and $(5,5)$.]{}[]{data-label="fig:priormoments"}](mt1 "fig:") ![[Prior Distribution (boxplot) of moments for the projected Pólya tree with $M=4$, $\alpha=1$, $\delta=1.1$, for varying ${\boldsymbol{\mu}}$. Mean $\nu_\theta$ (first column), and concentration $\varrho_\theta$ (second column). In each of the panels of top row $(0,1)$, $(1,0)$, $(0,-1)$ and $(-1,0)$. In each of the panels of bottom row $(0,0)$, $(1,1)$, $(2,2)$ and $(5,5)$.]{}[]{data-label="fig:priormoments"}](vt1 "fig:")

![[Posterior density estimates for simulated data with $n=500$. Across columns $\alpha=0.5$ and $\alpha=2$. Across rows ${\boldsymbol{\mu}}=(0,0)$, ${\boldsymbol{\mu}}=(1,1)$ and ${\boldsymbol{\mu}}=(2,2)$.]{}[]{data-label="fig:postsim1"}](postSimulatedData1_0_0M4a05 "fig:") ![[Posterior density estimates for simulated data with $n=500$. Across columns $\alpha=0.5$ and $\alpha=2$. Across rows ${\boldsymbol{\mu}}=(0,0)$, ${\boldsymbol{\mu}}=(1,1)$ and ${\boldsymbol{\mu}}=(2,2)$.]{}[]{data-label="fig:postsim1"}](postSimulatedData1_0_0M4a2 "fig:")

![[Posterior density estimates for simulated data with $n=500$. Across columns $\alpha=0.5$ and $\alpha=2$. Across rows ${\boldsymbol{\mu}}=(0,0)$, ${\boldsymbol{\mu}}=(1,1)$ and ${\boldsymbol{\mu}}=(2,2)$.]{}[]{data-label="fig:postsim1"}](postSimulatedData1_1_1M4a05 "fig:") ![[Posterior density estimates for simulated data with $n=500$. Across columns $\alpha=0.5$ and $\alpha=2$. Across rows ${\boldsymbol{\mu}}=(0,0)$, ${\boldsymbol{\mu}}=(1,1)$ and ${\boldsymbol{\mu}}=(2,2)$.]{}[]{data-label="fig:postsim1"}](postSimulatedData1_1_1M4a2 "fig:")

![[Posterior density estimates for simulated data with $n=500$. Across columns $\alpha=0.5$ and $\alpha=2$. Across rows ${\boldsymbol{\mu}}=(0,0)$, ${\boldsymbol{\mu}}=(1,1)$ and ${\boldsymbol{\mu}}=(2,2)$.]{}[]{data-label="fig:postsim1"}](postSimulatedData1_2_2M4a05 "fig:") ![[Posterior density estimates for simulated data with $n=500$. Across columns $\alpha=0.5$ and $\alpha=2$. Across rows ${\boldsymbol{\mu}}=(0,0)$, ${\boldsymbol{\mu}}=(1,1)$ and ${\boldsymbol{\mu}}=(2,2)$.]{}[]{data-label="fig:postsim1"}](postSimulatedData1_2_2M4a2 "fig:")

![[Posterior density estimates for the temporal activity of three animals from El Triunfo Reserve. Peccary (top left), tapir (top right) and deer (bottom).]{}[]{data-label="fig:postreal"}](postpecari_0_0M4a05 "fig:") ![[Posterior density estimates for the temporal activity of three animals from El Triunfo Reserve. Peccary (top left), tapir (top right) and deer (bottom).]{}[]{data-label="fig:postreal"}](posttapir_0_0M4a2 "fig:")

![[Posterior density estimates for the temporal activity of three animals from El Triunfo Reserve. Peccary (top left), tapir (top right) and deer (bottom).]{}[]{data-label="fig:postreal"}](postvenado_0_0M4a2)

![[Posterior distribution (boxplot) of the mean time of activity, $\nu_\theta$, for the three animals from El Triunfo Reserve.]{}[]{data-label="fig:postmeans"}](mtriunfo)
