<div>
<p> </p>
<h1> Bayesian nonparametrics in NIMBLE: Density estimation </h1>
<h2> Overview </h2>
<p>NIMBLE is a hierarchical modeling package that uses nearly the same language for model specification as the popular MCMC packages WinBUGS, OpenBUGS and JAGS, while making the modeling language extensible — you can add distributions and functions — and also allowing customization of the algorithms used to estimate the parameters of the model.</p>
<p>Recently, we added support for Markov chain Monte Carlo (MCMC) inference for Bayesian nonparametric (BNP) mixture models to NIMBLE. In particular, starting with version 0.6-11, NIMBLE provides functionality for fitting models involving Dirichlet process priors using either the Chinese Restaurant Process (CRP) or a truncated stick-breaking (SB) representation of the Dirichlet process prior.</p>
<p>In this post we illustrate NIMBLE’s BNP capabilities by showing how to use nonparametric mixture models with different kernels for density estimation. In a later post, we will take a parametric generalized linear mixed model and show how to switch to a nonparametric representation of the random effects that avoids the assumption of normally-distributed random effects.</p>
<p>For more detailed information on NIMBLE and Bayesian nonparametrics in NIMBLE, see the <a rel="nofollow" href="https://r-nimble.org/documentation" target="_blank">NIMBLE User Manual</a>.</p>
<h2> Basic density estimation using Dirichlet Process Mixture models </h2>
<p>NIMBLE provides the machinery for nonparametric density estimation by means of Dirichlet process mixture (DPM) models (Ferguson, 1974; Lo, 1984; Escobar, 1994; Escobar and West, 1995). For an independent and identically distributed sample <img alt="y_1, ldots, y_n" title="y_1, ldots, y_n" src="http://s0.wp.com/latex.php?latex=y_1%2C+%5Cldots%2C+y_n&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">, the model takes the form</p>
<p><img alt="  y_i mid theta_i sim p(y_i mid theta_i), quadquad theta_i mid G sim G, quadquad G mid  alpha, H sim mbox{DP}(alpha, H), quadquad i=1,ldots, n .  " title="  y_i mid theta_i sim p(y_i mid theta_i), quadquad theta_i mid G sim G, quadquad G mid  alpha, H sim mbox{DP}(alpha, H), quadquad i=1,ldots, n .  " src="http://s0.wp.com/latex.php?latex=++y_i+%5Cmid+%5Ctheta_i+%5Csim+p%28y_i+%5Cmid+%5Ctheta_i%29%2C+%5Cquad%5Cquad+%5Ctheta_i+%5Cmid+G+%5Csim+G%2C+%5Cquad%5Cquad+G+%5Cmid++%5Calpha%2C+H+%5Csim+%5Cmbox%7BDP%7D%28%5Calpha%2C+H%29%2C+%5Cquad%5Cquad+i%3D1%2C%5Cldots%2C+n+.++&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"></p>
<p>The NIMBLE implementation of this model is flexible and allows for mixtures of arbitrary kernels, <img alt="p(y_i mid theta)" title="p(y_i mid theta)" src="http://s0.wp.com/latex.php?latex=p%28y_i+%5Cmid+%5Ctheta%29&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">, which can be either conjugate or non-conjugate to the (also arbitrary) base measure <img alt="H" title="H" src="http://s0.wp.com/latex.php?latex=H&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">. In the case of conjugate kernel / base measure pairs, NIMBLE is able to detect the presence of the conjugacy and use it to improve the performance of the sampler.</p>
<p>To illustrate these capabilities, we consider the estimation of the probability density function of the waiting time between eruptions of the Old Faithful volcano data set available in R. </p>
<div>
<div>
<div>
<pre><span>data</span><span>(faithful)</span>
</pre>
</div>
</div>
</div>
<p>The observations <img alt="y_1, ldots, y_n" title="y_1, ldots, y_n" src="http://s0.wp.com/latex.php?latex=y_1%2C+%5Cldots%2C+y_n&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> correspond to the second column of the dataframe, and <img alt="n = 272" title="n = 272" src="http://s0.wp.com/latex.php?latex=n+%3D+272&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">.</p>
<h2> Fitting a location-scale mixture of Gaussian distributions using the CRP representation </h2>
<h3> Model specification</h3>
<p>We first consider a location-scale Dirichlet process mixture of normal distributionss fitted to the transformed data <img alt="y_i^{*} = log (y_i)" title="y_i^{*} = log (y_i)" src="http://s0.wp.com/latex.php?latex=y_i%5E%7B%2A%7D+%3D+%5Clog+%28y_i%29&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">:</p>
<p><img alt="  y^{*}_i mid mu_i, sigma^2_i sim mbox{N}(mu_i, sigma^2_i), quad (mu_i, sigma^2_i) mid G sim G, quad G mid alpha, H sim mbox{DP}(alpha, H), quad i=1,ldots, n,  " title="  y^{*}_i mid mu_i, sigma^2_i sim mbox{N}(mu_i, sigma^2_i), quad (mu_i, sigma^2_i) mid G sim G, quad G mid alpha, H sim mbox{DP}(alpha, H), quad i=1,ldots, n,  " src="http://s0.wp.com/latex.php?latex=++y%5E%7B%2A%7D_i+%5Cmid+%5Cmu_i%2C+%5Csigma%5E2_i+%5Csim+%5Cmbox%7BN%7D%28%5Cmu_i%2C+%5Csigma%5E2_i%29%2C+%5Cquad+%28%5Cmu_i%2C+%5Csigma%5E2_i%29+%5Cmid+G+%5Csim+G%2C+%5Cquad+G+%5Cmid+%5Calpha%2C+H+%5Csim+%5Cmbox%7BDP%7D%28%5Calpha%2C+H%29%2C+%5Cquad+i%3D1%2C%5Cldots%2C+n%2C++&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"></p>
<p>where <img alt="H" title="H" src="http://s0.wp.com/latex.php?latex=H&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> corresponds to a normal-inverse-gamma distribution. This model can be interpreted as providing a Bayesian version of kernel density estimation for <img alt="y^{*}_i" title="y^{*}_i" src="http://s0.wp.com/latex.php?latex=y%5E%7B%2A%7D_i&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> using Gaussian kernels and <em>adaptive bandwidths</em>. On the original scale of the data, this translates into an adaptive log-Gaussian kernel density estimate.</p>
<p>Introducing auxiliary variables <img alt="xi_1, ldots, xi_n" title="xi_1, ldots, xi_n" src="http://s0.wp.com/latex.php?latex=%5Cxi_1%2C+%5Cldots%2C+%5Cxi_n&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> that indicate which component of the mixture generates each observation, and integrating over the random measure <img alt="G" title="G" src="http://s0.wp.com/latex.php?latex=G&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">, we obtain the CRP representation of the model (Blackwell and MacQueen, 1973):</p>
<p><img alt="  y_i^{*} mid { tilde{mu}_k }, { tilde{sigma}_k^{2} } sim mbox{N}left( tilde{mu}_{xi_i}, tilde{sigma}^2_{xi_i} right), quadquad xi mid alpha sim mbox{CRP}(alpha), quadquad (tilde{mu}_k, tilde{sigma}_k^2) mid H sim H, quadquad i=1,ldots, n ,  " title="  y_i^{*} mid { tilde{mu}_k }, { tilde{sigma}_k^{2} } sim mbox{N}left( tilde{mu}_{xi_i}, tilde{sigma}^2_{xi_i} right), quadquad xi mid alpha sim mbox{CRP}(alpha), quadquad (tilde{mu}_k, tilde{sigma}_k^2) mid H sim H, quadquad i=1,ldots, n ,  " src="http://s0.wp.com/latex.php?latex=++y_i%5E%7B%2A%7D+%5Cmid+%5C%7B+%5Ctilde%7B%5Cmu%7D_k+%5C%7D%2C+%5C%7B+%5Ctilde%7B%5Csigma%7D_k%5E%7B2%7D+%5C%7D+%5Csim+%5Cmbox%7BN%7D%5Cleft%28+%5Ctilde%7B%5Cmu%7D_%7B%5Cxi_i%7D%2C+%5Ctilde%7B%5Csigma%7D%5E2_%7B%5Cxi_i%7D+%5Cright%29%2C+%5Cquad%5Cquad+%5Cxi+%5Cmid+%5Calpha+%5Csim+%5Cmbox%7BCRP%7D%28%5Calpha%29%2C+%5Cquad%5Cquad+%28%5Ctilde%7B%5Cmu%7D_k%2C+%5Ctilde%7B%5Csigma%7D_k%5E2%29+%5Cmid+H+%5Csim+H%2C+%5Cquad%5Cquad+i%3D1%2C%5Cldots%2C+n+%2C++&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"></p>
<p>where</p>
<p><img alt="  p(xi mid alpha) = frac{Gamma(alpha)}{Gamma(alpha + n)} alpha^{K(xi)} prod_k  Gammaleft(m_k(xi)right),  " title="  p(xi mid alpha) = frac{Gamma(alpha)}{Gamma(alpha + n)} alpha^{K(xi)} prod_k  Gammaleft(m_k(xi)right),  " src="http://s0.wp.com/latex.php?latex=++p%28%5Cxi+%5Cmid+%5Calpha%29+%3D+%5Cfrac%7B%5CGamma%28%5Calpha%29%7D%7B%5CGamma%28%5Calpha+%2B+n%29%7D+%5Calpha%5E%7BK%28%5Cxi%29%7D+%5Cprod_k++%5CGamma%5Cleft%28m_k%28%5Cxi%29%5Cright%29%2C++&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"></p>
<p><img alt="K(xi) le n" title="K(xi) le n" src="http://s0.wp.com/latex.php?latex=K%28%5Cxi%29+%5Cle+n&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> is the number of unique values in the vector <img alt="xi" title="xi" src="http://s0.wp.com/latex.php?latex=%5Cxi&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">, and <img alt="m_k(xi)" title="m_k(xi)" src="http://s0.wp.com/latex.php?latex=m_k%28%5Cxi%29&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> is the number of times the <img alt="k" title="k" src="http://s0.wp.com/latex.php?latex=k&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">-th unique value appears in <img alt="xi" title="xi" src="http://s0.wp.com/latex.php?latex=%5Cxi&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">. This specification makes it clear that each observation belongs to any of <em>at most</em> <img alt="n" title="n" src="http://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> normally distributed clusters, and that the CRP distribution corresponds to the prior distribution on the partition structure. </p>
<p>NIMBLE’s specification of this model is given by</p>
<div>
<div>
<div>
<pre><span>code</span> <span>&lt;-</span> <span>nimbleCode</span><span>({</span>
  <span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>n) {</span>
    <span>y[i]</span> <span>~</span> <span>dnorm</span><span>(mu[i],</span> <span>var</span> <span>= s2[i])</span>
    <span>mu[i]</span> <span>&lt;-</span> <span>muTilde[xi[i]]</span>
    <span>s2[i]</span> <span>&lt;-</span> <span>s2Tilde[xi[i]]</span>
  <span>}</span>
  <span>xi[</span><span>1</span><span>:</span><span>n]</span> <span>~</span> <span>dCRP</span><span>(alpha,</span> <span>size</span> <span>= n)</span>
  <span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>n) {</span>
    <span>muTilde[i]</span> <span>~</span> <span>dnorm</span><span>(</span><span>0</span><span>,</span> <span>var</span> <span>= s2Tilde[i])</span>
    <span>s2Tilde[i]</span> <span>~</span> <span>dinvgamma</span><span>(</span><span>2</span><span>,</span> <span>1</span><span>)</span>
  <span>}</span>
  <span>alpha</span> <span>~</span> <span>dgamma</span><span>(</span><span>1</span><span>,</span> <span>1</span><span>)</span>
<span>})</span>
</pre>
</div>
</div>
</div>
<p>Note that in the model code the length of the parameter vectors <tt>muTilde</tt> and <tt>s2Tilde</tt> has been set to <img alt="n" title="n" src="http://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">. We do this because the current implementation of NIMBLE requires that the length of vector of parameters be set in advance and does not allow for their number to change between iterations. Hence, if we are to ensure that the algorithm always performs as intended we need to work with the worst case scenario, i.e., the case where there are as many components as observations. While this ensures that the algorithm always works as intended, it is also somewhat inefficient, both in terms of memory requirements (when <img alt="n" title="n" src="http://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> is large a large number of unoccupied components need to be maintained) and in terms of computational burden (a large number of parameters that are not required for posterior inference need to be updated at every iteration). When we use a mixture of gamma distributions below, we will show a computational shortcut that improves the efficiency.</p>
<p>Note also that the value of <img alt="alpha" title="alpha" src="http://s0.wp.com/latex.php?latex=%5Calpha&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> controls the number of components we expect a priori, with larger values of <img alt="alpha" title="alpha" src="http://s0.wp.com/latex.php?latex=%5Calpha&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> corresponding to a larger number of components occupied by the data. Hence, by assigning a prior to <img alt="alpha" title="alpha" src="http://s0.wp.com/latex.php?latex=%5Calpha&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> we add flexibility to the model specification. The particular choice of a Gamma prior allows NIMBLE to use a data-augmentation scheme to efficiently sample from the corresponding full conditional distribution. Alternative prior specifications for <img alt="alpha" title="alpha" src="http://s0.wp.com/latex.php?latex=%5Calpha&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> are possible, in which case the default sampler for this parameter is an adaptive random-walk Metropolis-Hastings algorithm.</p>
<h3> Running the MCMC algorithm</h3>
<p>The following code sets up the data and constants, initializes the parameters, defines the model object, and builds and runs the MCMC algorithm. Because the specification is in terms of a Chinese restaurant process, the default sampler selected by NIMBLE is a collapsed Gibbs sampler (Neal, 2000).</p>
<div>
<div>
<div>
<pre><span>set.seed</span><span>(</span><span>1</span><span>)</span>
<span># Model Data</span>
<span>lFaithful</span> <span>&lt;-</span> <span>log</span><span>(faithful</span><span>$</span><span>waiting)</span>
<span>standlFaithful</span> <span>&lt;-</span> <span>(lFaithful</span> <span>-</span> <span>mean</span><span>(lFaithful))</span> <span>/</span> <span>sd</span><span>(lFaithful)</span>
<span>data</span> <span>&lt;-</span> <span>list</span><span>(</span><span>y</span> <span>= standlFaithful)</span>
<span># Model Constants</span>
<span>consts</span> <span>&lt;-</span> <span>list</span><span>(</span><span>n</span> <span>=</span> <span>length</span><span>(standlFaithful))</span>
<span># Parameter initialization</span>
<span>inits</span> <span>&lt;-</span> <span>list</span><span>(</span><span>xi</span> <span>=</span> <span>sample</span><span>(</span><span>1</span><span>:</span><span>10</span><span>,</span> <span>size</span><span>=consts</span><span>$</span><span>n,</span> <span>replace</span><span>=</span><span>TRUE</span><span>),</span>
              <span>muTilde</span> <span>=</span> <span>rnorm</span><span>(consts</span><span>$</span><span>n,</span> <span>0</span><span>,</span> <span>sd</span> <span>=</span> <span>sqrt</span><span>(</span><span>10</span><span>)),</span>
              <span>s2Tilde</span> <span>=</span> <span>rinvgamma</span><span>(consts</span><span>$</span><span>n,</span> <span>2</span><span>,</span> <span>1</span><span>),</span>
              <span>alpha</span> <span>=</span> <span>1</span><span>)</span>
<span># Model creation and compilation</span>
<span>rModel</span> <span>&lt;-</span> <span>nimbleModel</span><span>(code,</span> <span>data</span> <span>= data,</span> <span>inits</span> <span>= inits,</span> <span>constants</span> <span>= consts)</span>
</pre>
</div>
<div>
<pre>## defining model...
</pre>
</div>
<div>
<pre>## building model...
</pre>
</div>
<div>
<pre>## setting data and initial values...
</pre>
</div>
<div>
<pre>## running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
## checking model sizes and dimensions...
## model building finished.
</pre>
</div>
<div>
<pre><span>cModel</span> <span>&lt;-</span> <span>compileNimble</span><span>(rModel)</span>
</pre>
</div>
<div>
<pre>## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
## compilation finished.
</pre>
</div>
<div>
<pre><span># MCMC configuration, creation, and compilation</span>
<span>conf</span> <span>&lt;-</span> <span>configureMCMC</span><span>(rModel,</span> <span>monitors</span> <span>=</span> <span>c</span><span>(</span><span>&quot;xi&quot;</span><span>,</span> <span>&quot;muTilde&quot;</span><span>,</span> <span>&quot;s2Tilde&quot;</span><span>,</span> <span>&quot;alpha&quot;</span><span>))</span>
<span>mcmc</span> <span>&lt;-</span> <span>buildMCMC</span><span>(conf)</span>
<span>cmcmc</span> <span>&lt;-</span> <span>compileNimble</span><span>(mcmc,</span> <span>project</span> <span>= rModel)</span>
</pre>
</div>
<div>
<pre>## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
## compilation finished.
</pre>
</div>
<div>
<pre><span>samples</span> <span>&lt;-</span> <span>runMCMC</span><span>(cmcmc,</span> <span>niter</span> <span>=</span> <span>7000</span><span>,</span> <span>nburnin</span> <span>=</span> <span>2000</span><span>,</span> <span>setSeed</span> <span>=</span> <span>TRUE</span><span>)</span>
</pre>
</div>
<div>
<pre>## running chain 1...
</pre>
</div>
<div>
<pre>## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
</pre>
</div>
</div>
</div>
<p>We can extract the samples from the posterior distributions of the parameters and create trace plots, histograms, and any other summary of interest. For example, for the concentration parameter <img alt="alpha" title="alpha" src="http://s0.wp.com/latex.php?latex=%5Calpha&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> we have:</p>
<div>
<div>
<div>
<pre><span># Trace plot for the concentration parameter</span>
<span>ts.plot</span><span>(samples[ ,</span> <span>&quot;alpha&quot;</span><span>],</span> <span>xlab</span> <span>=</span> <span>&quot;iteration&quot;</span><span>,</span> <span>ylab</span> <span>=</span> <span>expression</span><span>(alpha))</span>
</pre>
</div>
</div>
<div><img src="https://i1.wp.com/r-nimble.org/wp-content/uploads/2018/12/MCMC-normaldensity-somegraphs-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
<div>
<div>
<pre><span># Posterior histogram</span>
<span>hist</span><span>(samples[ ,</span> <span>&quot;alpha&quot;</span><span>],</span> <span>xlab</span> <span>=</span> <span>expression</span><span>(alpha),</span> <span>main</span> <span>=</span> <span>&quot;&quot;</span><span>,</span> <span>ylab</span> <span>=</span> <span>&quot;Frequency&quot;</span><span>)</span>
</pre>
</div>
</div>
<div><img src="https://i2.wp.com/r-nimble.org/wp-content/uploads/2018/12/MCMC-normaldensity-somegraphs-2.png?w=456&ssl=1" data-recalc-dims="1"></div>
<div>
<div>
<pre><span>quantile</span><span>(samples[ ,</span> <span>&quot;alpha&quot;</span><span>],</span> <span>c</span><span>(</span><span>0.5</span><span>,</span> <span>0.025</span><span>,</span> <span>0.975</span><span>))</span>
</pre>
</div>
<div>
<pre>##        50%       2.5%      97.5% 
## 0.42365754 0.06021512 1.52299639
</pre>
</div>
</div>
</div>
<p>Under this model, the posterior predictive distribution for a new observation <img alt="tilde{y}" title="tilde{y}" src="http://s0.wp.com/latex.php?latex=%5Ctilde%7By%7D&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">, <img alt="p(tilde{y} mid y_1, ldots, y_n)" title="p(tilde{y} mid y_1, ldots, y_n)" src="http://s0.wp.com/latex.php?latex=p%28%5Ctilde%7By%7D+%5Cmid+y_1%2C+%5Cldots%2C+y_n%29&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">, is the optimal density estimator (under squared error loss). Samples for this estimator can be easily computed from the samples generated by our MCMC:</p>
<div>
<div>
<div>
<pre><span># posterior samples of the concentration parameter</span>
<span>alphaSamples</span> <span>&lt;-</span> <span>samples[ ,</span> <span>&quot;alpha&quot;</span><span>]</span>
<span># posterior samples of the cluster means</span>
<span>muTildeSamples</span> <span>&lt;-</span> <span>samples[ ,</span> <span>grep</span><span>(</span><span>'muTilde'</span><span>,</span> <span>colnames</span><span>(samples))]</span>
<span># posterior samples of the cluster variances</span>
<span>s2TildeSamples</span> <span>&lt;-</span> <span>samples[ ,</span> <span>grep</span><span>(</span><span>'s2Tilde'</span><span>,</span> <span>colnames</span><span>(samples))]</span>
<span># posterior samples of the cluster memberships</span>
<span>xiSamples</span> <span>&lt;-</span> <span>samples [ ,</span> <span>grep</span><span>(</span><span>'xi'</span><span>,</span> <span>colnames</span><span>(samples))]</span>

<span>standlGrid</span> <span>&lt;-</span> <span>seq</span><span>(</span><span>-</span><span>2.5</span><span>,</span> <span>2.5</span><span>,</span> <span>len</span> <span>=</span> <span>200</span><span>)</span> <span># standardized grid on log scale</span>

<span>densitySamplesStandl</span> <span>&lt;-</span> <span>matrix</span><span>(</span><span>0</span><span>,</span> <span>ncol</span> <span>=</span> <span>length</span><span>(standlGrid),</span> <span>nrow</span> <span>=</span> <span>nrow</span><span>(samples))</span>
<span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>nrow</span><span>(samples)){</span>
  <span>k</span> <span>&lt;-</span> <span>unique</span><span>(xiSamples[i, ])</span>
  <span>kNew</span> <span>&lt;-</span> <span>max</span><span>(k)</span> <span>+</span> <span>1</span>
  <span>mk</span> <span>&lt;-</span> <span>c</span><span>()</span>
  <span>li</span> <span>&lt;-</span> <span>1</span>
  <span>for</span><span>(l</span> <span>in</span> <span>1</span><span>:</span><span>length</span><span>(k)) {</span>
    <span>mk[li]</span> <span>&lt;-</span> <span>sum</span><span>(xiSamples[i, ]</span> <span>==</span> <span>k[li])</span>
    <span>li</span> <span>&lt;-</span> <span>li</span> <span>+</span> <span>1</span>
  <span>}</span>
  <span>alpha</span> <span>&lt;-</span> <span>alphaSamples[i]</span>

  <span>muK</span> <span>&lt;-</span>  <span>muTildeSamples[i, k]</span>
  <span>s2K</span> <span>&lt;-</span>  <span>s2TildeSamples[i, k]</span>
  <span>muKnew</span> <span>&lt;-</span>  <span>muTildeSamples[i, kNew]</span>
  <span>s2Knew</span> <span>&lt;-</span>  <span>s2TildeSamples[i, kNew]</span>

  <span>densitySamplesStandl[i, ]</span> <span>&lt;-</span> <span>sapply</span><span>(standlGrid,</span>
                <span>function</span><span>(</span><span>x</span><span>)(</span><span>sum</span><span>(mk</span> <span>*</span> <span>dnorm</span><span>(x, muK,</span> <span>sqrt</span><span>(s2K)))</span> <span>+</span>
                <span>alpha</span> <span>*</span> <span>dnorm</span><span>(x, muKnew,</span> <span>sqrt</span><span>(s2Knew)) )</span><span>/</span><span>(alpha</span><span>+</span><span>consts</span><span>$</span><span>n))</span>
<span>}</span>

<span>hist</span><span>(data</span><span>$</span><span>y,</span> <span>freq</span> <span>=</span> <span>FALSE</span><span>,</span> <span>xlim</span> <span>=</span> <span>c</span><span>(</span><span>-</span><span>2.5</span><span>,</span> <span>2.5</span><span>),</span> <span>ylim</span> <span>=</span> <span>c</span><span>(</span><span>0</span><span>,</span><span>0.75</span><span>),</span> <span>main</span> <span>=</span> <span>&quot;&quot;</span><span>,</span>
     <span>xlab</span> <span>=</span> <span>&quot;Waiting times on standardized log scale&quot;</span><span>)</span>
<span>## pointwise estimate of the density for standardized log grid</span>
<span>lines</span><span>(standlGrid,</span> <span>apply</span><span>(densitySamplesStandl,</span> <span>2</span><span>, mean),</span> <span>lwd</span> <span>=</span> <span>2</span><span>,</span> <span>col</span> <span>=</span> <span>'black'</span><span>)</span>
<span>lines</span><span>(standlGrid,</span> <span>apply</span><span>(densitySamplesStandl,</span> <span>2</span><span>, quantile,</span> <span>0.025</span><span>),</span> <span>lty</span> <span>=</span> <span>2</span><span>,</span> <span>col</span> <span>=</span> <span>'black'</span><span>)</span>
<span>lines</span><span>(standlGrid,</span> <span>apply</span><span>(densitySamplesStandl,</span> <span>2</span><span>, quantile,</span> <span>0.975</span><span>),</span> <span>lty</span> <span>=</span> <span>2</span><span>,</span> <span>col</span> <span>=</span> <span>'black'</span><span>)</span>
</pre>
</div>
</div>
<div><img src="https://i2.wp.com/r-nimble.org/wp-content/uploads/2018/12/MCMC-normaldensity-predictive-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
</div>
<p>Recall, however, that this is the density estimate for the logarithm of the waiting time. To obtain the density on the original scale we need to apply the appropriate transformation to the kernel.</p>
<div>
<div>
<div>
<pre><span>lgrid</span> <span>&lt;-</span> <span>standlGrid</span><span>*</span><span>sd</span><span>(lFaithful)</span> <span>+</span> <span>mean</span><span>(lFaithful)</span> <span># grid on log scale</span>
<span>densitySamplesl</span> <span>&lt;-</span> <span>densitySamplesStandl</span> <span>/</span> <span>sd</span><span>(lFaithful)</span> <span># density samples for grid on log scale</span>

<span>hist</span><span>(faithful</span><span>$</span><span>waiting,</span> <span>freq</span> <span>=</span> <span>FALSE</span><span>,</span> <span>xlim</span> <span>=</span> <span>c</span><span>(</span><span>40</span><span>,</span> <span>100</span><span>),</span> <span>ylim</span><span>=</span><span>c</span><span>(</span><span>0</span><span>,</span> <span>0.05</span><span>),</span>
     <span>main</span> <span>=</span> <span>&quot;&quot;</span><span>,</span> <span>xlab</span> <span>=</span> <span>&quot;Waiting times&quot;</span><span>)</span>
<span>lines</span><span>(</span><span>exp</span><span>(lgrid),</span> <span>apply</span><span>(densitySamplesl,</span> <span>2</span><span>, mean)</span><span>/</span><span>exp</span><span>(lgrid),</span> <span>lwd</span> <span>=</span> <span>2</span><span>,</span> <span>col</span> <span>=</span> <span>'black'</span><span>)</span>
<span>lines</span><span>(</span><span>exp</span><span>(lgrid),</span> <span>apply</span><span>(densitySamplesl,</span> <span>2</span><span>, quantile,</span> <span>0.025</span><span>)</span><span>/</span><span>exp</span><span>(lgrid),</span> <span>lty</span> <span>=</span> <span>2</span><span>,</span>
      <span>col</span> <span>=</span> <span>'black'</span><span>)</span>
<span>lines</span><span>(</span><span>exp</span><span>(lgrid),</span> <span>apply</span><span>(densitySamplesl,</span> <span>2</span><span>, quantile,</span> <span>0.975</span><span>)</span><span>/</span><span>exp</span><span>(lgrid),</span> <span>lty</span> <span>=</span> <span>2</span><span>,</span>
      <span>col</span> <span>=</span> <span>'black'</span><span>)</span>
</pre>
</div>
</div>
<div><img src="https://i1.wp.com/r-nimble.org/wp-content/uploads/2018/12/MCMC-normaldensity-predictive-originalScale-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
</div>
<p>In either case, there is clear evidence that the data has two components for the waiting times.</p>
<h3> Generating samples from the mixing distribution</h3>
<p>While samples from the posterior distribution of linear functionals of the mixing distribution <img alt="G" title="G" src="http://s0.wp.com/latex.php?latex=G&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> (such as the predictive distribution above) can be computed directly from the realizations of the collapsed sampler, inference for non-linear functionals of <img alt="G" title="G" src="http://s0.wp.com/latex.php?latex=G&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> requires that we first generate samples from the mixing distribution. In NIMBLE we can get posterior samples from the random measure <img alt="G" title="G" src="http://s0.wp.com/latex.php?latex=G&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">, using the <tt>getSamplesDPmeasure</tt> function. Note that, in order to get posterior samples from <img alt="G" title="G" src="http://s0.wp.com/latex.php?latex=G&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">, we need to monitor all the random variables involved in its computations, i.e., the membership variable, <tt>xi</tt>, the cluster parameters, <tt>muTilde</tt> and <tt>s2Tilde</tt>, and the concentration parameter, <tt>alpha</tt>.</p>
<p>The following code generates posterior samples from the random measure <img alt="G" title="G" src="http://s0.wp.com/latex.php?latex=G&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">. The <tt>cMCMC</tt> object includes the model and posterior samples from the parameters. The <tt>getSamplesDPmeasure</tt> function estimates a truncation level of <img alt="G" title="G" src="http://s0.wp.com/latex.php?latex=G&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">, namely <tt>truncG</tt>. The posterior samples are in a matrix with <img alt="(mbox{truncG} cdot (p+1))" title="(mbox{truncG} cdot (p+1))" src="http://s0.wp.com/latex.php?latex=%28%5Cmbox%7BtruncG%7D+%5Ccdot+%28p%2B1%29%29&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> columns, where <img alt="p" title="p" src="http://s0.wp.com/latex.php?latex=p&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> is the dimension of the vector of parameters with distribution <img alt="G" title="G" src="http://s0.wp.com/latex.php?latex=G&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> (in this example <img alt="p=2" title="p=2" src="http://s0.wp.com/latex.php?latex=p%3D2&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">).</p>
<div>
<div>
<div>
<pre><span>outputG</span> <span>&lt;-</span> <span>getSamplesDPmeasure</span><span>(cmcmc)</span>
</pre>
</div>
<div>
<pre>## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
</pre>
</div>
<div>
<pre>## compilation finished.
</pre>
</div>
<div>
<pre>## sampleDPmeasure: Approximating the random measure by a finite stick-breaking representation with and error smaller than 1e-10, leads to a truncation level of 33.
</pre>
</div>
<div>
<pre><span>if</span><span>(</span><span>packageVersion</span><span>(</span><span>'nimble'</span><span>)</span> <span>&lt;=</span> <span>'0.6-12'</span><span>)</span>
  <span>samplesG</span> <span>&lt;-</span> <span>outputG</span> <span>else</span> <span>samplesG</span> <span>&lt;-</span> <span>outputG</span><span>$</span><span>samples</span>
</pre>
</div>
</div>
</div>
<p>The following code computes posterior samples of <img alt="P(tilde{y} &gt; 70)" title="P(tilde{y} &gt; 70)" src="http://s0.wp.com/latex.php?latex=P%28%5Ctilde%7By%7D+%3E+70%29&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> using the posterior samples from the random measure <img alt="G" title="G" src="http://s0.wp.com/latex.php?latex=G&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">. Note that these samples are computed based on the transformed model and a value larger than 70 corresponds to a value larger than 0.03557236 on the above defined grid.</p>
<div>
<div>
<div>
<pre><span>if</span><span>(</span><span>packageVersion</span><span>(</span><span>'nimble'</span><span>)</span> <span>>=</span> <span>'0.6.13'</span><span>)</span>
     <span>truncG</span> <span>&lt;-</span> <span>outputG</span><span>$</span><span>trunc</span> <span># truncation level for G</span>

<span>weightIndex</span> <span>&lt;-</span> <span>grep</span><span>(</span><span>'weight'</span><span>,</span> <span>colnames</span><span>(samplesG))</span>
<span>muTildeIndex</span> <span>&lt;-</span> <span>grep</span><span>(</span><span>'muTilde'</span><span>,</span> <span>colnames</span><span>(samplesG))</span>
<span>s2TildeIndex</span> <span>&lt;-</span> <span>grep</span><span>(</span><span>'s2Tilde'</span><span>,</span> <span>colnames</span><span>(samplesG))</span>

<span>probY70</span> <span>&lt;-</span> <span>rep</span><span>(</span><span>0</span><span>,</span> <span>nrow</span><span>(samples))</span>  <span># posterior samples of P(y.tilde > 70)</span>
<span>for</span><span>(i</span> <span>in</span> <span>seq_len</span><span>(</span><span>nrow</span><span>(samples))) {</span>
  <span>probY70[i]</span> <span>&lt;-</span> <span>sum</span><span>(samplesG[i, weightIndex]</span> <span>*</span>
                    <span>pnorm</span><span>(</span><span>0.03557236</span><span>,</span> <span>mean</span> <span>= samplesG[i, muTildeIndex],</span>
                      <span>sd</span> <span>=</span> <span>sqrt</span><span>(samplesG[i, s2TildeIndex]),</span> <span>lower.tail</span> <span>=</span> <span>FALSE</span><span>))</span>
<span>}</span>

<span>hist</span><span>(probY70,</span>  <span>xlab</span> <span>=</span> <span>&quot;Probability&quot;</span><span>,</span> <span>ylab</span> <span>=</span> <span>&quot;P(yTilde > 70 | data)&quot;</span> <span>,</span> <span>main</span> <span>=</span> <span>&quot;&quot;</span> <span>)</span>
</pre>
</div>
</div>
<div><img src="https://i0.wp.com/r-nimble.org/wp-content/uploads/2018/12/G-normaldensity-probYtildeLarger70-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
</div>
<h2> Fitting a mixture of gamma distributions using the CRP representation</h2>
<p>NIMBLE is not restricted to using Gaussian kernels in DPM models. In the case of the Old Faithful data, an alternative to the mixture of Gaussian kernels on the logarithmic scale that we presented in the previous section is a (scale-and-shape) mixture of Gamma distributions on the <em>original</em> scale of the data. </p>
<h3> Model specification</h3>
<p>In this case, the model takes the form</p>
<p><img alt="  y_i mid { tilde{beta}_k }, { tilde{lambda}_k } sim mbox{Gamma}left( tilde{beta}_{xi_i}, tilde{lambda}_{xi_i} right), quadquad xi mid alpha sim mbox{CRP}(alpha), quadquad (tilde{beta}_k, tilde{lambda}_k) mid H sim H ,  " title="  y_i mid { tilde{beta}_k }, { tilde{lambda}_k } sim mbox{Gamma}left( tilde{beta}_{xi_i}, tilde{lambda}_{xi_i} right), quadquad xi mid alpha sim mbox{CRP}(alpha), quadquad (tilde{beta}_k, tilde{lambda}_k) mid H sim H ,  " src="http://s0.wp.com/latex.php?latex=++y_i+%5Cmid+%5C%7B+%5Ctilde%7B%5Cbeta%7D_k+%5C%7D%2C+%5C%7B+%5Ctilde%7B%5Clambda%7D_k+%5C%7D+%5Csim+%5Cmbox%7BGamma%7D%5Cleft%28+%5Ctilde%7B%5Cbeta%7D_%7B%5Cxi_i%7D%2C+%5Ctilde%7B%5Clambda%7D_%7B%5Cxi_i%7D+%5Cright%29%2C+%5Cquad%5Cquad+%5Cxi+%5Cmid+%5Calpha+%5Csim+%5Cmbox%7BCRP%7D%28%5Calpha%29%2C+%5Cquad%5Cquad+%28%5Ctilde%7B%5Cbeta%7D_k%2C+%5Ctilde%7B%5Clambda%7D_k%29+%5Cmid+H+%5Csim+H+%2C++&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"></p>
<p>where <img alt="H" title="H" src="http://s0.wp.com/latex.php?latex=H&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> corresponds to the product of two independent Gamma distributions. The following code provides the NIMBLE specification for the model:</p>
<div>
<div>
<div>
<pre><span>code</span> <span>&lt;-</span> <span>nimbleCode</span><span>({</span>
  <span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>n) {</span>
    <span>y[i]</span> <span>~</span> <span>dgamma</span><span>(</span><span>shape</span> <span>= beta[i],</span> <span>scale</span> <span>= lambda[i])</span>
    <span>beta[i]</span> <span>&lt;-</span> <span>betaTilde[xi[i]]</span>
    <span>lambda[i]</span> <span>&lt;-</span> <span>lambdaTilde[xi[i]]</span>
  <span>}</span>
  <span>xi[</span><span>1</span><span>:</span><span>n]</span> <span>~</span> <span>dCRP</span><span>(alpha,</span> <span>size</span> <span>= n)</span>
  <span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>50</span><span>) {</span> <span># only 50 cluster parameters</span>
    <span>betaTilde[i]</span> <span>~</span> <span>dgamma</span><span>(</span><span>shape</span> <span>=</span> <span>71</span><span>,</span> <span>scale</span> <span>=</span> <span>2</span><span>)</span>
    <span>lambdaTilde[i]</span> <span>~</span> <span>dgamma</span><span>(</span><span>shape</span> <span>=</span> <span>2</span><span>,</span> <span>scale</span> <span>=</span> <span>2</span><span>)</span>
  <span>}</span>
  <span>alpha</span> <span>~</span> <span>dgamma</span><span>(</span><span>1</span><span>,</span> <span>1</span><span>)</span>
<span>})</span>
</pre>
</div>
</div>
</div>
<p>Note that in this case the vectors <tt>betaTilde</tt> and <tt>lambdaTilde</tt> have length <img alt="50 ll n = 272" title="50 ll n = 272" src="http://s0.wp.com/latex.php?latex=50+%5Cll+n+%3D+272&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">. This is done to reduce the computational and storage burdens associated with the sampling algorithm. You could think about this approach as truncating the process, except that it can be thought of as an *exact* truncation. Indeed, under the CRP representation, using parameter vector(s) with a length that is shorter than the number of observations in the sample will lead to a proper algorithm as long as the number of components instatiated by the sampler is strictly lower than the length of the parameter vector(s) for every iteration of the sampler.</p>
<h3> Running the MCMC algorithm</h3>
<p>The following code sets up the model data and constants, initializes the parameters, defines the model object, and builds and runs the MCMC algorithm for the mixture of Gamma distributions. Note that, when building the MCMC, a warning message about the number of cluster parameters is generated. This is because the lengths of <tt>betaTilde</tt> and <tt>lambdaTilde</tt> are smaller than <img alt="n" title="n" src="http://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">. Also, note that no error message is generated during execution, which indicates that the number of clusters required never exceeded the maximum of 50.</p>
<div>
<div>
<div>
<pre><span>data</span> <span>&lt;-</span> <span>list</span><span>(</span><span>y</span> <span>= faithful</span><span>$</span><span>waiting)</span>
<span>set.seed</span><span>(</span><span>1</span><span>)</span>
<span>inits</span> <span>&lt;-</span> <span>list</span><span>(</span><span>xi</span> <span>=</span> <span>sample</span><span>(</span><span>1</span><span>:</span><span>10</span><span>,</span> <span>size</span><span>=consts</span><span>$</span><span>n,</span> <span>replace</span><span>=</span><span>TRUE</span><span>),</span>
              <span>betaTilde</span> <span>=</span> <span>rgamma</span><span>(</span><span>50</span><span>,</span> <span>shape</span> <span>=</span> <span>71</span><span>,</span> <span>scale</span> <span>=</span> <span>2</span><span>),</span>
              <span>lambdaTilde</span> <span>=</span> <span>rgamma</span><span>(</span><span>50</span><span>,</span> <span>shape</span> <span>=</span> <span>2</span><span>,</span> <span>scale</span> <span>=</span> <span>2</span><span>),</span>
              <span>alpha</span> <span>=</span> <span>1</span><span>)</span>
<span>rModel</span> <span>&lt;-</span> <span>nimbleModel</span><span>(code,</span> <span>data</span> <span>= data,</span> <span>inits</span> <span>= inits,</span> <span>constants</span> <span>= consts)</span>
</pre>
</div>
<div>
<pre>## defining model...
</pre>
</div>
<div>
<pre>## building model...
</pre>
</div>
<div>
<pre>## setting data and initial values...
</pre>
</div>
<div>
<pre>## running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
## checking model sizes and dimensions...
## model building finished.
</pre>
</div>
<div>
<pre><span>cModel</span> <span>&lt;-</span> <span>compileNimble</span><span>(rModel)</span>
</pre>
</div>
<div>
<pre>## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
## compilation finished.
</pre>
</div>
<div>
<pre><span>conf</span> <span>&lt;-</span> <span>configureMCMC</span><span>(rModel,</span> <span>monitors</span> <span>=</span> <span>c</span><span>(</span><span>&quot;xi&quot;</span><span>,</span> <span>&quot;betaTilde&quot;</span><span>,</span> <span>&quot;lambdaTilde&quot;</span><span>,</span> <span>&quot;alpha&quot;</span><span>))</span>
<span>mcmc</span> <span>&lt;-</span> <span>buildMCMC</span><span>(conf)</span>
</pre>
</div>
<div>
<pre>## Warning in samplerFunction(model = model, mvSaved = mvSaved, target = target, : sampler_CRP: The number of cluster parameters is less than the number of potential clusters. The MCMC is not strictly valid if ever it proposes more components than cluster parameters exist; NIMBLE will warn you if this occurs.
</pre>
</div>
<div>
<pre><span>cmcmc</span> <span>&lt;-</span> <span>compileNimble</span><span>(mcmc,</span> <span>project</span> <span>= rModel)</span>
</pre>
</div>
<div>
<pre>## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
## compilation finished.
</pre>
</div>
<div>
<pre><span>samples</span> <span>&lt;-</span> <span>runMCMC</span><span>(cmcmc,</span> <span>niter</span> <span>=</span> <span>7000</span><span>,</span> <span>nburnin</span> <span>=</span> <span>2000</span><span>,</span> <span>setSeed</span> <span>=</span> <span>TRUE</span><span>)</span>
</pre>
</div>
<div>
<pre>## running chain 1...
</pre>
</div>
<div>
<pre>## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
</pre>
</div>
</div>
</div>
<p>In this case we use the posterior samples of the parameters to construct a trace plot and estimate the posterior distribution of <img alt="alpha" title="alpha" src="http://s0.wp.com/latex.php?latex=%5Calpha&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">:</p>
<div>
<div>
<div>
<pre><span># Trace plot of the posterior samples of the concentration parameter</span>
<span>ts.plot</span><span>(samples[ ,</span> <span>'alpha'</span><span>],</span> <span>xlab</span> <span>=</span> <span>&quot;iteration&quot;</span><span>,</span> <span>ylab</span> <span>=</span> <span>expression</span><span>(alpha))</span>
</pre>
</div>
</div>
<div><img src="https://i0.wp.com/r-nimble.org/wp-content/uploads/2018/12/results-CRPalpha2-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
<div>
<div>
<pre><span># Histogram of the posterior samples for the concentration parameter </span>
<span>hist</span><span>(samples[ ,</span> <span>'alpha'</span><span>],</span> <span>xlab</span> <span>=</span> <span>expression</span><span>(alpha),</span> <span>ylab</span> <span>=</span> <span>&quot;Frequency&quot;</span><span>,</span> <span>main</span> <span>=</span> <span>&quot;&quot;</span><span>)</span>
</pre>
</div>
</div>
<div><img src="https://i1.wp.com/r-nimble.org/wp-content/uploads/2018/12/results-CRPalpha2-2.png?w=456&ssl=1" data-recalc-dims="1"></div>
</div>
<h3> Generating samples from the mixing distribution</h3>
<p>As before, we obtain samples from the posterior distribution of <img alt="G" title="G" src="http://s0.wp.com/latex.php?latex=G&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> using the <tt>getSamplesDPmeasure</tt> function. </p>
<div>
<div>
<div>
<pre><span>outputG</span> <span>&lt;-</span> <span>getSamplesDPmeasure</span><span>(cmcmc)</span>
</pre>
</div>
<div>
<pre>## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
</pre>
</div>
<div>
<pre>## compilation finished.
</pre>
</div>
<div>
<pre>## sampleDPmeasure: Approximating the random measure by a finite stick-breaking representation with and error smaller than 1e-10, leads to a truncation level of 28.
</pre>
</div>
<div>
<pre><span>if</span><span>(</span><span>packageVersion</span><span>(</span><span>'nimble'</span><span>)</span> <span>&lt;=</span> <span>'0.6-12'</span><span>)</span>
  <span>samplesG</span> <span>&lt;-</span> <span>outputG</span> <span>else</span> <span>samplesG</span> <span>&lt;-</span> <span>outputG</span><span>$</span><span>samples</span>
</pre>
</div>
</div>
</div>
<p>We use these samples to create an estimate of the density of the data along with a pointwise 95% credible band:</p>
<div>
<div>
<div>
<pre><span>if</span><span>(</span><span>packageVersion</span><span>(</span><span>'nimble'</span><span>)</span> <span>>=</span> <span>'0.6.13'</span><span>)</span>
     <span>truncG</span> <span>&lt;-</span> <span>outputG</span><span>$</span><span>trunc</span> <span># truncation level for G</span>
<span>grid</span> <span>&lt;-</span> <span>seq</span><span>(</span><span>40</span><span>,</span> <span>100</span><span>,</span> <span>len</span> <span>=</span> <span>200</span><span>)</span>

<span>weightSamples</span> <span>&lt;-</span> <span>samplesG[ ,</span> <span>grep</span><span>(</span><span>'weight'</span><span>,</span> <span>colnames</span><span>(samplesG))]</span>
<span>betaTildeSamples</span> <span>&lt;-</span> <span>samplesG[ ,</span> <span>grep</span><span>(</span><span>'betaTilde'</span><span>,</span> <span>colnames</span><span>(samplesG))]</span>
<span>lambdaTildeSamples</span> <span>&lt;-</span> <span>samplesG[ ,</span> <span>grep</span><span>(</span><span>'lambdaTilde'</span><span>,</span> <span>colnames</span><span>(samplesG))]</span>

<span>densitySamples</span> <span>&lt;-</span> <span>matrix</span><span>(</span><span>0</span><span>,</span> <span>ncol</span> <span>=</span> <span>length</span><span>(grid),</span> <span>nrow</span> <span>=</span> <span>nrow</span><span>(samples))</span>
<span>for</span><span>(iter</span> <span>in</span> <span>seq_len</span><span>(</span><span>nrow</span><span>(samples))) {</span>
  <span>densitySamples[iter, ]</span> <span>&lt;-</span> <span>sapply</span><span>(grid,</span> <span>function</span><span>(</span><span>x</span><span>)</span>
    <span>sum</span><span>( weightSamples[iter, ]</span> <span>*</span> <span>dgamma</span><span>(x,</span> <span>shape</span> <span>= betaTildeSamples[iter, ],</span>
                <span>scale</span> <span>= lambdaTildeSamples[iter, ])))</span>
<span>}</span>


<span>hist</span><span>(faithful</span><span>$</span><span>waiting,</span> <span>freq</span> <span>=</span> <span>FALSE</span><span>,</span> <span>xlim</span> <span>=</span> <span>c</span><span>(</span><span>40</span><span>,</span><span>100</span><span>),</span> <span>ylim</span> <span>=</span> <span>c</span><span>(</span><span>0</span><span>,</span> <span>.05</span><span>),</span> <span>main</span> <span>=</span> <span>&quot;&quot;</span><span>,</span>
   <span>ylab</span> <span>=</span> <span>&quot;&quot;</span><span>,</span> <span>xlab</span> <span>=</span> <span>&quot;Waiting times&quot;</span><span>)</span>
<span>lines</span><span>(grid,</span> <span>apply</span><span>(densitySamples,</span> <span>2</span><span>, mean),</span> <span>lwd</span> <span>=</span> <span>2</span><span>,</span> <span>col</span> <span>=</span> <span>'black'</span><span>)</span>
<span>lines</span><span>(grid,</span> <span>apply</span><span>(densitySamples,</span> <span>2</span><span>, quantile,</span> <span>0.025</span><span>),</span> <span>lwd</span> <span>=</span> <span>2</span><span>,</span> <span>lty</span> <span>=</span> <span>2</span><span>,</span> <span>col</span> <span>=</span> <span>'black'</span><span>)</span>
<span>lines</span><span>(grid,</span> <span>apply</span><span>(densitySamples,</span> <span>2</span><span>, quantile,</span> <span>0.975</span><span>),</span> <span>lwd</span> <span>=</span> <span>2</span><span>,</span> <span>lty</span> <span>=</span> <span>2</span><span>,</span> <span>col</span> <span>=</span> <span>'black'</span><span>)</span>
</pre>
</div>
</div>
<div><img src="https://i2.wp.com/r-nimble.org/wp-content/uploads/2018/12/results-normaldensity-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
</div>
<p>Again, we see that the density of the data is bimodal, and looks very similar to the one we obtained before.</p>
<h2> Fitting a DP mixture of Gammas using a stick-breaking representation</h2>
<h3> Model specification</h3>
<p>An alternative representation of the Dirichlet process mixture uses the stick-breaking representation of the random distribution <img alt="G" title="G" src="http://s0.wp.com/latex.php?latex=G&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> (Sethuraman, 1994). NIMBLE allows us to specify an approximation that involves a truncation of the Dirichlet process to a finite number of atoms, <img alt="L" title="L" src="http://s0.wp.com/latex.php?latex=L&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">. The resulting model therefore reduces to a finite mixture with <img alt="L" title="L" src="http://s0.wp.com/latex.php?latex=L&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> components and a very particular prior on the weights of the mixture components.</p>
<p>Introducing auxiliary variables, <img alt="z_1, ldots, z_n" title="z_1, ldots, z_n" src="http://s0.wp.com/latex.php?latex=z_1%2C+%5Cldots%2C+z_n&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">, that indicate which component generated each observation, the corresponding model for the mixture of Gamma densities discussed in the previous section takes the form</p>
<p><img alt="  y_i mid { {beta}_k^{star} }, { {lambda}_k^{star} }, z_i sim mbox{Gamma}left( {beta}_{z_i}^{star}, {lambda}_{z_i}^{star} right), quadquad boldsymbol{z} mid boldsymbol{w} sim mbox{Discrete}(boldsymbol{w}), quadquad ({beta}_k^{star}, {lambda}_k^{star}) mid H sim H ,  " title="  y_i mid { {beta}_k^{star} }, { {lambda}_k^{star} }, z_i sim mbox{Gamma}left( {beta}_{z_i}^{star}, {lambda}_{z_i}^{star} right), quadquad boldsymbol{z} mid boldsymbol{w} sim mbox{Discrete}(boldsymbol{w}), quadquad ({beta}_k^{star}, {lambda}_k^{star}) mid H sim H ,  " src="http://s0.wp.com/latex.php?latex=++y_i+%5Cmid+%5C%7B+%7B%5Cbeta%7D_k%5E%7B%5Cstar%7D+%5C%7D%2C+%5C%7B+%7B%5Clambda%7D_k%5E%7B%5Cstar%7D+%5C%7D%2C+z_i+%5Csim+%5Cmbox%7BGamma%7D%5Cleft%28+%7B%5Cbeta%7D_%7Bz_i%7D%5E%7B%5Cstar%7D%2C+%7B%5Clambda%7D_%7Bz_i%7D%5E%7B%5Cstar%7D+%5Cright%29%2C+%5Cquad%5Cquad+%5Cboldsymbol%7Bz%7D+%5Cmid+%5Cboldsymbol%7Bw%7D+%5Csim+%5Cmbox%7BDiscrete%7D%28%5Cboldsymbol%7Bw%7D%29%2C+%5Cquad%5Cquad+%28%7B%5Cbeta%7D_k%5E%7B%5Cstar%7D%2C+%7B%5Clambda%7D_k%5E%7B%5Cstar%7D%29+%5Cmid+H+%5Csim+H+%2C++&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"></p>
<p>where <img alt="H" title="H" src="http://s0.wp.com/latex.php?latex=H&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> is again the product of two independent Gamma distributions,</p>
<p><img alt="  w_1=v_1, quadquad w_l=v_lprod_{m=1}^{l-1}(1-v_m), quad l=2, ldots, L-1,quadquad w_L=prod_{m=1}^{L-1}(1-v_m)  " title="  w_1=v_1, quadquad w_l=v_lprod_{m=1}^{l-1}(1-v_m), quad l=2, ldots, L-1,quadquad w_L=prod_{m=1}^{L-1}(1-v_m)  " src="http://s0.wp.com/latex.php?latex=++w_1%3Dv_1%2C+%5Cquad%5Cquad+w_l%3Dv_l%5Cprod_%7Bm%3D1%7D%5E%7Bl-1%7D%281-v_m%29%2C+%5Cquad+l%3D2%2C+%5Cldots%2C+L-1%2C%5Cquad%5Cquad+w_L%3D%5Cprod_%7Bm%3D1%7D%5E%7BL-1%7D%281-v_m%29++&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"></p>
<p>with <img alt="v_l mid alphasim mbox{Beta}(1, alpha), l=1, ldots, L-1" title="v_l mid alphasim mbox{Beta}(1, alpha), l=1, ldots, L-1" src="http://s0.wp.com/latex.php?latex=v_l+%5Cmid+%5Calpha%5Csim+%5Cmbox%7BBeta%7D%281%2C+%5Calpha%29%2C+l%3D1%2C+%5Cldots%2C+L-1&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">. The following code provides the NIMBLE specification for the model:</p>
<div>
<div>
<div>
<pre><span>code</span> <span>&lt;-</span> <span>nimbleCode</span><span>(</span>
  <span>{</span>
    <span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>n) {</span>
      <span>y[i]</span> <span>~</span> <span>dgamma</span><span>(</span><span>shape</span> <span>= beta[i],</span> <span>scale</span> <span>= lambda[i])</span>
      <span>beta[i]</span> <span>&lt;-</span> <span>betaStar[z[i]]</span>
      <span>lambda[i]</span> <span>&lt;-</span> <span>lambdaStar[z[i]]</span>
      <span>z[i]</span> <span>~</span> <span>dcat</span><span>(w[</span><span>1</span><span>:</span><span>Trunc])</span>
    <span>}</span>
    <span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>(Trunc</span><span>-</span><span>1</span><span>)) {</span> <span># stick-breaking variables</span>
      <span>v[i]</span> <span>~</span> <span>dbeta</span><span>(</span><span>1</span><span>, alpha)</span>
    <span>}</span>
    <span>w[</span><span>1</span><span>:</span><span>Trunc]</span> <span>&lt;-</span> <span>stick_breaking</span><span>(v[</span><span>1</span><span>:</span><span>(Trunc</span><span>-</span><span>1</span><span>)])</span> <span># stick-breaking weights</span>
    <span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>Trunc) {</span>
      <span>betaStar[i]</span> <span>~</span> <span>dgamma</span><span>(</span><span>shape</span> <span>=</span> <span>71</span><span>,</span> <span>scale</span> <span>=</span> <span>2</span><span>)</span>
      <span>lambdaStar[i]</span> <span>~</span> <span>dgamma</span><span>(</span><span>shape</span> <span>=</span> <span>2</span><span>,</span> <span>scale</span> <span>=</span> <span>2</span><span>)</span>
    <span>}</span>
    <span>alpha</span> <span>~</span> <span>dgamma</span><span>(</span><span>1</span><span>,</span> <span>1</span><span>)</span>
  <span>}</span>
<span>)</span>
</pre>
</div>
</div>
</div>
<p>Note that the truncation level <img alt="L" title="L" src="http://s0.wp.com/latex.php?latex=L&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> of <img alt="G" title="G" src="http://s0.wp.com/latex.php?latex=G&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> has been set to a value <tt>Trunc</tt>, which is to be defined in the <tt>constants</tt> argument of the <tt>nimbleModel</tt> function.</p>
<h3> Running the MCMC algorithm</h3>
<p>The following code sets up the model data and constants, initializes the parameters, defines the model object, and builds and runs the MCMC algorithm for the mixture of Gamma distributions. When a stick-breaking representation is used, a blocked Gibbs sampler is assigned (Ishwaran, 2001; Ishwaran and James, 2002).</p>
<div>
<div>
<div>
<pre><span>data</span> <span>&lt;-</span> <span>list</span><span>(</span><span>y</span> <span>= faithful</span><span>$</span><span>waiting)</span>
<span>set.seed</span><span>(</span><span>1</span><span>)</span>
<span>consts</span> <span>&lt;-</span> <span>list</span><span>(</span><span>n</span> <span>=</span> <span>length</span><span>(faithful</span><span>$</span><span>waiting),</span> <span>Trunc</span> <span>=</span> <span>50</span><span>)</span>
<span>inits</span> <span>&lt;-</span> <span>list</span><span>(</span><span>betaStar</span> <span>=</span> <span>rgamma</span><span>(consts</span><span>$</span><span>Trunc,</span> <span>shape</span> <span>=</span> <span>71</span><span>,</span> <span>scale</span> <span>=</span> <span>2</span><span>),</span>
              <span>lambdaStar</span> <span>=</span> <span>rgamma</span><span>(consts</span><span>$</span><span>Trunc,</span> <span>shape</span> <span>=</span> <span>2</span><span>,</span> <span>scale</span> <span>=</span> <span>2</span><span>),</span>
              <span>v</span> <span>=</span> <span>rbeta</span><span>(consts</span><span>$</span><span>Trunc</span><span>-</span><span>1</span><span>,</span> <span>1</span><span>,</span> <span>1</span><span>),</span>
              <span>z</span> <span>=</span> <span>sample</span><span>(</span><span>1</span><span>:</span><span>10</span><span>,</span> <span>size</span> <span>= consts</span><span>$</span><span>n,</span> <span>replace</span> <span>=</span> <span>TRUE</span><span>),</span>
              <span>alpha</span> <span>=</span> <span>1</span><span>)</span>

<span>rModel</span> <span>&lt;-</span> <span>nimbleModel</span><span>(code,</span> <span>data</span> <span>= data,</span> <span>inits</span> <span>= inits,</span> <span>constants</span> <span>= consts)</span>
</pre>
</div>
<div>
<pre>## defining model...
</pre>
</div>
<div>
<pre>## building model...
</pre>
</div>
<div>
<pre>## setting data and initial values...
</pre>
</div>
<div>
<pre>## running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
## checking model sizes and dimensions...
## model building finished.
</pre>
</div>
<div>
<pre><span>cModel</span> <span>&lt;-</span> <span>compileNimble</span><span>(rModel)</span>
</pre>
</div>
<div>
<pre>## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
## compilation finished.
</pre>
</div>
<div>
<pre><span>conf</span> <span>&lt;-</span> <span>configureMCMC</span><span>(rModel,</span> <span>monitors</span> <span>=</span> <span>c</span><span>(</span><span>&quot;w&quot;</span><span>,</span> <span>&quot;betaStar&quot;</span><span>,</span> <span>&quot;lambdaStar&quot;</span><span>,</span> <span>'z'</span><span>,</span> <span>'alpha'</span><span>))</span>
<span>mcmc</span> <span>&lt;-</span> <span>buildMCMC</span><span>(conf)</span>
<span>cmcmc</span> <span>&lt;-</span> <span>compileNimble</span><span>(mcmc,</span> <span>project</span> <span>= rModel)</span>
</pre>
</div>
<div>
<pre>## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
## compilation finished.
</pre>
</div>
<div>
<pre><span>samples</span> <span>&lt;-</span> <span>runMCMC</span><span>(cmcmc,</span> <span>niter</span> <span>=</span> <span>24000</span><span>,</span> <span>nburnin</span> <span>=</span> <span>4000</span><span>,</span> <span>setSeed</span> <span>=</span> <span>TRUE</span><span>)</span>
</pre>
</div>
<div>
<pre>## running chain 1...
</pre>
</div>
<div>
<pre>## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
</pre>
</div>
</div>
</div>
<p>Using the stick-breaking approximation automatically provides an approximation, <img alt="G_L" title="G_L" src="http://s0.wp.com/latex.php?latex=G_L&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">, of the random distribution <img alt="G" title="G" src="http://s0.wp.com/latex.php?latex=G&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex">. The following code computes posterior samples of <img alt="G_L" title="G_L" src="http://s0.wp.com/latex.php?latex=G_L&bg=ffffff&%23038;fg=000000&%23038;s=0" class="latex"> using posterior samples from the <tt>samples</tt> object, and from them, a density estimate for the data.</p>
<div>
<div>
<div>
<pre><span>betaStarSamples</span> <span>&lt;-</span> <span>samples[ ,</span> <span>grep</span><span>(</span><span>'betaStar'</span><span>,</span> <span>colnames</span><span>(samples))]</span>
<span>lambdaStarSamples</span> <span>&lt;-</span> <span>samples[ ,</span> <span>grep</span><span>(</span><span>'lambdaStar'</span><span>,</span> <span>colnames</span><span>(samples))]</span>
<span>weightSamples</span> <span>&lt;-</span> <span>samples[ ,</span> <span>grep</span><span>(</span><span>'w'</span><span>,</span> <span>colnames</span><span>(samples))]</span>

<span>grid</span> <span>&lt;-</span> <span>seq</span><span>(</span><span>40</span><span>,</span> <span>100</span><span>,</span> <span>len</span> <span>=</span> <span>200</span><span>)</span>

<span>densitySamples</span> <span>&lt;-</span> <span>matrix</span><span>(</span><span>0</span><span>,</span> <span>ncol</span> <span>=</span> <span>length</span><span>(grid),</span> <span>nrow</span> <span>=</span> <span>nrow</span><span>(samples))</span>
<span>for</span><span>(i</span> <span>in</span> <span>1</span><span>:</span><span>nrow</span><span>(samples)) {</span>
  <span>densitySamples[i, ]</span> <span>&lt;-</span> <span>sapply</span><span>(grid,</span> <span>function</span><span>(</span><span>x</span><span>)</span>
    <span>sum</span><span>(weightSamples[i, ]</span> <span>*</span> <span>dgamma</span><span>(x,</span> <span>shape</span> <span>= betaStarSamples[i, ],</span>
                                    <span>scale</span> <span>= lambdaStarSamples[i, ])))</span>
<span>}</span>

<span>hist</span><span>(faithful</span><span>$</span><span>waiting,</span> <span>freq</span> <span>=</span> <span>FALSE</span><span>,</span>  <span>xlab</span> <span>=</span> <span>&quot;Waiting times&quot;</span><span>,</span> <span>ylim</span><span>=</span><span>c</span><span>(</span><span>0</span><span>,</span><span>0.05</span><span>),</span>
     <span>main</span> <span>=</span> <span>''</span><span>)</span>
<span>lines</span><span>(grid,</span> <span>apply</span><span>(densitySamples,</span> <span>2</span><span>, mean),</span> <span>lwd</span> <span>=</span> <span>2</span><span>,</span> <span>col</span> <span>=</span> <span>'black'</span><span>)</span>
<span>lines</span><span>(grid,</span> <span>apply</span><span>(densitySamples,</span> <span>2</span><span>, quantile,</span> <span>0.025</span><span>),</span> <span>lwd</span> <span>=</span> <span>2</span><span>,</span> <span>lty</span> <span>=</span> <span>2</span><span>,</span> <span>col</span> <span>=</span> <span>'black'</span><span>)</span>
<span>lines</span><span>(grid,</span> <span>apply</span><span>(densitySamples,</span> <span>2</span><span>, quantile,</span> <span>0.975</span><span>),</span> <span>lwd</span> <span>=</span> <span>2</span><span>,</span> <span>lty</span> <span>=</span> <span>2</span><span>,</span> <span>col</span> <span>=</span> <span>'black'</span><span>)</span>
</pre>
</div>
</div>
<div><img src="https://i0.wp.com/r-nimble.org/wp-content/uploads/2018/12/MCMC-gammadensity-bnp-SB-G-1.png?w=456&ssl=1" data-recalc-dims="1"></div>
</div>
<p>As expected, this estimate looks identical to the one we obtained through the CRP representation of the process.</p>
<h2> More information and future development </h2>
<p>Please see our <a rel="nofollow" href="https://r-nimble.org/documentation" target="_blank">User Manual</a> for more details.</p>
<p>We’re in the midst of improvements to the existing BNP functionality as well as adding additional Bayesian nonparametric models, such as hierarchical Dirichlet processes and Pitman-Yor processes, so please add yourself to our <a rel="nofollow" href="https://groups.google.com/forum/#!forum/nimble-announce" target="_blank">announcement</a> or <a rel="nofollow" href="https://groups.google.com/forum/#!forum/nimble-users" target="_blank">user support/discussion</a> Google groups. </p>
<h3> References </h3>
<p>Blackwell, D. and MacQueen, J. 1973. Ferguson distributions via Polya urn schemes. The Annals of Statistics 1:353-355.</p>
<p>Ferguson, T.S. 1974. Prior distribution on the spaces of probability measures. Annals of Statistics 2:615-629.</p>
<p>Lo, A.Y. 1984. On a class of Bayesian nonparametric estimates I: Density estimates. The Annals of Statistics 12:351-357.</p>
<p>Escobar, M.D. 1994. Estimating normal means with a Dirichlet process prior. Journal of the American Statistical Association 89:268-277.</p>
<p>Escobar, M.D. and West, M. 1995. Bayesian density estimation and inference using mixtures. Journal of the American Statistical Association 90:577-588.</p>
<p>Ishwaran, H. and James, L.F. 2001. Gibbs sampling methods for stick-breaking priors. Journal of the American Statistical Association 96: 161-173.</p>
<p>Ishwaran, H. and James, L.F. 2002. Approximate Dirichlet process computing in finite normal mixtures: smoothing and prior information. Journal of Computational and Graphical Statistics 11:508-532.</p>
<p>Neal, R. 2000. Markov chain sampling methods for Dirichlet process mixture models. Journal of Computational and Graphical Statistics 9:249-265.</p>
<p>Sethuraman, J. 1994. A constructive definition of Dirichlet prior. Statistica Sinica 2: 639-650.</p>
<p><br>
</p>
<p> <img width="1" src="https://r-nimble.org/?feed-stats-post-id=898" height="1"></p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://r-nimble.org/bayesian-nonparametric-models-in-nimble-part-1-density-estimation"> R – NIMBLE</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/_YJJvUNfW7E" height="1"></a></div></div>