<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/fitting-a-distribution-in-stan-from-scratch/"></div></div><div><a data-url="https://www.r-bloggers.com/fitting-a-distribution-in-stan-from-scratch/" data-counturl="https://www.r-bloggers.com/fitting-a-distribution-in-stan-from-scratch/" data-count="horizontal" data-text="Fitting a distribution in Stan from scratch" 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.magesblog.com/2016/09/fitting-distribution-in-stan-from.html"> mages' blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Last week the French National Institute of Health and Medical Research (<a target="_blank" href="http://www.inserm.fr/">Inserm</a>) organised with the <a target="_blank" href="http://stan.fit/">Stan Group</a> a training programme on <i><a target="_blank" href="http://mc-stan.org/events/">Bayesian Inference with Stan for Pharmacometrics</a></i> in Paris. </p>
<p><a target="_blank" href="http://syclik.com/">Daniel Lee</a> and <a target="_blank" href="http://betanalpha.github.io/">Michael Betancourt</a>, who run the course over three days, are not only members of Stan’s development team, but also excellent teachers. Both were supported by <a target="_blank" href="https://about.me/ericnovik">Eric Novik</a>, who gave an <i><a target="_blank" href="http://pages.dataiku.com/hubfs/PDF/Meetup-Presentation/Intro_to_stan_by_Eric_Novik.pdf">Introduction to Stan</a></i> talk at the <a target="_blank" href="http://www.meetup.com/Dataiku-User-Group-France/">Paris Dataiku User Group</a> last week as well.</p>
<table align="center" cellpadding="0" cellspacing="0">
<tbody>
<tr>
<td><a target="_blank" href="https://i1.wp.com/3.bp.blogspot.com/-BzjFDDBPrds/V-gt6x8ZLEI/AAAAAAAACmk/8y-ZgmkOvhMCmVQNf_rDMy33GddjoTSrQCLcB/s1600/IMG_0650.jpg?ssl=1"><img src="https://i1.wp.com/3.bp.blogspot.com/-BzjFDDBPrds/V-gt6x8ZLEI/AAAAAAAACmk/8y-ZgmkOvhMCmVQNf_rDMy33GddjoTSrQCLcB/s400/IMG_0650.jpg?resize=400%2C195&ssl=1" border="0" data-recalc-dims="1"></a></td>
</tr>
<tr>
<td>Eric Kramer (Dataiku), Daniel Lee, Eric Novik &amp; Michael Betancourt (Stan Group)</td>
</tr>
</tbody>
</table>
<p>I have been playing around with Stan <a target="_blank" href="http://www.magesblog.com/search/label/Stan">on and off for some time</a>, but as Eric pointed out to me, Stan is not that kind of girl(boy?). Indeed, having spent three days working with Stan has revitalised my relationship. Getting down to the basics has been really helpful and I shall remember, Stan is not drawing samples from a distribution. Instead, it is calculating the joint distribution function (in log space), and evaluating the probability distribution function (in log space).</p>
<p>Thus, here is a little example of fitting a set of random numbers in R to a Normal distribution with Stan. Yet, instead of using the built-in functions for the Normal distribution, I define the log probability function by hand, which I will use in the model block as well, and even generate a random sample, starting with a uniform distribution. However, I do use pre-defined distributions for the priors.</p>
<p>Why do I want to do this? This will be a template for the day when I have to use a distribution, which is not predefined in Stan, e.g. the <a target="_blank" href="https://cran.r-project.org/web/packages/actuar/">actuar</a> package has some <a target="_blank" href="https://cran.r-project.org/web/packages/actuar/vignettes/lossdist.pdf">interesting candidates</a>.</p>
<p><br>
<h3>Testing</h3>
<p>I start off by generating fake data, a sample of 100 random numbers drawn from a Normal distribution with a mean of 4 and a standard deviation of 2. Note, the sample mean of the 100 figures is 4.2 and not 4.<br>
<table align="center" cellpadding="0" cellspacing="0">
<tbody>
<tr>
<td><a target="_blank" href="https://i0.wp.com/3.bp.blogspot.com/-cClakE2zyV4/V-gkE9D9H_I/AAAAAAAACmQ/jDqr6VDWndIS7fwQP6xi68dRGujdTDHTgCLcB/s1600/Histogram.png?ssl=1"><img src="https://i2.wp.com/3.bp.blogspot.com/-cClakE2zyV4/V-gkE9D9H_I/AAAAAAAACmQ/jDqr6VDWndIS7fwQP6xi68dRGujdTDHTgCLcB/s400/Histogram.png?resize=400%2C343&ssl=1" border="0" data-recalc-dims="1"></a></td>
</tr>
<tr>
<td>Histogram of 100 random numbers drawn from N(4,2).</td>
</tr>
</tbody>
</table>
<p>I then use the Stan script to fit the data, i.e. to find the the parameters \(\mu\) and \(\sigma\), assuming that the data was generated by a Gaussian process. <br><br>
<table align="center" cellpadding="0" cellspacing="0">
<tbody>
<tr>
<td><a target="_blank" href="https://i1.wp.com/4.bp.blogspot.com/-ijf4oA1noxM/V-gkE7lecpI/AAAAAAAACmI/jGm5Jy8LqTQQeM_3Dbx89N_AIgLhC85ggCLcB/s1600/Traceplot.png?ssl=1"><img src="https://i0.wp.com/4.bp.blogspot.com/-ijf4oA1noxM/V-gkE7lecpI/AAAAAAAACmI/jGm5Jy8LqTQQeM_3Dbx89N_AIgLhC85ggCLcB/s400/Traceplot.png?resize=400%2C203&ssl=1" border="0" data-recalc-dims="1"></a></td>
</tr>
<tr>
<td>Traceplot of 4 chains, including warm-up phase</td>
</tr>
</tbody>
</table>
<table align="center" cellpadding="0" cellspacing="0">
<tbody>
<tr>
<td><a target="_blank" href="https://i0.wp.com/2.bp.blogspot.com/-A4O9Bvv1qD4/V-gkEwQhmXI/AAAAAAAACmM/YMzo2hSFzOADc9s1CZkEfeJp1T9xUs6cgCLcB/s1600/StanOutput.png?ssl=1"><img src="https://i1.wp.com/2.bp.blogspot.com/-A4O9Bvv1qD4/V-gkEwQhmXI/AAAAAAAACmM/YMzo2hSFzOADc9s1CZkEfeJp1T9xUs6cgCLcB/s400/StanOutput.png?resize=400%2C203&ssl=1" border="0" data-recalc-dims="1"></a></td>
</tr>
<tr>
<td>Histograms of posterior parameter and predictive samples</td>
</tr>
</tbody>
</table>
<table align="center" cellpadding="0" cellspacing="0">
<tbody>
<tr>
<td><a target="_blank" href="https://i2.wp.com/3.bp.blogspot.com/-iZoTbE_Gjpc/V-g7eW3zIaI/AAAAAAAACm0/DU4YdlOVSHYc3YH3BsKed-GLGV4ZrONMgCLcB/s1600/PPC.png?ssl=1"><img src="https://i0.wp.com/3.bp.blogspot.com/-iZoTbE_Gjpc/V-g7eW3zIaI/AAAAAAAACm0/DU4YdlOVSHYc3YH3BsKed-GLGV4ZrONMgCLcB/s400/PPC.png?resize=400%2C400&ssl=1" border="0" data-recalc-dims="1"></a></td>
</tr>
<tr>
<td>Comparison of the emperical distributions</td>
</tr>
</tbody>
</table>
<p>The posterior parameter distributions include both \(\mu\) and \(\sigma\) in the 95% credible interval. The distribution of posterior predictive check (<code>y_ppc</code>) is wider, taking into account the uncertainty of the parameters. The interquartile range and mean of my initial fake data and the sample of the posterior predictive distribution look very similar. That’s good, my model generates data, which looks like the original data.<br>
<h3>Bayesian Mixer Meetup</h3>
<p>Btw, tonight we have the 4th <a target="_blank" href="http://www.meetup.com/Bayesian-Mixer-London/">Bayesian Mixer Meetup</a> in London. <br>
<h3>Session Info</h3>
<pre><code>R version 3.3.1 (2016-06-21)<br>Platform: x86_64-apple-darwin13.4.0 (64-bit)<br>Running under: OS X 10.12 (Sierra)<br><br>locale:<br>[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8<br><br>attached base packages:<br>[1] stats graphics grDevices utils datasets methods base     <br><br>other attached packages:<br>[1] MASS_7.3-45 rstan_2.12.1 StanHeaders_2.12.0 ggplot2_2.1.0     <br><br>loaded via a namespace (and not attached):<br> [1] Rcpp_0.12.7      codetools_0.2-14 digest_0.6.10    grid_3.3.1      <br> [5] plyr_1.8.4       gtable_0.2.0     stats4_3.3.1     scales_0.4.0    <br> [9] labeling_0.3     tools_3.3.1      munsell_0.4.3    inline_0.3.14   <br>[13] colorspace_1.2-6 gridExtra_2.2.1</code></pre>
<div>This post was originally published on <a target="_blank" href="http://www.magesblog.com/">mages’ blog</a>.</div>
<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/fitting-a-distribution-in-stan-from-scratch/"></div></div><div><a data-url="https://www.r-bloggers.com/fitting-a-distribution-in-stan-from-scratch/" data-counturl="https://www.r-bloggers.com/fitting-a-distribution-in-stan-from-scratch/" data-count="vertical" data-text="Fitting a distribution in Stan from scratch" 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.magesblog.com/2016/09/fitting-distribution-in-stan-from.html"> mages' blog</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/PbW9FPLTqI4"></a></div>