-
Notifications
You must be signed in to change notification settings - Fork 1
/
a-gibbs-sampler-in-crystal.html
48 lines (48 loc) · 14 KB
/
a-gibbs-sampler-in-crystal.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
<!doctype html><html lang=en-uk><head><script data-goatcounter=https://ruivieira-dev.goatcounter.com/count async src=//gc.zgo.at/count.js></script><script src=https://unpkg.com/@alpinejs/intersect@3.x.x/dist/cdn.min.js></script><script src=https://unpkg.com/alpinejs@3.x.x/dist/cdn.min.js></script><script type=module src=https://ruivieira.dev/js/deeplinks/deeplinks.js></script><link rel=preload href=https://ruivieira.dev/lib/fonts/fa-brands-400.woff2 as=font type=font/woff2 crossorigin=anonymous><link rel=preload href=https://ruivieira.dev/lib/fonts/fa-regular-400.woff2 as=font type=font/woff2 crossorigin=anonymous><link rel=preload href=https://ruivieira.dev/lib/fonts/fa-solid-900.woff2 as=font type=font/woff2 crossorigin=anonymous><link rel=preload href=https://ruivieira.dev/fonts/firacode/FiraCode-Regular.woff2 as=font type=font/woff2 crossorigin=anonymous><link rel=preload href=https://ruivieira.dev/fonts/vollkorn/Vollkorn-Regular.woff2 as=font type=font/woff2 crossorigin=anonymous><link rel=stylesheet href=https://ruivieira.dev/css/kbd.css type=text/css><meta charset=utf-8><meta http-equiv=X-UA-Compatible content="IE=edge"><title>A Gibbs Sampler in Crystal · Rui Vieira</title>
<link rel=canonical href=https://ruivieira.dev/a-gibbs-sampler-in-crystal.html><meta name=viewport content="width=device-width,initial-scale=1"><meta name=robots content="all,follow"><meta name=googlebot content="index,follow,snippet,archive"><meta property="og:title" content="A Gibbs Sampler in Crystal"><meta property="og:description" content="Recently, I’ve been following with interest the development of the Crystal language.
Crystal is a statically typed language with a syntax resembling Ruby’s. The main features which drawn me to it were its simple boilerplate-free syntax (which is ideal for quick prototyping), tied with the ability to compile directly to native code along with a dead simple way of creating bindings to existing C code.
These features make it quite attractive, in my opinion, for scientific computing."><meta property="og:type" content="article"><meta property="og:url" content="https://ruivieira.dev/a-gibbs-sampler-in-crystal.html"><meta property="article:section" content="posts"><meta property="article:modified_time" content="2023-10-01T20:46:29+01:00"><meta name=twitter:card content="summary"><meta name=twitter:title content="A Gibbs Sampler in Crystal"><meta name=twitter:description content="Recently, I’ve been following with interest the development of the Crystal language.
Crystal is a statically typed language with a syntax resembling Ruby’s. The main features which drawn me to it were its simple boilerplate-free syntax (which is ideal for quick prototyping), tied with the ability to compile directly to native code along with a dead simple way of creating bindings to existing C code.
These features make it quite attractive, in my opinion, for scientific computing."><link rel=stylesheet href=https://ruivieira.dev/css/styles.css><!--[if lt IE 9]><script src=https://oss.maxcdn.com/html5shiv/3.7.2/html5shiv.min.js></script><script src=https://oss.maxcdn.com/respond/1.4.2/respond.min.js></script><![endif]--><link rel=icon type=image/png href=https://ruivieira.dev/images/favicon.ico></head><body class="max-width mx-auto px3 ltr" x-data="{currentHeading: undefined}"><div class="content index py4"><div id=header-post><a id=menu-icon href=#><i class="fas fa-eye fa-lg"></i></a>
<a id=menu-icon-tablet href=#><i class="fas fa-eye fa-lg"></i></a>
<a id=top-icon-tablet href=# onclick='$("html, body").animate({scrollTop:0},"fast")' style=display:none aria-label="Top of Page"><i class="fas fa-chevron-up fa-lg"></i></a>
<span id=menu><span id=nav><ul><li><a href=https://ruivieira.dev/>Home</a></li><li><a href=https://ruivieira.dev/blog/>Blog</a></li><li><a href=https://ruivieira.dev/draw/>Drawings</a></li><li><a href=https://ruivieira.dev/map/>All pages</a></li><li><a href=https://ruivieira.dev/search.html>Search</a></li></ul></span><br><div id=share style=display:none></div><div id=toc><h4>Contents</h4><nav id=TableOfContents><ul></ul></nav><h4>Related</h4><nav><ul><li class="header-post toc"><span class=backlink-count>1</span>
<a href=https://ruivieira.dev/a-simple-python-benchmark-exercise.html>A simple Python benchmark exercise</a></li><li class="header-post toc"><span class=backlink-count>1</span>
<a href=https://ruivieira.dev/mcmc-performance-on-substrate-vm.html>MCMC performance on Substrate VM</a></li></ul></nav></div></span></div><article class=post itemscope itemtype=http://schema.org/BlogPosting><header><h1 class=posttitle itemprop="name headline">A Gibbs Sampler in Crystal</h1><div class=meta><div class=postdate>Updated <time datetime="2023-10-01 20:46:29 +0100 BST" itemprop=datePublished>2023-10-01</time>
<span class=commit-hash>(<a href=https://ruivieira.dev/log/index.html#e23fe25>e23fe25</a>)</span></div></div></header><div class=content itemprop=articleBody><p>Recently, I’ve been following with interest the development of the <a href=http://crystal-lang.org/>Crystal</a> language.</p><p>Crystal is a statically typed language with a syntax resembling Ruby’s. The main features which drawn me to it were its simple boilerplate-free syntax (which is ideal for quick prototyping), tied with the ability to compile directly to native code along with a dead simple way of creating bindings to existing <a href>C</a> code.</p><p>These features make it quite attractive, in my opinion, for scientific computing. To test it against more popular languages, I’ve decided to run the Gibbs sampling examples created in <a href=https://darrenjw.wordpress.com/2011/07/16/gibbs-sampler-in-various-languages-revisited/>Darren Wilkinson’s blog</a>.</p><p>I recommend reading this post, and in fact, if you are interested in Mathematics and scientific computing in general, I <em>strongly</em> recommend you follow the <a href=https://darrenjw.wordpress.com/>blog</a>.</p><p>As explained in the linked post, I will make a Gibbs sampler for</p><p>$$
f\left(x,y\right)=kx^2\exp\left\lbrace-xy^2-y^2+2y-4x\right\rbrace
$$</p><p>with</p><p>$$
\begin{aligned}
x|y &\sim Ga\left(3,y^2+4\right) \\
y|x &\sim N\left(\frac{1}{1+x},\frac{1}{2\left(1+x\right)}\right)
\end{aligned}
$$</p><p>The original examples were ran again, without any code alterations. I’ve just added the Crystal version.</p><p>This implementation uses a <a href=https://github.com/ruivieira/crystal-gsl>very simple wrapper</a> I wrote to the famous <a href=http://www.gnu.org/software/gsl/>GNU Scientific Library</a> (GSL).</p><div class=highlight><pre tabindex=0 style=background-color:#fff;-moz-tab-size:4;-o-tab-size:4;tab-size:4><code class=language-crystal data-lang=crystal><span style=display:flex><span><span style=font-weight:700>require</span> <span style=color:#b84>"../libs/gsl/statistics.cr"</span>
</span></span><span style=display:flex><span>
</span></span><span style=display:flex><span><span style=font-weight:700>require</span> <span style=color:#b84>"math"</span>
</span></span><span style=display:flex><span>
</span></span><span style=display:flex><span><span style=font-weight:700>def</span> <span style=color:#900;font-weight:700>gibbs</span>(n : <span style=color:#999>Int</span> <span style=font-weight:700>=</span> <span style=color:#099>50000</span>, thin : <span style=color:#999>Int</span> <span style=font-weight:700>=</span> <span style=color:#099>1000</span>)
</span></span><span style=display:flex><span>
</span></span><span style=display:flex><span> x <span style=font-weight:700>=</span> <span style=color:#099>0.0</span>
</span></span><span style=display:flex><span> y <span style=font-weight:700>=</span> <span style=color:#099>0.0</span>
</span></span><span style=display:flex><span>
</span></span><span style=display:flex><span> <span style=color:#999>puts</span> <span style=color:#b84>"Iter x y"</span>
</span></span><span style=display:flex><span>
</span></span><span style=display:flex><span> (<span style=color:#099>0</span><span style=font-weight:700>..</span>n)<span style=font-weight:700>.</span>each <span style=font-weight:700>do</span> <span style=font-weight:700>|</span>i<span style=font-weight:700>|</span>
</span></span><span style=display:flex><span> (<span style=color:#099>0</span><span style=font-weight:700>..</span>thin)<span style=font-weight:700>.</span>each <span style=font-weight:700>do</span> <span style=font-weight:700>|</span>j<span style=font-weight:700>|</span>
</span></span><span style=display:flex><span> x <span style=font-weight:700>=</span> Statistics<span style=font-weight:700>::</span>Gamma<span style=font-weight:700>.</span>sample(<span style=color:#099>3.0</span>, y\<span style=font-weight:700>*</span>y<span style=font-weight:700>+</span><span style=color:#099>4.0</span>)
</span></span><span style=display:flex><span> y <span style=font-weight:700>=</span> Statistics<span style=font-weight:700>::</span>Normal<span style=font-weight:700>.</span>sample(<span style=color:#099>1.0</span><span style=font-weight:700>/</span>(x<span style=font-weight:700>+</span><span style=color:#099>1.0</span>), <span style=color:#099>1.0</span><span style=font-weight:700>/</span>Math<span style=font-weight:700>.</span>sqrt(<span style=color:#099>2.0</span>\<span style=font-weight:700>*</span>x<span style=font-weight:700>+</span><span style=color:#099>2.0</span>))
</span></span><span style=display:flex><span> <span style=font-weight:700>end</span>
</span></span><span style=display:flex><span>
</span></span><span style=display:flex><span> <span style=color:#999>puts</span> <span style=color:#b84>"</span><span style=color:#b84>#{</span>i<span style=color:#b84>}</span><span style=color:#b84> </span><span style=color:#b84>#{</span>x<span style=color:#b84>}</span><span style=color:#b84> </span><span style=color:#b84>#{</span>y<span style=color:#b84>}</span><span style=color:#b84>"</span>
</span></span><span style=display:flex><span>
</span></span><span style=display:flex><span> <span style=font-weight:700>end</span>
</span></span><span style=display:flex><span>
</span></span><span style=display:flex><span><span style=font-weight:700>end</span>
</span></span><span style=display:flex><span>
</span></span><span style=display:flex><span>gibbs
</span></span></code></pre></div><p>(As you can see, the Crystal code is quite similar to the Python one).</p><p>To make sure it’s a fair comparison, I ran it in compiled (and optimised mode build using</p><div class=highlight><pre tabindex=0 style=background-color:#fff;-moz-tab-size:4;-o-tab-size:4;tab-size:4><code class=language-shell data-lang=shell><span style=display:flex><span>$ crystal build gibbs.cr --release
</span></span><span style=display:flex><span>$ <span style=color:#999>time</span> ./gibbs > gibbs_crystal.csv
</span></span></code></pre></div><p>Looking at the results, you can see that they are consistent with the other implementations:</p><figure><img src=https://ruivieira.dev/site/images/gibbs_crystal.png alt=gibbs_crystal.png loading=lazy></figure><p>The timings for each of the different versions<sup id=fnref:1><a href=#fn:1 class=footnote-ref role=doc-noteref>1</a></sup> were</p><table><thead><tr><th>Language</th><th>Time (s)</th></tr></thead><tbody><tr><td>R</td><td>364.8</td></tr><tr><td>Python</td><td>144.0</td></tr><tr><td>Scala</td><td>9.896</td></tr><tr><td>Crystal</td><td>5.171</td></tr><tr><td>C</td><td>5.038</td></tr></tbody></table><p>So there you have it. A Ruby-like language which can easily compete with <a href>C</a> performance-wise.</p><p>I sincerely hope that Crystal gets some traction in the scientific community.</p><p>That of course won’t depend solely on its merits but rather on an active community along with a strong library ecosystem.</p><p>This is lacking at the moment, simply because it is relatively new language with the specs and standard library still being finalised.</p><div class=footnotes role=doc-endnotes><hr><ol><li id=fn:1><p>Ran in a 1.7 GHz Intel Core i7 Macbook Air. <a href=#fnref:1 class=footnote-backref role=doc-backlink>↩︎</a></p></li></ol></div></div></article><div id=footer-post-container><div id=footer-post><div id=nav-footer style=display:none><ul><li><a href=https://ruivieira.dev/>Home</a></li><li><a href=https://ruivieira.dev/blog/>Blog</a></li><li><a href=https://ruivieira.dev/draw/>Drawings</a></li><li><a href=https://ruivieira.dev/map/>All pages</a></li><li><a href=https://ruivieira.dev/search.html>Search</a></li></ul></div><div id=toc-footer style=display:none><nav id=TableOfContents></nav></div><div id=share-footer style=display:none></div><div id=actions-footer><a id=menu-toggle class=icon href=# onclick='return $("#nav-footer").toggle(),!1' aria-label=Menu><i class="fas fa-bars fa-lg" aria-hidden=true></i> Menu</a>
<a id=toc-toggle class=icon href=# onclick='return $("#toc-footer").toggle(),!1' aria-label=TOC><i class="fas fa-list fa-lg" aria-hidden=true></i> TOC</a>
<a id=share-toggle class=icon href=# onclick='return $("#share-footer").toggle(),!1' aria-label=Share><i class="fas fa-share-alt fa-lg" aria-hidden=true></i> share</a>
<a id=top style=display:none class=icon href=# onclick='$("html, body").animate({scrollTop:0},"fast")' aria-label="Top of Page"><i class="fas fa-chevron-up fa-lg" aria-hidden=true></i> Top</a></div></div></div><footer id=footer><div class=footer-left>Copyright © 2024 Rui Vieira</div><div class=footer-right><nav><ul><li><a href=https://ruivieira.dev/>Home</a></li><li><a href=https://ruivieira.dev/blog/>Blog</a></li><li><a href=https://ruivieira.dev/draw/>Drawings</a></li><li><a href=https://ruivieira.dev/map/>All pages</a></li><li><a href=https://ruivieira.dev/search.html>Search</a></li></ul></nav></div></footer></div></body><link rel=stylesheet href=https://ruivieira.dev/css/fa.min.css><script src=https://ruivieira.dev/js/jquery-3.6.0.min.js></script><script src=https://ruivieira.dev/js/mark.min.js></script><script src=https://ruivieira.dev/js/main.js></script><script>MathJax={tex:{inlineMath:[["$","$"],["\\(","\\)"]]},svg:{fontCache:"global"}}</script><script type=text/javascript id=MathJax-script async src=https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js></script></html>