-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathControlFlowDemo.html
466 lines (464 loc) · 37.3 KB
/
ControlFlowDemo.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
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml" lang="" xml:lang="">
<head>
<meta charset="utf-8" />
<meta name="generator" content="pandoc" />
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes"/>
<link rel="shortcut icon" href="images/favicon.ico" type="image/x-icon" />
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.3.0/css/all.min.css">
<link rel="stylesheet" href="https://cdn.jsdelivr.net/gh/jpswalsh/academicons@1/css/academicons.min.css">
<link href="https://cdn.jsdelivr.net/npm/[email protected]/dist/css/bootstrap.min.css" rel="stylesheet">
<script src="https://cdn.jsdelivr.net/npm/[email protected]/dist/js/bootstrap.bundle.min.js"></script>
<title>Discussion of program transformations and laziness</title>
<style>
html {
line-height: 1.5;
font-family: Georgia, serif;
font-size: 20px;
color: #1a1a1a;
background-color: #fdfdfd;
}
body {
margin: 0 auto;
max-width: 56em;
padding-left: 50px;
padding-right: 50px;
padding-top: 50px;
padding-bottom: 50px;
hyphens: auto;
overflow-wrap: break-word;
text-rendering: optimizeLegibility;
font-kerning: normal;
}
@media (max-width: 600px) {
body {
font-size: 0.9em;
padding: 1em;
}
h1 {
font-size: 1.8em;
}
}
@media print {
body {
background-color: transparent;
color: black;
font-size: 12pt;
}
p, h2, h3 {
orphans: 3;
widows: 3;
}
h2, h3, h4 {
page-break-after: avoid;
}
}
p {
margin: 1em 0;
}
a {
color: #1a1a1a;
}
a:visited {
color: #1a1a1a;
}
img {
max-width: 100%;
}
h1, h2, h3, h4, h5, h6 {
margin-top: 1.4em;
}
h5, h6 {
font-size: 1em;
font-style: italic;
}
h6 {
font-weight: normal;
}
ol, ul {
padding-left: 1.7em;
}
li > ol, li > ul {
margin-top: 0;
}
blockquote {
margin: 1em 0 1em 1.7em;
padding-left: 1em;
border-left: 2px solid #e6e6e6;
color: #606060;
}
code {
font-family: Menlo, Monaco, 'Lucida Console', Consolas, monospace;
font-size: 85%;
margin: 0;
}
pre {
margin: 1em 0;
overflow: auto;
}
pre code {
padding: 0;
overflow: visible;
overflow-wrap: normal;
}
.sourceCode {
background-color: transparent;
overflow: visible;
}
hr {
background-color: #1a1a1a;
border: none;
height: 1px;
margin: 1em 0;
}
table {
margin: 1em 0;
border-collapse: collapse;
width: 100%;
overflow-x: auto;
display: block;
font-variant-numeric: lining-nums tabular-nums;
}
table caption {
margin-bottom: 0.75em;
}
tbody {
margin-top: 0.5em;
border-top: 1px solid #1a1a1a;
border-bottom: 1px solid #1a1a1a;
}
th {
border-top: 1px solid #1a1a1a;
padding: 0.25em 0.5em 0.25em 0.5em;
}
td {
padding: 0.125em 0.5em 0.25em 0.5em;
}
header {
margin-bottom: 4em;
text-align: center;
}
#TOC li {
list-style: none;
}
#TOC ul {
padding-left: 1.3em;
}
#TOC > ul {
padding-left: 0;
}
#TOC a:not(:hover) {
text-decoration: none;
}
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
span.underline{text-decoration: underline;}
div.column{display: inline-block; vertical-align: top; width: 50%;}
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
ul.task-list{list-style: none;}
pre > code.sourceCode { white-space: pre; position: relative; }
pre > code.sourceCode > span { line-height: 1.25; }
pre > code.sourceCode > span:empty { height: 1.2em; }
.sourceCode { overflow: visible; }
code.sourceCode > span { color: inherit; text-decoration: inherit; }
div.sourceCode { margin: 1em 0; }
pre.sourceCode { margin: 0; }
@media screen {
div.sourceCode { overflow: auto; }
}
@media print {
pre > code.sourceCode { white-space: pre-wrap; }
pre > code.sourceCode > span { display: inline-block; text-indent: -5em; padding-left: 5em; }
}
pre.numberSource code
{ counter-reset: source-line 0; }
pre.numberSource code > span
{ position: relative; left: -4em; counter-increment: source-line; }
pre.numberSource code > span > a:first-child::before
{ content: counter(source-line);
position: relative; left: -1em; text-align: right; vertical-align: baseline;
border: none; display: inline-block;
-webkit-touch-callout: none; -webkit-user-select: none;
-khtml-user-select: none; -moz-user-select: none;
-ms-user-select: none; user-select: none;
padding: 0 4px; width: 4em;
color: #aaaaaa;
}
pre.numberSource { margin-left: 3em; border-left: 1px solid #aaaaaa; padding-left: 4px; }
div.sourceCode
{ background-color: #f8f8f8; }
@media screen {
pre > code.sourceCode > span > a:first-child::before { text-decoration: underline; }
}
code span.al { color: #ef2929; } /* Alert */
code span.an { color: #8f5902; font-weight: bold; font-style: italic; } /* Annotation */
code span.at { color: #204a87; } /* Attribute */
code span.bn { color: #0000cf; } /* BaseN */
code span.cf { color: #204a87; font-weight: bold; } /* ControlFlow */
code span.ch { color: #4e9a06; } /* Char */
code span.cn { color: #8f5902; } /* Constant */
code span.co { color: #8f5902; font-style: italic; } /* Comment */
code span.cv { color: #8f5902; font-weight: bold; font-style: italic; } /* CommentVar */
code span.do { color: #8f5902; font-weight: bold; font-style: italic; } /* Documentation */
code span.dt { color: #204a87; } /* DataType */
code span.dv { color: #0000cf; } /* DecVal */
code span.er { color: #a40000; font-weight: bold; } /* Error */
code span.ex { } /* Extension */
code span.fl { color: #0000cf; } /* Float */
code span.fu { color: #204a87; font-weight: bold; } /* Function */
code span.im { } /* Import */
code span.in { color: #8f5902; font-weight: bold; font-style: italic; } /* Information */
code span.kw { color: #204a87; font-weight: bold; } /* Keyword */
code span.op { color: #ce5c00; font-weight: bold; } /* Operator */
code span.ot { color: #8f5902; } /* Other */
code span.pp { color: #8f5902; font-style: italic; } /* Preprocessor */
code span.sc { color: #ce5c00; font-weight: bold; } /* SpecialChar */
code span.ss { color: #4e9a06; } /* SpecialString */
code span.st { color: #4e9a06; } /* String */
code span.va { color: #000000; } /* Variable */
code span.vs { color: #4e9a06; } /* VerbatimString */
code span.wa { color: #8f5902; font-weight: bold; font-style: italic; } /* Warning */
#title-logo {
margin-right: 30px;
}</style>
<!--[if lt IE 9]>
<script src="//cdnjs.cloudflare.com/ajax/libs/html5shiv/3.7.3/html5shiv-printshiv.min.js"></script>
<![endif]-->
</head>
<body>
<nav class="navbar navbar-expand-lg navbar-light bg-light">
<div class="container-fluid">
<a class="navbar-brand" href="index.html"><img src="images/logo.png" width="45" height="45" class="align-top d-inline-block" alt="LazyPPL"></a>
<button class="navbar-toggler" type="button" data-bs-toggle="collapse" data-bs-target="#navbarSupportedContent" aria-controls="navbarSupportedContent" aria-expanded="false" aria-label="Toggle navigation">
<span class="navbar-toggler-icon"></span>
</button>
<div class="collapse navbar-collapse" id="navbarSupportedContent">
<ul class="mb-2 navbar-nav me-auto mb-lg-0">
<li class="nav-item dropdown">
<a class="nav-link dropdown-toggle" href="#" id="navbarDropdown" role="button" data-bs-toggle="dropdown" aria-expanded="false">
Regression examples
</a>
<ul class="dropdown-menu" aria-labelledby="navbarDropdown">
<li><a class="dropdown-item" href="RegressionDemo.html">Linear and piecewise linear</a></li>
<li><a class="dropdown-item" href="GaussianProcessDemo.html">Gaussian process regression</a></li>
<li><a class="dropdown-item" href="WienerDemo.html">Wiener process regression</a></li>
<li><a class="dropdown-item" href="ProgramInductionDemo.html">Program induction</a></li>
</ul>
</li>
<li class="nav-item dropdown">
<a class="nav-link dropdown-toggle" href="#" id="navbarDropdown" role="button" data-bs-toggle="dropdown" aria-expanded="false">
Other examples
</a>
<ul class="dropdown-menu" aria-labelledby="navbarDropdown">
<li><a class="dropdown-item" href="MinimalDemo.html">Minimal demonstration</a></li>
<li><a class="dropdown-item" href="ClusteringDemo.html">Non-parametric clustering</a></li>
<li><a class="dropdown-item" href="IrmDemo.html">Infinite relational model</a></li>
<li><a class="dropdown-item" href="AdditiveClusteringDemo.html">Indian Buffet Process</a></li>
<li><a class="dropdown-item" href="ControlFlowDemo.html">Program transformations</a></li>
<li><a class="dropdown-item" href="GraphDemo.html">Graph inference</a></li>
<li><a class="dropdown-item" href="Physics.html">Inference in a physics model</a></li>
</ul>
</li>
<li class="nav-item"><a href="https://dl.acm.org/doi/abs/10.1145/3571239" class="nav-link"><i class="ai ai-acm"></i>
</a></li>
<li class="nav-item"><a href="https://arxiv.org/abs/2212.07250" class="nav-link"><i class="ai ai-arxiv"></i>
</a></li>
<li class="nav-item"><a href="https://github.com/lazyppl-team/lazyppl" class="nav-link"><i class="fa-brands fa-github"></i></a></li>
</ul>
</div>
</div>
</nav> <header id="title-block-header">
<h1 class="title">Discussion of program transformations and
laziness</h1>
</header>
<details class="code-details">
<summary>
Extensions and imports for this Literate Haskell file
</summary>
<div class="sourceCode" id="cb1"><pre
class="sourceCode haskell literate"><code class="sourceCode haskell"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a><span class="kw">module</span> <span class="dt">ControlFlowDemo</span> <span class="kw">where</span></span>
<span id="cb1-2"><a href="#cb1-2" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb1-3"><a href="#cb1-3" aria-hidden="true" tabindex="-1"></a><span class="kw">import</span> <span class="dt">LazyPPL</span></span>
<span id="cb1-4"><a href="#cb1-4" aria-hidden="true" tabindex="-1"></a><span class="kw">import</span> <span class="dt">LazyPPL.SingleSite</span></span>
<span id="cb1-5"><a href="#cb1-5" aria-hidden="true" tabindex="-1"></a><span class="kw">import</span> <span class="dt">LazyPPL.Distributions</span></span>
<span id="cb1-6"><a href="#cb1-6" aria-hidden="true" tabindex="-1"></a><span class="kw">import</span> <span class="dt">Control.Monad</span></span>
<span id="cb1-7"><a href="#cb1-7" aria-hidden="true" tabindex="-1"></a><span class="kw">import</span> <span class="dt">Graphics.Matplotlib</span></span></code></pre></div>
</details>
<p><br></p>
We consider a simple but slightly contrived model to illustrate how
control flow transformations and laziness can be advantageous. Our first
prior on pairs of booleans.
<div class="sourceCode" id="cb2"><pre
class="sourceCode haskell literate"><code class="sourceCode haskell"><span id="cb2-1"><a href="#cb2-1" aria-hidden="true" tabindex="-1"></a><span class="ot">prior1 ::</span> <span class="dt">Prob</span> (<span class="dt">Bool</span>,<span class="dt">Bool</span>)</span>
<span id="cb2-2"><a href="#cb2-2" aria-hidden="true" tabindex="-1"></a>prior1 <span class="ot">=</span> </span>
<span id="cb2-3"><a href="#cb2-3" aria-hidden="true" tabindex="-1"></a> <span class="kw">do</span> x <span class="ot"><-</span> bernoulli <span class="fl">0.5</span></span>
<span id="cb2-4"><a href="#cb2-4" aria-hidden="true" tabindex="-1"></a> y <span class="ot"><-</span> <span class="kw">if</span> x <span class="kw">then</span> bernoulli <span class="fl">0.4</span> <span class="kw">else</span> bernoulli <span class="fl">0.7</span></span>
<span id="cb2-5"><a href="#cb2-5" aria-hidden="true" tabindex="-1"></a> <span class="fu">return</span> (x,y)</span></code></pre></div>
Our model conditions on the two booleans being equal. This means that if
<code>x</code> changes, <code>y</code> has to change too, and vice
versa.
<div class="sourceCode" id="cb3"><pre
class="sourceCode haskell literate"><code class="sourceCode haskell"><span id="cb3-1"><a href="#cb3-1" aria-hidden="true" tabindex="-1"></a><span class="ot">model ::</span> <span class="dt">Prob</span> (<span class="dt">Bool</span>,<span class="dt">Bool</span>) <span class="ot">-></span> <span class="dt">Meas</span> <span class="dt">Bool</span></span>
<span id="cb3-2"><a href="#cb3-2" aria-hidden="true" tabindex="-1"></a>model prior <span class="ot">=</span></span>
<span id="cb3-3"><a href="#cb3-3" aria-hidden="true" tabindex="-1"></a> <span class="kw">do</span> (x,y) <span class="ot"><-</span> sample prior</span>
<span id="cb3-4"><a href="#cb3-4" aria-hidden="true" tabindex="-1"></a> score (<span class="kw">if</span> x<span class="op">==</span>y <span class="kw">then</span> <span class="dv">1</span> <span class="kw">else</span> <span class="dv">0</span>)</span>
<span id="cb3-5"><a href="#cb3-5" aria-hidden="true" tabindex="-1"></a> <span class="fu">return</span> x</span></code></pre></div>
<p>The posterior probability of <code>(model prior1)</code> returning
true is 4/7.</p>
<p>Since <code>x</code> and <code>y</code> are correlated, it is not a
good idea to use a naive Metropolis Hastings simulation. But this is an
extreme example, and less extreme examples might appear in more complex
probabilistic programs that do benefit from Metropolis Hastings
simulation, where the correlations are not known in advance.</p>
Here is a second implementation of the prior. It is mathematically
equivalent, and computationally equivalent on each run, because of
laziness. With Metropolis-Hastings, it is asymptotically the same (of
course, as long as the proposal is irreducible – see below). But in
practice this prior works better with Metropolis-Hastings because it
explicitly uses different sites for the two different conditional
branches.
<div class="sourceCode" id="cb4"><pre
class="sourceCode haskell literate"><code class="sourceCode haskell"><span id="cb4-1"><a href="#cb4-1" aria-hidden="true" tabindex="-1"></a><span class="ot">prior2 ::</span> <span class="dt">Prob</span> (<span class="dt">Bool</span>,<span class="dt">Bool</span>)</span>
<span id="cb4-2"><a href="#cb4-2" aria-hidden="true" tabindex="-1"></a>prior2 <span class="ot">=</span></span>
<span id="cb4-3"><a href="#cb4-3" aria-hidden="true" tabindex="-1"></a> <span class="kw">do</span> x <span class="ot"><-</span> bernoulli <span class="fl">0.5</span></span>
<span id="cb4-4"><a href="#cb4-4" aria-hidden="true" tabindex="-1"></a> ytrue <span class="ot"><-</span> bernoulli <span class="fl">0.4</span></span>
<span id="cb4-5"><a href="#cb4-5" aria-hidden="true" tabindex="-1"></a> yfalse <span class="ot"><-</span> bernoulli <span class="fl">0.7</span></span>
<span id="cb4-6"><a href="#cb4-6" aria-hidden="true" tabindex="-1"></a> <span class="fu">return</span> (<span class="kw">if</span> x <span class="kw">then</span> (x,ytrue) <span class="kw">else</span> (x,yfalse)) </span></code></pre></div>
<p>Of course, we could solve this particular program analytically. The
point is rather that we have a universal program transformation
manipulation that could be applied automatically.</p>
The following chart shows that the expected sample size (the ratio
between the variance of a single true sample and the variance among the
generated samples) is much better with this second implementation of the
prior, although Metropolis Hastings is still far from optimal for this
contrived problem. <img src="images/controlflow-essA.png" />
<details class="code-details">
<summary>
(Code for precision of estimators)
</summary>
<div class="sourceCode" id="cb5"><pre
class="sourceCode haskell literate"><code class="sourceCode haskell"><span id="cb5-1"><a href="#cb5-1" aria-hidden="true" tabindex="-1"></a><span class="co">-- Take a sampler (m) returning reals, regarded as an estimator,</span></span>
<span id="cb5-2"><a href="#cb5-2" aria-hidden="true" tabindex="-1"></a><span class="co">-- and calculate the mean and precision for the estimator.</span></span>
<span id="cb5-3"><a href="#cb5-3" aria-hidden="true" tabindex="-1"></a><span class="co">-- k is the number of samples to use in the estimator.</span></span>
<span id="cb5-4"><a href="#cb5-4" aria-hidden="true" tabindex="-1"></a><span class="co">-- n is the number of runs to average over to approximate the mean, variance and ess.</span></span>
<span id="cb5-5"><a href="#cb5-5" aria-hidden="true" tabindex="-1"></a><span class="ot">meanPrec ::</span> <span class="dt">IO</span> [(<span class="dt">Double</span>,a)] <span class="ot">-></span> <span class="dt">Int</span> <span class="ot">-></span> <span class="dt">Int</span> <span class="ot">-></span> <span class="dt">IO</span> (<span class="dt">Double</span>,<span class="dt">Double</span>)</span>
<span id="cb5-6"><a href="#cb5-6" aria-hidden="true" tabindex="-1"></a>meanPrec m n k <span class="ot">=</span></span>
<span id="cb5-7"><a href="#cb5-7" aria-hidden="true" tabindex="-1"></a> <span class="kw">do</span> xwss <span class="ot"><-</span> replicateM n m</span>
<span id="cb5-8"><a href="#cb5-8" aria-hidden="true" tabindex="-1"></a> <span class="kw">let</span> as <span class="ot">=</span> <span class="fu">map</span> ((\x <span class="ot">-></span> x <span class="op">/</span> (<span class="fu">fromIntegral</span> k))<span class="op">.</span> <span class="fu">sum</span> <span class="op">.</span> <span class="fu">map</span> <span class="fu">fst</span> <span class="op">.</span> <span class="fu">take</span> k) xwss</span>
<span id="cb5-9"><a href="#cb5-9" aria-hidden="true" tabindex="-1"></a> <span class="kw">let</span> mean <span class="ot">=</span> (<span class="fu">sum</span> as) <span class="op">/</span> (<span class="fu">fromIntegral</span> n) <span class="co">-- sample mean</span></span>
<span id="cb5-10"><a href="#cb5-10" aria-hidden="true" tabindex="-1"></a> <span class="kw">let</span> var <span class="ot">=</span> (<span class="fu">sum</span> (<span class="fu">map</span> (\x <span class="ot">-></span> (x <span class="op">-</span> mean)<span class="op">^</span><span class="dv">2</span>) as)) <span class="op">/</span> (<span class="fu">fromIntegral</span> n) <span class="co">-- sample variance</span></span>
<span id="cb5-11"><a href="#cb5-11" aria-hidden="true" tabindex="-1"></a> <span class="fu">return</span> (mean,<span class="dv">1</span><span class="op">/</span>var)</span></code></pre></div>
</details>
<details class="code-details">
<summary>
(Plotting code)
</summary>
<div class="sourceCode" id="cb6"><pre
class="sourceCode haskell literate"><code class="sourceCode haskell"><span id="cb6-1"><a href="#cb6-1" aria-hidden="true" tabindex="-1"></a>plotESSA <span class="ot">=</span> </span>
<span id="cb6-2"><a href="#cb6-2" aria-hidden="true" tabindex="-1"></a> <span class="kw">do</span> <span class="fu">putStrLn</span> <span class="st">"Plotting images/controlflow-essA.png..."</span></span>
<span id="cb6-3"><a href="#cb6-3" aria-hidden="true" tabindex="-1"></a> <span class="kw">let</span> xs <span class="ot">=</span> [<span class="dv">1</span>,<span class="dv">1000</span>,<span class="dv">10000</span>]</span>
<span id="cb6-4"><a href="#cb6-4" aria-hidden="true" tabindex="-1"></a> <span class="co">-- A rejection sampler where the scores are known to be 0 or 1.</span></span>
<span id="cb6-5"><a href="#cb6-5" aria-hidden="true" tabindex="-1"></a> <span class="kw">let</span> rejectionsampler m <span class="ot">=</span> <span class="fu">fmap</span> (<span class="fu">filter</span> (\(_,w) <span class="ot">-></span> w<span class="op">></span><span class="fl">0.00001</span>)) <span class="op">$</span> weightedsamples <span class="op">$</span> m</span>
<span id="cb6-6"><a href="#cb6-6" aria-hidden="true" tabindex="-1"></a> <span class="kw">let</span> truevar <span class="ot">=</span> (<span class="fl">4.0</span><span class="op">/</span><span class="fl">7.0</span>) <span class="op">-</span> (<span class="fl">4.0</span><span class="op">/</span><span class="fl">7.0</span>)<span class="op">^</span><span class="dv">2</span> <span class="co">-- variance for one true sample</span></span>
<span id="cb6-7"><a href="#cb6-7" aria-hidden="true" tabindex="-1"></a> <span class="co">-- Tests to plot</span></span>
<span id="cb6-8"><a href="#cb6-8" aria-hidden="true" tabindex="-1"></a> as <span class="ot"><-</span> forM xs <span class="op">$</span> <span class="fu">fmap</span> (\(_,prec)<span class="ot">-></span>truevar<span class="op">*</span>prec) <span class="op">.</span> meanPrec (mh <span class="fl">0.3</span> <span class="op">$</span> <span class="fu">fmap</span> fromBool <span class="op">$</span> model prior1) <span class="dv">100</span></span>
<span id="cb6-9"><a href="#cb6-9" aria-hidden="true" tabindex="-1"></a> cs <span class="ot"><-</span> forM xs <span class="op">$</span> <span class="fu">fmap</span> (\(_,prec)<span class="ot">-></span>truevar<span class="op">*</span>prec) <span class="op">.</span> meanPrec (mh <span class="fl">0.33</span> <span class="op">$</span> <span class="fu">fmap</span> fromBool <span class="op">$</span> model prior2) <span class="dv">100</span></span>
<span id="cb6-10"><a href="#cb6-10" aria-hidden="true" tabindex="-1"></a> es <span class="ot"><-</span> forM xs <span class="op">$</span> <span class="fu">fmap</span> (\(_,prec)<span class="ot">-></span>truevar<span class="op">*</span>prec) <span class="op">.</span> meanPrec (mh <span class="dv">1</span> <span class="op">$</span> <span class="fu">fmap</span> fromBool <span class="op">$</span> model prior1) <span class="dv">100</span></span>
<span id="cb6-11"><a href="#cb6-11" aria-hidden="true" tabindex="-1"></a> fs <span class="ot"><-</span> forM xs <span class="op">$</span> <span class="fu">fmap</span> (\(_,prec)<span class="ot">-></span>truevar<span class="op">*</span>prec) <span class="op">.</span> meanPrec (rejectionsampler <span class="op">$</span> <span class="fu">fmap</span> fromBool <span class="op">$</span> model prior1) <span class="dv">100</span></span>
<span id="cb6-12"><a href="#cb6-12" aria-hidden="true" tabindex="-1"></a> file <span class="st">"images/controlflow-essA.png"</span> <span class="op">$</span></span>
<span id="cb6-13"><a href="#cb6-13" aria-hidden="true" tabindex="-1"></a> plot xs as <span class="op">@@</span> [o2 <span class="st">"label"</span> <span class="st">"mh 0.3 with prior1"</span>, o2 <span class="st">"color"</span> <span class="st">"tomato"</span>] <span class="op">%</span></span>
<span id="cb6-14"><a href="#cb6-14" aria-hidden="true" tabindex="-1"></a> plot xs cs <span class="op">@@</span> [o2 <span class="st">"label"</span> <span class="st">"mh 0.3 with prior2"</span>, o2 <span class="st">"color"</span> <span class="st">"firebrick"</span>] <span class="op">%</span></span>
<span id="cb6-15"><a href="#cb6-15" aria-hidden="true" tabindex="-1"></a> plot xs es <span class="op">@@</span> [o2 <span class="st">"label"</span> <span class="st">"all sites mh"</span>, o2 <span class="st">"color"</span> <span class="st">"forestgreen"</span>] <span class="op">%</span></span>
<span id="cb6-16"><a href="#cb6-16" aria-hidden="true" tabindex="-1"></a> plot xs fs <span class="op">@@</span> [o2 <span class="st">"label"</span> <span class="st">"rejection sampling"</span>, o2 <span class="st">"color"</span> <span class="st">"black"</span>] <span class="op">%</span></span>
<span id="cb6-17"><a href="#cb6-17" aria-hidden="true" tabindex="-1"></a> ylim (<span class="dv">1</span><span class="ot">::</span><span class="dt">Double</span>) (<span class="dv">500</span><span class="ot">::</span><span class="dt">Double</span>) <span class="op">%</span></span>
<span id="cb6-18"><a href="#cb6-18" aria-hidden="true" tabindex="-1"></a> xlabel <span class="st">"Samples"</span> <span class="op">%</span> ylabel <span class="st">"Effective sample size"</span> <span class="op">%</span> legend <span class="op">@@</span> [o2 <span class="st">"loc"</span> <span class="st">"upper right"</span>]</span>
<span id="cb6-19"><a href="#cb6-19" aria-hidden="true" tabindex="-1"></a> <span class="fu">putStrLn</span> <span class="st">"done."</span></span></code></pre></div>
</details>
<h2 id="irreducibility">Irreducibility</h2>
This example also illustrates the issue with potential irreducibility in
single-site Metropolis Hastings (<code>mh1</code>): the first
implementation of the prior is especially bad, because <em>all</em>
proposals will be rejected except perhaps one. This irreducibility
arises from the hard constraints (scores of <code>0</code>) and is
avoided to some extent with a softer constraint:
<div class="sourceCode" id="cb7"><pre
class="sourceCode haskell literate"><code class="sourceCode haskell"><span id="cb7-1"><a href="#cb7-1" aria-hidden="true" tabindex="-1"></a><span class="ot">softModel ::</span> <span class="dt">Prob</span> (<span class="dt">Bool</span>,<span class="dt">Bool</span>) <span class="ot">-></span> <span class="dt">Meas</span> <span class="dt">Bool</span></span>
<span id="cb7-2"><a href="#cb7-2" aria-hidden="true" tabindex="-1"></a>softModel prior <span class="ot">=</span></span>
<span id="cb7-3"><a href="#cb7-3" aria-hidden="true" tabindex="-1"></a> <span class="kw">do</span> (x,y) <span class="ot"><-</span> sample prior</span>
<span id="cb7-4"><a href="#cb7-4" aria-hidden="true" tabindex="-1"></a> score (<span class="kw">if</span> x<span class="op">==</span>y <span class="kw">then</span> <span class="fl">0.9</span> <span class="kw">else</span> <span class="fl">0.1</span>)</span>
<span id="cb7-5"><a href="#cb7-5" aria-hidden="true" tabindex="-1"></a> <span class="fu">return</span> x</span></code></pre></div>
<img src="images/controlflow-essB.png" />
<details class="code-details">
<summary>
(Plotting code)
</summary>
<div class="sourceCode" id="cb8"><pre
class="sourceCode haskell literate"><code class="sourceCode haskell"><span id="cb8-1"><a href="#cb8-1" aria-hidden="true" tabindex="-1"></a>fromBool b <span class="ot">=</span> <span class="kw">if</span> b <span class="kw">then</span> <span class="fl">1.0</span> <span class="kw">else</span> <span class="fl">0.0</span></span>
<span id="cb8-2"><a href="#cb8-2" aria-hidden="true" tabindex="-1"></a>plotESSB <span class="ot">=</span> </span>
<span id="cb8-3"><a href="#cb8-3" aria-hidden="true" tabindex="-1"></a> <span class="kw">do</span> <span class="fu">putStrLn</span> <span class="st">"Plotting images/controlflow-essB.png..."</span></span>
<span id="cb8-4"><a href="#cb8-4" aria-hidden="true" tabindex="-1"></a> <span class="kw">let</span> xs <span class="ot">=</span> [<span class="dv">1</span>,<span class="dv">2</span>,<span class="dv">3</span>,<span class="dv">4</span>,<span class="dv">5</span>,<span class="dv">6</span>,<span class="dv">7</span>,<span class="dv">8</span>,<span class="dv">9</span>,<span class="dv">10</span>,<span class="dv">11</span>,<span class="dv">12</span>,<span class="dv">13</span>,<span class="dv">14</span>,<span class="dv">15</span>,<span class="dv">20</span>,<span class="dv">30</span>,<span class="dv">40</span>,<span class="dv">50</span>,<span class="dv">60</span>,<span class="dv">70</span>,<span class="dv">80</span>,<span class="dv">90</span>,<span class="dv">100</span>]</span>
<span id="cb8-5"><a href="#cb8-5" aria-hidden="true" tabindex="-1"></a> <span class="kw">let</span> truevar <span class="ot">=</span> (<span class="fl">4.0</span><span class="op">/</span><span class="fl">7.0</span>) <span class="op">-</span> (<span class="fl">4.0</span><span class="op">/</span><span class="fl">7.0</span>)<span class="op">^</span><span class="dv">2</span> <span class="co">-- variance for one true sample</span></span>
<span id="cb8-6"><a href="#cb8-6" aria-hidden="true" tabindex="-1"></a> <span class="co">-- Tests to plot</span></span>
<span id="cb8-7"><a href="#cb8-7" aria-hidden="true" tabindex="-1"></a> as <span class="ot"><-</span> forM xs <span class="op">$</span> <span class="fu">fmap</span> (\(_,prec)<span class="ot">-></span>truevar<span class="op">*</span>prec) <span class="op">.</span> meanPrec (mh1 <span class="op">$</span> <span class="fu">fmap</span> fromBool <span class="op">$</span> model prior1) <span class="dv">1000</span> </span>
<span id="cb8-8"><a href="#cb8-8" aria-hidden="true" tabindex="-1"></a> bs <span class="ot"><-</span> forM xs <span class="op">$</span> <span class="fu">fmap</span> (\(_,prec)<span class="ot">-></span>truevar<span class="op">*</span>prec) <span class="op">.</span> meanPrec (mh1 <span class="op">$</span> <span class="fu">fmap</span> fromBool <span class="op">$</span> softModel prior1) <span class="dv">1000</span></span>
<span id="cb8-9"><a href="#cb8-9" aria-hidden="true" tabindex="-1"></a> file <span class="st">"images/controlflow-essB.png"</span> <span class="op">$</span></span>
<span id="cb8-10"><a href="#cb8-10" aria-hidden="true" tabindex="-1"></a> plot xs as <span class="op">@@</span> [o2 <span class="st">"label"</span> <span class="st">"mh1 with hard constraint"</span>, o2 <span class="st">"color"</span> <span class="st">"mediumblue"</span>] <span class="op">%</span></span>
<span id="cb8-11"><a href="#cb8-11" aria-hidden="true" tabindex="-1"></a> plot xs bs <span class="op">@@</span> [o2 <span class="st">"label"</span> <span class="st">"mh1 with soft constraint"</span>, o2 <span class="st">"color"</span> <span class="st">"deepskyblue"</span>] <span class="op">%</span></span>
<span id="cb8-12"><a href="#cb8-12" aria-hidden="true" tabindex="-1"></a> xlabel <span class="st">"Samples"</span> <span class="op">%</span> ylabel <span class="st">"Effective sample size"</span> <span class="op">%</span> legend <span class="op">@@</span> [o2 <span class="st">"loc"</span> <span class="st">"upper right"</span>]</span>
<span id="cb8-13"><a href="#cb8-13" aria-hidden="true" tabindex="-1"></a> <span class="fu">putStrLn</span> <span class="st">"done."</span></span></code></pre></div>
</details>
<h2
id="note-about-different-metropolis-hastings-kernels-in-general">Note
about different Metropolis-Hastings kernels in general</h2>
These are contrived examples and the reader should not conclude that
all-sites Metropolis-Hastings outperforms the more complex
Metropolis-Hastings algorithms of LazyPPL in general. To emphasise this,
we plot the precision in the intercept for a linear regression example
(similar to <a href="Regression.html">here</a>).
<div class="sourceCode" id="cb9"><pre
class="sourceCode haskell literate"><code class="sourceCode haskell"><span id="cb9-1"><a href="#cb9-1" aria-hidden="true" tabindex="-1"></a><span class="ot">linreg ::</span> <span class="dt">Meas</span> (<span class="dt">Double</span>,<span class="dt">Double</span>)</span>
<span id="cb9-2"><a href="#cb9-2" aria-hidden="true" tabindex="-1"></a>linreg <span class="ot">=</span></span>
<span id="cb9-3"><a href="#cb9-3" aria-hidden="true" tabindex="-1"></a> <span class="kw">do</span></span>
<span id="cb9-4"><a href="#cb9-4" aria-hidden="true" tabindex="-1"></a> a <span class="ot"><-</span> sample <span class="op">$</span> normal <span class="dv">0</span> <span class="dv">3</span></span>
<span id="cb9-5"><a href="#cb9-5" aria-hidden="true" tabindex="-1"></a> b <span class="ot"><-</span> sample <span class="op">$</span> normal <span class="dv">0</span> <span class="dv">3</span></span>
<span id="cb9-6"><a href="#cb9-6" aria-hidden="true" tabindex="-1"></a> <span class="kw">let</span> dataset <span class="ot">=</span> [(<span class="dv">0</span>,<span class="dv">5</span>), (<span class="fl">0.1</span>,<span class="dv">10</span>), (<span class="fl">0.2</span>,<span class="dv">0</span>)]</span>
<span id="cb9-7"><a href="#cb9-7" aria-hidden="true" tabindex="-1"></a> <span class="fu">mapM</span> (\(x, y) <span class="ot">-></span> score <span class="op">$</span> normalPdf (a<span class="op">*</span>x <span class="op">+</span> b) <span class="fl">0.5</span> y) dataset</span>
<span id="cb9-8"><a href="#cb9-8" aria-hidden="true" tabindex="-1"></a> <span class="fu">return</span> (a,b)</span></code></pre></div>
<p><img src="images/controlflow-prec.png" /></p>
<p>(Here we plot precision, which is proportional to effective sample
size, because a general effective sample size analysis isn’t written
yet.) <br></p>
<details class="code-details">
<summary>
(Plotting code)
</summary>
<div class="sourceCode" id="cb10"><pre
class="sourceCode haskell literate"><code class="sourceCode haskell"><span id="cb10-1"><a href="#cb10-1" aria-hidden="true" tabindex="-1"></a>plotPrec <span class="ot">=</span> </span>
<span id="cb10-2"><a href="#cb10-2" aria-hidden="true" tabindex="-1"></a> <span class="kw">do</span> <span class="fu">putStrLn</span> <span class="st">"Plotting images/controlflow-prec.png..."</span></span>
<span id="cb10-3"><a href="#cb10-3" aria-hidden="true" tabindex="-1"></a> <span class="kw">let</span> xs <span class="ot">=</span> [<span class="dv">1</span>,<span class="dv">100</span>,<span class="dv">200</span>,<span class="dv">400</span>,<span class="dv">700</span>,<span class="dv">1000</span>]</span>
<span id="cb10-4"><a href="#cb10-4" aria-hidden="true" tabindex="-1"></a> <span class="co">-- Tests to plot</span></span>
<span id="cb10-5"><a href="#cb10-5" aria-hidden="true" tabindex="-1"></a> as <span class="ot"><-</span> forM xs <span class="op">$</span> <span class="fu">fmap</span> <span class="fu">snd</span> <span class="op">.</span> meanPrec (mh <span class="fl">0.5</span> <span class="op">$</span> <span class="fu">fmap</span> <span class="fu">snd</span> <span class="op">$</span> linreg) <span class="dv">1000</span></span>
<span id="cb10-6"><a href="#cb10-6" aria-hidden="true" tabindex="-1"></a> bs <span class="ot"><-</span> forM xs <span class="op">$</span> <span class="fu">fmap</span> <span class="fu">snd</span> <span class="op">.</span> meanPrec (mh1 <span class="op">$</span> <span class="fu">fmap</span> <span class="fu">snd</span> <span class="op">$</span> linreg) <span class="dv">1000</span></span>
<span id="cb10-7"><a href="#cb10-7" aria-hidden="true" tabindex="-1"></a> es <span class="ot"><-</span> forM xs <span class="op">$</span> <span class="fu">fmap</span> <span class="fu">snd</span> <span class="op">.</span> meanPrec (mh <span class="dv">1</span> <span class="op">$</span> <span class="fu">fmap</span> <span class="fu">snd</span> <span class="op">$</span> linreg) <span class="dv">1000</span></span>
<span id="cb10-8"><a href="#cb10-8" aria-hidden="true" tabindex="-1"></a> file <span class="st">"images/controlflow-prec.png"</span> <span class="op">$</span></span>
<span id="cb10-9"><a href="#cb10-9" aria-hidden="true" tabindex="-1"></a> plot xs bs <span class="op">@@</span> [o2 <span class="st">"label"</span> <span class="st">"mh1"</span>, o2 <span class="st">"color"</span> <span class="st">"mediumblue"</span>] <span class="op">%</span></span>
<span id="cb10-10"><a href="#cb10-10" aria-hidden="true" tabindex="-1"></a> plot xs as <span class="op">@@</span> [o2 <span class="st">"label"</span> <span class="st">"mh 0.5"</span>, o2 <span class="st">"color"</span> <span class="st">"firebrick"</span>] <span class="op">%</span></span>
<span id="cb10-11"><a href="#cb10-11" aria-hidden="true" tabindex="-1"></a> plot xs es <span class="op">@@</span> [o2 <span class="st">"label"</span> <span class="st">"all sites mh"</span>, o2 <span class="st">"color"</span> <span class="st">"forestgreen"</span>] <span class="op">%</span></span>
<span id="cb10-12"><a href="#cb10-12" aria-hidden="true" tabindex="-1"></a> xlabel <span class="st">"Samples"</span> <span class="op">%</span> ylabel <span class="st">"Precision for mean b"</span> <span class="op">%</span> legend <span class="op">@@</span> [o2 <span class="st">"loc"</span> <span class="st">"upper right"</span>]</span>
<span id="cb10-13"><a href="#cb10-13" aria-hidden="true" tabindex="-1"></a> <span class="fu">putStrLn</span> <span class="st">"done."</span></span></code></pre></div>
</details>
<details class="code-details">
<summary>
(Epilogue code)
</summary>
<div class="sourceCode" id="cb11"><pre
class="sourceCode haskell literate"><code class="sourceCode haskell"><span id="cb11-1"><a href="#cb11-1" aria-hidden="true" tabindex="-1"></a>main <span class="ot">=</span> <span class="kw">do</span></span>
<span id="cb11-2"><a href="#cb11-2" aria-hidden="true" tabindex="-1"></a> plotESSA</span>
<span id="cb11-3"><a href="#cb11-3" aria-hidden="true" tabindex="-1"></a> plotESSB</span>
<span id="cb11-4"><a href="#cb11-4" aria-hidden="true" tabindex="-1"></a> plotPrec</span></code></pre></div>
</details>
<br>
<hr>
<small><i> Generated by Pandoc from Literate Haskell. Full source on <a href="https://github.com/lazyppl-team/lazyppl"><i class="fa-brands fa-github"></i></a>.</i></small>
</body>
</html>