Spatially-dependent Indian Buffet Processes¹¹1This version: September 3, 2024

Let $Z$ be a binary feature matrix, with $n$ rows for subjects and an unbounded number of columns for features, where each element $z_{ik}\in\{0,1\}$ indicates whether the $i$th subject possess the $k$th feature (i.e., $z_{ik}=1$) or not. In the Indian buffet metaphor, the rows of $Z$ represent customers, and the columns represent dishes, with $Z$ being constructed through a sequential sampling process. The first customer enters the restaurant and selects dishes according to ${\rm Po}(\alpha)$ with a parameter $\alpha$. For the $i$th customer, they sample each dish that has been previously taken with a probability of $m_{k}/i$, where $m_{k}$ denotes the number of prior customers who have chosen dish $k$. Additionally, the $i$th customer selects new dishes based on ${\rm Po}(\alpha/i)$. The selection of dishes by each customer can be represented in a binary feature allocation matrix $Z$. Such process is called Indian Buffet Process (IBP) (Griffiths and Ghahramani (2005, 2011)) and is denoted by ${\rm IBP}(\alpha)$.

As shown in Griffiths and Ghahramani (2011), the IBP can also be defined as the limit of a Beta-Bernoulli model. Here we consider the following latent variable model with finite $K$ features:

with a parameter $\alpha$, where the $b_{k}$ are generated independently and each $z_{ik}$ is independent of all other assignments under given $b_{k}$. Griffiths and Ghahramani (2011) showed that the distribution of $Z$ under $K\to\infty$ is equivalent to the one obtained by IBP. The IBP has several desirable properties. In particular, the expectation of the number of features possessed by each object under $K\to\infty$ is $\alpha$, which indicates the sparsity of feature selection. As an alternative representation of the Beta-Bernoulli model of (1), Teh et al. (2007) introduced the stick-breaking representation of (1) under finite $K$, given by

where $v_{j}\sim{\rm Beta}(\alpha,1)$. The above representation guarantees that the sequence $b_{1},\ldots,b_{K}$ is strictly decreasing. We leverage this representation to introduce spatial correlation in the subsequent section.

To take account of spatial information, let $z_{k}(s)$ is a binary random variable at location $s$ and consider the following model:

where $\sigma(x)=1/(1\!+\!e^{-x})$ is a logistic function and $u_{j}(s)$ follows a Gaussian process, independently for $j=1,\ldots,K$. We call Spatially-dependent Indian buffet process (SIBP) of the process (2). The process (2) can be regarded as an extension of the logistic stick-breaking process (Ren et al., 2011). For $n$ locations denoted by $s_{1},\ldots,s_{n}$, the joint model for $z_{ik}\equiv z_{k}(s_{i})\ (i=1,\ldots,n)$ is described as

Here $u_{ik}$ is a continuous latent variable and the joint distribution is

independently for $k=1,\ldots,K$, where $\mu$ is a location parameter, $1_{n}$ is an $n$-dimensional vector of $1$, $\tau$ is a precision parameter and $Q(\psi)$ is a correlation matrix parametrized by $\psi$. One example of $Q(\psi)$ is that its $(i,i^{\prime})$-element is $\rho(d_{ii^{\prime}};\psi)$, where $d_{ii^{\prime}}$ denotes some distance between $i$th and $i^{\prime}$th locations. Standard choices are the Matérn function, $\rho(d;\psi)=\{2^{\kappa-1}\Gamma(\kappa)\}^{-1}(d/\phi)^{\kappa}K_{\kappa}(d/\phi)$ with $\psi=(\phi,\kappa)$ and the exponential function, $\rho(d;\psi)=\exp(-d/\phi)$, where $\phi>0$ is a spatial range parameter and $\kappa>0$ is a smoothness parameter. Here $K_{\kappa}(\cdot)$ denotes the modified Bessel function of the third kind of order $\kappa$. We call Spatially-dependent Indian buffet process (SIBP) of the model (3) with spatially correlated latent variable $u_{ik}$. It should be noted that when the spatial range parameter is very large, the covariance matrix $Q(\psi)$ is degenerate and it holds that $u_{1k},\ldots,u_{nk}$ are identical with probability one and follow $N(\mu,\tau^{-1})$. This corresponds to the stick-breaking representation of IBP with different distributional assumptions for $v_{k}$.

We present some theoretical properties of SIBP. First, we prove that the expected number of non-zero entries is finite in the limiting case.

The Proposition 1 ensures that the number of sampled features is finite, namely, $P(c_{i}<\infty)=1$. Furthermore, it can be seen that the expectation of the total number of non-zero elements is $n\delta_{1}/(1-\delta_{1})$. Further, when $\mu$ and $\tau$ are chosen such that $2\delta_{2}/(1-\delta_{2})=\delta_{1}/(1-\delta_{1})$, the limit of $E[c_{i}]$ and ${\rm Var}(c_{i})$ are the same as the standard IBP. It should also be noted that the spatial range parameter $\psi$ (controlling the correlation given a fixed distance) in the Gaussian process does not appear in the formula. This indicates that the basic performance of SIBP depends only on the parameters, $\mu$ and $\tau$, of the marginal distributions.

We next investigate how the spatial correlation is introduced in the generated binary matrices through the latent Gaussian process. To this end, we evaluate the joint probability that two binary variables are both 1, as given in the following proposition:

Proposition 2 indicates that the probability that the two subjects possess the $j$th feature simultaneously is dictated by the spatial correlation $\rho_{ij}$ and the spatial range parameter $\psi$.

We next consider the number of (non-null) common features, defined as $K^{\ast}=\sum_{k=1}^{K}I(c_{k}>0)$, where $c_{k}=\sum_{i=1}^{n}z_{ik}$. Then, it holds that

where the expectation is taken with respect to the joint distribution of $u_{k}\sim N(\mu 1_{n},\tau^{-1}Q(\psi))$. Although the above expectation cannot be expressed in a closed form, we simulated the expectations by Monte Carlo integration under various choices of $(\mu,\tau,\psi)$. Specifically, we consider the exponential covariance function, $Q(\psi)_{ij}=\exp(\|s_{i}-s_{j}\|/\psi)$, and generated the location information from the uniform distribution on $[0,1]\times[0,1]$. Regarding the parameter values, we consider all the combinations of $\mu\in\{-1,0,1\}$, $\tau\in\{0.5,1\}$ and $\psi\in\{0.2,0.5,1\}$. In Figure 1, we show $E[K^{\ast}]$ as a function of the sample size. The results indicate that the parameter $\mu$ significantly controls the number of features, that is, larger value of $\mu$ leads to larger number of features in the prior distribution. On the other hand, $\tau$ and $\psi$ do not affect much on the number of features.

Let $x_{1},\ldots,x_{n}$ be an observed data, and $s_{i}\ (i=1,\ldots,n)$ be a location information (e.g., longitude and latitude) associated with $x_{i}$. We consider the following latent variable model:

where $f(x_{i}\mid\{z_{ik}\}_{k=1}^{K};\{\theta_{k}\}_{k=1}^{K})$ is a density (or probability) function of $x_{i}$. Note that $\theta_{k}$ is related to the distribution of $x_{i}$, and $(u_{1k},\ldots,u_{nk})$ are independent for $k=1,\ldots,K$. The unknown parameters in the latent variable model (4) are $\theta_{k}$, $\tau$, and $\psi$, for which we assign prior distributions for Bayesian inference. In particular, we employ a conditionally conjugate prior, $\tau\sim{\rm Ga}(a_{\tau},b_{\tau})$, $\mu\sim{\rm N}(m_{\mu},S_{\mu})$ and gamma priors for each element of $\psi$. Regarding the prior for $\theta_{k}$, we use a class of repulsive prior (e.g. Petralia et al., 2012; Xie and Xu, 2020) instead of assuming independent priors for $\theta_{1},\ldots,\theta_{K}$ to prevent redundant expression of the latent factors. Let $\pi(\theta_{k})$ be a marginal prior of $\theta_{k}$. Then, the repulsive prior for $\theta_{1},\ldots,\theta_{K}$ is given by

where $Z_{K}$ is a normalizing constant and $g:\mathcal{R}_{+}\to[0,1]$ is a strictly increasing function with $g(0)=0$. The repulsive prior eliminates the prior mass from the region where $\theta_{k}$ and $\theta_{k^{\prime}}$ are close. For the specific choice of $g(\cdot)$, we use $g(x)=I(x>\delta)(1-\delta/x)$ with a tuning constant $\delta>0$. Since the results would not be sensitive to the choice of $\delta$ as long as $\delta$ is small, we set $\delta=10^{-3}$ throughout the numerical studies in this paper.

For a general latent variable model with SIBP (4), we here provide a Gibbs sampler to simulate the posterior distributions by using a novel combination of binomial expansion and the Pólya-gamma data augmentation (Polson et al., 2013). Define $Z$, $U$ and $\Theta$ be sets of $z_{ik}$, $u_{ik}$ and $\theta_{k}$, respectively. Under the prior for these parameters, the joint posterior distribution $\pi(Z,U,\Theta,\tau,\psi|D)$ of the latent variables ($Z$ and $U$) and unknown parameters ($\Theta,\tau$ and $\psi$) given a set of observed data $D=\{x_{1},\ldots,x_{n}\}$ can be expressed as

We use Gibbs sampling to generate random samples from the above posterior distribution. The details are described as follows:

• -

  (Update of $u_{k}$)  For each $i$ and $k$, the full conditional distribution of $u_{k}$ is proportional to

  		$\displaystyle\phi(u_{k};\mu 1_{n},\tau^{-1}Q(\psi))\prod_{i=1}^{n}\prod_{j=k}^{K}\sigma(u_{ik})^{z_{ij}}\bigg{\{}1-\sigma(u_{ik})\prod_{h=1,h\neq k}^{j}\sigma(u_{ih})\bigg{\}}^{1-z_{ij}}$
  	$\displaystyle=$	$\displaystyle\ \phi(u_{k};\mu 1_{n},\tau^{-1}Q(\psi))\prod_{i=1}^{n}\prod_{j=k}^{K}\frac{(e^{u_{ik}})^{z_{ij}}(1+C_{ikj}e^{u_{ik}})^{1-z_{ij}}}{1+e^{u_{ik}}},$

  where $C_{ikj}=1-\prod_{h=1,h\neq k}^{j}\sigma(u_{ih})\in(0,1)$. By introducing latent binary variables, $s_{ikj}$ for $j=k,\ldots,K$ with $P(s_{ikj}=1)=C_{ikj}/(1+C_{ikj})$, the likelihood term can be augmented as

  $\displaystyle\prod_{j=k}^{K}\frac{(e^{u_{iik}})^{z_{ij}}(e^{u_{ik}})^{(1-z_{ij})s_{ikj}}}{1+e^{u_{ik}}}=\frac{(e^{u_{ik}})^{\sum_{j=k}^{K}\{z_{ij}+(1-z_{ij})s_{ikj}\}}}{(1+e^{u_{ik}})^{K-k+1}}.$

  Note that the full conditional probability being $s_{ikj}=1$ is $C_{ikj}e^{u_{ik}}/(C_{ikj}e^{u_{ik}}+e^{u_{ik}z_{ij}})$. Then, applying the Pólya-gamma data augmentation (Polson et al., 2013), the above function can be expressed as

  $$\frac{(e^{u_{ik}})^{\sum_{j=k}^{K}\{z_{ij}+(1-z_{ij})s_{ikj}\}}}{(1+e^{u_{ik}})^{K-k+1}}\propto\exp(\kappa_{ik}u_{ik})\int_{0}^{\infty}\exp\left(-\frac{1}{2}\omega_{ik}u_{ik}^{2}\right)f_{\rm PG}(\omega_{ik}|K-k+1,0)d\omega_{ik},$$

  where $f_{\rm PG}(\cdot|b,c)$ is the density function of ${\rm PG}(b,c)$ and

  $$\kappa_{ik}=\sum_{j=k}^{K}\{z_{ij}+(1-z_{ij})s_{ikj}\}-\frac{1}{2}(K-k+1).$$

  Note that the full conditional distribution of $\omega_{ik}$ is ${\rm PG}(K-k+1,u_{ik})$. Using the augmentation, the full conditional distribution of $u_{k}$ is $N(A_{k}^{-1}B_{k},A_{k}^{-1})$, where

  $$A_{k}={\rm diag}(\omega_{1k},\ldots,\omega_{nk})+\tau Q(\psi)^{-1},\ \ \ B_{k}=(\kappa_{1k},\ldots,\kappa_{nk})^{\top}+\tau Q(\psi)^{-1}1_{n}\mu.$$

• -

  (Update of $z_{i1},\ldots,z_{iK}$)  For each $i$, the full conditional distribution of $(z_{i1},\ldots,z_{iK})$ is proportional to

  $$f\big{(}x_{i}\mid\{z_{ik}\}_{k=1}^{K};\{\theta_{k}\}_{k=1}^{K}\big{)}\bigg{\{}\prod_{k=1}^{K}\sigma(u_{ik})\bigg{\}}^{z_{ik}}\bigg{\{}1-\prod_{k=1}^{K}\sigma(u_{ik})\bigg{\}}^{1-z_{ik}}.$$

  When $K$ is not large (e.g., $K<8$), we can compute the above (unnormalized) probability on $(z_{i1},\ldots,z_{iK})\in\{0,1\}^{K}$, and update the binary vector $(z_{i1},\ldots,z_{iK})$ simultaneously. When $K$ is large, one can update $z_{ik}$ separately by computing the full conditional probability of $z_{ik}=1$.

• -

  (Update of $\theta_{k}$)  For each $k$, the full conditional of $\theta_{k}$ is proportional to

  $\displaystyle\pi(\theta_{k})\prod_{i=1}^{n}f\big{(}x_{i}\mid\{z_{ij}\}_{j=1}^{K};\{\theta_{j}\}_{j=1}^{K}\big{)}\min_{1\leq k^{\prime}\leq K,k^{\prime}\neq k}g(\|\theta_{k}-\theta_{k^{\prime}}\|),$

  and the specific sampling algorithms depend on the distributional assumption for $x_{i}$. A typical strategy is to generate a proposal $\theta_{k}^{\ast}$ from a distribution proportional to $\pi(\theta_{k})\prod_{i=1}^{n}f\big{(}x_{i}\mid\{z_{ij}\}_{j=1}^{K};\{\theta_{j}\}_{j=1}^{K}\big{)}$ and accept it with probability

  $$\min_{1\leq k^{\prime}\leq K,k^{\prime}\neq k}g(\|\theta_{k}^{\ast}-\theta_{k^{\prime}}\|)/\min_{1\leq k^{\prime}\leq K,k^{\prime}\neq k}g(\|\theta_{k}^{\dagger}-\theta_{k^{\prime}}\|),$$

  for the current value $\theta_{k}^{\dagger}$

• -

  (Update of $\tau$)  The full conditional of $\tau$ is ${\rm Ga}(\tilde{a}_{\tau},\tilde{b}_{\tau})$, where $\tilde{a}_{\tau}=a_{\tau}+nK/2$ and $\tilde{b}_{\tau}=b_{\tau}+\frac{1}{2}\sum_{k=1}^{K}(u_{k}-\mu 1_{n})^{\top}Q(\psi)^{-1}(u_{k}-\mu 1_{n})/2$.

• -

  (Update of $\mu$)  The full conditional of $\mu$ is $N(A_{\mu}^{-1}B_{\mu},A_{\mu}^{-1})$, where $A_{\mu}=K\tau 1_{n}^{\top}Q(\psi)^{-1}1_{n}+S_{\mu}^{-1}$ and $B_{\mu}=\tau 1_{n}^{\top}Q(\psi)^{-1}\sum_{k=1}^{K}u_{k}+S_{\mu}^{-1}m_{\mu}$.

• -

  (Update of $\psi$)  The full conditional of $\psi$ is proportional to

  $$\pi(\psi)|Q(\psi)|^{-K/2}\exp\left\{-\frac{\tau}{2}\sum_{k=1}^{K}(u_{k}-\mu 1_{n})^{\top}Q(\psi)^{-1}(u_{k}-\mu 1_{n})\right\},$$

  for which we use a random-walk Metropolis-Hastings algorithm.