<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/an-attempt-at-ep-abc-from-scratch-nothing-more-except-for-a-few-bugs/"></div></div><div><a data-url="https://www.r-bloggers.com/an-attempt-at-ep-abc-from-scratch-nothing-more-except-for-a-few-bugs/" data-counturl="https://www.r-bloggers.com/an-attempt-at-ep-abc-from-scratch-nothing-more-except-for-a-few-bugs/" data-count="horizontal" data-text="an attempt at EP-ABC from scratch, nothing more… [except for a few bugs]" 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="https://xianblog.wordpress.com/2016/10/19/an-attempt-at-ep-abc-from-scratch-nothing-more-except-for-a-few-bugs/"> R – Xi'an's Og</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><strong>F</strong>ollowing a request from one of the reviewers of our chapter <a target="_blank" href="https://xianblog.wordpress.com/2015/03/27/likelihood-free-model-choice/">Likelihood-free model choice</a>, I tried to run <a target="_blank" href="https://xianblog.wordpress.com/2011/08/24/expectation-propagation-and-abc/">EP-ABC</a> on a toy problem and to compare it with the outcome of a random forest ABC. Literally starting from scratch, namely from the description found in <a target="_blank" href="http://www.tandfonline.com/doi/abs/10.1080/01621459.2013.864178#.VaYBz7zh4-o">Simon and Nicolas’ JASA paper</a>.  To run my test, I chose as my modelled data an exponential Exp(λ) sample of size 20, with a Gaussian N(0,1) prior on the log parameter (here replicated 100 times):</p>
<pre>
n=20
trueda=matrix(rexp(100*n),ncol=n)
#prior values
er0=0
Q0=1
</pre>
<p>Then I applied the formulas found in the paper for approximating the evidence, first defining the log normalising constant for an unnormalised Gaussian density as on the top of p.318</p>
<pre>
Psi=function(er,Q){
  -.5*(log(Q/2/pi)-Q*er*er)}
</pre>
<p>[which exhibited a typo in the <a target="_blank" href="http://www.tandfonline.com/doi/abs/10.1080/01621459.2013.864178#.VaYBz7zh4-o">paper</a>, as Simon Barthelmé figured out after emails exchanges that the right formulation was the inverse]</p>
<pre>
Psi=function(er,Q){
  -.5*(log(Q/2/pi)-er*er/Q)}
</pre>
<p>and iterating the site updates as described in Sections 2.2, 3.1 and Algorithms 1 and 3:</p>
<pre>
ep=function(dat,M=1e4,eps,Nep=10){
 n=length(dat)
 #global and ith natural parameters for Gaussian approximations
 #initialisation:
 Qi=rep(0,n)
 Q=Q0+sum(Qi)
 eri=rep(0,n)
 er=er0+sum(eri)
 #constants Ci in (2.6)
 normz=rep(1,n)
 for (t in 1:Nep){
  for (i in sample(1:n)){
  #site i update
   Qm=Q-Qi[i] #m for minus i
   erm=er-eri[i]
   #ABC sample
   thetas=rnorm(M,me=erm/Qm,sd=1/sqrt(Qm))
   dati=rexp(M,rate=exp(thetas))
#update by formula (3.3)
   Mh=sum((abs(dati-dat[i])&lt;eps))
   mu=mean(thetas[abs(dati-dat[i])&lt;eps])
  sig=var(thetas[abs(dati-dat[i])&lt;eps])
  if (Mh>1e3){
#enough proposals accepted
   Qi[i]=1/sig-Qm
   eri[i]=mu/sig-erm
   Q=Qm+Qi[i]
   er=erm+eri[i]
#update of Ci as in formula (2.7)
#with correction bottom of p.319
  normz[i]=log(Mh/M/2/eps)- Psi(er,Q)+Psi(erm,Qm)
  }}}
#approximation of the log evidence as on p.318
return(sum(normz)+Psi(er,Q)-Psi(er0,Q0)) }
</pre>
<p>except that I made an R beginner’s mistake (!) when calling the normal simulation to be</p>
<pre>
thetas=rnorm(M,me=erm/Qm)/sqrt(Qm)
</pre>
<p>to be compared with the genuine evidence under a conjugate Exp(1) prior <em>[values of the evidence should differ but should not differ that much]</em></p>
<pre>
ev=function(dat){
 n=length(dat)
 lgamma(n+1)-(n+1)*log(1+sum(dat))
 }
</pre>
<p>After correcting for those bugs and too small sample sizes, thanks to Simon kind-hearted help!, running this code results in minor discrepancies between both quantities:</p>
<pre>
evid=ovid=rep(0,100)
M=1e4 #number of ABC samples
propep=.1 #tolerance factor
for (t in 1:100){
 #tolerance scaled by data
 epsi=propep*sd(trueda[t,])
 evid[t]=ep(trueda[t,],M,epsi)
 ovid[t]=ev(trueda[t,])
 }
plot(evid,ovid,cex=.6,pch=19)
</pre>
<p>as shown by the comparison below between the evidence and the EP-ABC approximations to the evidence, called <em>epabence</em> (the dashed line is the diagonal):</p>
<p><a target="_blank" href="https://i1.wp.com/xianblog.files.wordpress.com/2015/07/epabc2.jpg?ssl=1"><img alt="epabc" class="wp-image-28145" src="https://i1.wp.com/xianblog.files.wordpress.com/2015/07/epabc2.jpg?resize=383%2C383&ssl=1" data-recalc-dims="1"></a>Obviously, this short (if lengthy at my scale) experiment is not intended to conclude that EP-ABC work or does not work. (Simon also provided additional comparison with the genuine evidence, that is under the right prior, and with the zero Monte-Carlo version of the EP-ABC evidence, achieving a high concentration over the diagonal.) I just conclude that the method does require a certain amount of calibration to become stable. Calibrations steps that are clearly spelled out in the <a target="_blank" href="http://www.tandfonline.com/doi/abs/10.1080/01621459.2013.864178#.VaYBz7zh4-o">paper</a> (adaptive choice of M, stabilisation by quasi-random samples, stabilisation by stochastic optimisation steps and importance sampling), but also imply that EP-ABC is also “far from routine use because it make take days to tune on real problems” (pp.315-316). Thinking about the reasons for such discrepancy over the past days, one possible explanation I see for the impact of the calibration parameters is that the validation of EP if any occurs at a numerical rather than Monte Carlo level, hence that the Monte Carlo error must be really small to avoid interfering with the numerical aspects.</p>
<p>Filed under: <a target="_blank" href="https://xianblog.wordpress.com/category/statistics/r-statistics/">R</a>, <a target="_blank" href="https://xianblog.wordpress.com/category/statistics/">Statistics</a>, <a target="_blank" href="https://xianblog.wordpress.com/category/university-life/">University life</a> Tagged: <a target="_blank" href="https://xianblog.wordpress.com/tag/abc/">ABC</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/computer-experiment/">computer experiment</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/expectation-propagation/">expectation-propagation</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/gaussian-approximation/">Gaussian approximation</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/implementation/">implementation</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/monte-carlo-statistical-methods/">Monte Carlo Statistical Methods</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/r/">R</a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/xianblog.wordpress.com/28035/"><img alt="" src="http://feeds.wordpress.com/1.0/comments/xianblog.wordpress.com/28035/" border="0"></a> <img alt="" src="https://i1.wp.com/pixel.wp.com/b.gif?resize=1%2C1&ssl=1" border="0" data-recalc-dims="1"></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/an-attempt-at-ep-abc-from-scratch-nothing-more-except-for-a-few-bugs/"></div></div><div><a data-url="https://www.r-bloggers.com/an-attempt-at-ep-abc-from-scratch-nothing-more-except-for-a-few-bugs/" data-counturl="https://www.r-bloggers.com/an-attempt-at-ep-abc-from-scratch-nothing-more-except-for-a-few-bugs/" data-count="vertical" data-text="an attempt at EP-ABC from scratch, nothing more… [except for a few bugs]" 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="https://xianblog.wordpress.com/2016/10/19/an-attempt-at-ep-abc-from-scratch-nothing-more-except-for-a-few-bugs/"> R – Xi'an's Og</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/pzabxgzFiw0"></a></div>