<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/cute-gibbs-sampling-for-rounded-observations/"></div></div><div><a data-url="https://www.r-bloggers.com/cute-gibbs-sampling-for-rounded-observations/" data-counturl="https://www.r-bloggers.com/cute-gibbs-sampling-for-rounded-observations/" data-count="horizontal" data-text="Cute Gibbs sampling for rounded observations" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.quantumforest.com/2016/11/cute-gibbs-sampling-for-rounded-observations/"> rblogs – Quantum Forest</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>I was attending a course of Bayesian Statistics where this problem showed up:</p>
<blockquote><p>There is a number of individuals, say 12, who take a pass/fail test 15 times. For each individual we have recorded the number of passes, which can go from 0 to 15. Because of confidentiality issues, we are presented with rounded-to-the-closest-multiple-of-3 data (\(\mathbf{R}\)). We are interested on estimating \(\theta\) of the Binomial distribution behind the data.</p>
<p>Rounding is probabilistic, with probability 2/3 if you are one count away from a multiple of 3 and probability 1/3 if the count is you are two counts away. Multiples of 3 are not rounded.</p>
</blockquote>
<p>We can use Gibbs sampling to alternate between sampling the posterior for the unrounded \(\mathbf{Y}\) and \(\theta\). In the case of \(\mathbf{Y}\) I used:</p>

<pre># Possible values that were rounded to R
possible &lt;- function(rounded) {
  if(rounded == 0) {
    options &lt;- c(0, 1, 2)
  } else {
    options &lt;- c(rounded - 2, rounded - 1, rounded,
                 rounded + 1, rounded + 2)
  }
  return(options)
}

# Probability mass function of numbers rounding to R
# given theta
prior_y &lt;- function(options, theta) {
  p &lt;- dbinom(options, 15, prob = theta)
  return(p)
}

# Likelihood of rounding 
like_round3 &lt;- function(options) {
  if(length(options) == 3) {
    like &lt;- c(1, 2/3, 1/3) }
  else {
    like &lt;- c(1/3, 2/3, 1, 2/3, 1/3)
    }
  return(like)
}

# Estimating posterior mass function and drawing a
# random value of it
posterior_sample_y &lt;- function(R, theta) {
  po &lt;- possible(R)
  pr &lt;- prior_y(po, theta)
  li &lt;- like_round3(po)
  post &lt;- li*pr/sum(li*pr)
  samp &lt;- sample(po, 1, prob = post)
  return(samp)
}</pre>

<p>While for \(\theta\) we are assuming a vague \(\mbox{Beta}(\alpha, \beta)\), with \(\alpha\) and \(\beta\) equal to 1, as prior density function for \(\theta\), so the posterior density is a \(\mbox{Beta}(\alpha + \sum Y_i, \beta + 12*15 – \sum Y_i)\).</p>

<pre>## Function to sample from the posterior Pr(theta | Y, R)
posterior_sample_theta &lt;- function(alpha, beta, Y) {
  theta &lt;- rbeta(1, alpha + sum(Y), beta + 12*15 - sum(Y))
  return(theta)
}</pre>

<p>I then implemented the sampler as:</p>

<pre>## Data
R &lt;- c(0, 0, 3, 9, 3, 0, 6, 3, 0, 6, 0, 3)
nsim &lt;- 10000
burnin &lt;- 1000
alpha &lt;- 1
beta &lt;- 1
store &lt;- matrix(0, nrow = nsim, ncol = length(R) + 1)

starting.values &lt;- c(R, 0.1)

## Sampling
store[1,] &lt;- starting.values
for(draw in 2:nsim){
  current &lt;- store[draw - 1,]
  for(obs in 1:length(R)) {
    y &lt;- posterior_sample_y(R[obs], current[length(R) + 1])
    # Jump or not still missing
    current[obs] &lt;- y
  }
  theta &lt;- posterior_sample_theta(alpha, beta, current[1:length(R)])
  # Jump or not still missing
  current[length(R) + 1] &lt;- theta
  
  store[draw,] &lt;- current
}</pre>

<p>And plotted the results as:</p>

<pre>plot((burnin+1):nsim, store[(burnin+1):nsim,13], type = 'l')

library(ggplot2)

ggplot(data.frame(theta = store[(burnin+1):nsim,13]), aes(x = theta)) +
  geom_density(fill = 'blue', alpha = 0.5)</pre>

<div class="wp-caption"><a target="_blank" rel="nofollow" href="https://i2.wp.com/www.quantumforest.com/wp-content/uploads/2016/11/theta-density.png"><img sizes="(max-width: 1239px) 100vw, 1239px" alt="Posterior density for Binomials's theta." class="wp-image-2703" src="https://i2.wp.com/www.quantumforest.com/wp-content/uploads/2016/11/theta-density.png?resize=456%2C321" data-recalc-dims="1"></a>
<p class="wp-caption-text">Posterior density for <img title="\theta" alt="\theta" class="latex" src="https://s0.wp.com/latex.php?latex=%5Ctheta&bg=ffffff&fg=000&s=0">.</p>
</div>
<pre>multiple_plot &lt;- data.frame(Y = matrix(store[(burnin+1):nsim, 1:12], 
                                       nrow = (nsim - burnin)*12,
                                       ncol = 1))
multiple_plot$obs &lt;- factor(rep(1:12, each = (nsim - burnin)))

ggplot(multiple_plot, aes(x = Y)) + geom_histogram() + facet_grid(~obs)</pre>

<div class="wp-caption"><a target="_blank" rel="nofollow" href="https://i2.wp.com/www.quantumforest.com/wp-content/uploads/2016/11/y-mass.png"><img sizes="(max-width: 1239px) 100vw, 1239px" alt="Posterior mass for each rounded observation." class="wp-image-2702" src="https://i2.wp.com/www.quantumforest.com/wp-content/uploads/2016/11/y-mass.png?resize=456%2C321" data-recalc-dims="1"></a>
<p class="wp-caption-text">Posterior mass for each rounded observation.</p>
</div>
<p>I thought it was a nice, cute example of simultaneously estimating a latent variable and, based on that, estimating the parameter behind it.</p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="https://www.r-bloggers.com/cute-gibbs-sampling-for-rounded-observations/"></div></div><div><a data-url="https://www.r-bloggers.com/cute-gibbs-sampling-for-rounded-observations/" data-counturl="https://www.r-bloggers.com/cute-gibbs-sampling-for-rounded-observations/" data-count="vertical" data-text="Cute Gibbs sampling for rounded observations" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.quantumforest.com/2016/11/cute-gibbs-sampling-for-rounded-observations/"> rblogs – Quantum Forest</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 title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/N1W9xIf3evM"></a></div>