# A Massively Parallel SMC Sampler for Decision Trees +
## 4 Bayesian Decision Trees
### Sequential Monte Carlo  
The general idea behind SMC samplers is to perform sampling and resampling in sequence for \(K\) iterations, for \(k=0,1,\ldots,K-1\), each of which will generate and update \(N\) independent, weighted samples, \(\mathbf{t}_{k}^{i}\), for \(i=0,1,2,\ldots,N-1\). The samples are first initialized from the initial proposal and weighted as follows:  
\[\mathbf{t}_{0}^{i} \sim q_{0}(\cdot), \tag{4a}\] \[\mathbf{w}_{0}^{i} =\frac{\pi(\mathbf{t}_{0}^{i}|\mathbf{Y},\mathbf{x})}{q_{0}( \mathbf{t}_{0}^{i})}. \tag{4b}\]  
Then, during each SMC iteration, for \(k=1,2,\ldots,K-1\), the samples are updated from the proposal and weighted as follows:  
\[\mathbf{t}_{k}^{i} \sim q(\cdot|\mathbf{t}_{k-1}^{i}), \tag{5a}\] \[\mathbf{w}_{k}^{i} =\mathbf{w}_{k-1}^{i}\frac{\pi(\mathbf{t}_{k}^{i}|\mathbf{Y}, \mathbf{x})}{\pi(\mathbf{t}_{k-1}^{i}|\mathbf{Y},\mathbf{x})}\frac{q\Big{(} \mathbf{t}_{k-1}^{i}|\mathbf{t}_{k}^{i}\Big{)}}{q\Big{(}\mathbf{t}_{k}^{i}| \mathbf{t}_{k-1}^{i}\Big{)}}. \tag{5b}\]  
We note that Equation (5) is often referred to as the IS step. The weights, \(\mathbf{w}_{k}^{i}\), for \(i=0,1,\ldots,N-1\), are then normalized such that \(\sum_{i}\tilde{\mathbf{w}}_{k}^{i}=1\), i.e., 100% of the total probability space. The normalized weights are then computed as follows:  
\[\tilde{\mathbf{w}}_{k}^{i}=\frac{\mathbf{w}_{k}^{i}}{\sum_{j=0}^{N-1}\mathbf{ w}_{k}^{j}},\ \ \ \forall i=0,1,2,\ldots,N-1. \tag{6}\]  
Since the samples are generated from \(q(\cdot|\cdot)\), and not directly from the posterior distribution of interest, this sampling strategy may suffer from a numerical error, called degeneracy, which could make most of the weights equal to 0, leading to poor estimations of the state. This problem is commonly tackled by using resampling when the effective sample size (ESS) of the weights,  
\[N_{eff}=\frac{1}{\sum_{i=0}^{N-1}\big{(}\tilde{\mathbf{w}}_{k}^{i}\big{)}^{2}}, \tag{7}\]  
falls below a user-defined threshold, usually set at \(\frac{N}{2}\). Resampling overwrites the samples with low weights with copies of the samples with high weights. Several variants of resampling can be found in the literature [41]. In this work, we use multinomial resampling, which performs three steps. In the first step, we compute \(\mathbf{ncopies}\in\mathbb{N}_{0}^{N}\), where \(\mathbb{N}_{0}\) is the set of non-negative integers. The \(i\)-th element, \(\mathbf{ncopies}^{i}\), indicates how many copies of the \(i\)-th sample, \(\mathbf{t}_{k}^{i}\), are necessary to keep. To do this, multinomial resampling first computes \(\mathbf{cdf}\in\mathbb{R}^{N}\), the cumulative sum (or prefix sum; the terms are used interchangeably) of \(N\tilde{\mathbf{w}}\), as follows:  
\[\mathbf{cdf}^{i}=N\sum_{j=0}^{i-1}\tilde{\mathbf{w}}_{k}^{i},\ \ \forall i=0,1,2,\ldots,N-1. \tag{8}\]  
Then, each term \(\mathbf{ncopies}^{i}\) is computed by  
\[\mathbf{ncopies}^{i}=\left\lceil\mathbf{cdf}^{i}+N\tilde{\mathbf{w}}_{k}^{i}- u\right\rceil-\left\lceil\mathbf{cdf}^{i}-u\right\rceil, \tag{9}\]  
where \(u\sim\)Uniform\([0,1]\), and \(\left\lceil\cdot\right\rceil\) is the ceiling operator. From (8) and (9), it is relativelystraightforward to infer that  
\[0\leq\mathbf{ncopies}^{i} \leq N, \tag{10a}\] \[\sum_{i=0}^{N-1}\mathbf{ncopies}^{i} =N. \tag{10b}\]  
The second step is redistribution, which creates a new population of samples by duplicating each sample, \(\mathbf{t}_{k}^{i}\), a total of \(\mathbf{ncopies}^{i}\) times, meaning that those samples for which \(\mathbf{ncopies}^{i}=0\) will be deleted. Algorithm 3 illustrates a possible implementation of sequential redistribution (S-R). ```
1:t,ncopies,\(N\)
2:t_new
3:\(i\gets 0\)
4:for\(j\gets 0\); \(j\leq N\); \(j\gets j+1\)do
5:forcopy\(\gets 0\); \(copy<\mathbf{ncopies}^{j}\); \(copy\gets copy+1\)do
6:\(\mathbf{t}_{new}^{i}\leftarrow\mathbf{t}^{i}\)
7:\(i\gets i+1\)
8:endfor
9:endfor
```  
**Algorithm 3** Sequential Redistribution (S-R)  
In the third and final step of resampling, each weight, \(\mathbf{w}^{i}\), is reset to \(1/N\).