-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathcode.html
2478 lines (2090 loc) · 875 KB
/
code.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
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
<title></title>
<script type="text/javascript">
window.onload = function() {
var imgs = document.getElementsByTagName('img'), i, img;
for (i = 0; i < imgs.length; i++) {
img = imgs[i];
// center an image if it is the only element of its parent
if (img.parentElement.childElementCount === 1)
img.parentElement.style.textAlign = 'center';
}
};
</script>
<!-- Styles for R syntax highlighter -->
<style type="text/css">
pre .operator,
pre .paren {
color: rgb(104, 118, 135)
}
pre .literal {
color: #990073
}
pre .number {
color: #099;
}
pre .comment {
color: #998;
font-style: italic
}
pre .keyword {
color: #900;
font-weight: bold
}
pre .identifier {
color: rgb(0, 0, 0);
}
pre .string {
color: #d14;
}
</style>
<!-- R syntax highlighter -->
<script type="text/javascript">
var hljs=new function(){function m(p){return p.replace(/&/gm,"&").replace(/</gm,"<")}function f(r,q,p){return RegExp(q,"m"+(r.cI?"i":"")+(p?"g":""))}function b(r){for(var p=0;p<r.childNodes.length;p++){var q=r.childNodes[p];if(q.nodeName=="CODE"){return q}if(!(q.nodeType==3&&q.nodeValue.match(/\s+/))){break}}}function h(t,s){var p="";for(var r=0;r<t.childNodes.length;r++){if(t.childNodes[r].nodeType==3){var q=t.childNodes[r].nodeValue;if(s){q=q.replace(/\n/g,"")}p+=q}else{if(t.childNodes[r].nodeName=="BR"){p+="\n"}else{p+=h(t.childNodes[r])}}}if(/MSIE [678]/.test(navigator.userAgent)){p=p.replace(/\r/g,"\n")}return p}function a(s){var r=s.className.split(/\s+/);r=r.concat(s.parentNode.className.split(/\s+/));for(var q=0;q<r.length;q++){var p=r[q].replace(/^language-/,"");if(e[p]){return p}}}function c(q){var p=[];(function(s,t){for(var r=0;r<s.childNodes.length;r++){if(s.childNodes[r].nodeType==3){t+=s.childNodes[r].nodeValue.length}else{if(s.childNodes[r].nodeName=="BR"){t+=1}else{if(s.childNodes[r].nodeType==1){p.push({event:"start",offset:t,node:s.childNodes[r]});t=arguments.callee(s.childNodes[r],t);p.push({event:"stop",offset:t,node:s.childNodes[r]})}}}}return t})(q,0);return p}function k(y,w,x){var q=0;var z="";var s=[];function u(){if(y.length&&w.length){if(y[0].offset!=w[0].offset){return(y[0].offset<w[0].offset)?y:w}else{return w[0].event=="start"?y:w}}else{return y.length?y:w}}function t(D){var A="<"+D.nodeName.toLowerCase();for(var B=0;B<D.attributes.length;B++){var C=D.attributes[B];A+=" "+C.nodeName.toLowerCase();if(C.value!==undefined&&C.value!==false&&C.value!==null){A+='="'+m(C.value)+'"'}}return A+">"}while(y.length||w.length){var v=u().splice(0,1)[0];z+=m(x.substr(q,v.offset-q));q=v.offset;if(v.event=="start"){z+=t(v.node);s.push(v.node)}else{if(v.event=="stop"){var p,r=s.length;do{r--;p=s[r];z+=("</"+p.nodeName.toLowerCase()+">")}while(p!=v.node);s.splice(r,1);while(r<s.length){z+=t(s[r]);r++}}}}return z+m(x.substr(q))}function j(){function q(x,y,v){if(x.compiled){return}var u;var s=[];if(x.k){x.lR=f(y,x.l||hljs.IR,true);for(var w in x.k){if(!x.k.hasOwnProperty(w)){continue}if(x.k[w] instanceof Object){u=x.k[w]}else{u=x.k;w="keyword"}for(var r in u){if(!u.hasOwnProperty(r)){continue}x.k[r]=[w,u[r]];s.push(r)}}}if(!v){if(x.bWK){x.b="\\b("+s.join("|")+")\\s"}x.bR=f(y,x.b?x.b:"\\B|\\b");if(!x.e&&!x.eW){x.e="\\B|\\b"}if(x.e){x.eR=f(y,x.e)}}if(x.i){x.iR=f(y,x.i)}if(x.r===undefined){x.r=1}if(!x.c){x.c=[]}x.compiled=true;for(var t=0;t<x.c.length;t++){if(x.c[t]=="self"){x.c[t]=x}q(x.c[t],y,false)}if(x.starts){q(x.starts,y,false)}}for(var p in e){if(!e.hasOwnProperty(p)){continue}q(e[p].dM,e[p],true)}}function d(B,C){if(!j.called){j();j.called=true}function q(r,M){for(var L=0;L<M.c.length;L++){if((M.c[L].bR.exec(r)||[null])[0]==r){return M.c[L]}}}function v(L,r){if(D[L].e&&D[L].eR.test(r)){return 1}if(D[L].eW){var M=v(L-1,r);return M?M+1:0}return 0}function w(r,L){return L.i&&L.iR.test(r)}function K(N,O){var M=[];for(var L=0;L<N.c.length;L++){M.push(N.c[L].b)}var r=D.length-1;do{if(D[r].e){M.push(D[r].e)}r--}while(D[r+1].eW);if(N.i){M.push(N.i)}return f(O,M.join("|"),true)}function p(M,L){var N=D[D.length-1];if(!N.t){N.t=K(N,E)}N.t.lastIndex=L;var r=N.t.exec(M);return r?[M.substr(L,r.index-L),r[0],false]:[M.substr(L),"",true]}function z(N,r){var L=E.cI?r[0].toLowerCase():r[0];var M=N.k[L];if(M&&M instanceof Array){return M}return false}function F(L,P){L=m(L);if(!P.k){return L}var r="";var O=0;P.lR.lastIndex=0;var M=P.lR.exec(L);while(M){r+=L.substr(O,M.index-O);var N=z(P,M);if(N){x+=N[1];r+='<span class="'+N[0]+'">'+M[0]+"</span>"}else{r+=M[0]}O=P.lR.lastIndex;M=P.lR.exec(L)}return r+L.substr(O,L.length-O)}function J(L,M){if(M.sL&&e[M.sL]){var r=d(M.sL,L);x+=r.keyword_count;return r.value}else{return F(L,M)}}function I(M,r){var L=M.cN?'<span class="'+M.cN+'">':"";if(M.rB){y+=L;M.buffer=""}else{if(M.eB){y+=m(r)+L;M.buffer=""}else{y+=L;M.buffer=r}}D.push(M);A+=M.r}function G(N,M,Q){var R=D[D.length-1];if(Q){y+=J(R.buffer+N,R);return false}var P=q(M,R);if(P){y+=J(R.buffer+N,R);I(P,M);return P.rB}var L=v(D.length-1,M);if(L){var O=R.cN?"</span>":"";if(R.rE){y+=J(R.buffer+N,R)+O}else{if(R.eE){y+=J(R.buffer+N,R)+O+m(M)}else{y+=J(R.buffer+N+M,R)+O}}while(L>1){O=D[D.length-2].cN?"</span>":"";y+=O;L--;D.length--}var r=D[D.length-1];D.length--;D[D.length-1].buffer="";if(r.starts){I(r.starts,"")}return R.rE}if(w(M,R)){throw"Illegal"}}var E=e[B];var D=[E.dM];var A=0;var x=0;var y="";try{var s,u=0;E.dM.buffer="";do{s=p(C,u);var t=G(s[0],s[1],s[2]);u+=s[0].length;if(!t){u+=s[1].length}}while(!s[2]);if(D.length>1){throw"Illegal"}return{r:A,keyword_count:x,value:y}}catch(H){if(H=="Illegal"){return{r:0,keyword_count:0,value:m(C)}}else{throw H}}}function g(t){var p={keyword_count:0,r:0,value:m(t)};var r=p;for(var q in e){if(!e.hasOwnProperty(q)){continue}var s=d(q,t);s.language=q;if(s.keyword_count+s.r>r.keyword_count+r.r){r=s}if(s.keyword_count+s.r>p.keyword_count+p.r){r=p;p=s}}if(r.language){p.second_best=r}return p}function i(r,q,p){if(q){r=r.replace(/^((<[^>]+>|\t)+)/gm,function(t,w,v,u){return w.replace(/\t/g,q)})}if(p){r=r.replace(/\n/g,"<br>")}return r}function n(t,w,r){var x=h(t,r);var v=a(t);var y,s;if(v){y=d(v,x)}else{return}var q=c(t);if(q.length){s=document.createElement("pre");s.innerHTML=y.value;y.value=k(q,c(s),x)}y.value=i(y.value,w,r);var u=t.className;if(!u.match("(\\s|^)(language-)?"+v+"(\\s|$)")){u=u?(u+" "+v):v}if(/MSIE [678]/.test(navigator.userAgent)&&t.tagName=="CODE"&&t.parentNode.tagName=="PRE"){s=t.parentNode;var p=document.createElement("div");p.innerHTML="<pre><code>"+y.value+"</code></pre>";t=p.firstChild.firstChild;p.firstChild.cN=s.cN;s.parentNode.replaceChild(p.firstChild,s)}else{t.innerHTML=y.value}t.className=u;t.result={language:v,kw:y.keyword_count,re:y.r};if(y.second_best){t.second_best={language:y.second_best.language,kw:y.second_best.keyword_count,re:y.second_best.r}}}function o(){if(o.called){return}o.called=true;var r=document.getElementsByTagName("pre");for(var p=0;p<r.length;p++){var q=b(r[p]);if(q){n(q,hljs.tabReplace)}}}function l(){if(window.addEventListener){window.addEventListener("DOMContentLoaded",o,false);window.addEventListener("load",o,false)}else{if(window.attachEvent){window.attachEvent("onload",o)}else{window.onload=o}}}var e={};this.LANGUAGES=e;this.highlight=d;this.highlightAuto=g;this.fixMarkup=i;this.highlightBlock=n;this.initHighlighting=o;this.initHighlightingOnLoad=l;this.IR="[a-zA-Z][a-zA-Z0-9_]*";this.UIR="[a-zA-Z_][a-zA-Z0-9_]*";this.NR="\\b\\d+(\\.\\d+)?";this.CNR="\\b(0[xX][a-fA-F0-9]+|(\\d+(\\.\\d*)?|\\.\\d+)([eE][-+]?\\d+)?)";this.BNR="\\b(0b[01]+)";this.RSR="!|!=|!==|%|%=|&|&&|&=|\\*|\\*=|\\+|\\+=|,|\\.|-|-=|/|/=|:|;|<|<<|<<=|<=|=|==|===|>|>=|>>|>>=|>>>|>>>=|\\?|\\[|\\{|\\(|\\^|\\^=|\\||\\|=|\\|\\||~";this.ER="(?![\\s\\S])";this.BE={b:"\\\\.",r:0};this.ASM={cN:"string",b:"'",e:"'",i:"\\n",c:[this.BE],r:0};this.QSM={cN:"string",b:'"',e:'"',i:"\\n",c:[this.BE],r:0};this.CLCM={cN:"comment",b:"//",e:"$"};this.CBLCLM={cN:"comment",b:"/\\*",e:"\\*/"};this.HCM={cN:"comment",b:"#",e:"$"};this.NM={cN:"number",b:this.NR,r:0};this.CNM={cN:"number",b:this.CNR,r:0};this.BNM={cN:"number",b:this.BNR,r:0};this.inherit=function(r,s){var p={};for(var q in r){p[q]=r[q]}if(s){for(var q in s){p[q]=s[q]}}return p}}();hljs.LANGUAGES.cpp=function(){var a={keyword:{"false":1,"int":1,"float":1,"while":1,"private":1,"char":1,"catch":1,"export":1,virtual:1,operator:2,sizeof:2,dynamic_cast:2,typedef:2,const_cast:2,"const":1,struct:1,"for":1,static_cast:2,union:1,namespace:1,unsigned:1,"long":1,"throw":1,"volatile":2,"static":1,"protected":1,bool:1,template:1,mutable:1,"if":1,"public":1,friend:2,"do":1,"return":1,"goto":1,auto:1,"void":2,"enum":1,"else":1,"break":1,"new":1,extern:1,using:1,"true":1,"class":1,asm:1,"case":1,typeid:1,"short":1,reinterpret_cast:2,"default":1,"double":1,register:1,explicit:1,signed:1,typename:1,"try":1,"this":1,"switch":1,"continue":1,wchar_t:1,inline:1,"delete":1,alignof:1,char16_t:1,char32_t:1,constexpr:1,decltype:1,noexcept:1,nullptr:1,static_assert:1,thread_local:1,restrict:1,_Bool:1,complex:1},built_in:{std:1,string:1,cin:1,cout:1,cerr:1,clog:1,stringstream:1,istringstream:1,ostringstream:1,auto_ptr:1,deque:1,list:1,queue:1,stack:1,vector:1,map:1,set:1,bitset:1,multiset:1,multimap:1,unordered_set:1,unordered_map:1,unordered_multiset:1,unordered_multimap:1,array:1,shared_ptr:1}};return{dM:{k:a,i:"</",c:[hljs.CLCM,hljs.CBLCLM,hljs.QSM,{cN:"string",b:"'\\\\?.",e:"'",i:"."},{cN:"number",b:"\\b(\\d+(\\.\\d*)?|\\.\\d+)(u|U|l|L|ul|UL|f|F)"},hljs.CNM,{cN:"preprocessor",b:"#",e:"$"},{cN:"stl_container",b:"\\b(deque|list|queue|stack|vector|map|set|bitset|multiset|multimap|unordered_map|unordered_set|unordered_multiset|unordered_multimap|array)\\s*<",e:">",k:a,r:10,c:["self"]}]}}}();hljs.LANGUAGES.r={dM:{c:[hljs.HCM,{cN:"number",b:"\\b0[xX][0-9a-fA-F]+[Li]?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+(?:[eE][+\\-]?\\d*)?L\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+\\.(?!\\d)(?:i\\b)?",e:hljs.IMMEDIATE_RE,r:1},{cN:"number",b:"\\b\\d+(?:\\.\\d*)?(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\.\\d+(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"keyword",b:"(?:tryCatch|library|setGeneric|setGroupGeneric)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\.",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\d+(?![\\w.])",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\b(?:function)",e:hljs.IMMEDIATE_RE,r:2},{cN:"keyword",b:"(?:if|in|break|next|repeat|else|for|return|switch|while|try|stop|warning|require|attach|detach|source|setMethod|setClass)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"literal",b:"(?:NA|NA_integer_|NA_real_|NA_character_|NA_complex_)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"literal",b:"(?:NULL|TRUE|FALSE|T|F|Inf|NaN)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"identifier",b:"[a-zA-Z.][a-zA-Z0-9._]*\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"<\\-(?!\\s*\\d)",e:hljs.IMMEDIATE_RE,r:2},{cN:"operator",b:"\\->|<\\-",e:hljs.IMMEDIATE_RE,r:1},{cN:"operator",b:"%%|~",e:hljs.IMMEDIATE_RE},{cN:"operator",b:">=|<=|==|!=|\\|\\||&&|=|\\+|\\-|\\*|/|\\^|>|<|!|&|\\||\\$|:",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"%",e:"%",i:"\\n",r:1},{cN:"identifier",b:"`",e:"`",r:0},{cN:"string",b:'"',e:'"',c:[hljs.BE],r:0},{cN:"string",b:"'",e:"'",c:[hljs.BE],r:0},{cN:"paren",b:"[[({\\])}]",e:hljs.IMMEDIATE_RE,r:0}]}};
hljs.initHighlightingOnLoad();
</script>
<style type="text/css">
body, td {
font-family: sans-serif;
background-color: white;
font-size: 13px;
}
body {
max-width: 800px;
margin: auto;
padding: 1em;
line-height: 20px;
}
tt, code, pre {
font-family: 'DejaVu Sans Mono', 'Droid Sans Mono', 'Lucida Console', Consolas, Monaco, monospace;
}
h1 {
font-size:2.2em;
}
h2 {
font-size:1.8em;
}
h3 {
font-size:1.4em;
}
h4 {
font-size:1.0em;
}
h5 {
font-size:0.9em;
}
h6 {
font-size:0.8em;
}
a:visited {
color: rgb(50%, 0%, 50%);
}
pre, img {
max-width: 100%;
}
pre {
overflow-x: auto;
}
pre code {
display: block; padding: 0.5em;
}
code {
font-size: 92%;
border: 1px solid #ccc;
}
code[class] {
background-color: #F8F8F8;
}
table, td, th {
border: none;
}
blockquote {
color:#666666;
margin:0;
padding-left: 1em;
border-left: 0.5em #EEE solid;
}
hr {
height: 0px;
border-bottom: none;
border-top-width: thin;
border-top-style: dotted;
border-top-color: #999999;
}
@media print {
* {
background: transparent !important;
color: black !important;
filter:none !important;
-ms-filter: none !important;
}
body {
font-size:12pt;
max-width:100%;
}
a, a:visited {
text-decoration: underline;
}
hr {
visibility: hidden;
page-break-before: always;
}
pre, blockquote {
padding-right: 1em;
page-break-inside: avoid;
}
tr, img {
page-break-inside: avoid;
}
img {
max-width: 100% !important;
}
@page :left {
margin: 15mm 20mm 15mm 10mm;
}
@page :right {
margin: 15mm 10mm 15mm 20mm;
}
p, h2, h3 {
orphans: 3; widows: 3;
}
h2, h3 {
page-break-after: avoid;
}
}
</style>
</head>
<body>
<pre><code class="r">### IMPLEMENTING APPROXIMATE BAYESIAN INFERENCE USINGN ADAPTIVE QUADRATURE: THE AGHQ PACKAGE ###
# This script reproduces the results in that paper
# See https://github.com/awstringer1/aghq-software-paper-code for README
# Alex Stringer
# 2021/05/11
#
# IMPORTANT NOTES:
# - Use CRAN version 2.1 of aghq.
# - The script depends on external resources only through the following two lines (1024, 1025):
# cameroonBorderLL = raster::getData("GADM", country=c('CMR'), level=2)
# nigeriaBorderLL = raster::getData("GADM", country=c('NGA'), level=2)
# I have observed a small number of times that the server where this pulls data from has been down.
# If this happens, it usually comes back up in a day or so.
# - The compilation of the TMB templates creates a ton of output from the compiler. Each
# template can take a few minutes at least to compile.
#
# TIMING:
# - The most recent testing on a private Ubuntu server yielded run times as follows:
# - Example 2: <10s
# - Example 4.1: 6m24s
# - Example 4.2: 5m43s
# - Example 5.1: 47m50s
# - Example 5.2: 10m44s
# - Example 6.1: 11s
#
# - The only operation that takes a long time is the MCMC within Example 5.1.
# To skip Example 5.1, set dofast = TRUE
dofast = FALSE
#
# VARIABLES TO SET:
# Example 5.2:
# Set the resolution for the spatial interpolations.
# The results shown in the paper use:
# reslist <- list(nrow = 200,ncol = 400)
# but this takes a couple hours. Here I set:
reslist <- list(nrow = 50,ncol = 100)
# which is faster but produces a blockier map.
# So some plots show up properly in the spinning:
par(mar = c(5,5,0,0))
## Check for missing packages ----
# The astro example requires ipoptr, and IPOPT (https://coin-or.github.io/Ipopt/INSTALL.html).
# This is a laborious installation.
# If you want to not run the astro example set doastro = FALSE
doastro <- TRUE
if (!('ipoptr' %in% installed.packages()[ ,'Package'])) {
warning("No installation of ipoptr found. Skipping the astro example.\n")
doastro <- FALSE
}
# Obtain the disease example data
# If can't obtain, then skip that example
dodisease <- TRUE
if (dodisease) {
# Get the Tomato disease data
if (require('EpiILMCT')) {
data("tswv", package = "EpiILMCT")
} else {
download.file("https://github.com/waleedalmutiry/EpiILMCT/blob/master/data/tswv.RData?raw=true",destfile = file.path(globalpath,"tomato.RData"))
load(file.path(globalpath,'tomato.RData'))
}
if (!exists('tswv')) {
warning('You asked to do the disease example but the data cannot be obtained. This is either because you do not have the EpiILMCT package installed or because something is preventing downloading the data from github on https using download.file(). In any event, skipping the disease example.\n')
dodisease <- FALSE
}
}
</code></pre>
<pre><code>## Loading required package: EpiILMCT
</code></pre>
<pre><code>## Loading required package: coda
</code></pre>
<pre><code>## Loading required package: parallel
</code></pre>
<pre><code class="r"># See if INLA is installed.
# If INLA is not installed then cannot compare to it in the loaloa example so skip it.
doloaloa <- TRUE
if (!('INLA' %in% installed.packages()[ ,'Package'])) {
warning("No installation of INLA found. Skipping the loaloa example.\n")
doloaloa <- FALSE
}
# But, don't do loaloa if the dofast flag is TRUE
if (dofast) doloaloa <- FALSE
## Install and load packages ----
# Check if installed
neededpackages <- c(
'aghq',
'TMB',
'tmbstan',
'parallel',
'glmmTMB',
'geostatsp',
'PrevMap',
'geoR',
'trustOptim',
'numDeriv'
)
missingpackages <- setdiff(neededpackages,installed.packages()[ ,'Package'])
if (length(missingpackages)) stop(paste0("Missing the following packages, please install: ",missingpackages,"\n"))
# Code for users to run interactively, if they like
if (FALSE) {
install.packages(missingpackages)
}
library(aghq)
library(TMB)
precompile()
</code></pre>
<pre><code>## Removing: libTMB.cpp libTMB.o libTMBdbg.cpp libTMBomp.cpp
</code></pre>
<pre><code>## Precompilation sources generated
</code></pre>
<pre><code class="r">library(tmbstan)
</code></pre>
<pre><code>## Loading required package: rstan
</code></pre>
<pre><code>## Loading required package: StanHeaders
</code></pre>
<pre><code>## Loading required package: ggplot2
</code></pre>
<pre><code>## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
</code></pre>
<pre><code>## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
</code></pre>
<pre><code>##
## Attaching package: 'rstan'
</code></pre>
<pre><code>## The following object is masked from 'package:coda':
##
## traceplot
</code></pre>
<pre><code class="r">library(parallel)
options(mc.cores = parallel::detectCores())
library(Matrix)
library(glmmTMB)
library(geostatsp)
</code></pre>
<pre><code>## Loading required package: raster
</code></pre>
<pre><code>## Loading required package: sp
</code></pre>
<pre><code>##
## Attaching package: 'raster'
</code></pre>
<pre><code>## The following object is masked from 'package:rstan':
##
## extract
</code></pre>
<pre><code class="r">library(PrevMap)
</code></pre>
<pre><code>## Loading required package: maxLik
</code></pre>
<pre><code>## Loading required package: miscTools
</code></pre>
<pre><code>##
## Please cite the 'maxLik' package as:
## Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
##
## If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
## https://r-forge.r-project.org/projects/maxlik/
</code></pre>
<pre><code>##
## Attaching package: 'maxLik'
</code></pre>
<pre><code>## The following object is masked from 'package:raster':
##
## maxValue
</code></pre>
<pre><code>## Loading required package: pdist
</code></pre>
<pre><code class="r">library(geoR)
</code></pre>
<pre><code>## Warning: no DISPLAY variable so Tk is not available
</code></pre>
<pre><code>## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.8-1 (built on 2020-02-08) is now loaded
## --------------------------------------------------------------
</code></pre>
<pre><code>##
## Attaching package: 'geoR'
</code></pre>
<pre><code>## The following objects are masked from 'package:geostatsp':
##
## matern, variog
</code></pre>
<pre><code class="r">## Set up the directory structure ----
# Each example gets its own directory within the tempdir()
globalpath <- tempdir()
## Example 2: Basic Use ----
plotpath <- file.path(globalpath,"basic-use")
if (!dir.exists(plotpath)) dir.create(plotpath)
set.seed(84343124)
y <- rpois(10,5) # True lambda = 5, n = 10
# Define the log posterior, log(pi(theta,y)) here
logpithetay <- function(theta,y) {
sum(y) * theta - (length(y) + 1) * exp(theta) - sum(lgamma(y+1)) + theta
}
objfunc <- function(x) logpithetay(x,y)
objfuncgrad <- function(x) numDeriv::grad(objfunc,x)
objfunchess <- function(x) numDeriv::hessian(objfunc,x)
# Now create the list to pass to aghq()
funlist <- list(
fn = objfunc,
gr = objfuncgrad,
he = objfunchess
)
# AGHQ with k = 3
# Use theta = 0 as a starting value
thequadrature <- aghq::aghq(ff = funlist,k = 3,startingvalue = 0)
summary(thequadrature)
</code></pre>
<pre><code>## AGHQ on a 1 dimensional posterior with 3 quadrature points
##
## The posterior mode is: 1.493925
##
## The log of the normalizing constant/marginal likelihood is: -23.32123
##
## The posterior Hessian at the mode is:
## [,1]
## [1,] 49
##
## The covariance matrix used for the quadrature is...
## [,1]
## [1,] 0.02040816
##
## ...and its Cholesky is:
## [,1]
## [1,] 0.1428571
##
## Here are some moments and quantiles for theta:
##
## mean median mode sd 2.5% 97.5%
## theta1 1.483742 1.482532 1.493925 0.1424943 1.204135 1.762909
</code></pre>
<pre><code class="r">plot(thequadrature)
</code></pre>
<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////isF19AAAACXBIWXMAAAsSAAALEgHS3X78AAAgAElEQVR4nO2dCXgVRbaAexDm4ZjH5GlUQMeZhzrPBUeR55iQDcISwIR9CbJvgoAKCIMiAuIIAmEVhcQl6IgsYRN1BMJOwg4uLEnYZN+EBEPIdpPU3CXr7b7Jrb5Vdaq66v8+Q+zue84hP7l9u7vqlIbIkNAvjkvmatoQrf9YbVANbdKYZzrHfVRLmzPTvi0ubowWMD8uQmtQfmwd7f4/1/+d1sq+6+759i93xY3W7o+Na/381Lg/aeMrx3VsmK9pI+K6ao/FddaaxcU10EY6t5ZmsR/UTOvG9C+LQT+NlPgEQoEIc1PTTv5+2CztdE3tCEqd0uEhTbt+RbvLvme+9hpCP2iB5cfW1ezUH5mLPtCG2/+3nnbpV3+tRtCUqwg9ox2oHNexIVv7fTHaYY/QT3MS6zqsJIv9oBHaXIZ/VSwSJBB/4fnGMfchu/idtfx6x9dxiA+w75mhjUbocGXx21zffKCNsH99QDuHbkx4TNP8T3gSXxuhZHuEgdrLy+z86NxamgUp8ZA4xL9W84Foh/iR2gT7r3ip+K+1x/PQeLu24itXip3Hlonfqv0lG+3U7kbrRiSi9CDtA7vRvc7jyo51bCgTH6uNQmjL4hPOraVZkBIPiUP8Uk17zyF+hnZ353q/0y66xBc8rP01WLNru+6yVEE8CtcadPTT5qON2h+69/DTdqAm2ovHHceVHevYUCY+4+5a/Xvecc8159bSLEiJh8Qh/hdN2+QQf7uL32NfdPtLgks8Ot2qzpNxxuKzXn34rkaL7d989dwf73rmC4SW3XfHmkriHRvKxKOjrfzv6ZDq2lqaBSnxfJL5YzpCe7UXMF92dxaValgjsfj0mndM++hxLQHvVT89QKUY5kgsHm0I9fdr9Anea24230GnGNbILF5qlHhJUeIlRYmXFCVeUpR4SVHiPfLbpn8t2ZMPXQUtlHgPbI1uPTHuozEhA9OgK6GDEm/IpU5Dz7q+2/fC2FzYWuigxBuxI3BP+f98EXYOrhJqKPEGrGl5o+L/Hg48BlUJPZR4Peva5lTecCbwOEwlFFHidexuftt908mgaxCV0ESJd+dioIHk3S0L2FdCFSXeDVvLg0ab48ewLoQySrwb73xgvP3FjWzroI0SX5mDLxQb78gMvMm2Esoo8ZWwhXq8Zv92CMtCqKPEVyJ2nud9MbvY1UEfJb4il0ILPe+8EFbErhLqKPEV6b+9qr1TMQdmco0SX4Efu1S5Ozcwm1EhDFDiKxBdzSPYz//Jpg4WKPHl7BxUzQGFwRlMCmGBEl9O5JnqjkicwKIOJijxZeys/kK9OCSTQSFMUOLLiDpd/THLptCvgw1KfCmHentxUGGQ7pGtoCjxpbz4szdHLVpAuw5GKPElnI3y6rCcwCru7YmEEl/CmA3eHTdxDd06WKHEu8gO9vA41p3LkXQLYYUS72LRQm+P7OXVZwHuUeJdhNzy9shdQ2nWwQwl3snOEd4fG2qJ7kdKvJNeR7w/Ns7rswLPKPEObjTHOPhWOK0yWKLEO5i3GOfowQeqP4Z7lHgHIVgjLPa8TKsOhijxdvYPxjs+OKf6Y3hHibczAnP87MyldOpgiSnxNwzeGQUWnx+E+YLL3t3X5xpc8VG56HjjGrUiL+sCJRCqiD1rsIfSRV+iUQdTcMVr2ShiXG7+5M66QAmEKmJP52qHXLnz1SwadTDFhHh/G0JF9+kCJRCqiDkZrbBfcjuMQh1swRZ/Cj17EqGj9XWBEghVxJyP4/Bf0+co+TrYgis+7EE//0i0JSBWFyiBUEXMaXWj+mPc+f4t8nWwBf9TfV7qHpSyWR8ogUQ5AFxsb+JFNtwLAe4weR2fuV4XKMHXUoCYs8TMq17eR7oOxpgUn1K77Nvk9520GUuqJMY0M/WUdfvrpOtgjO937s4kOWnTi0Q57DlT9URJTxQFejlUi1fwxbvmkuj6ggwcSKAaAGKXm3vdq7vJ1sEaXPHHnqjRYLX9E57uZaKKb2Zy6vMOwd/rccWHzyjYUjfZOuLP6W5BeklRkNjv9bji/YsQWvN4gWXEzzX1md7B8P0k62AOrvhHHE8wOw+3jPiWv5l95eY3SdbBHFzxiXUibqCM5xpbRPzVdqZfamtCsA72YH+qv7TG/mmoIHGc+3Yxxcf70M+oH8bIXP4gNgJHTPFRV82/9muhO+LILf43nGHV7uSEkyoDArnFL9c9ZMSh03lSdQAgt/ieJ3x59eIPSdUBgNTiC4J9evmvLxCqAwKpxW95w7fXtxB4+qTU4kel+Pb6mYlk6oBAavFNfOxnk9qXSBkgyCzed29NxG1kLrP4WJ/fqUeK+1BeZvGtTD+gKSVJ3N62Eou/6Xv/qnxxJ1ZILH6FT7ftXHS66HsMGCQW3z/V9xgff+x7DBjkFV9M4nn6RbMjt8CRV/zBYSSihIm65qy84t/7hkSUN7aRiAKAvOJbEFlSartuKJIgSCs+sw2RMAWiXtBJKz6RUFOLToJ2RZFW/EBCQyXjEsjEYY204gMJxTnTg1Agxsgq/ugAUpGCxXxCJ6v42SYnyep5RcwWCbKKf+E6qUjfiDm8XlLxeeHEQt1qQSwUSyQVv4lg16qWQg65lFT8OIJ3Wsnc+2WNpOLD8snF2vcKuVjskFP8ryTbTxf6Ni0DCDnFL59DMlqXCySjMUJO8YMPk4y2CGtFG06QUzzZxkUnRGzxJ6X4kz3JxhOxAZaU4kk/URskYFMUKcV3O0c23tJ5ZOOxQEbxRaSvv66a6XwOjIzifyC+YGCojXRE6sgoftZK0hFH7iEdkToyio8m9ki2lHVTSUekjoTiKQyM/a018ZC0kVB8CoWG483yyMeki4Ti//kt+ZhvCjehRkLxkT73Q9CzcSL5mHSRT3xuMwpBb/vSGxUE+cRvpdJnvsVtGlEpIp/4SRvFiUoR+cTT+d3cJtp6FdKJz6EzGprKJweaSCd+09t04rYSbJC1dOLf0q+HTIQp39OJSwvpxDfPoRN35z/oxKWFbOJv05rwlN+UUmBKyCY+aRKtyK3FOsnLJv6tLbQiv/tvWpGpYEJ8cZbRoFJBxNM6xQt3kscVnzfx0VpazUcm6eaeiSGe2inefpIX60oeV3yfDskZtozdMf3dd4ghfhPFp2iRQp3kccUH5Dr/KHzIfYcY4t+mdBXvQKwreVzxT691/rGtkfsOMcTTfIi2Xajb9bji99V7qvugmEb1dYunCyE+h+Zj89wIisGJg/2p3pYUPz0+ST+OXAjxWwl2QNHTkkh3XEaYvI7PXO++RQjxkzfQjD4xiWZ0wpgUn1K77NvPGzu5V4RpRK1u0Yy+mdKDPypIdeeO8qU2xZsE5MEXn+n8etN9swjik8fSjd88l258kuCKP/ZEjQarEcrTvUwE8e9RGFFfEZFG1+OKD59RsKVusqDi2+rep8jy/RS68UmCK96/CKE1jxcIKd5GezWJLIGm0OGKf2SX/Uvn4UKK3/ca7QxNxVmTCld8Yp2IGyjjucYiio9dTTvDaHEWGcb+VH9pTTZCBYm6xZcEEN/+Gu0Ma6fTzkAMia7ji0Kop7geTT0FKSQS//NQ+jlChVmnRCLxHy6hn2PYD/RzkEEi8d3P08+xbD79HGSQSDyL7uIXuzJIQgR5xJ/oyyJLsChtbeUR/9knLLL0T2eRhQDyiGejhM0/LwLII57NAiIn+zBJ4zvSiGf1sUuUBWqkEb+c0YVWd0EWqJFG/HBGt1YWLGWTx1ekER9WyCbPT8R7otNBFvGZLzBKVBTKKJGPyCL+22msMkXfYJXJJ2QRP24nq0zvr2OVySdkER/BrK94CuUx3ISQRDyltoZG5Isxd1IS8duozpasTEsh+hlLIv5dhk0LJlDrr0QSScS3obA4gSc2vMMul3nkEG8LZ5gsqw3DZKaRQ/yBV1lmCxdh+UE5xM9dwTLbyH0ss5lEDvGdL7PMtnI2y2wmkUJ8Mf2pFBW50olpOnNIIT5d146RLiECjLiUQvynn7HNNzCNbT4zSCG+33G2+RI+ZpvPDFKIZ3uKR+hkX8YJTSCD+EtdWGcUYMSlDOJXzGWdsetF1hmxkUH8qwdZZ5y/jHVGbGQQH85onGU5h0awzoiNBOJ/a8s8ZSHt/lq+I4H4799ln/MFyh31fEcC8W9tZZ/zve/Y58RDAvEtqS085Znt49nnxMP64vMgBj8yHNxpEuuLTwZZD645s+HcJrG+eJgJDm8wm8BhEuuLh5nSxG7KlkksLx5oEmNGFEha77G8+J+Bpi2Hcd7j0vLiWfSzNGLYjzB5vcXy4nucg8m7dAFMXm+xvHioR+MXugMl9hKri/+lF1Rm1sN+MLG6+C8WQWXufQoqs1dYXfxLR6Eyxy+GyuwVVhcPN8Q9dRBUZq+wuPhrHcBSM56+g4vFxa+Ohcvd6Qpc7uqxuPjRe+Byz14Jl7t6LC4ecgXA/SPhclePKfE3svXbuBR/KxIwOdM2HNjgio/KRccb16gVqZtwzqX4jZMhs7NsvIONofhRyR4fLWnZKGJcbv7kzu47uBT/9mbI7FPWQ2avBkPxExrWHZpk3MjFLt7fvqfoPvcdXIpvBdpybivD5nrYeHirPxEbEtD361z9Du0UevYkQkfru+/gUXx+M9D0XI+49CA+++sh9/5fyD36nkFhD/r5R6ItAbrrYx7FQ7eVZddAFx9D8XNa/neLufbf60MPGOzMS92DUvTnTh7FQzeSZtcyGx9D8T1XZdm/5iPbWk8vyyz/3JK6wknEixSq85GoDNj8374Hm78qDMU3dHzJe7CKl6XULvv2QJyTYP4GHhRBT11ktiyGCQzE16yp1XSA17SLw7f6H4ZDV8BqIRwTGP7GR5sIxKF4+PYErzBvyeA1Vr5XD9+QJHEOdAUeMRDvt7OhC6Pj00px38Gf+GL4FkRXdDc4ucFA/PqMIy6Mjo/U7qzrxH0Hf+LTGPezNILfHpfGb/XpqGDVZ8ZPNIcMNQ7En/iPE6ArQGjQMegKPGEo/u0/FL7/xPMvGb5g23TjQPyJ52GYK9wg3+owFH/PafTgnkzdc5gq4U88/CkeoV96QlfgCUPxd2ccqFd0yx8rEHfiz/SArsABD//6DDEU/9LTD8+4Fog30Zc78f9aCF2Bg14cnG8MMRRfmLjUdul9vPEj3IkfDDaVoiLcTquw7g0cPq6keLimNMRQ/Pogx/2bV7AC8Sb+Mifrg/B6kjcU3+Dz1LS0tHSsQLyJX868ZbUx3S5AV2CMofj2JgLxJn74D9AVuFjwFXQFxhiKn7ERPxBv4kM56UHz8xDoCowxFB/+X/We9PCQxiOcib/eDrqCEnidO2kovoqHNB7hTPzqmdAVlNLxKnQFhni4nCu+jhuIM/H8LPM5JxG6AkMMxZ+P8GuYFnIaKxBn4sO4WdiX09UqDMV3fSO/YdFUvOkAfInnaJhjIUxrzeowFH+/DTVEtruxAvElfh1HrWTbYZ82WWAo/qlku/hDT2AF4kv867ugKygndjV0BUYYit8c0COgfwDeNBS+xIcDNkRwZ/9r0BUYYfyp/tqnkxditgLlSvzNNtAVVIDPJams+XSOr7lLUSAd86vBSPyh7o/6PRqDea+bK/FjU6ArqMgMj1MQATEQv/WPb+1I2zHBfxtWIK7EQ/Y80rOXxy5IBuKfd3V4XxqEFYgn8b/xdIpHyMbjSd5A/J2unlY5d2EF4kn8d1yd4vk8yRvNli35s7ZuT1XwJH4MV6d4+0l+DXQFegzE31HycO73WIF4Es/TVbyDfRxeyRuI9y8FKxBH4rm6infA4+16K17Hr5sKXYE70fzdrrei+NG7oStwJ3YVdAU6rCien2fxpRzk75m8BcXf4OdZfCmF/A28s6D41TOgK9DTgbtFCywo/pX90BXomQveh8kdC4oP4bDH2E/cja63nvircOsPeYaDRkxuWE/8Mk4mzVWm23noCtywnviXfoauwIiFn0NX4Ib1xAdzMS/enfS+0BW4YTnx5/hrpuykCXQBblhOfEIcdAXG9DkOXUFlLCe+9wnoCoxZzFnHO8uJ5+0ttZTzXaErqIzVxKfy2mwIBXPSqaEEq4lf8CV0BZ54+RB0BZWwmvhOuiUweWElN60anFhMPIfPP0vh7GmxxcTvfRW6As+E50NXUBGLiZ8KvNJcVYzDm5pEGYuJb8XxAs5Jb0NXUBFrib8dAV1BFeTArnTrhrXEr58EXUFVtL4JXUEFrCV+DMeLuSL0Pk8zqawlPpSzuVOVOTgMuoIKWEr8FTPdl9lRxNNzBEuJXzIfuoKqifkFuoJyLCW+r279S774NB66gnJMiC/OMhrcxIH4YrweHuw5x9GjWVzxeRMfraXVfGSS7vYjB+IPD4KuoDqC+Rnyjyu+T4fkDFvG7hjdY28OxM9aAV1Bdby6B7qCMnDFB+Q6/yh8yH0HB+LbcthqpjLfTYGuoAxc8U+7erZta+S+A158Lle3RA3Jbg5dQRm44vfVe6r7oJhG9XUTE+HFb5gIXUH18HPXFvtTvS0pfnp8kr71ALz40cnQFVTPDG5aY5i8js9c774FXnwod40w9PxovDY7ACbFp5Q3wVvSwskD0LdLeRu/bEgxN3dtrXPn7pNPgAvwigGp0BWUYB3xXXmbiGzI8jnQFZRgGfE2DpsIGpDRFrqCEnDFp5XivgNa/M6xsPm9pflt6Apc4IqP1O6s68R9B7T48Vtg83vLu99BV+AC+61+yFDj7dDiw7gefFPOAU56HWKL3zbdeDuw+EudQNN7Dy8XdFb5cMfTGIeq6cfHaBGriO9yATQ9BstnQ1fgxCLiC8S4mHOQ2Rq6AicWEb95PGR2PFpmQVfgwCLiXxfgyVwp07l4QmcR8RwNZquWI9A3t51YQ/yJXoDJsQnioRuONcTPWQqYHJtX9kFXgKwivnUmYHJs1vMwRswS4m+2gsttgjwelhy1hPjls+Bym6EzB0MHLCG+F2d9YqsjYSF0BdYQb+Nu+YdquBYFXYE1xG/7B1hqkzS/BV2BJcSPEui2nYvpK6ErsIR4Lm6IYJHWB7oCK4j/eQBUZvM0AZ/8YQHx/+S4m6UnxoGPELSA+PAcqMzm2QXec1d88ec6AyX2haJA6LWyxBc//19AiX1i6EHgAsQXHynUA5pSNkwALkB48df56v/vLQXQKyoIL/5TTteZq44+wKOshRff7gpMXl9ZOxU2v+jif2sJktZ3csJh84su/st5IGkJ0P00aHrRxXfiYEyDOVbALkcmuPhbPK9FUjXZsG35BBe/lI+JaKbochYyu+DiO5+DyEqGZbGQ2cUWnyXuO739vb4pZHaxxS/hpYeUKbpCfq4XW3wHYWbFG5HoobkIE4QWnynq3RsXOZATK4QW/9lH7HOSpCfg/Xqhxbe5xj4nSdZNhsstsvgrYj6RLScfcPUkkcXPW8w8JWFeOgSWWmTxzTheM9w7tr8Ollpg8SdjWGckThHcXBCBxb/zDeuM5Bm3CSqzwOKbCNK9tioO94PKLK743cMZJ6RCGFQXc3HFD+Nn1UYfmL0EKLGw4vM56QLtI5eh7kUIK37le2zz0aId0HMmYcW3F3gIRkVWTIPJK6r4y7ws6uMreUCnLFHFz/yKaTqKDN8FklZU8U1ymaajyMFBIGkFFZ8yjGU2uoRlQ2QVVPxA6OnlBPngM4isYorPCmeYjDYZzSGyiik+bhHDZNTpfQQgqZjiw4V/El+RHRCdkEyJv2HwcYSl+AMwH4SpEQzQtwtXfFQuOt64Rq3Iy+47WIoftJ9dLhbMTWCfE1e8lo0ixuXmT9b1GGMoPhN07hEFMgEmzpoQ729DqOg+9x0Mxc/7lFkqRgw4wDwltvhT6NmTCB2t776DnfjiJgK2sqyag/2Zp8QVH/agn38k2hKgm+LLTnzSWFaZ2BFxnXVG/E/1eal7UMpm3WZ24jvCNo+hwhLm8ydNXsdnri/79toBJ+37EqqoOk51ZJSIJfmBrJfKNCk+pXbZtxvHOfkbq17Co8BGJNNkCus1K4S7c5cVDN33mQpXWzBOaEJ8cZbRj56V+LmfsMnDmgGMb0rhis+b+GgtreYjk/LddzASXxhomREYlTncg20+XPF9OiRn2DJ2x+guPBmJXzGFSRoAos4wTYcrPsD1C1f4kPsORuKbMr/gZcWmkUzT4Yp/eq3zj22N3HewEb/1FRZZYAi7wTIbrvh99Z7qPiimUX3dRxE24qN+YZEFhmVMz2LYn+ptSfHT45P0y6YxEX+oF4MkUBQGsRx1KdZ1fHeIQUrMiJ/LMJlQ4o91oZ8DkPzAPHbJhBLf20KDqo2Yz7Bvn0ji09tTTwFLbqDuvhg1RBLf22JD7fTMXcgslUDiUzvQzgBOLruzvEDie/xAOwM8H85nlUkc8T92p5yAB/IDWV3LiyO+XSrlBFywmFWHF2HEb7fY7BkPFAYzegolivjiZkIvRuE934xik0cU8cugl91mRuQpJmkEEZ8bmEUzPE8c6MYkjSDip8XTjM4X/XawyCKG+EthrIedA3KRyV9WDPF9mPwS8MK0OAZJhBCf3JNebA7Jb8JgEJYI4m0hl6jF5pINg+nnEEF87DxqoTklhn63SwHEn2kq0Sc7F5dCqD+YF0B8FNwaXWAsepd2Bv7Ff/4GpcA8Uxx5lHIG7sVftE67YhxONNOPYCcK7+KL21li6Rl85lN+Psu7+HhpHs64Udya7kcbzsWnN7XA4nLmOBdMdWUyvsXnhx+nEFUQEofQjM63+JGfUwgqDINXUAzOtfjVYAtwckFO6Al6wXkWnx4Ctf4mJ6SH0vsBcCz+VrDEJ3gXa+g9luRXfFGX7whHFJC3Z9CKzK/48dT+zgJR1OUbSpG5FZ8wgGw8QckOo3Qfh1fxG9pKe+emMpcC6bRB41T8vnBphlNXx9Em12iE5VP8kaCrBKMJzu6wmxSicik+PdAia4STYXNzCm9/PIpPC2Qzi0gY1rcgb55D8YcDLbgEhW9saJZBOiR/4ncHq/d5HdtDdQv9+Qh34r9pSuVDrOgcCkojG5A38Qs6gKymzj+nm2wjGo8v8QUvj5JuDL23ZLQluoQ2V+IvtkjwPYhlsY0cTHDAMU/i1wdZvGWprywLSycWix/xeaN7WGpVeBocDye2FBM34vcFf0GmEEuT/2b7i2QicSL+1ugOknS18pW9wR8VkYjDhfjiZYE0B5Rai4KpTfcSCMOD+F3N37xFqAopONuj1y8+B4EXf6hD37OEapCGPa1G+HpmhBa/s13PY4QqkIqkFkN9G3QPKj7vy6YvU5wzYG12duy0yYf1lQEXFT48Omgm0zX2rMbJ14Omm24LBbWo8OnpoX23WnJBcJbkLYtu86m5Xx6IRYUL97wdGpOYg5dYYcy1Ra0jZ5u4k8t6UeGiwws6B4/cwHCBNeuTuaxf4IAvMEdhs1xU+OzaCW2Chy4hdM9RUZFjC3sGd3o/yftVDpgsKlxw/N+zhzQL7Rm7VY2Wp8iVb97pGBI5ctGWc17c1KW5qHDW8eRV89/sHRHWYuC0VcfUzBgmZB9YMqln0/CWAybFrdt7zvMHKZPX8Znr3be4i98X3jSq//h5iSnn1ZgaAPJObP9y1ujercPCwiO6DBo9eWac260+k+JTapd9+/VLTh7v5EudCnrYfj15YMu6ZW4frXy/c5d1ysk0oiPCFLQhducuIcHnWhQMIXbnTokXC2J37pR4sSB2506JFwtid+6UeLEgdudOiRcLYnfulHixIDYCR4kXCyVeUpR4SVHiJUWJlxQlXlLIie8XV5kOvcjSI4pwwC6dCAds341wwLY9CQfsWq6nHynxF928x9UjXHT0XwkHDH2OcMC/tSQc8KGuhAPWL9cTT0q8jqaE4x0dRjjg8g8JB5xEtj0RQp29HzvpHRWdKPHEUOJ1SUigxPuOEu9AiaeDEu8zSrwDJd53mIj/inC8zH8TDphGuqvedtL9m1bqRjb6SEUn1MQr+EaJlxQlXlKUeElR4iVFiZcUJV5SlHhJUeIlhYb4v5esmLS/0b19ifS3Kg248lH/ZqkkAyJ0pA7JeGea+zc5TjLgmkf/GE1gudXyH1yZE/Litw/WXFXb6q3JaT+RYMDLdVKKpj9JMCBChc/VrvJQzHgNv8yfQOAhRVnA63d+ndnF9wXVy39w5U7Ii58+9E5X1UlPIZT8V4IB10QgVFAjk1xAhKa9SEB8Wbxdz9p/tATWRy0LuLeu/bf1OZ/jlf/gyp3QeKuv66o6vjtCGbVI9C0tCZh9A6EtDQjEKw2Ijj15ioD4snifRb34v9FEFkgt/Rvf9+X56LE+Ryv/wZU7oSh++iD7v3+NRGu7umWn5LX3ryIQrzRgUdD2KyTFT6+5NuutvxMMiBZote/9lUTAkh9cuROav/Ex9n9dNQn+xqOMLo8lEwhXFnDWCERU/KJIhPLuILEAcEnA7U+k5cc943u4sh9cuROK4pMaIbT7EYIBC/4+jFCbxJKAvfz8/qD57SIW7/tWCOXXJPgeN2kkQsW1fP6VL//BlTuhJj4l01Z/W2GPSQQDJjbOs0MwoP0Lwd/4lMyCepuL3wknGPC7R368Pf9hn6OV/uAqOqEmvvZ6tP+ZP5G5ji8J+Ibm4Ca5gIioeHu83U8HRBJZYac04Ly/+Dc97HO00h9cRSfqzp2kKPGSosRLihIvKUq8pCjxkqLES4oSLylKvKQo8ZKixEuKEi8pSrykKPGSosRLihIvKUq8pCjxkqLES4oSLylKvBs10ZGGbhvKKZ9sKTxKvBtViK8w2VJ8lPjKRGsNDz3xxkN/3ojQlr/d3+2CY0P+vAf8nk+vNNlSfJR4N+y/8dqi4onB6NeAHbZxEY4NF2v/nPvSEMfOukq8ZbGLv6cY/dQQfdoDoXy/YvuGvEsoZ3SMY6cSb13s4p9EjvP85ICGdq7bNxRNafz/4Uq8xXF9uLP/t7A/QkVpjg0rnjqPPlfiLR+l/VMAAABuSURBVE7NvBLxZwN+sE0NcWyYF1l0rmmUY6cSb106PnDIJR59/dj/RJx2bLgcfv/zqxqtRkq8QnyUeElR4iVFiZcUJV5SlHhJUeIlRYmXFCVeUpR4SVHiJUWJlxQlXlKUeElR4iVFiZcUJV5S/gPxGG32hE2rTwAAAABJRU5ErkJggg==" alt="plot of chunk unnamed-chunk-1"/><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAAC+lBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxsdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7///+HlbSeAAAACXBIWXMAAAsSAAALEgHS3X78AAAaQElEQVR4nO2dCXQUVbrHyzExKo9hNDMCysxD8LmAkKdINghZMC5BIhASEJAlAoM8EZeBh5DgSgL6FJQZEoGIojKIAjIq2ooEE8ISEZRVJyAiBgQ6CCQk6aTvOa+XdCfpqiR9q27VV3Xv9/MknVO3+Opv/3qpulV1r0TYUDA2z7S8IkmTpHFPSpm/k7Kf6D007++h0ssLXMvy8p6QwhflJUhd/as+JF3n+7O9NN31e4zUPe8KyY2/Ie86aVbeIkmampcm3Zg3VIrPy+sqPepZ6qvuWileGm7k/yMtYyVW4gsYFdKBs5L078umvCQdDpH2kgPPpHaRpNMnpCtdLYukaYR8I0X6V/1Kus73Z0fpc9fvF6V00kHa1qxeb6mUXJAuc5Itrn851vOqkF70LPVVd600VXrFoP89VRQIIv7nvrdl/Im4xH8V2m50fnu3+HBXy3zpMUK+ayL+pHTJ966HjbdnNohPl2a1ID6MkCLXv5wg/XWVi92epb7qBMWbAbf4aSHXDnKLf1Sa7XqL+8Svl26qJrNc+pwnTjg96w6R4s8Se7T0hFf8F5dd8q1PvH+d3tL2RvEvStMJ2fTGD56lvuoExZsBt/h3Jel5t/j50lVDO10iHfeKr71euiFGcuk77bVFyL87SaE3XCH9+bRL/A23d5WkUcQn3r9OtDTye794+1Wh4x649OpfPUt91QmKNwNu8Uck1xvYJb5yWLsb3xz+lwKveHL4zva35DUVT05l3nzFfz1ywv0dL0lX9nqmTi5+1Z8uXesXT/bd2eHq1APepb7qBMWbm4rdhwjZLt0b5OpXndM1jZGoEO8851QqVKA5CwCHQi6d9/ebpILg1t5zra5hDIVWfHVW91AppFt2jaxQAaNExvJpvw7tIpYGt+7ZxC36hjESWvFjUovsDntJxjhZoQJGiRBDoBUfftHzUNdFVqiARRzEKGjF91rnedgcIStUwCIOYhS04nd06pmemRHReaesUAGjRIghUO/VO2z5ufk2h7xQAZM8iEGoPI6v2CgrVKA1CmIkKsUXh/n/3DDRw+1TWEVCVFNnt/9aVla2p7S0tMTmZsPqBgryjzdfVXvP3W9lHobLDvAQHakt37/1o3fycmZNm5g2KCk2KSYp2vVfUkpa2mjXu3D6DBc5HhZ7T78vXcpcfAMTJjAqhLTMmT0b8rIn3jcgrl/C8Mkzc/JWffTV12XlFWpKoXhLcNS2ePqg/v3vf2TBW1/ut7OoSCv+oI/ABhSvD+e2vJrZv/+IZ1eVnmVal1Z8snR5Rw+BDSieOc69eWNj7/7f1QfrdChO/VE/abLychTPloOvDYmduOKQbvWpxW/OVV6O4tlR/cmU6LErf9F1G7hzZzYurh3ZL7tU6ZIHpqB4U+Hckhnz/A9GbAnFm4jTC6KnfWPQtlC8adgzPvHtasO2huJNwpcpD5QauT0Ubwo+T5pyxNgtongTsD354Z+N3iaKB+dIxqjDxm8VxQNTlZ2wA2K7KB6WjZErdO+rUQTFQ3JmzIQzQJtG8YD8K9IGtm0UD0bV5HGA92CieCj2Rq+B3DyKB2JF/I+g20fxINROeVh2v7GxoHgIfh1YAB0BxQOwL2ordAQUD8CmWNivdw8o3nBWJau6A4IxKN5oFqcZd7VFK6B4g5k3QY+r5OlB8caS9QjMORkZKN5Q5jwOncAHijeSuY9BJ/CD4g1kwVSTfM4TFG8k+WPN4x3FG8eaofIho+BA8UZRmFQFHaEpKN4g9scwGciCGSjeGE5GmaB/viko3hCqE7e1vZKhoHhDeHAVdIJAULwRvPgUdAIZKN4APk+th44gA8Xrz9FotiOVMQHF6051/HfQERRA8bozdQV0AiVQvN6sMecTg+J15khsJXQERVC8vjgS9kJHUAbF60v2YugELYDidaV4sIlOwTcDxevJ+agT0BFaAsXryaS10AlaBMXryKdjoBO0DIrXj98izXCvVAugeP2Y+C/oBK2A4nXji9HQCVoDxetFZeQp6AitgeL14m9vQydoFRSvE3tSoBO0DorXh/p4gIGJaUDx+rDkeegEbYDideHXaODRzNoExevC+C+gE7QFiteDbRnQCdoExetA/YBj0BHaBMXrwLLnoBO0DYpnz2+RF6EjtA2KZ8/fQMcjDxIUz5yyO6ETBAOKZ85wQ2eMVAuKZ81WE1920wR14hVOOKJ4L87Eo9ARgoJW/MlxfZ481Tv0etltAijey/szoRMEB634lGFrBv0x3/lKQmADivfgMOMt0UrQim9/lvx0eQ2p/ENgA4r3sOQl6ARBQiu+23rylnSQ7O4a2IDi3VRaoe/GA634D8I6XvPazdO6LQxsQPFucpZBJwgW6r36U9svkMLsj2XLUbyLs9HmmIUgCFQex1dsDFyC4l3MWQ2dIGhUii8O8//52QwPtw5lFcm6nOpv1ntj5WjvuTtZ6mHwgwzSWJwnN0AnCB4V4p3nlF7X+FFPygdAJ6CAVnx1VvdQKaRbtuxaQhRPpn8CnYACWvFjUovsDntJxrjABhT/i6w308zQig/3dlDUdQlsQPGPyo50zAyt+F7rPA+bIwIbhBd/Ih46ARW04nd06pmemRHReWdgg/Din/gIOgEV1Hv1Dlt+br5NPq2O6OJPxUEnoAOvwGHErHXQCehA8WyoiLVOp50HFM+G596BTkAJimdCZZRlTss1gOKZsDAPOgEtKJ4FtZHV0BFoQfEsWDEfOgE1KJ4BzpjfoCNQg+IZsMEi19I3BcUzYGA5dAJ6ULx2tmVCJ1ABitfO8APQCVSA4jVTNhg6gRpQvGb+ZzN0AjWgeK3Y+0MnUAWK10rOSugEqkDxGqntWwsdQRUoXiPvzINOoA4Ur5EBZ6ATqAPFa2PrZOgEKkHx2sjYD51AJSheE8fuhU6gFhSviVnWupi+CSheCxej6qEjqAXFa2H5K9AJVIPitdDfelfe+EDxGiiaAp1APSheAyP3QSdQD4pXT/ld0Ak0gOLV8+z70Ak0gOJV4+grv1ncOqB41bz/NHQCLaB41dz1C3QCLaB4tRwcDp1AEyheLdM3QSfQBIpXSVW0xYbACADFq6TAut30HlC8SuLt0Am0geLVsXssdAKNoHh1/LUEOoFGULwqzveDTqAVFK+K1xdDJ9AKildFvEWmFWwZFK+Gb8ZDJ9AMilfDFKvv2qF4VVRaftcOxauiYBF0Au2geBUkWrzXzg2Kp2ffKOgEDEDx9DxmyUFvAkDx1FRHWvuErBcUT81qi46B0RwUT809lr7WzgeKp+VoKnQCJqB4Wuauh07ABBRPSb2lb6NoBMVTYvsbdAI2oHhKRlhxqGoFUDwd9kToBIxA8XQsXgqdgBEono64c9AJGIHiqfj2QegErEDxVEzn4fyMB1Xi7Qqfd0KIr+Xi/IwHWvHfDyg7fMelIXE/BTYIIX7tM9AJmEErPmqOI2V2dXX23YENQoi//0foBMygFd+hhvzlPCH14YENIog/aeVhrgKgFR//Nhm5npDPIwIbRBD/8lvQCdhBK/5Ij9vuDblrYKdtgQ0iiO9XCZ2AHfR79VuX5y75RD4BjwDid1lxKtGWUHkcX7ExcIkA4qd9BZ2AISrFF4f5/9w0w8OtaawimZXaKG4O4gmLnrufbB7uHsMijplZx89BPFEl3nlO6ZXP/0f90CPQCVhCK746q3uoFNItuyawgXvxp5OhEzCFVvyY1CK7w16SMS6wgXvxrxZAJ2AKrfjwi56Hui6BDdyLH8DLmXgvtOJ7rfM8bBau527faOgEbKEVv6NTz/TMjIjOOwMbeBc/0wadgC3Ue/UOW35uvk1+bTnn4usjLTvDnDJ4BU5w2GZCJ2AMig+O0RaecEoRFB8U5wdAJ2ANig+KFRwMd9QcFB8Ud/0KnYA1KD4YjvFxT3xTUHww5LwHnYA5KD4YYi5CJ2AOig+CrydCJ2APig+CR3m65qoBFN82Dq6uuWoAxbfNR3OgE+gAim+bEd9DJ9ABFN8mv8VDJ9ADFN8my1+DTqAHKL5Nkk9BJ9ADFN8WP/HXXesGxbcFh921blB8W8Ty113rBsW3wa6HoBPog4L4q10/b1MX4lX844XQCfRBQXxIww8dnIqv47G71g2Kb51PZ0En0AkU3zqj90Mn0AkF8Zfu3r3b/bObqhCf4i9wd3WtDwXx4T6oCvEpfuXL0An0Ag/nWiXlBHQCvUDxrVGeAp1AN5TE70rv3q57xjd0hbgU/zJ9f4ZVUBD/5e+f2nJwy+wOdCN0cyk+7gJ0At1QEN/X+zJ/N4qqEI/i93M8kpeC+Mu9L/OqK6kK8Sh+1mfQCfSjhQ4cN2GyltbgULwzsg46gn4odeDs9XIZVSEOxRc+Dp1ARxTEd/BBVYhD8Q/tgk6gI3gc3yLVsdAJ9KQl8fJxyduAP/Fr5kEn0BMl8UtTCFl1bQ7dME/8iR8im3CJJxTEv9z1A0IufnzzfKpC3Ik/w9fYtYEoiL/eO3jhd92pCnEn/h/LoBPoioL4K70jU9f+B1Uh7sQnnoVOoCsK4iO9Ew1tv52qEG/iy9KhE+iLgvg3btrr+n2o52KqQryJf/pD6AT6orRXv+D3PZJ7X0l5lSFv4qOpD2itheJxvP2TvPXHKQtxJn7bw9AJdAZ77pSZWgKdQGdQvCK10dAJ9AbFK/Lh09AJ9AbFKzK8DDqB3qB4Jc4mQifQHRSvxOtLoBPoDopXYuAZ6AS6g+IVODIUOoH+oHgFnvsAOoH+oHgFomUz5/IHipez/a/QCQwAxcuZuhU6gQGgeBm10ZwOe9MMFC9j/TPQCYwAxctIOwydwAhQfCD2gdAJDAHFB7LkdegEhoDiA+H86lof6sQrHO/wIv6HDOgExqBOfDv5Il7EZ30EncAYaMW3D3MjhclGTeBEvDPKAR3BGGjF7+0z8vCJE1eckI3/xon4zY9BJzAI6o/6unk9Cjn+qB9PN5CrdVHxHb8/aiq34iv7QScwCjU7d/UvjZIv5EP8yv+DTmAUKo/jKzYGLuFD/D3cjl0biErxxY179UU5HvqksYoEyLHB0AkMQ3vP3RGbh7sVPv4txwt8TjWmhArxznNK56u5+KiProZOYBi04quzuodKId2yZVel8SC+hPdbZJtAK35MapHdYS/JGBfYwIP4yTugExgHrfhw77yLdV0CGzgQXxUDncBAaMX3Wud52BwR2MCB+HfpBnizNrTid3TqmZ6ZEdF5Z2ADB+Lv/QU6gYFQ79U7bPm5+Tb5KSzriz92H3QCI8ErcPw8vwY6gZGgeB/OKAFunGoExfv46hHoBIaC4n2Mp5xvzeKg+AbOC3Mm3guKb2D5q9AJjAXFN5DE/+gnzUDxXg4Jcjm9HxTvZcan0AkMBsV7cETSzcBjfVC8hw+zoRMYDYr3kPojdAKjQfFuyu+BTmA4KN5Nzj+hExgOinfhFGFguwBQPBHoTskmoHgXo/ZBJzAeFE+IPQk6AQAonpCFb0AnAADFExJbCZ0AABRPtgp0/0wjKJ6MFWUQjGag+Ip46AQgoPhFfM8T3xIoPvYCdAIQhBdfNBU6AQzCix/9LXQCGEQXf1rEXjs3oot/8S3oBEAILt4ZdRE6AhCCi//sSegEUAgufsi/oRNAIbb4n1KgE4AhtvjZG6ATgCG0+Jq+ddARwBBa/Ds50AngEFp80inoBHCILP4b2fCcAiGy+PGl0AkAEVj86UToBJAILH7eu9AJIBFXvCOyFjoCJOKKf+9Z6ASgiCs+6SR0AlCEFf/1eOgEsAgrfrSQV9M3Iqr443dBJwBGVPGzPoROAIyg4iujRBveLBBBxf/j79AJoBFTfH20iHdGN0NM8eufgk4AjpjiBwozaXSLCCl+WyZ0AniEFJ92ADoBPCKK//5+6AQmQETxk4qhE5gAAcWXi3qDbDMEFD9T3LsomiCe+LOxTugIZkA88c+/DZ3AFAgnvjJSPhO2iAgn/pV86ATmgF58hef32cDFFhFf3bcaOoI5oBW//+bfdf3A9fzJ/plFxC9ZCJ3AJNCKj5tfu6ljkWXF10ZVQUcwCbTiO9QTsvamWquKX/YidAKzQCu+21bXr6EPW1R8baSY45cqQCv+vfYJZ4i9z23WFL98PnQC00C9V//LWtebpva9GYHLrSAe3/CNqDyOr9gYuMQK4l/Hb3g/KsUXh/n/LM3zEJPOKpJu1ETiLr0f7T13B1Z7SBjJIo6uLF4EncBEqBDvPKd0esv8H/VVkdhp1wit+Oqs7qFSSLds2WSs5he/YCl0AjNBK35MapHdYS/JkA0YZXrxZ6PwtFwTaMWHe4f5rusS2GB68bPFmyq8NWjF91rnedgcEdhgdvHlcXjhTVNoxe/o1DM9MyOi887ABrOLn2KDTmAuqPfqHbb83Hyb/OvS5OIPiTtAuTKiXIEzTPCRT2QIIv6rsdAJzIYY4p0DjkFHMBtiiF+ZBZ3AdAghvhJPx8oQQnz2m9AJzIcI4n+Mx74bGSKIT98OncCECCD+C9Mmg4R/8bUxYg9T3QL8i5+/GDqBKeFe/NE4cScVbA3uxQ+RnUdE3PAufu0j0AlMCufif4s8Bx3BpHAufsp66ARmhW/xWzKgE5gWrsVXReEhfEtwLf7xd6ATmBeexRcNg05gYjgWfwE/6FuBY/GTPoBOYGb4Fb9BdpcX0gRuxZdHYddNa/Aq3plSAh3B3PAqfsEL0AlMDqfiS+4VfSbJtuBT/OlInF+sDbgUXz+oEDqC6eFS/Fwc1qxNeBT/YQZeR98mHIo/0A9vmGob/sSfiT4MHcEKcCe+9i7csQsG7sRPKIBOYA14E//cHOgEFoEz8QUP4g59cPAlfsOgWugIVoEr8YWJeCAXLDyJ397PDh3BOnAk/mu8H5oCfsR/HV0OnMBScCN+O3qnghfxhf3wc54KTsR/mIj7dXTwIf71+yoBt25JeBDvnD0BZx2hhQPxVSPwilp6rC/+2ID3gbZsaSwv/suoPTAbtjgWF+984b4zENu1PtYWfyIlF0/DqsPS4tdHbTN+o5xgYfEVYx86b/Q2+cG64t+Pks1hjwSPVcUfSZ2G979rwZriK7MSdhm4OR6xovi6ZZErcWdeI9YTX786OqfKoG1xjNXE170bM6fCkC1xjrXEn3816lnUzgQriT84PSYPP+QZYRnxF1Ykp2/SdxNCYQ3xNf8aFffKKR03IB4qxDvPKR1L6Sf+3Huj+s09oFd1UaEVX53VPVQK6ZZdE9igj/j60pzkOxd8r0dpwaEVPya1yO6wl2TIBoplL/5i0YLBsVPXYc+sLtCKD7/oeajrEtjAVHzlzvzJcQlPrMOvdd2gFd9rnedhc0RgAyPxJ4uXz0yNufOxN7/DeQJ1hVb8jk490zMzIjrLpvHTJt6+b9Ob8x4e1L9/2sxlxXgxlQFQ79U7bPm5+Tb5dex04h32owd22tYse2nO1BHJsf0H3Ddp7usf78Wvc+NQeRxf0XgRxO48D7Ejmq9x3LN0UY6LrBkzZkyZOHFUWlpyUkxiv4TYpNikIQ9MnPHcwhXrC789flFLfkQlKsUXh/n/bBA/fG7zNY4vXe3hM5vNVlRaWvpDWdmvdnxLmwZmPXcFBYwKIYbArOcOxVsLZj13KN5aMOu5Q/HWglnPHYq3Fsx67lC8tWDWc4firQWznjsUby3wOF5QULygoHhBQfGCguIFBcULCjvxY/OakzqKLSNSGBccNoRxwcHDGRe85wHGBdMa9YxlJf54gPe8ToxDD7qBccF+fRgXvHUg44Jd0hgX7NyoJ5+VeBkDGNfbN4VxwX8uZlwwezPjgkNPMy7Y1AmKZwaKl22EBSheOyjeDYrXBxSvGRTvBsVrxxDx7zCuV/Ex44IHv2ZcsPBnxgXXyK5s1EhTJ7qJR8wNihcUFC8oKF5QULygoHhBQfGCguIFBcULih7i7zjofdwZ8ccHq1kWXNO9QzyToQ59BQnZ255lvR8TO0QzGZXPV3Bt998PYjBPduMT53fCXnzhQ5I3taPT2qrBWQwLlrcvrs+9hWFBQur6hLW6KmW9HitrZjM4SeEvePry9RXDxmuu1/jENTphLz538uXe1LaehBTdwLDg2gRCan+nfdxyf0FC5o1kIN5fb+t/u57aMoYFt3d0vVv7aK7X+MQ1OtHjo76jN3V+OiH2UBaTiDQUvHCGkE1dGdTzFST7byljIN5fb3nKyP8cdJhhwQt/Wnls0JOaqzU+cY1OdBSfm+l6/UssBjzq6P9KXncNkymEGwrWRxWeYCk+N2TduafuYFiQvCaF/ZHJ+J4NT1yjEz3f8RmuV1cIw3c8sQ+7sYhBOX/Bl6YSpuKXJBNSfamdXcHCmw/W5PXWXs7/xDU60VG8LYKQkm4MC9beMaWWRTl/wVHt2l0htdvKrN4ndxJSE8LwMy77UUKcoZrf8o1PXKMT3cQXVzg6b64bkc2w4Hu3VbtgWND1i+E7vriittMXzqfjGBb8qNvuykXXa67me+KaOtFNfNhGsrP3dWyO4xsKzpTcnGVXkDAV76pX0is8+SjLggv/0mHAd5qr+Z64pk6w505QULygoHhBQfGCguIFBcULCooXFBQvKCheUFC8oKB4QUHxgoLiBQXFCwqKFxQULygoXlBQvKCgeEFB8YKC4gMIIXt7BCxopPFmS8uD4gNoRXyTmy2tD4pvziCpx66bZ3b582eEbLr1muE/uxfULLy2Xd9DzW62tD4oPgDXO15a4syKIafCtzhmJLgXHA/79uLESe7GjiieW1zir3aSPT3IshGE1LRzuhZU/0KqHstwN6J4fnGJv4W4v+fnhvdwcdq1oP6Z226PQ/Gc4925c/38Yxwh9QfdC1b3PEZWoHjOCaluEH80/BvHC7HuBQuT638akOJuRPH8cv+1u7ziyfob/5Bw2L2gPO6avu9HfEBQPGJ9ULygoHhBQfGCguIFBcULCooXFBQvKCheUFC8oKB4QUHxgoLiBQXFCwqKFxQULygoXlD+HxH6wceZI6wFAAAAAElFTkSuQmCC" alt="plot of chunk unnamed-chunk-1"/></p>
<pre><code class="r"># The posterior
thequadrature$normalized_posterior$nodesandweights
</code></pre>
<pre><code>## theta1 weights logpost logpost_normalized
## 1 1.246489 0.2674745 -23.67784 -0.3566038
## 2 1.493925 0.2387265 -22.29426 1.0269677
## 3 1.741361 0.2674745 -23.92603 -0.6047982
</code></pre>
<pre><code class="r"># The log normalization constant:
thequadrature$normalized_posterior$lognormconst
</code></pre>
<pre><code>## [1] -23.32123
</code></pre>
<pre><code class="r"># Compare to the truth:
lgamma(1 + sum(y)) - (1 + sum(y)) * log(length(y) + 1) - sum(lgamma(y+1))
</code></pre>
<pre><code>## [1] -23.31954
</code></pre>
<pre><code class="r"># Quite accurate with only n = 10 and k = 3; this example is very simple.
# The mode found by the optimization:
thequadrature$optresults$mode
</code></pre>
<pre><code>## [1] 1.493925
</code></pre>
<pre><code class="r"># The true mode:
log((sum(y) + 1)/(length(y) + 1))
</code></pre>
<pre><code>## [1] 1.493925
</code></pre>
<pre><code class="r"># Compute the pdf for theta
transformation <- list(totheta = log,fromtheta = exp)
pdfwithlambda <- compute_pdf_and_cdf(
thequadrature,
transformation = transformation
)[[1]]
head(pdfwithlambda,n = 2)
</code></pre>
<pre><code>## theta pdf cdf transparam pdf_transparam
## 1 0.9990534 0.008604132 0.000000e+00 2.715710 0.003168281
## 2 1.0000441 0.008809832 8.728201e-06 2.718402 0.003240813
</code></pre>
<pre><code class="r">lambdapostsamps <- sample_marginal(thequadrature,1e04,transformation = transformation)[[1]]
# Plot along with the true posterior
# pdf(file = file.path(plotpath,'lambda-post-plot.pdf'))
with(pdfwithlambda,{
hist(lambdapostsamps,breaks = 50,freq = FALSE,main = "",xlab = expression(lambda))
lines(transparam,pdf_transparam)
lines(transparam,dgamma(transparam,1+sum(y),1+length(y)),lty='dashed')
})
</code></pre>
<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////isF19AAAACXBIWXMAAAsSAAALEgHS3X78AAAgAElEQVR4nO3deVxUVcPAcd+net+eMkXFjVxatFxQYNiGfYhNEEUFVFQwFzDUFlwTF3AXFbdIFs0VyzQ1FCUXSFHABcHhtFhG+bRYmUuPUSII552NkmFELnPuWe49vz8YPhe5nXO/Mfvc0wLyZFkL0gPgkYnDyzQOL9M4vEzj8DKNw8s0Di/TOLxM4/AyjcPLNA4v0zi8TOPwMo3DyzSpw1dsSte3o5b0UOhK6vAXXBL0vfw76aHQleThXwX6vDl8vTi8TOPwMo3DyzQOL9M4vEzj8DKNw8s0icHXlhv63rCBwz8kicGf7R2mr3OlfgOHf0gSgy+INji7/6XfwOEfEoeXaRKHr1jYQ+EUlnyRwxsnafia9S4zQsvUWXH9E8s4fP2kDH99wKpq/VX9xSle7hy+XhKG/49L4T+38dtbl5EeG11JF/4Xl8/hA3funB3/Q3pwVCVZ+FveF7QbdPDqzTuLvM+63yY9OpqSLPy43boNF17dlQzOx4137/t77uAawqOjKanC94rWbzjZIzjH8Dj+bPxisoOjKonC5z9zQ78hwe/vJ3C+dnc6S3Z0NCVR+EH9TDxzVzGg912yw6MoacJvD+nsrFQq7doqez34lG112jzSA6QmScKr7fPaFANw1u4gWFbvufpa3xzSI6QlScIvnwK08KNTgRE8PNvmB9JDpCQpwpf2P6eDLwUN4OEIxX3CQ6QkKcInzgAa+KJiYAL+es9zhIdISRKEv9T/ghZeddAUPEzcRXaEtCRB+KQ3Nc5tVo4FJuHvuBzjT+DB5sHfuoN8GKjSwtsXauEVRabh4fKwdKJDpCSB8F+ryr91euxxr+8f/U+JpIHfHKl1blMIHgJ/x8WFv1ojGN5lfnXwvMrKhEBxRmN2Gngf3XPzbYqN4JUfH9f3U8KcTNLDpCCB8Bb34HN/QFhjKc5ozK4g+ugrWubXnjGGt4wcryvwpZD24doi5f30rUB4711wVBaEJ+zEGY3ZFURPSNUon1Y0+Iu31N/mg9UhOYNXvp+Tk+P+G+nBEk0g/HfW9gMfH+DXmdaXuQom9FdrdCetfTj8GHDYx/8QAP4cXlCFW5LScqrEGAqKCny1j+XAVnVj8MBnQzCHb9Zv3f4E8TBQVdAtV+/bKPzmKO+jHL45FTz597eZvroU76AakVkd6K6xvZT0CPgy29wSDo+mrVsR7ci8Yn01tslTHwEP4hP5Vb3g36i9Y+pUgZTA9x+vYXU59Sj4IoejMzm8kCoX9HyixeM9Eu4Z/4AO+Ish0ZrHciMfdRsPQOhuWx8OL6CoIWduVd8qGjnO+Ad0wE9PqXuX7SPg3x++oDeHF5Cl/umu+92Mf0AFfI3ydBPhgeI0hxeSzce6i5MNnrmjAv50XEE0WLq0KfCzl/LbeCGd79xvxMSRdlYXjH9ABfzrZzXwzmeaAp/v6i/vt+IIvVdffTwjKeN4dYPtNMDXuNQWRGcHNumqHvh5BMj6fZdSehx/ehosiD78cdPgN/TYkUh6wCSTEvxbBf98du6R8KUW11xID5hkEoKvVdbAgpHHmwgPuuVcIz1ikkkI/uJrEBb02ttUeOcJpAdMNAnBL/hE8zizY1Ov6oG/Q3XGDdJjJpeE4D3uQZjs0HT46Ue2bSA9ZnJJB/7qcM2XTyc0Hf5M5B9epAdNLunAp+yEUMC9euD/m/vduf8lPWpiSQd+oOYG++Q7QuCX7ic9ZoJJBP5GOXAtLy8fLAj+ykhI7ZsHRU8i8H5hbtZhYcMsFULgoddfC/KJDptgEoH3UQ/XPIDPmr9AEPyKfRdjiA6bYJKBty3TagqDvzIKujR4L5FMkgr8gaFAODz0vJtXQXTc5JIK/IzVAGyfKhR+yUGioyaZVOA9zgAQtUMo/BevwnMyva6XCLy3k4bSRi0UHrpVLztAdODEkgi8XQwA6iTBt/Ew/vjXo4gOnFgSge++DTTnzh08PwV6NnwfmRySCLzFRQCyhcHbvjFbk9Xst04QHTmppAF/3VINzngJg++4ZpMm77mLYkmOnFjSgN/9khokzRYIn6/9mjEpy02W9+ulAR/trAaDDzcHvkSR5XCU5NBJJQ14Vx812Cfwzp0eHgRmhL9GcuikkgT81QgftV5TOHzSxNiJBIdOLEnAb8nwUacXNQ++QMHv3JkTUfjIb3zUzsXNgwd28nybtSTgXaHPaY9mXtWD8f5wMMGxk0oK8FeioM+Gmc2FT30Zjv2C3OBJJQX4TVugT6m6ufBZne7vX0Fu8KSSAvyYctj8e/Ugy6agQoaLF0gB3hVCtyXNhx8kyzXJJABfPgbCl5c1H36SF6z6kdjoSSUB+K0ZELbPbT58bOiP10KJjZ5UEoAf9xWEVmbcxsdu3gQ9K4kNn1ASgHeH0Kw7d7E/hcL5x4kNn1Dsw/8UDs2Eh55VP8vuRp59+A/XQzjaPPg5n5IaPLnYh59aAu96mwefPwte+JnU+AnFPrznfZg31zz4ag/44RpS4ycU2/BFXZSOFkpl735PmAUPw368PYDE+AnGNvzJWJA6VavX2jz4ze9BP5l9VJ55+JjNoOy8ufA/DId/kBg/wZoDf9PEJ0yJwSvPg/S3zIHf45Went4zNT19s6z+5gXCB9+FX9v/64mABveBScGXOgAwYZs58CmuycnJXq8nB9rJ6kyXAuFbVMBXZt+9l9jguW1S8B9oHJ0umgWvXckk7TUw1JPDN/LPK6BFNYQ1HYx/QAp+VjIAa8y6c6eDL3YEyb05fCP/vBwqvoHwcyvjH5CCDzih0zMXHrifKXyBwz88zy4tLQJgnuVq4x+QgrcFIFuNAH5mMgji8I1V+eVZWJDbYDMh+MhAAJwvIoA/EMrhm1DDtWUJwfu9DYqVKK7qy2yA6g0SUyCV2WvLZsfoUi1BNSIhnez/AdgSgwIehBwMdCIxBVKZ/czdjWJdiRkoRiO0k1al4FQuEvikOUGDvyMxB0KxvbbssS4GRfPhT3kH7S4jMQdCsb22bIotuJSBBh4oAvmdu4dH2dqykwPArihE8OOV1343dV0m0dheW1YVBaZvQAS/qce1WUUkJkEmtteWtY4FfvmI4IvbXTuykMQkyMT02rI/qWKBng8BPGh/ucKHwCQIxfTasvsmxxrQUMD32gJzCEyCUEy/A2dWSuzmbGTwbg3usUo5puF9j8YG5yCDD3SAsJjALMjEMvx9z5OxCnRX9UGDv69xwz8LQrEMr558cvxAhPCrtsDw7/FPg0wsw2dsO4nyzl1Q0SjNLvFPg0wsw0/8Ai38NZfaP3/BPw0ysQzvUXPSDil8NMA/CVIxDF8RAHO6IoX/cB3M/w/2eZCJYfj8OXCNPVL464PgHrl8eJJh+NUH4LwhSOGhR9VvwdjnQSaG4Yf/CBHfuYPTC6Af9nmQiWF4N4gc/sgi2awvzS78b0Pg5VcRw/8hl793luGPLIWb3kYMD33+hPNxT4RMbMIvCw8P7+sZ/pwTaviEo9BPHp+UZxNeVVBQ4JZTMOlN1PD5s+GiIzhnQixG4TVO2tflkF7Vv+ikVDq3UvbrnohzKqRiFv54ILikRgpvlaf54lEADsticSJm4dfNBJMy0cPHrefwgsIOH70NuJ5DD797FJjhj3MqpGIW3v2cWoH2Nl4Hr9npO71xToVUrMKXKUDxGhHggc+n2R1xToVUrMIfCdEqiQAfn3RYKYdPUrEKv3oOyBUF/uNQfudOSLjhx2VqH8iLAF9me3jibZxzIRSr8C7FBV6iwIOgLRO9cM6FUIzCa+7bac9eLAZ84psTI7/BORkyMQqfPQR8clQc+BzviVvTcU6GTIzCJ83VIYkBD3pPvHkW52TIxCj82F1le8SC9wvDORVSMQrvUnxoqFjwM1zh1Zs4Z0MkNuG97MHSBLHgdzwPP9iAczZEYhPeaSgYuU8s+MMda34ehnM2RGITvs9c3amLxYHf9u/+ymeUmjpI+fxnbMJ33aU3EgV+k/MsoDvPRshVnHPCHJvwbYqP7hMPPthb/w2Hb0JY4WtbgxnrxIMPVajBag7ftLDCX+kEAnJFhI98HygucfgmhRV+d0/gCESEf+ctELmLwzcprPCz7EChmPBF7iBlJodvUljhAzwMViLBA8eSshIO36Rwwtd6qNbkiwo/frv2Gw7fhHDCfxulcjsnKnzaZLBhP4dvSjjhP1rrZTivoVjw51207/Pg8MYVNtyEEz7+lHK4uPDAufisK4dvUMuGm3DCB95RAZHhJ20Cgzj8A7V6UluLJ580/gFOeA8oOvx72jXNOPw/feY46ttffnnqlwYngMQI/0MEfE5s+GInDl+/+8utTxG+qs9adbeN2PDA7SwYyuHr9YXLVLLwCScK605oKR786+8CZTCHr1dN8piGGzHCD765r7/o8JljwWuuHN64BxYV/nKPrsnrUY2okYp6aN8Y84xS+aTo8KX2YI8zhzfugUWFL6brGttgRXkROva6huSUHwDPiA4PvM7wO3dNCctVvQ4+9Y0z4zHAT1sHBnH4ByO4qLAO/vXUtKkY4D8cBfpuxzAnUjG1qLAO3u/k1DQM8GoF8JiEYU6kYmpRYR28AswsxAAPfD4d5IxhTqRialFhLXyh9ulaHPDxSSEpGOZEKqYWFdbCb47FBH8gjN+r/yeyiwpr4WeuX7sQC3yZbcjVKxgmRSimFhXWwgcffTUTCzwI8vvWE8OkCMXc43gFUBbjgU+0uxp8A8OsyMQa/AV3MB/PbTw40uXqqv0YZkUm1uAzx2tNsMCDtlevS3cVOtbg5yadVeOC75aLYU6kYg0+9FD4IVzw9othwe8YpkUk1uAd1PbY/uIDAuG63RimRSTG4Eudi5UAF3yIQ41ask/XMwa/d9T5VHzwURdrVRimRSTG4Bcv1pFggk9bBW9hmBaRGIOP+CgbI/ylgRgmRSjG4J1KFBjhr3pW3UvFMC8SsQU/1eGMN074maehC4Z5kcgUfNyZGsH7wQMfOSx9Ck74TxLh6HIMEyOQKfh51p1eM/ECXKPhgQ+Ye+IYTvgKX7jjIwwTI5Dpq/orq90tx2bdFbAfPPCKTD0JLnjoW4FhWkQyDV+RNan9y+7t9jR9P3jgu53fjhd+4SePHhWbmYJf6/eM77pvICx5tun7wQNvdSIQL/yZGfAdab4NxxT86H13NF/vweqPm74fLPCZL6+diRe+yhNmvothZvgzBW+t/VLZRdB+sMAnukdjvqqHg377KRTDzPDXEP7xx1s8rk3YKduxwEcOParGDL9mD5yJYWb4M/UXP6gZ+8EC71rHig9esi/PMfXMnfWErbjha10xzItEDeFbnrbWJ2g/OOCvu/kl4IaHEd9VDhF/avhrCP/Jrc/0CdoPDvijr/b7CDv85k3Qo0r8uWHP5FX9V7Bq3xZhs8UBv2Jx10vY4a+OgNMLxJ8b9kzBz3/q/oo+zjGC9oMDfsQHkwB2eOhaW9rgA2MSyBR8u29hl7O3OwjaDw54D93HpHHDx5aKPzMCmYJve6u4c80fFoL2gwH+TtDSOjSc8B+tgr/cF31y2DMFH2Pz4srrymBB+8EAnx/vF0UA/lYQXHRU9MlhzxT8/b0fVF9b8V9B+8EAv26vNYmreqiqPDVb9Mlhj6EncKK+dicCH593z0v0yWHPFPwnLtrnb14XtB8M8B61RO7cwU/nwHzRJ4c9U/AvbP/y8uXLXwnaj/jwd/0gGfh7XmJPjUSm4EOasR/x4c9Pm5NDBB4OvFEjvT95U/Arjwnfj/jwqZluZP7iYfLeGjexZ4c9U/Be/9e5L30v0kzKG04Iviwahv8g9vRwZwqezhdpVEXvEYKvdYFpO8WeHu5MP5yrFXzSH9Hhq70hoTt3mgeSX/0luQ9PmoL/4ZWW1pfdvxW0H1HhfysvLz8cWV6+nRD8Tgme4tIUfPjb96xrlvkK2o+o8J5hYWEO9t3D+uGFd46drW/7EPiJ1M5yaQq+YzW0htVtBe1HVHjt6WtHrwkDcXjhu6Z8qOvdKI+q3etEnB+JTMH3O6OBL+kjaD+iwzsnLsQOf1R/eSxqZv5vzXlyg+ZMwedaRliOszwoaD9iw6sdlueQgz8+D0rtHAkm79Vffy8x9ftGfsnEJ2nFhj+o1SAGf9dLxNmRSeCrc79O8oi/ZvuYa4O7/GLDJ80rJQgPgyR3UlsT8CUjerbsOdL0G44GBu+Jssq4vSzI+Adiw49dGEsSfu1uGCfiBAnUEP7T1nPzL+fPszhp6p9b3ISfPV0Faxq8IU9seOW0NSThPx8HA6R1ksuG8M67dBcfmDz5S9cr8P5+zZ0AK+MfiAxfphhwgiQ8VNauEPDZYQZqCP9v/Ukg/nra1D9Pe8pP83W79TzjH4gMfyQkERCFjwafvyfiDPFn4tOyhssGS8TrurxD8yXtQIPtIsMnzwFk4T/CsZImzhrCP2Z4ce5/G/mtB9aWvXdLV4qYfw8qMCGtmCz87UAR50eihvAWdTXyWw+sLbs3XJfTQhEGV5cKuI/fThYe+vy5R1LX9Uy8y1YFFO7nCMMvPvxDuIhTxB4Ta8uqjgbpzmRKCn6k5rYsN+aWs5SexWFibVnV2hlpBOHft/L39/ez8O9v9Z2Ik8QcE2vLqmK2AYLwO0O0X4OPgFAJnfmMibVlVR47ycOviAeB8oUns7asSmGvJg6f7w06lYg4ScwxsbasS4ArIA4PnIqtJXSuQybWlrWOiKYAfnKaT4SIk8QcE4/ju7+TSwH8+2NC94o4ScwxAd/OgEYWXm0r43v1D01U+NapNMCDQQEXz4o4S7yxAP8fyyAq4JNszjfnbK90xgL8/o5vUwF/ptMVVaV408QbC/DxHd6nAh50KFtg8g1pLMYCfKBzGR3w1ml/8b94o8SEd1cBOuB9R4o3SdwxAP9diC0l8KF2te//JNo88cYA/N4hfWmBjyjdlibaPPHGAPxsbxda4DMW/yyVk9czAB/gQ8ttfCjwhqNFmyfe6Iev9aiiBv7K4F/Emibu6If/ehykBz5ti1jTxB398LtG/koP/I/DpHJdTz98nP09euCh591R0lhXnH54bxVFV/Uw8UjmO2LNFGvUw993WkwT/MXYm9L4QA318OA1SBN8rdLUx0kYjHr4Le9RBQ9jL1b/IdJUsUY9fOwwuuBzEr6QxHKz1MM7RdIFX+kBXaRwbU87fGWfNHrgbazt7e3b9Gvfe7k4k8UZ7fBnB3xDD/zzhzVfVrx9cLcE3l9PO/yG3ZAy+AJXkMfh/04M+Eu+vr6dXTVfLGmCB+6nN4WhnyzuaIY/HAdAj5kPONMBH79kuhP6yeKOcvjCbluogz/hs7cL+snijnL4jC4XqIMHTjkdTJzGmbEoh399zIPOlMDHLRqKfrK4oxze+zSF8Nl+/F7934kE309NITzoN2JdDfrp4o1u+KOdztAIH6V6sxD9dPFGN/zKLvWcaYHf1COP+RXl6YYfrKISPs+yMgr9dPFGN7zdbjrh++Shny3mqIZ/w6m+MzXw/lPQzxZzVMOPCqIUPsKlZgT6+WKNanhXG1rhp+cP/RH9hHFGNbzVFFrhz76+nfGTHVIN//Q+WuFrXW7uQD9hnNEMv6t7Ga3wcPop9PPFGs3wC9yNnCmCPzel+lf0M8aYYPjbuq8NFt8TA94/kF74WuXPg9HPGGMC4b/o868X9kNY2eDXxIBvP5xeeDg71++/6KeML4HwXiur8jqdwQNf/fSbFMOXRK/fhXzKGBMIb1ED4YHeVVjgz3SMoxgeut74BvmUMSYQvof25cjQKVjgN8RRDZ+QjXzGOBMIv7fVKzfhLUd7HPARW6mG/3J0eTHyOeNL6L36awcqIKza2+DlaBHgux+mGh56Xg5FPmd8Ne9x/ANryxpCD/9lZ8rhV33wyh3Uk8ZX8+AfWFs201dXnwWoRlTX7IHUwh9so52yR4eXX/wB9ayxRe8zdx47qIXfq3+52C1/xOeoZ40teteW9ayiHX7hgpGygceytuxnxzUdsj++lHL4QqU3u2fAonFtWdfxmhyfG6+gHB4E+HmhmzXmaFxbVvfW2v7xYDLt8Bv6uf+Gbtp4o3FtWR18qwL64UstU5hdbZbGtWW18JesAP3w4CV2X6ehcW1ZLfye0SzA+w//s8HdXEai8XG8Fn5GMgvwI+0+XIlu3lijFb5tHhPw03d4ops31iiFP9wWMAF/cuCrAN3EcUYp/GR7NuA/H1x0Dd3EcUYpfGQ8E/Du4X6OMdq2oZs9piiFD/6ECfie6Yd6p6bn5Oxh75V5OuEL+gE24D8GUYs1Pz3N4VGkAuMcmYH/ONj+PIdHkwr0mMMMPHCe9g6HR5PqfPtsduCXzOFX9YhSlfYC7MBfUHB4RKmOBjEEDyKWzObwSFKNTmAJft+g/qc4PIqUXQ6zBA+cI9ZzeBS91A0wBb80ZjWHR1G7ALbgi23yOTyKeixhCx6MW+WPbvaYohG+w3HG4HM8LNHNHlMUwr/bCjAGDzy6FSGbPqYohLe1Yg4+uXs0suljij74cru+zMHnt7uKavq4og8+28GNOfjTthmopo8r+uArveoWlmQIPsRtPmMrDdMHnzeHQfjQ6UMZO4U9dfCHgvNYhL/qHY7oAGCKOvhxjpUswsMRkxAdAExRBH97j6bdL9rvsWER/twE8w8iziiCP+AVFxc3xVYV9ySL8NAf3Df/MOKLJvhZ2mMYcuhvVrbgD/kx9eZ62uBP2AJG4WuVbuYfRnxRBp/1yghW4eE2t9NmHwN8UQb/ZvAGZuGrHFj6GB1l8E5O55mFh+vTzD4G+KIMfos7YBf+TyeGzpJAGfyqWQzDw6Re7JwEiy74NYOyWYa/8zLyM/qKFlXwn6psH2RlB/5E13Bdfbv8bP6hxBNV8HNjotiEz/Yu0HWsHTPnvaMK3m/MFkbh/Q0bbKc3PBMcnVEFX2p3iXH4EVabzD+WWKIK/nBQPVYG4cdMfZGRl2pogrefsZp5+HNdPzT/YOKIIviMl5yLmIe/lJJk/sHEEUXwkQO86rMyCV+lZGPlUZrg/RIlAA9392PiVl4g/OW6jH+A4ja+xykpwNc+n2L2ocCQQPiAFv/upMv4B+bD3971nBErm/DwSEcWHssLvaqf9Jrp7WbD17i95SsNeOheaOaxwJFQ+JMPudNqNnzeTOVkicB/FWjmscARNXfuxhbazJIIPJzcYOld+qJmbdminRMlAO/oqF181OsxjzlmHg7RM3tt2SzdabtjVIvMHcngLRKA7/WR7mKkzWBzD4fYmX1Vf6dc10ozX5zYddP/gHTg1R38zD2uYkfJ2rJfjcxIlxA82NrZrMOBIUrWlp11IuCGlOCLutP+8iwFa8tmp6en2y23S4+VEvwrre40+4BgiYK1ZV2TNQVHJftICX5wBOW38hSsLav7OLyiGEjh4dw/8JUdzzX7iOCIgrVlVeDg+A/DgMTg4Y/uVK8+SsHasiowetuoHZKDh+kzmn1IMETBU7aqc4qLNmXSg7/eluYrexrgT25d/RaQFvzJdn369Gn/RO8+HT5DdIBRRwM8AJ55EoM/4a35on5xBBh3HtEBRh0F8I4njupeiZccPCgJ2sThH17b3JiN0oQHhX0jOPzDKulU2l8tUXhgbUXrHTzy8HGOq98AUoU/0GoqogOMOvLwUOV6SrLwwLHnGURHGHEUwDsMBtKFH3dE+QuiQ4w20vA/joadd0sZ/vwR3ypExxhppOHfOP5rGyBp+GOuryM6xkgjDH/NBybUrUQiUXg47RUal68gDP/r5T9dvCQOXzV8AIWLGJC+qocpqXWnp5cqPIS33b5EdJjRRRZ+dXWV8q704WGakrqzYRGFB8Pg1pVQDvARXrR9ap4o/KAv77v8IQd4OPWtgLuIjjSiSMJXrIE7lkFZwFfv2DuUrofzZG/jqzR/8LKA17RtBFUfmycIn3cXpq2BEocPW75HX+mIpRE0yeOHr7ilb59l2LB2w8LDw+vWEJYkvHVEnK6pfhccFw39VTNvSlakxA9v76/v6ZEFE+dpTwBbpyhNeMMLEYUd/F06WXfw9X+Bkhfo8cMbrtk/7RefZ69+UFHS8AXaFRhAsltRdAGiI25mxOBBVHzItnqKMoBfuEYxWubw20HU2KD6ijKA/8AureNOREfczAjBpw4BkV1z6yvKAB7s6z/whaXHdX2O6Mg3MzLw+f2LgL2fkaIc4EGhbWgPD+3d/Dd9EB35ZkYGfuFOcLjTHCNFWcAD21UbJw84C4BalvCaebuExBspygT+Pf9pG233yRJ+10EApk2Lkin8zrJJMceV8ZfkB3/UpgDscVPLFh6AEnApVuWB6Mg3M/zwHoosUGSbC2QMD8CcpMzWmYgOffPCD++SCcr8U4G84YtHDfaYOfgqooPfnIjcuXtjKpA5PADv2sASr+XkzpaCDf6Etb229j3bgnWBZUD28Jp79adupLocQnT8BYcNPjNBM9uyiDfAM5nOOhQOn+ey48ZbAxqcRwpPeOHTNZRP2+mPOIeHf83fCL8ZFU7kbCl44TVlPXZC/43c4Z2Lde3OKg4bXopIQUAY4U/6lADwkc3ThiMgc/hSq/G6hj3VfVRYz25BMTEbEFE0LXzwE/trjsE2u9w6NLnDOxs2OLw2Midne2ifqe6IKJoWPninUwAscSsAHF7XP/BbATg3cPO5BIvoIkQYTQkT/M9/aG7jS0aHlQAOr7+sBw+ORNp+pLowxXXxFUQejwwL/P0NHt9nJmQ7LnkQjcPr08EDUHxRlXuoImuk53I8H7DEAu+3ugpu93fOrofG4fUZ4AFQ/ZDoEgv/2u/RqrviFe2aRg9Z4g9N4i9N8lkN1PzCyZcD1PXROLy+v+G7xMTERMTEOL2wakeMg+f03WqVYBsBibw0SU22wu3d9PR4G/eRCUZoHF7f3/D9c3Rt72o7CYC8daMVbRNzbjXHtEmJvDTJhbm2yaujezvMTnbj8LoeDu9ouFRmaHkKzUQAAAWeSURBVH7qbRe20/3YokHuw5fn/Ngc2Ecl3tIk363z1b7i7DzZZuJxzXT4X7y+psEDoD6Q3TXy2d7uwwY4PdveqpcyYPjEmDXCfRfGGDJ+clC0pUlqJ73/e1X+XPeOKaW66XB4fU2FB9rr/uy0tzfmbG5lN2bZggmBCmuLVxPeO/b5f+E36YYOGw72lp5KQ8Yudvrbj5xpxj7iLU2SFe+nmnniXt0HZzi8PgHwdRvcF2WcAHud7J36ZC0fFuznoejaqZ+b/9AxMXGK3/XHOiXZ8E8b3B+s28cSM+GbvjRJbdqpCu0lh9dlDrzHRv2luvv4CHdF36HjO/WwtV6cHO2gaBXg5Wrn4Bhkrxw7KW7eouQUm8LiL8qv3fodPby+BxYVvqS/0nl1bf1/sbeFof9L0Oc4yPDNk4ZLZZDRBvcAow0qX6MNvt5GGwZ4Gb7593z9ZbC7YcNT8frLIS6GDS3f1l+GORs2tJqpv4xwMGxoPU1/GakwbGj7pv5yvI1hg+VU/eUka8OGDpP1l1N7GzZ0mqS/jHvZsMEqWn85q4dhQ9fx+st5zxs2dB9r+OY5w+ULkXXfGC57jjJ8U7ePXiO0X+eMG9XGcJD/x/nZjp0CQlwtHq877M8+27p16xd7dnkiWP8rIWjgH1hU2AC/sqT+v/ip7mZovuFy1QqjDauNNyQvM9qwxnjD2iWGb+YZLtctNtqw3njDhoV1+0jTX6YkGm9IqNuQqr98t25Dwkb95cYFhg2J7+ovU+sGlJiiv0yr27CwbkPdf3/RBqM5LFpvNMIl64w2LF1r9CvL1hhtWJ5s9CsrVtfbsH7tilUbly5dujF9w3jDfy7jJyTwPOZD9cwdj7FQPXPHYyxUz9zxGAvVM3c8xkL1zB2PsVA9c8djLFTP3PEYiz+Ol2kcXqZxeJnG4WUah5dpHF6m4Yeve1kSafFi7HT+RhF2ukiM+ScJP00mfvh+Ikw8vVeqCDt1WfbofyO4wFki7HTUJsEM+OFF+ZiAz30RdjpWjJMTxYtx9uqUPYJ/hcM/PA6PNA6PPg6PNA6PNA6PPg6PNA6PtPdF2akYbwDNqhBhp3nXRNhp8WXBv8KfuZNpHF6mcXiZxuFlGoeXaRxepnF4mcbhZRqHl2m44dO6WQy6jn63n7VCvkvPFi1aBKPe6VUfC9evEe9zte4UGEKfCcYM/9VTX98MRn+qzvuOTz76HwmsS3lFxV3UO7XOvDcP9YsVVRUVFaXOQj/chBl+hx+Ee9Avtbd8FHL4ypao96ipUAFhdbkIOw4VvLwJ9tv4+59HxaPe5xd9y5HDf2lhb+GD+iTiW4JHPT/oW8Q71XQiWvCvYIff/7TVz4h3WeNy6hfk8Gf9r1TNckC806THP74z1wnxTjUHwE74C8j479Xfy7BFvMfkqRA9vLa7/7qJdodpAZrbkMeQn5n49ADhv4MZPm0HhL8/hvhT1mNatnyqRctCtDstyNP8P/oE4lfkc/w1O338DtqdQhi3XfjvYIbf3++32uWu6PeL/i8+t/0XNQv9Ee+0qnNu7UIvxDuFsHszzm+N+6p+rlWHQBHu3YhwVb+qc/shyN8tU2RjGfAf1Du92KUZv8SfuZNpHF6mcXiZxuFlGoeXaRxepnF4mcbhZRqHl2kcXqZxeJnG4WUah5dpHF6mcXiZxuFlGoeXaRxepnF4mcbhG22s5X7SQxApDt94m3uSHoFIcfjGq2pzkfQQxInDN15lu1mkhyBOHL7xtrd8jvQQxInDN55zSruzpMcgShy+0Uos/oyJIz0IUeLwjTYhDuY9K8YJkonH4RvrdssrsKbzadLDECMO31jrAjVf3pxCehhixOFlGoeXaRxepnF4mcbhZRqHl2kcXqZxeJnG4WUah5dpHF6mcXiZxuFlGoeXaRxepnF4mcbhZdr/AxbmTPP/b3gSAAAAAElFTkSuQmCC" alt="plot of chunk unnamed-chunk-1"/></p>
<pre><code class="r"># dev.off()
# Check if the posterior integrates to 1, by computing the "moment" of "1":
compute_moment(thequadrature$normalized_posterior,
ff = function(x) 1)
</code></pre>
<pre><code>## [1] 1
</code></pre>
<pre><code class="r"># Posterior mean for theta:
compute_moment(thequadrature$normalized_posterior,
ff = function(x) x)
</code></pre>
<pre><code>## [1] 1.483742
</code></pre>
<pre><code class="r"># Posterior mean for lambda = exp(theta)
compute_moment(thequadrature$normalized_posterior,
ff = function(x) exp(x))
</code></pre>
<pre><code>## [1] 4.454407
</code></pre>
<pre><code class="r"># Compare to the truth:
(sum(y) + 1)/(length(y) + 1)
</code></pre>
<pre><code>## [1] 4.454545
</code></pre>
<pre><code class="r"># Quantiles
compute_quantiles(
thequadrature,
q = c(.01,.25,.50,.75,.99),
transformation = transformation
)[[1]]
</code></pre>
<pre><code>## 1% 25% 50% 75% 99%
## 3.166469 4.000544 4.404081 4.848323 6.149735
</code></pre>
<pre><code class="r"># The truth:
qgamma(c(.01,.25,.50,.75,.99),1+sum(y),1+length(y))
</code></pre>
<pre><code>## [1] 3.108896 4.010430 4.424279 4.865683 6.067076
</code></pre>
<pre><code class="r">#### END EXAMPLE 2 ####
## Example 4.1: Infectious Disease Modelling ----
if (dodisease) {
set.seed(8097968)
# use temp dirs
plotpath <- file.path(globalpath,"disease")
if (!dir.exists(plotpath)) dir.create(plotpath)
datapath <- plotpath
# the TMB template is part of the package. move it to a temp dir
# for compiling since this generates a bunch of new files
file.copy(system.file('extsrc/02_disease.cpp',package='aghq'),globalpath)
# Compile TMB template-- only need to do once
compile(file.path(globalpath,"02_disease.cpp"))
dyn.load(dynlib(file.path(globalpath,"02_disease")))
# Create the functions
dat <- tswv$tswvsir
dat$epidat <- dat$epidat[order(dat$epidat[ ,4]), ]
I <- dat$epidat[ ,4]
R <- dat$epidat[ ,2]
infected <- !is.infinite(I)
datlist <- list(
D = as.matrix(dist(dat$location[dat$epidat[ ,1], ])),
I = I,
R = R,
infected = as.numeric(infected[infected])
)
ff <- MakeADFun(data = datlist,
parameters = list(theta1 = 0,theta2 = 0),
DLL = "02_disease",
ADreport = FALSE,
silent = TRUE)
## Inference ----
# AGHQ
tm <- Sys.time()
quadmod <- aghq(ff,9,c(0,0),control = default_control(negate = TRUE))
aghqtime <- difftime(Sys.time(),tm,units='secs')
# STAN
stanmod <- tmbstan(
ff,
chains = 4,
cores = 4,
iter = 1e04,
warmup = 1e03,
init = quadmod$optresults$mode,
seed = 124698,
algorithm = "NUTS"
)
# save(stanmod,file = file.path(datapath,"disease-stanmod-20210503.RData"))
# load(file.path(globalpath,"data/disease-stanmod-20210405.RData"))
## Summarize ----
# pdf(file = file.path(plotpath,"stanmod-trace.pdf"),width=7,height=7)
# traceplot(stanmod,window = c(9000,10000))
# dev.off()
# Run time
max(get_elapsed_time(stanmod)[,2])
# Number of iterations
as.numeric(aghqtime) * stanmod@sim$iter / max(get_elapsed_time(stanmod)[,2])
stansamps <- as.data.frame(stanmod)
stansamps$alpha <- exp(stansamps$`par[1]`)
stansamps$beta <- exp(stansamps$`par[2]`)
posttrans <- list(totheta = log,fromtheta = exp)
quaddens <- compute_pdf_and_cdf(quadmod,posttrans)
quaddensalpha <- quaddens[[1]]
quaddensbeta <- quaddens[[2]]
# alpha
# pdf(file.path(plotpath,"alpha-postplot.pdf"),width=7,height=7)
hist(stansamps$alpha,freq=FALSE,breaks=50,main = "",xlab = "",cex.lab=1.5,cex.axis = 1.5)
with(quaddensalpha,lines(transparam,pdf_transparam))
# dev.off()
# beta
# pdf(file.path(plotpath,"beta-postplot.pdf"),width=7,height=7)
hist(stansamps$beta,freq=FALSE,breaks=50,main = "",xlab = "",cex.lab=1.5,cex.axis = 1.5)
with(quaddensbeta,lines(transparam,pdf_transparam))
# dev.off()
# summary stats
moms <- compute_moment(quadmod,function(x) exp(x))
getks <- function(x,y) {
suppressWarnings(capture.output(ks <- ks.test(x,y)))
unname(ks$statistic)
}
M <- nrow(stansamps)
quadsamps <- sample_marginal(quadmod,M)
summstats <- data.frame(
stat = c('mean','sd','q2.5','q97.5','KS'),
alphamcmc = c(
mean(stansamps$alpha),
sd(stansamps$alpha),
quantile(stansamps$alpha,c(.025,.975)),
NA
),
alphaaghq = c(
moms[1],
sqrt( compute_moment(quadmod$normalized_posterior,function(x) ( (exp(x)[1] - moms[1])^2 )) ),
exp(compute_quantiles(quadmod$marginals[[1]])),
getks(stansamps$`par[1]`,quadsamps[[1]])
),
betamcmc = c(
mean(stansamps$beta),
sd(stansamps$beta),
quantile(stansamps$beta,c(.025,.975)),
NA
),
betaaghq = c(
moms[2],
sqrt( compute_moment(quadmod$normalized_posterior,function(x) ( (exp(x)[2] - moms[2])^2 )) ),
exp(compute_quantiles(quadmod$marginals[[2]])),
getks(stansamps$`par[2]`,quadsamps[[2]])
)
)
# readr::write_csv(summstats,file.path(plotpath,"summstattable.csv"))
knitr::kable(summstats,digits = 3)
# Joint moment
compute_moment(
quadmod,
function(x) exp(x)[1] * 2^(-exp(x)[2])
)
mean(stansamps$alpha * 2^(-stansamps$beta))
#### END EXAMPLE 4.1 ####
}
</code></pre>
<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////isF19AAAACXBIWXMAAAsSAAALEgHS3X78AAAgAElEQVR4nO3deUCUZeLA8Te3326ZJ56YoGBqS3IoIocDjIAIZGkKeYui4q0ghqIQluCBmCeLiHmsmmaYmZYZeeSKhKIhb4em5mqh7qZUaroq8v7mnQuYGXLe4zne93m+fwwPvrPvPM/z2VCG4R2GoxEZg3oCNDRReEKj8IRG4QmNwhMahSc0Ck9oFJ7jLhR/X4V6DtAjHv7eUr8xcyb4pvyKeiKQIx2+LGDjI92H6g/9DqOeCtwIhz8U9JNx9Gv/rUhnAjuy4U9ofzOPHw7ahXAm0CMa/qr/L7U+ux96EtlM4Ecy/KOQs3U+v+5XiWgmCCIZPmO1ebg1jy/BNy9vB8IJwYxg+HN9H5uGj1zT9fUcnt4R5ZQgRjB81A4nN2Mv+rH6ij1KfVBPC1Lkwn+ccHAGa+y0EZ7NmEbhVV5V71s24M96e6GeGKSIhd+SydmAZ9e1QT0xSJEKX+V/xyY826QC9dTgRCr8zoWcbfguiainBidS4YMq64H3CbqJem5QIhT+yAyuPvidGagnByVC4Qdd4OqDf+T3APXsYEQm/JUB/K1teG7pdtTTgxGZ8PM/4W/rgb8Zgnp6MCIS/pG//ln6euC50eWoJwgh8uALYmJ6vxTDF1QP/PFpqOcIIfLgU/KKgncX8U2qB57rfQ/1JMFHIPy2YxqDcYIteM/S0tKZmbqb0q8foZ4qyEiET8n8E/jG0dHR/dvqbqLdvkQ9VZCRCO9T8mfw/I32kO5m8hHUUwUZgfDZkeyT4LNmU3jVlRK9+onwp3pQeNWV0uX0E+HZlz+i8GprQiD7ZPicKRRebQXOtgP+TA8Kr7acNtsBzw74iMKrqwsvbbMHPmcahVdXWa/ZBX/am8LX6cG1+0DmAa3QJLvg2ahPKbyp2/kaB4Zhmvvn3wU2H9D9t3+KffDZsym8sUp3xnn4jPnTR7gynr89+e54tjHXTvgSfwpvLLbxXuPoaMs4IJOB0OCf7IRng8YcQT1ZkAmA71jzivNFLgCmAqMHQZy98KkhR1DPFmQC4B3nmocrHQFMBUaFqXbDF3Y6gnq2IBMAH93mtHF0zjUGyGTAN+uE3fCsYyHq2YJMAPwV5wYByVm5y+YGP+18FdyMgBZQZT+893LUswWZkG/nKlKdGD6nNKX+YuHloZz98IOU+mXNrgQ+gXPnQslF5X4Tz63bKAB+UjfU0wWZ8Kdsb5QpV35QhQD4yaGXUM8XYALgEwp0N592YpintN8Amw/QHmk4IfAp/0A9YYAJgGemcdzxBg4Jq5PbNLoMbEIgK5rFCYHfMwD1hAEmED6k7Q3dqLL9EGATAtmCA5wQ+CNBKv69WYHwDoYncVI7gJkN4Prw/zoRAD/3KOoZg0sgfGvDVQOWNzL/6bYwfR5Zcs9M/m6H87cC4A/NQz1lcAmBj67gxvo+1I2qg7tbHty0Sb45gWqf/v+1AuD/F4x4xgATAs8wjV2YJI4r6cdstDyoBPjEYv5WADwX9cuTzqnYBMBXHN0wZ7DHaI4b+ly21UElwGv4r1aC4LM/QD1nYIl5zR1r47eIFQB/8xX9ByHwZycinjO45HqxpQLgCwxfp4TAV2sQzxlcBMFPM/xUWQg8N+Qy0ikDjCB4jeHN5QTBr7f6V6xaEgB/vHaWB3GHv195IapSn50vrzbAXxyJeuKgEgDflKmV5UHc4YP9unTw09dICDznj3rioBIAfz+LcdtpyvIg7vBadsT7Bk9PQfCx51HPHFCC/o7XaOs9hD+891kx8JvzUM8cUILgUxQMXxzMioG/Ohz1zAElCP7W5XoPYQ+/NkEUPNcb8cRBRcq3c9q4reLgx6r0L3li4H1Pi4NX61/ypMAH+bPi4C+PQDxzQJEC7zVRJLxa/5InBd4lVyz8yB/RzhxQpMC3OCEWfv1mtDMHFCHwj5uaPYXCfz8e7dQBRQj82Xai4asD0U4dUITA5/5dNDw36AbSqQOKEPhRvsLhBwyboy9k4BwVXrmeEHh/rXD4nunv68t8OT0d9fzljwz4G6+Jgd9k+Fjmsykd9QLkjwz4PVkS4Nne69JRL0D+yIBP/pcU+ElvpKNegPyRAR9yTwp83sB01AuQPyLgHwZzUuCL3dNRr0D+iIA/mSAJnu2chnoF8kcE/Jr3pMGHqPA3qYiAH3FZGnx8FOoVyB8R8P6cNPgsD9QrkD8S4G8OlAi/yRn1EuSPBPhPFkmF//t11GuQPRLg0w5LhQ//EPUaZI8E+IjfpcLHvYF6DbJHALz+lRTS4FNDUS9C9giAP8+/j4o0+HSt6i51SAD8lnWcZPiEEtSrkDsC4Ked4STD71yFehVyRwC8lr/KmUT4K8NQr0Lu1A//oA9/KxFefb9Po374kwn8rVT4QWp7Ckf98P/Yxt9Khc/6CPU6ZE798GMu8LdS4b9MQb0OmVM/vKaav5UKfy8c9TpkTvXwd/vpP0iF54Kq0K5D7lQPf8zwNVoy/KRytOuQO9XDZxt+sCYZflM+2nXInerhh/6k/yAJPm9UYeG7kYV8LOLlyJbq4Y0XHpcEP9c3MTHh+US+YKSLkTF1w1+5VBZ6SZ/52kei4Cfrbnp/xQ/rv8SjwlI1/D3H6CC3aH1/kwwfv4nC2whL+LtBbEKOpadY+NWJFN5GuMKHH5EL/khfCm8jXOF7WHmKhTecisJbhCn88RD54MMPUXjrMIXPmyYf/KzVFN46TOFn5Fp5iobfFE/hrcMUPvyofPAlGgpvHabwPaw9RcOz3cspvFV4wgf0kRN+4H4KbxWe8O7TrD3Fw6ctpvBW4QnvslZO+PdHUHir8IRvdVhO+K97UXir8IRvasNTPDzr8zWFtwxL+J9byAs/bBeFtwxL+P0d5IXPTKfwlmEJn/mSvPB7B1N4y7CEj/G14SkBvrwnhbcMS/iAIHnhWf9SCm8RjvC3w+WGH7uNwluEI/yXyXLDZ//J+2krLDXDr9omN/yBVyi8RTjCx7Jyw7PeFN4iHOGD7sgOr1HNm9CpGP5++F3Z4Sd6oV6VXKkYviRJfvjVnVCvSq5UDL9um/zwh1qjXpVcqRg+/jv54dmmqFclVyqGD64CAN/iJuplyZR64R/14QDAd/wc9bpkSr3w7GQQ8O6LUK9LptQLvyUfBHzvwajXJVPqhU8oBQGv1aBel0ypFz7sPhD41/6LemHypFr46iAOBHyLTt3DDEU9Rr1ESakW/sfRQOB91pl+R8PvEeolSkq18AUrwcCbf+OewuvDDj71KBh48zU2KLw+7OD7/wYIPvwwha8VdvD8D86BwCespfC1wg3+P4M4QPD5Uyh8rXCD/yyDAwRf1IfC1wo3+MWfcIDg2e4Uvla4wQ+5xoGCjzhE4WvCDV7/okgw8PrLnlF4YxjB/7xkyZK3uy7hb4HAb5hE4WvCCH5fdH7+nMh8XSuAwBcHU/iacIJPYtn5WbzNMSDwxufuKLw+zOCjPwYIH1VI4c1hBt/zLED42asovDm84Mt6sQDhN00kEP7Btfs2/xwv+D1DQMIXBxIGfztf48AwTHP//LtWx/CCX5QOEt7w3B058JXujPPwGfOnj3BlPH+zPIgX/OgdQOGjPicKPrbxXuPoaMs4y4N4wfufAgr/xkqi4DsmmoeLXCwP4gVvfJUMKHj9WxaQA+841zxc6Wh5ECv4g/3BwpcEEgUf3ea0cXTONcbyIFbwq5PAwuu/opADf8W5QUByVu6yucFPO1+1PIgV/NR8wPD9D5IEz1WkOjF8TmkVVsewgg89Bhg++R2i4HXduVBy0fqbeA4zeNMroIHBbxlPGrwinrmbZnonGmDw/BtSEQSvlGfuXp8IGp7tUU4QvGKeudOuAA7/ymcEwSvmmbuXDgCHn7OcIHjbz9yV5umLzZZzVpLa174cOPw/xxEEb/uZu+936ZuySs5ZSWqP+Q1JwMGf0hAEr5Rn7laa31cUHLzuX3fkwCvlmbsZ/SDAD/iUHHilPHP3yigI8CnZBMFzynjmzjMBAvzWOLLg6w0f+OpuSRDgTwVQeH34wF8KhQHP9qDw+vCBL5gABX6gFzHwx2tneRAf+LRFUODndVYZfE69V25sytTK8iA+8P13QIHf7qgyeObpfpt/t3nX+1mM205TlgfxgQ/cBwW+tKnK4MvmvcA8M+gDmz9119T/5lvYwN8cBAeebfwA9VIlZevv+DNzXZnGow5Y/z/6T95nERv4z9+GBN/6G9RLlVQ9/7j7lwfDtJxcbvGnty7Xex5s4JfthQTvshn1UiVlC/7X7TGNmFbjRzRssNHu82ADP/wqJHj3GaiXKikr+J9y+v4f4zzzaJXuP3D/tnafBxt4DQcJ3jcE9VIlZf2veubFlFPG8fx2dp8HF/h74bDg/YKqUC9WSlbwGd/VjB/Z/w9XXOC/mg0NfiKLerFSsoJP+NHw8dSbgs6DC/y6rdDg129CvVgp1YWvqqpijlbxPVzcUNB5cIGf+C00+NNTUS9WSnXhaz8pGyToPLjAax9Bg7+v6PeSrwufkZHBTMjQt/SCoPNgAl/Vh4MG/0j7EPVyJWT1d7xW3D9ZMIH/diJE+OlnUC9XQir7efz2XIjwW9ajXq6E6sI38+C05gSdBxP4pBKI8N/Go16uhOrCt+yucPi+f0CEfxyMerkSUtmX+mAOIjwXdg/1esVXD/z5c9XCzoMH/JURHEz42cWoFyw+K/hHCyM5bgrD+N0UdB484D/if3UTHvzONagXLD4r+AVMCFfGvLy2cZKg8+ABn36Igwl/aTTqBYvPCt41Qrd/f7nOjXtB0HmQwy/v4efn17yn7sYVGny1gt9M3gr+mSyO0/biuJxnBJ0HOfysXazxvSMmQYPnIm2/LFUJWcF3jOUq/jKf41Lt/1k8Hxbw/wqFCN/a18/veTfdVxi/l04iXruYrOATn0kLeOrs/fymgwSdBwv4DVMgwjctY1nDJTRnHES8djFZwf8WwTyVxp1jnL6zdfd6wwJ+1hrI8If7qQee4279quMvvCPsPFjAR30BGd7wlgVqgRcVFvCGd32FCR/ypWrgH2WGafRNEHQeHOBPaqDDz/iHauBnM26hYXyvCToPDvBb46DDr5+iGvh2S0WdBwf4ednQ4fVvN6oO+BbXRZ0HB/hB+6HD658yUgf8gEOizoMDvPdZ+PD8WxaoA/5ckKjnoTCAP2PUgwo/N1st8CMDmQ4aBb4CZ9auXSMRwG8boxZ4pb70atau9EwE8KW+aoEXGQbwQ/YggGd7lKkH/ucztwWfBwP4XmUo4IfsVgv8F24Mc/zLwGPCzoMefqcviwJ+4QKVwBf99cWlzPFLbn89Ieg86OGXv44E/qNolcCHd/6DY45z93v2E3Qe9PBTFyCBP+utEvhmGRwPz61sKeg86OGj3kcCzwacVAe8Y6YBPtP+69/woYfvdhoN/Pgt6oB/tYv+S/3tTq8IOg9y+EQ3kwhk+FVJ6oD/ruGLq5i3Mto/K+z6fcjhx4Qigj8crg547qQ/fz0MbYmw8yCHj5yACJ7trhJ4jrt54rTVO0k+KeTw3otRwYePUwu8mJDDO29HBZ/UXw3wdxaHd3q2U/higS+yRQ5f3X4XKvjNPVUAv78V08Knv09LptUnws6DGv58N2Twp5yVD/9Dk44fPNZ9rN7t2vSioPOght8Rggyeff4A2rWLqi786Iama15923CMoPOghk8ajg7efR3atYuqLrzbq+bhwG6CzoMavu8MdPB9E9CuXVR14RvMMw/T/iLoPIjhqwNnoYMfFYF07eKyuKRphnm4RNg3eojhz8chhJ8h7IsjHqkEfnsOSvhelUgXLyoL+DFfmBqnKPikr1DCj1Lg93P1X71aUfBh91HCL3oL6eJFVZd3Xe0EnQct/OMgDiX87iiUixeXOp6r/348UviDGoGXg8QgdcBvyUMLP+57lKsXlTrgZ5xGC7/B/jfowyV1wIc8QAv/jbDLh+CQKuAfaTm08I+Vd4lLVcCXTUUK39fFu6mXN59nLsJdEJYq4N/diBTeP4+dmqcfbZ+LcBeEpQr4id+ghl8/mcIjKLgKNXyxhsLD734ohxqe7XGWwkOvmH9PBcTwQ3ZTeOit3skhh1+cSuGhN5p/XShi+M/6U3jo6X9Eghie9aLwMEuLiYkZ6Ki7iemCGD78CwoPseCioqJVY3U3RV0Rwycvp/AQ07LGS4ezbojhd4yk8BDj4UOP4QD/tTeFhxgPb3hfCtTwbMBJCg8vHXxhJB7wkzZQeHjp4Jcn4wGfN5XCw0sHH7cVD/ivNBQeXjp4v1I84NmeZRQeWlr2bE8WE/iROyk8tLTsB8NwgX8nmcJDS8umLcYF/lhfCg8tLTtgPy7wbI9yCg8rLdu9HBv4wXspPKy0J4JZbOAz36TwsNKum4YP/MEoCg8r7eR8fOBZLwoPK62mGCP4V1ZQeEgF9zDtPg7wC+IpPKR8onGC3x9E4SHVdSFO8GwXCg8px4+xgg+YjGQXxKRw+KblWMHH9UeyC2JSNvytFmYHLOCz3VHsgqiUDb/PFS/47e1R7IKolA0/zwszeLfLKLZBTMqGDw3EDD4CyTaISSj8g2v3bf45EviHWi1m8BNHI9gGUQmBv52vcWAYprl//l2rY0jgv0rADX5uAIJtEJUA+Ep3xnn4jPnTR7gynlbvS4cEPrsAO/jRFxDsg5gEwMc23mscHW0ZZ3kQCfzga9jBb1bK+9MIgO+YaB4ucrE8iAQ+kMMO/ufXEeyDmATAO9Y8Eb3S0fIgCvgfxuIHz/k/hr8RYhIAH93mtHF0zjXG8iAK+I0bMISfVAZ/I8QkAP6Kc4OA5KzcZXODn3a+ankQBXzcOQzhP1gOfyPEJOTbuYpUJ/17ljilVVgdQwGvqcYOPqHyYnilvj/g74egBD6Bc+dCyUXrb+I5JPDXoznc4Be6hIc3DwvnC4S+H8IS/pTtjbI68tvD9HVLl2tGdrdrJXbwC8azbNxm/VALfT+EJQA+oUB382knhnlK+43VQQT/xU8/jSV8/kS1wTPTOO54A4eE1cltGl22PIgAPrAKS/jSXmqED2l7QzeqbD/E8iBM+OqPduna5K278cQPntUUqRDewfAkTmoHy4Mw4StfSNT1qkZ38yyG8EkrVAjf2vDms8sbWR6ECh/K72zs9jpM+MAXRKsNPrqCG+v7UDeqDu5ueRA+fM+vMYUv91AbPMM0dmGSOK6kH2P1PmvQ4YsD6zLhA8++ul9d8BVHN8wZ7DGa44Y+l211EDp8zgxs4ZfOVRe8Ofae9Z9Bhx+zFVv4Y1qVwtsKOrz+r3g84dlepRQeRDy84a94TOEnrafwIOLhVydgDL99FIUHEQ8/ajvG8Ge7U3gQ8fA9yjCGZ6M+o/AA0sEf62PJhBX84nkUHkA6+OxkrOF139BRePnTwcfsxhqe9T1F4eVPB+9Zjjf89DUUXv4qQw9GWTHhBV8QTeHlrzJ0wduYw7OeFF7+KkMjPscdfqg3vP0QlSLhQ8yXNcQWPrcjvP0QlSLhew3DHv50M3j7ISpFwndajT082+oKvA0RkyLhHb7CH/7vq+FtiJiUCP+zgw0m3OA1EfA2RExKhH+vkwLgtf0q4e2IiJQIP76XEuDXbIW3IyJSInyvUCXA/zQY3o6ISIHw54cpAp7T2vx9clxSIPyqDcqAX7ob2paISIHwUReVAf/jMGhbIiLlwd8NrVQGPBdk++qveKQ8+I8zlQKf+RGsPRGR8uAnlikF/ofhsPZERIqDr/bnFAHfo7Cw0H2/7qawCNLOCEtx8KcnKwPeKS4urmeI7ibOCcsrnykOfsGnyoDn73gggh8FYfn9vOLgA+8rB571KaHwMvXTIE5B8MlLKbxM5WxWEvyhPhRepiJ/URI8qzlG4WXp176couAz5lN4Wdq+ilMUfLEvhZelaP5C+QqCZyP3U3gZ+qMPf6sk+JzJFF5i38bHx/ftpbuJj1UQfJlHOYWXVsHkAwe0Gw/o2qQgeDb2XQovrYIUttRbv5efKQl+z6sUXlo6+JUzlQfP+vSm8JLSwUd8okD4tM4UXlIFKad6sgqEL25G4SVVkLIsSYnwbJti0FsjJiXB9zmoSHivCaC3RkwKgk8MYBUJH+SP49d6BcFHpSsUfsV60HsjIgXBuxxXKPx1HN91UjnwuV1NO6k0+LuTMXyhrXLghw1WLPw3GL7AXjHw1V2SFQvPRVq//TbqFAN/PCpFufB70wDvjvAUAz8uS8Hwj/1svH8T2pQCfyewQMHw3BrsvqNTCvy7OYqGvxNQDXZ/BKcU+D6/Khqem78X7P4ITiHwZ2M5ZcPfCAO6P8JTCPzkEwqH5yYdA7pBglMG/G0Np3T4S/1BbpDwlAGfk6dgeM3PlXxDDlf+BnKPBIY7/I4YvpaDYmL8lQrfspcfn6eDX4eDYDZJTLjDT1134MCBRYN1N4OVCu9w0vAxatfsj8Fskpiwh9+j27I+n+luxigdfl8Ihbc/Hn5/X1YN8OyAERTe7nj41/+pDvgDLhTe7nTwx3xZdcCz3d8Cs0liUgD8xFVqgZ/y0mMwuyQi/OGLvc6qBX72iC1gdklE+MMnZhp2TQ3wH/hhc7FD7OF3eho3WA3wH297E8w2CQ97+Ni3jLumCvjq0Mtg9klwuMPHd/nauGsqgA990c/dQf/8rUcKmO2yP9zh/cw4KoAPXsuyw3L40a5ZYLbL/jCHv+W427RrKoEv9uA/pfBP6I2oPaZdUwk8u3wUS+Gf1OWwqaqDZ/tup/BPaliJCuEPdy+l8H/e8VGcCuHZpaMo/J9WFVShSng2Mp/C2+5r/k18CqfHFxYOUCN8kce7FN5mvRJ1TWg3MzGxpRrh2X96U3ibafndCeeFOqsSnh3RW87tEhO+8MtH8jukUvj3O38u536JCFv4Lz31+6VS+F1T/S/KuWHCwxY+bIt+h9QKP+tcb7S/XoErfFq8YYdUCp/XJSa47WD9L4sgulgGpvAf+hq3VaXwS4cWFaVGHC/SFSznxtkfnvAar8+NO6RW+NG6m1kj+KFWzo2zPyzhq1utM+2QmuHZ+Akshee4h/l5xl5uY941VcOzY+MpPMdd91puaLhPQ/OuqRueHTe6nMJfjzRsRn7gmZrdVzk8O3NAsOSNExV+8Nt7FrPkwLNvt/hd+u6LCDv4nT2OsSTBs569kTyHhxv8tu5HWbLgtecD9kvff8FhBp/rc6zu7hMAz90envxAuoDA8IJPD/7KYvfVD98iLCzMrbm/7jYQ5gVvcYKPGDG0zHL31Q+vv+OnAdPOsGG3pDvYHUbwp1sssN59QuDZ8nSvfELhd/pobOw+KfC6yUa3hvkGdbjAV8TMvBxpY/fJgWfZ1s1beel/pdKvC/jX5+AB/79lgcXmZ+5YYuGL3guJfI8fpRTIxFJ/yOH3d/Xz7dTIxdfPz5PCs+yeaJ+Mk0TAv5eS6pGg35oPKTx/e2JujwFD35eJpf5Qw5dHOKWWGtZN4YuMgz0BXWILAL8kTyj8g2v3bf65QPicMD7/Ts3aubxpWjeFN8Gzs7I+mOnvMy7nWGnpDYFA9iYE/na+xoFhmOb++dbvkisQfsL+w8tf9xy4upQdRuFN1cCHh8TFxcVGdm/fpmPQl5WCdtbeBMBXujPOw2fMnz7ClfG0+jokAL7q0r5Frm595xSU82uk8OZRDXzfd4yD8gXNnm/euGlr5y7dvHvZT2VHAuBjG5veUOdoyzjLg0+G//3fp/w6tmvVtFHjFs93brbftEYKbx7ZgGdXDdPdnPk4Z96YiEYvvtQz5LVRU95If+fCL5KvlycAvmOiebjIxTwsM7xMbszrphfMrVnClxSs8fPz6+nl5da5UZNGjXU1bdGm4cCRU+al87Wclm6s18um0aRuplH6M+ZRm8mmkX+kaTTjRRt3dIw3jYLCTaM3Otu4Y/s40yg0xDSa72rjjh1iTaOIYPMfdjSPnk0zjVxHmkavaMyHO5lHz6WYRl2GmUaD/MyHa+bYONk0cosxjYb4mA+3ad++XdvWLR2aNWng0Fy3p7p9bfJs8zbPt2/f0dXVtVPztm3bOj7P19rd3b2HN1/XsH7GIt7Oq9N6AfCOc83DlY6W8FlLd5nas3t3YeHOhITEefPmLcjMzE42P1qaeZSeYxotW2Ia5abZuOOCtaZR9mLTaJ2tO761xjR6J9N8x1Qbd3x7lWm0wnzHvFQbo4yVptHKhX96x8wVptGqt20cTltnGi16xzRa85aNmaWa77gk2zRau8DGHdNyzXdctnbFiqWZmZmpc+LnzUswNDzeVPR402hink5mfV5NAuCj25w2js65xtj/P6NhmQD4K84NApKzcpfNDX7a+Sq4GdGgJOTbuYpUJ4bPKQ2/d8WmCUzgEzh3LpRctP4mnqa85HrKlqawKDyhUXhCo/CERuEJjcITGoUnNApPaBSe0Cg8oSGC35WTB6F182E8Sq0fywJttawCiOBfW/jkhcqwVe4wHiVv4EwoD+MhqwAi+Ak/wHiUuy/DeBQu8yCUh5H3KkkUXnoU3v4ovPAovN1R+Pqj8NKj8PZH4YWnCvj9UK7qVwX+l075jv8bysO8J+vZ6DN3hEbhCY3CExqFJzQKT2gQ4K++3s5lbKWNT82jCv0v6DAJYB6G4/7X5Litu4B5FLCLKenb1qHPFzbuIjTw8Nc6tJ4zs5HHHatPa0ZHmehpuiRd6qneh9GVxRy3cRdAjwJ0MScbdFyY2ZXZIX0x4OHTGpRy3F4m3+rTmtEGRvov49X7MLfmBzIGEou7AHoUoIsJbfwfjrvXtY30xYCHdwrgb501Vp/WjOY0Bvgwl7VaNwOJxV0APQrQxTQZxo8WMFckLwY4/B0mmf8wzMHy01oHBrvEtW0RVgTkYfh26kks7gLoUYAu5l6q/lqnscx/JS8GOPwFJov/MJ15YPFprQMeTK+MOW0afAriYfgPBhKLuwB6FOCL4bhTTYKlLwY4fAmTy3+Yz+BI9swAAAHISURBVFyz+LTWgdFJDznuhsMLIB6G/2AgsbgLoEcBvpjqjc+1PC99McDhLzLL+A/TmXsWn1oc4LgERsLF/Op9GP6DgcTqAYE8ijFgi/nGlwm/IsNigMPfZfTXTBrR3PJTiwMct5opA/AwfAYSqwcE8ijGQC0m92+d99m4i/DA/6u+vf7nyK7+Vp+aR+VDS/hR8lNSrrVR78NwZhKLu4B5FLCL2flUzG2bdxEcePjUp8/xz2rkW31qHv3RMKKa4yrbBgN5GM4Mb3EXMI8CdDGPXbrXcxfBgYevcHZdtaSVp+4/gE0vbKr9ac1oGaN9J82pEQvkYTgzfJ0/A/YoIBdzhvE3XMruN8mLgfFcfbRjxzG/6gYrmZW1P601eterkVPMj4AepuZv39p/Bu5RAC6mgDH2s+TF0J/OERqFJzQKT2gUntAoPKFReEKj8IRG4QmNwhMahSc0Ck9oFJ7QKDyhUXhCo/CERuEJjcITGoUnNApPaBSe0Cg8oVF4QqPwhEbhCY3CExqFJzQKT2gUntAoPKFReEKj8IRG4QmNwhMahSc0Ck9oFJ7QKDyhUXhCo/CERuEJjcITGoUnNApPaBSe0Cg8oVF4Qvt/yXf+7E3lBjMAAAAASUVORK5CYII=" alt="plot of chunk unnamed-chunk-1"/><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////isF19AAAACXBIWXMAAAsSAAALEgHS3X78AAAgAElEQVR4nO3dd1wUR/8H8LU8aYoa7AZREBVRUZqAHHIqgi1REayxYcRE0QdsWEA0Aio2rGhsj3lssYRoEnlii4+iBiVG3BQT9WcaeTSxJGrsOr+dA/TY4+Bmb3dnZnc+f8wt7M3uzPf9Ori73bvlAIsuw+EeAAueMHidhsHrNAxep2HwOg2D12kYvE7D4HUaBq/TMHidhsHrNAxep2HwOg2D12kYvLQ8WL8G5iTucUgOg5eWS17JQiZH4B6H5DB4abnUhxeSw+D1Fgav0zB4nYbB6zQMXqdh8DoNg9dpGLxOw+B1Ggav0zB4nYbB6zK/J/rW9A1Zco7B6yv/DtrzQx/+yKigjxi8njJj9P3CP/U7WnbBPRbJYfC2587OHTDjpoLi//G7HH/EPSipYfC253BgvJCglnC56MldaPDfuEclMQze9hweI1jvbTkeLhc/q9/xNu5RSQyDtz0Q/qzfsni4/Ozl3LD9uIclLQze9kD4SQk7SsJfC7yHe1ySwuBtjwB/yCdfBA/WpuIel6Sgw185e0eBcdAQAb7XJl4M/yT4Cu6BSQkCfNwuodnXhOMqGL9WbDwk5/CYD7vwFvDgkwmYxyUpCPBcLAA5FR3jlk2pW/WyYgMiOIfHdMkqBR4YaXzII8J3qgcnecOpv2IDIjiH+3XnS4PfOw3vuCQFEd5xqmkxsZEyoyE7h5tmlQr/NOgW3oFJCSJ8nRTT4qKqCg2H6Gxx4UuFB+uWYR2XpKDARxaAEf4PhaWnIV6KDYjgRPSyAn838CnWgUkJCjzHObhwEwHIDec2KDciYvNX6zGlwz9p3KBFQEBAH7zDQwwCfMGRdQl9PYcCMKDKQuUGRG6WT7QC/9h/X5jwgxHv8BAj5Z07/q7sw6AhQdkm+K19Dwh5vwQ8H/RfPcDrM8fGmY7O8ene8OBsVEn4+ZP0Ab8x49nioxum3JRrQMRm2LlC+PlDYLuxJHyetz7g3Z532xpqSms6D1XYnr9CQBnwfMROXcCbPeKLf7PR3pEQnvUry4RfP0IX8JbRPHzotTLh8z3ztQz/0++Ft3/+Kl6jdfjLfUGZ8PyQ97ULv60ux4V8D5fiLLppHT51dznwm4doFv4g5z5rWLXaPwM9whvulQZ/qC58YtsZwp9ro1n4EPe7AJyuEgZ0CP/dMFAa/L5Q2H4F4fnB3rgHiRQE+FfTYLuB+0SH8MmfgvLgNzjjHiRSEOALD8k+bdf0lv7gDQ9AefD51XEPEikI8KEuV+EN/0KU7uDPD4dt2fB8ve9wDxMlCPBfVHQYCD89sKqCV6DO4FP2wrYc+JbzcQ8TJSgv53I7VUuEt9tcOJ3Bh5g+NFEOfHAo7mGiBPGdu8IDsk/O7BKv0DT85cJzS8uBN3a7jnmcKGFv2dqQpVtMN+XBZ2zFPE6UMHgbEn4qBqZnOfA/DMU9UIQw+PJzMzxrdLaQmHLgQRBF51wy+PKzbUnWFCg7pTz4sV/iHqrtYfDlZ/BFG+H3UnQ2CoMvN48NwEb4W2G4x2p7GHy5OT7RNvjmAQHV/QJCcA/XxjD4cjPjkG3wfjwfu4qaE7AYfLnpcN9m+M1DGbxmcqU3sBn+rDeD10w2ZdoOz3c8yuC1kkGXEOAnL2LwGsmT9gABfld/Bq+R5I0FCPD53gxeI0nbAxDg+c6BuAdsYxh8OQn7C6DAT/XAPWAbw+DLzt+dYGs7/C4H0xHcWbjHXW4YfNn5LBm2tsPnV4UHcLPbYh52+WHwZSd6cIKQCJvh+RpHTcu4x11uGHzZaTb6AyG9bYdvuJTBayDX26RAx7dsh/cYyuA1kN39UeF9/Rm8BhKbjArv1z6XwdMfw3pk+LfXMXjqc7XPRmT41WMZPPXZnoEOf9LA4KnPO2fR4XnvfAZPewxPJMD338XgKc/VCCABfm4ig6c8H2RIgd/3OoOnPGO+kgLPezN4yhP8RBJ8p6MMnupcewNIgo9bzuCpTtYCafAbYhg81ZnwhTT4U0EMnuoYH0qD573zGTzFudMVSISPymLwFOcAPGNSEvyc2Qye4iQfBBLhP4pk8BQn7A6QCJ/vzeDpzSPTGfWS4Pn2PrhHX24YvLUcHpQnBP3UK9iMdMc9+nLD4K2lt2e0kGaS4DMa4h59uWHw1lJ4YnWEJPhDjrhHX24YvLW4LpUOz1fFPfpyw+Ct5Ac/e+BrXMM9/vLC4K1kYz974F/7D+7xlxcGbyWjEuyBbzYH9/jLCyr8g9/ulfp7zcEbFtsD79UL9/jLCwr8rbUGR47jXg1ce8dindbgb/ZcYg98y1YHhOTgnkUZQYC/0ZpzHjR+xrjBrlybP8UrtQb/aYpd8FWajYqPjyf5SnQI8MMc9hQtHakVLV6pNfjEQ3bBO8Qv58k+5RIBvnH8s8U0F/FKrcF3uW0f/Lp3tANff+qzxYz64pUag3/cAdgHf9KoHfjIusVX3jjvGiVeqTH4r8bYCQ9PrtcK/E/OFdtPSc9cMDWksvPP4pUag1/1b3vhux7SDDwoSGzIwTRMKrBYpzH4Ny/ZCz9huXbghdy+kHvR8kU80Bx8ELAXHj670xC81WgL/mpvu+GPG7UIvzHD4jcb7R0JSdkzz2540j86KQ3e7Xm3TT6muCTLNCAC8kVCwKCEhM52wod9rkF4bT/iZ03yWLF2rZed8P9cpUF4y2gKfq2PYGfvI351rIbgf/q98PbPX8VrNAWfPFAG+KOhmoHfVpfjQr6HS3EW3TQFPzRVBnjh2Z1G4A9y7rOGVasN37TTOHzQJ3LAd8zRCHyI+10ATleBF87VOLzbOTng31mnEfhX02C7gftE6/BTPXk54JdN1Ah8nRTYPm3X9JbG4Yf0kgX+QHeNwIe6XIU3/AtRGofvGC8LPE/0Z2YR4L+o6DBwv3C7qoJXoKbhmy2XBz6I5M/Morycy+1ULRHebnPhtAz/1Hm9PPAjWuCeShlBfOfurql9cmaXeIWG4C+0lQl+fiPcUykj7C1bi2zpIRP8ntq4p1JGGLxFxsfIBJ/vgHsqZYTBWyQkSSZ43uER7rlYD4MX536nWXLB1/ka92Ssh8GLkztBNvjGm3FPxnoYvDjLtssG7zEJ92Ssh8GLM+hH2eB9uuKejPUweHECgWzwfkG4J2M9DF6Ua2/ICD/A4kwlYsLgRcmeIyP83E9xT8dqGLwoyQdkhM9OxT0dq2HwovS8ISP8b/1wT8dqGLwoHYCM8HBrhIbBl8ylobLCd72Ne0LWwuBLZvtyWeEnn8Q9IWth8CUTf0pW+M2rcU/IWhh8yXS8Lys8/w7uCVkLgy+Rh0YgK3zhZS5IDIMvkTOxQFZ4EPIU95SshMGXyJr3gbzwIy7inpKVMPgSiT4P5IXP2I17SlbC4EskGP5llhP+8EzcU7ISBm+e26YD6HLCX++Ne05WwuDN89/psJUTHgRjnpK1MHjzLMiCrazwPf/CPCcrYfDmiTJ9Zaes8NOPYZ6TlTB48wTdgEmQE377CtyTKj0M3ixXmwWGCXGUE/67GNyzKj0M3iwfG3aYvOSEf9wR96xKD4M3S1Kk/PAgZEkUzADCjswzeLN0G6sA/EjjoeNCwgg74ZbBP89TQ7wC8Et9TsHlHgye2FwYrgT8500ZPOHZukIJ+Ov1GDzhiTstM3xd0ze6V2bwhKfTfZnhC9t/fM7gic4jI1AE/uW1DJ7onIlVBr5qAoMnOu9tUga+egSDJzojv1MG3tGXwROd4CcKwXufY/AE504YUAi+RzaDJzjwtCtl4OOWMniCs/BDpeBXjGXwBKdfgVLwn4UxeIIDv6FKGfhTPgye3PweARSDN5xi8KTmjyWxeXl5byoDP3yrFuCvnC3tAvK0w4906xAZGVlDGfjU2VTDx8HLUuxrwnEVjJbfykw7fLTvCcHIQxn4nYMpgF/5u9W7xgKQU9ExbtmUulUvi1dSD9+KVw7+y0AK4LnK4f8q/VM/EL5TvSvC0g2n/uKVtMNHdVYQnvehAP7sdDfupYid90q5qwDvONW0mNhIvJJ2+JCxSsJ3P0g+vJAzU105hyHZ4utqQPjCq02CRVXFfWiHb7lUSfi4TCrghRzz5Lha75wredfIAjDC/6Gw9DTES9yBdvj6HysJvyqeCvibW6KqcrXfGvxKxQ0l7spxDi7cRAByw7kN4j6Uwz+sn60k/IHu5MP/srLLPzjnfx55DMD1wHrmawqOrEvo6zkUgAFVFlpsh3L4Lz0Uhed9yYfnOPdpp4uWZzQotQ9/1/J3lMOvDlEWPuAM8fAp3z5ffvTA5u1QDh/dV1n4QbuIh4/7v8Lb02V9XdPGDIvfbJRrRFhiiFYWPjmNbPjHjx9zRx7DPJz7Shm93J5322T6tIiPS7Iy41Mnt7orDL9lBNnwnFnK+o59rT3iDycpDH/KQDZ8SkoKNyrFlPkXkLZDN/y8jxWG533Ihhdi5CVth2743leVhu/ShXR4iaEbPggoDT/Wn2j4Gp7A+CxI26Ea/pcBisNneBANX8urDPgc84hXUg2/a7Hi8PuciIYvM9XNn/OLV1INPyVHcfhzNemA//685ZUV7qVzHtuLI15JNXyXu4rD8zV/xD3LkrGAfzSnGwBjOC7gmuWdDdb/79MM/zgEKA/f6DDuaZaMBfwsrhM4y/VY4TDR8s7TtAl/dowK8K1W4p5myVjAu3YFILnS/8BIN8s7X79sdTv0wid1btEyNLSu0vDtY3FPtGQs4F9KB8DYDoCVLyFth174iJyIvcJfYqXhmzqGCjH8hnu6xbGAbzwMFFSaAUBi6cfirYVmeK98FeAbeMJ2wDe4p1scC/j4l5LaV8i/t7Z6BNJ2KIbfb+DVgA/KIRv+z65chSRwnmv4bWl3txqK4TPeVgV+yAay4YWncDcF/gOI37JNMXzMKlXgk6eSDi8pFMMHH1UF/v0IsuEfpYYaTBmFtB2K4U3PupSHP+hHNvwkzqMzfOER2gdpO/TCh4erA3/IL59o+AbzJW2HXvh28SrBR+4lGr7m/yRth174JutVgp++gGj4XockbYde+NpHVILfNIpo+PMdTknZDrXw92vnqAR/MoRo+DeDuUYGHZ16dbypWvC8N9HwejvnbpG/avAdj5EMLzHUwvfrphr82xsIh//1zC3k7VALHxShGvySBKLhD3pwXM5/g4+ibYdW+IJ+6sHv600y/PEX3OdzOZc8XjiBtB1a4XcuVg/+nB/J8GFN/wZcDrjnG460HVrhJx1XD573P0swfI0UAOFBRi2k7dAK3/G+ivD9sgiGr59aCJ9ar7S7Ww2l8A86AhXhk+YTDP9GM9Of+ltNXkfaDqXwX0xQE35zNMHw377ivpSbneL0suUXFZcVSuEX7VIT/pSBYHhwKhB+Ns6Yi7YdSuEjf1UTnvcmGR6Aaye+/BN1O5TCtweqwof1JhpeSuiE/3mAuvDjQoiFvz03rMnLTcLmIp5kSyn8tmXqwq9oQyr8J7W5mn49/WpxtT9F2w6N8JsSfIcmJDRXEf5gI0Lhf6jWeOcT4fbpbtfqF5G2QyN8523uWz/4wFFFeL42ofBDXyn+zqtvXhmOtB0q4U8FChZ11ISv9yXuSRenJLzHG88We7dC2g6V8OvfUhve3eKrRHClJHzF6c8WkyohbYdK+HEr1YYPSMY96eKIvtI05dniPLQXelTCd8hRG77bANyTLo6u4dvyasP3b4d70sURwQ8/WJyR2ocPiFAdfoDfE9yzLor1b6/WPrz7u+rDR57HPeuilORdbR6k7dAIX3jlKXXhE7fhnnVR9PxeffV89eG3JOCedVF0DP+/OurDG/o0jBFSgHvuuobf1VR9+KZrWmVnZ/c/jnvuuoaP88MA/1FX4XYUg8cZYycc8PErGTze3AnvjAN+TSyDx5uDyVjgj4YyeLyZtR8LPPyUPIPHlf4BAQE1/KpggTfmMHhsMfL8WX++Ghb4t9czeGwR4LcOwQS/dBKDxxYBfuJSTPD7e1AJ/+C3e6X+njr4TkcxwQvP7miDv7XW4Mhx3KuBa+9YrKMNPl94co0J3nCKMvgbrTnnQeNnjBvsyrWx+IgVbfC7+2ODH/k+ZfDDHPYULR2pFS1eSRv89PnY4BdOowy+cfyzxTQX8Ura4MMPYYP/tBdl8PWnPlvMqC9eSRu88C8eF/w5H8rgI+sWfwzkvGuUeCVl8B/3wQfPB46gC/4n54rtp6RnLpgaUtn5Z/FKyuCTUzHCD+9FFzwoSGxoOv22YZLlqUOUwXf/DCP8/CDK4IXcvpB70fJFPKAO3vRRClzwnzSnD95q6IL374UT/lx9auE3Zlj8ZqO9I1EzzefghOdf+xx3AaTCuz3vtjvKlHaz5RqRGqn7H6zwnmtxF0CGR/ydS6akr5drRGqkOo8VvtNE3AXQ6f/4H+rhhe/XE3cFdAq/2h0v/Ki2uCugP/h2AUJqvogZ3q/Ul8SqRm/wRqHu59o6YIYfjP/1HAJ8jnnEK2mCz+qLG37OUtx1QIGvXta3JtAEP30+bvhdb+KuAwr8vXTOY3txxCtpgu/8OW7448G464D2P95g/SKEFMHne/PY4XvcxF0IJPhpmoDfPgg//Expl26WMUjw1y9bXUURfPwy/PAfz8VdCB2+nDOcwA0/OGOH4cCBw3dxFkJ/8HnteNzwbV+Prhsd7ZeNsxD6g183Gj/8Zj7sEB+/D2ch9Ac/cgMJ8BOWMXhVY+R9viQBfuMoBq9qjMc68CTA5wYxeFVjXDiZCHjeO5/Bqxlj3ywy4PvuYfBqxtjmHBnws1MYvJpp1+u5CFb4rCgGr2bc5hMCn+/D4NVMzaOEwPPtYxm8erlfw0wEL/yoKAavXg42IgZ+hYHBq5dJXsTAH2vC4NVLcAgx8HwDBq9afulnJAfeA+unzvQF/956guA7T8JZCn3B9ykgCP7NHjhLoSv4Bx0AQfBxaNfrljm6gj+QSBJ8vM8tjLXQFXz8CaLgB+/HWAv9wF+9dMn7h0tBBMGnzMJYDv3ANwozNAgLe5kc+JiZfgeE3MBTDv3A+/GTF5YUwQwf3NspLj6+M6arDOsJPvAEUfCGzEHC7mZuxVMOHcEfN/CEwadPY/DKx29+AmnwB7syeOXj1+0/pMHzXgxe+fj6iEXww7+ezeAVT/MY8uBnpTB4xVNnO3nweyMYvNJ5WvUcefDnvBi80vmiroUIfng+7BCDVzjTmpMIPy2dwSucIF8S4Xf3Z/DK5ptoSxEC4PO9Gbyymf0pkfB854kMXtEE3ycTPiGSwSuZHwYDMuF3+TJ4JZOaRSh8fkMGr2QM9wiF51uswFMRfcB/PxiQCt99LJ6S6AN+zh5i4cdY/2JoRaMP+KD7xMLPbIanJLqA50cAcuGDL2CpiR7ge7t4hYaGkgr/9mosRdEDvJ/n2VJFyIBfEYmlKHqAbzGodBEy4LcGPsZRFHT4K2dLu1oeyfB1thANPz4XR1EQ4ON2Cc2+JhxXwfi1xUqC4R9WtSJCCPzHKTiqggDPxQKQU9ExbtmUulUvi1cSDP/xayTDD/fvUzcqKuqY2lVBhO9U74qwdMOpv3glwfD9PUmG7zH7uN/+47Mz1a4KIrzjVNNiYiPxSnLhb4ZaEyEDPp2PW8ankw5fp/D/0aKq4pXkwq9+j3T4nf0Jh48sACP8HwpLT0O8xCvJhTfeJB3+XFvC4TnOwYWbCEBuOLdBvJJY+G+LT8EgF57vvZdo+IIj6xL6eg4FYECVhRYriYWffJB8+EWTiYZ/Fr6USySSCv8o8An58CcNdMCXFlLhs2YD8uH59nMogd+YYfGbjfaORJm88QsN8BOGUALv9rzb4QRTemC/PHKp+akXoAF+jy8l8GaP+F8OmDJ5lVwjkjEP3J2aBwQEkA/PN1K9fJr+H38n2DO/DBGC4IMnq10bVPgHv90r9fdkwnvEliVCEPzoMLVrgwJ/a63BkeO4VwPXWh6RJxPe8TAl8POaqF0bBPgbrTnnQeNnjBvsyrX5U7ySSPgva5cpQhB8uv/3KtcGAX6Yw56ipSO1osUriYSPaUMN/Oh5KtcGAb5x/LPFNBfxSsLgH65bs2ZNRpMO1MAv7ahyhRDg6099tphRX7ySMPgC70WLFvXqQw98ZuQv6lYIAT6y7pdFS+ddo8QrSYPvzvP5ngepgZ85eGQ/4W/URfUqhAD/k3PF9lPSMxdMDans/LN4JYHwi8ccpQZ+UI+prsnJg1V8+xPl5VxBYkMOpmFSgcU6AuH9j1IEn8R3PMRnEgov5PaF3IulnVZPIPz7/Xiq4OdOJRneasiD7/wpXfAn/Rm8/SnonhXO0wXPh37G4O1OQffXd9AGvyiewdudAiO8bDRd8HneDN7uFDi/Tx083yuJwdubXNPhGcrg14UxeHvzZgCF8PkuqepVSJvwlzt1pxCeD7c46KlctAk/Yi+V8Mlt1SuRJuG/61VAJXym603VaqRJ+Mh8SuEjVqpWIy3CnxwKKIWf1UG1ImkRPvRHWuHnjjqtVpE0B/+Tk1MtDw83KuFXzTw85MaNGw/UqJPm4L/p1+Ykz++gEn5807BXO4X5J6pRJ+3Be87iaYWPjedTpvCbp6lRJ83BH66dTzN8Xpt8Bi8p3brwNMPzo1YyeCnJjhxAN/zBYAYvIX8HnqAcnu+RxuDRM3nLN7TD7whm8Gj53+iYPo1iomiH59uocnkiDcGf6LXZfdUHU6iHn+qnRrW0BP/WSEFpBfXwmxtdVaFaWoLvCRk0AB81RYVqaQh+f52j2oCf1kGFh7yG4EO78RqB/2SC8tXSDnzmkLe0Ag86W3waWfZoBj4v7Jh24I8qf9alVuD/CPjthHbgQZ98pQumEfiH4ceBVuDTa/n4tKzeNl3ZimkCPs/g5BEa6q0R+GQ4jyEz3lG2ZpqA39churhgWoE/0XyUsjWjHP7rGJi2zfO1Bs/H+ShbOcrht8VmZ2enNRj/vGBagV/9YpuAgADvkUpVjnZ4oWCbfWPMCqYV+BU9/IXhHRqkVOXoh9/qfSxWi/ADJk9g8FazLWmL9395bcLnt9/N4K1l2xCfo+KCaQWe/8zrNIO3knjXk5YF0wo8v7gvg7fILk/hSW+jV6aVVjCtwPNvTmbw4mSm819FDR9YesG0Av+VT7hS9aMY/vPAFGsF0wo8v6v6FniNr/Py149e+NGe260XTCvwO1xeeyc+Pra7/PWjFf5h19bHyyqYVuC7b/LP5U8zeJibeXl52727lVMwrcDzq4JOM3hT0sOGtXWKaqIXeH5F4BEGDzNvbJuUczYUTCvw/Fqv6j4wcXJWkSr4VBcPDw/Xyh2P21YwrcDz71XLgjdGOWtJySVGv9kB02czv7dHmI/NBdMK/FF/n4U44TFeYnTMmGQhTindOu9AKZhW4Dvk9e2fhw0e5yVGxwgFO7vYIfhDxIJpBZ7nF7b9Ny54PJcYvXcDZuRHn4zyjG0poWBagec/79rgmox1Jf4So0NCwsLCOrzYqntmvrSCaQWe5z2D5v8tW12Jv8TooIPbx/qG1dtuR8G0Ah/8x3Lf5Msy2ZN9idEn/CrnloMy8+wrmFbg64eFdWntWCdPltISeYnREfDtitb1OjZv7N+/tf0F0wp84aa9Y4Le/dr+GhN2idGLeXl5Jz9sNjmynU/PV0alLlq0yJHBl5yHx5qVY4Ncu2z/zb5KE3CJ0W6hMAP/vnzio+UN3OrXa9Cs8pI9X8leMK3AtxAeDYsWNsmIDIpM/ejCI6lVx/WW7V83frx4OqpHqDGoXf3I8EBv7yqtAsL7xbgcyFeqYFqBL2xrCo+VDl51mjeo7+TeYvCYKe8uXnhK+GuZ99VTWN0/8kx5WAaANPiNGc8Wz64xZXgTY0dDaGiPKJgRMTE9a8HUbGr6b93QtfFrrzVwbubaooVLoxau7u5NXFr5+AeHVvUNCO7Y1Sl2UmJycqMAmBfgG3TJL5naKtNgWyMetnUK37mLhq3bINi26gvbdj1g2zEUtj0NsB3oC9tRrWE7ualpQy6mtpmpfdlsBy8nwdZhMmxrjodt/dGwbTwUtu79YevVC7aGcNiGG2HbNwC2w7xgO64FbJMKd+BmaluY7UD5eURHRw/r+48mjRrUcaxYqVKlypUqNHRt1qJFnRdeFFIpOCIi8s3wgGAhHWauKZH3pMG7Pe9WBJ9+Bv7wEL7Z8uMlISdSYaZPmZu6aEnatCVLlqxcPhveb2USbFcnmjrNMLWm36xJtFxeDduZmbCdtQK27y6Fbepi2M5bANsF82C7OBW2S+fAdkUybDPNd2O+afMdFLWm3SSvhO3s5bBNyYDt3IWmec2H7aI02GakwLZoHqZCrrY6dlzzSFywMDV15ugp0+Li4kaMGRMTEzNwCPx44eg1pqMd/1pTHLsf8Sx0Rq7/8SyURa7DsiyURa7DsiyURa7DsiyURa7DsiyURa7DsiyURa7DsiyURa7DsiyURa7DsiyURa7DsiyURa7DsiyUhb1lq9MweJ1GWfjv09ZIz4zV0vsWHs6UmDkZdnSebkffwsOyEvMu2sk4ysJv6mfHTFovk953fjs7dvz6BDs6N7Oj77gIOzr7/oFEoyx81hI7Ovew41lkgT3fFjXrsB2d7fmc0755dnSOuI50dwZvGQZvdxg8Uhi8KQweKQwehsGXFwZvGQZvdy6csqPzB4+l972bZceOj/5iR+etdvS9fMKOzjsfIN2dvXOn0zB4nYbB6zQMXqdh8DqNIvA/92vgMuLGsx//nu7rEHBIWt8ToY71evG27/p+tRyr20Lqm9ulnmPHgxJ3DPddA+G1WYnOSOUS9bW9XErA/9aoTsI/q3reLvrxccgLI+e5/8O2V3aivl9UdkuZWa/qJZv3nc49L4NoW0h9T1VsPCe1ObdN0o6FPO3EIcCbd0Yql6gvQrmUgE+qmAfAHm5t0Y8fclsAuOYYIaVvr2q/A3C+QqxtO74+I5gzK4NoW0h9OztcBeBu87qSdixkcSWb4b9slIgAAAMxSURBVEWdkcol6otQLiXgG7aHrbOh6MdubrDds05K3xYBsH0tzLYdXzYaPczKINoWUt9qA2E7i/tJSmcAvn5xus3wos5I5RL1RSiXAvC3uSnwZqBj4Y9PXxgnuS8Y4PgXAL9UmmzzBrY/L4N4Wyh97ybuhzfDuN8ldAbgQduel1H+1Jt1RiqXeMcI5VIA/gJnuub9OK7wPcQ/uLQ0Twf/3VL6gm+dfNYsb9rS9vdQzcog3hZK38KcrhYiZccAJNS+8qNEeKRyiXeMUC4F4HO5THgzgyv8Qq5zXD3nxDQfbo2EvuDBWHgif+pTm3duVgbxtlD6wjzdUKXW91J2DI5VzAJS4ZHKJd4xQrkUgL/ILYA347i7ph/PcE43heeqXjVsOeYi6gter/rvG1eXVRpr887NyiDeFkpfIV/7c2G2/YcXd77lMgJIhkcql6gvSrkUgL/DmT5dOfjVwh+vcKZLaqRx36H3zeXS4M07FWz+6I5ZGUTbQuoLQOaLTT+2taOo8/yKaStWzObeWGHz3wuzzkjlEvVFKZcSz+qdTMekXQMLf3ry8hh4k8LZ9Hm7kn2zuQ/gzULO5rdwzPFKbgut7/YKUbds7SfuPIsrymYJndHKVbIvSrmUgE+sfB6AI89eP79VW3hqfK+Fs4S+VytGCu3j9g42H5o3xxONA6XvExcvW3uVsmMYqX/q0cpVsi9KuZSAL3B2XTqvdps7AGx02wjA+ZpOs+Z6VNwnpe9szrgwzYvbYPO+i8pg6vx8W8h9z3CBcabY/J0v5juGkQKPXi5RX4RyKfNefWT9xsOFpyggg4NfiHehj1OtsC+k9d0Z9GqdzghvmReVobDzs20h991V/Nf6V0k7BtLgJZRL1Nf2crGjczoNg9dpGLxOw+B1Ggav0zB4nYbB6zQMXqdh8DoNg9dpGLxOw+B1Ggav0zB4nYbB6zQMXqdh8DoNg9dpGLxOw+B1Ggav0zB4nYbB6zQMXqdh8DoNg9dpGLxOw+B1Ggav0zB4nYbB6zQMXqdh8DoNg9dpGLxOw+B1Ggav0zB4nYbB6zQMXqdh8DoNg9dpGLxOw+B1Ggav0zB4neb/AQLFfCMWn0AaAAAAAElFTkSuQmCC" alt="plot of chunk unnamed-chunk-1"/></p>
<pre><code>## [1] 0.004809063
</code></pre>
<pre><code class="r">if (doastro) {
## Example 4.2: Galactic Mass Estimation ----
set.seed(563478)
plotpath <- file.path(globalpath,"astro")
if (!dir.exists(plotpath)) dir.create(plotpath)
# the TMB template is part of the package. move it to a temp dir
# for compiling since this generates a bunch of new files
file.copy(system.file('extsrc/01_astro.cpp',package='aghq'),globalpath)
# Compile TMB template-- only need to do once
compile(file.path(globalpath,"01_astro.cpp"))
data("gcdatalist",package = 'aghq')
dyn.load(dynlib(file.path(globalpath,"01_astro")))
# Function and its derivatives
ff <- MakeADFun(data = gcdatalist,
parameters = list(theta1 = 0,
theta2 = 0,
theta3 = 0,
theta4 = 0
),
DLL = "01_astro",
ADreport = FALSE,
silent = TRUE)
# Nonlinear constraints and their jacobian
Es <- MakeADFun(data = gcdatalist,
parameters = list(theta1 = 0,
theta2 = 0,
theta3 = 0,
theta4 = 0
),
DLL = "01_astro",
ADreport = TRUE,
silent = TRUE)
## Parameter transformations ##
parambounds <- list(
Psi0 = c(1,200),
gamma = c(.3,.7),
alpha = c(3.0,3.7),
beta = c(-.5,1)
)
get_psi0 <- function(theta) {
# theta = -log( (Psi0 - 1) / (200 - 1) )
(parambounds$Psi0[2] - parambounds$Psi0[1]) *
exp(-exp(theta)) + parambounds$Psi0[1]
}
get_theta1 <- function(Psi0) log(
-log(
(Psi0 - parambounds$Psi0[1]) / (parambounds$Psi0[2] - parambounds$Psi0[1])
)
)
get_gamma <- function(theta) {
# theta = -log( (gamma - .3) / (.7 - .3) )
(parambounds$gamma[2] - parambounds$gamma[1]) *
exp(-exp(theta)) + parambounds$gamma[1]
}
get_theta2 <- function(gamma) log(
-log(
(gamma - parambounds$gamma[1]) / (parambounds$gamma[2] - parambounds$gamma[1])
)
)
get_alpha <- function(theta) {
# theta = log(alpha - 3)
exp(theta) + parambounds$alpha[1]
}
get_theta3 <- function(alpha) log(alpha - parambounds$alpha[1])
get_beta <- function(theta) {
# theta = -log( (beta - (-.5)) / (1 - (-.5)) )
(parambounds$beta[2] - parambounds$beta[1]) *
exp(-exp(theta)) + parambounds$beta[1]
}
get_theta4 <- function(beta) log(
-log(
(beta - parambounds$beta[1]) / (parambounds$beta[2] - parambounds$beta[1])
)
)
## Optimization using IPOPT ##
# The template returns the NEGATIVE log posterior
# So leave these as negatives. IPOPT will minimize.
ipopt_objective <- function(theta) ff$fn(theta)
ipopt_objective_gradient <- function(theta) ff$gr(theta)
ipopt_objective_hessian <- function(theta) {
H <- forceSymmetric(ff$he(theta))
H <- as(H,"dsTMatrix")
H
}
ipopt_objective_hessian_x <- function(theta,obj_factor,hessian_lambda)
obj_factor * ipopt_objective_hessian(theta)@x
ipopt_objective_hessian_structure <- function(theta) {
H <- ipopt_objective_hessian(theta)
H <- as(forceSymmetric(H),'dsTMatrix')
forStruct = cbind(H@i+1, H@j+1)
tapply(forStruct[,1], forStruct[,2], c)
}
# Box constraints, to improve stability of optimization
lowerbounds <- c(
get_theta1(parambounds$Psi0[2] - .001),
get_theta2(parambounds$gamma[2] - .001),
get_theta3(parambounds$alpha[1] + .001),
get_theta4(parambounds$beta[2] - .001)
)
upperbounds <- c(
get_theta1(parambounds$Psi0[1] + 1),
get_theta2(parambounds$gamma[1] + .01),
get_theta3(parambounds$alpha[2] - .01),
get_theta4(parambounds$beta[1] + .01)
)
thetastart <- (lowerbounds + upperbounds)/2 # Start in the middle
# Nonlinear constraints, specified as a function
ipopt_nonlinear_constraints <- function(theta) Es$fn(theta)
ipopt_nonlinear_constraints_jacobian <- function(theta) {
J <- Es$gr(theta)
as(J,"dgTMatrix")
}
ipopt_nonlinear_constraints_jacobian_x <- function(theta)
ipopt_nonlinear_constraints_jacobian(theta)@x
ipopt_nonlinear_constraints_jacobian_structure <- function(theta) {
J <- ipopt_nonlinear_constraints_jacobian(theta)
J <- as(J,'dgTMatrix')
forStruct = cbind(J@i+1, J@j+1)
tapply(forStruct[,2], forStruct[,1], c)
}
nonlinear_lowerbound <- rep(0,nrow(gcdatalist$y)+2)
nonlinear_upperbound <- rep(Inf,nrow(gcdatalist$y)+2)
stopifnot(all(ipopt_nonlinear_constraints(thetastart) > 0))
tm <- Sys.time()
ipopt_result <- ipoptr::ipoptr(
x0 = thetastart,
eval_f = ipopt_objective,
eval_grad_f = ipopt_objective_gradient,
eval_h = ipopt_objective_hessian_x,
eval_h_structure = ipopt_objective_hessian_structure(thetastart),
eval_g = ipopt_nonlinear_constraints,
eval_jac_g = ipopt_nonlinear_constraints_jacobian_x,
eval_jac_g_structure = ipopt_nonlinear_constraints_jacobian_structure(thetastart),
lb = lowerbounds,
ub = upperbounds,
constraint_lb = nonlinear_lowerbound,
constraint_ub = nonlinear_upperbound,
opts = list(obj_scaling_factor = 1,
tol = 1e-03)
)
optruntime <- difftime(Sys.time(),tm,units = 'secs')
cat('Run time for mass model optimization:',optruntime,'seconds.\n')
## AGHQ ----
# Create the optimization template
useropt <- list(
ff = list(
fn = function(theta) -1*ff$fn(theta),
gr = function(theta) -1*ff$gr(theta),
he = function(theta) -1*ff$he(theta)
),
mode = ipopt_result$solution,
hessian = ff$he(ipopt_result$solution)
)
# Do the quadrature
tm <- Sys.time()
astroquad <- aghq::aghq(ff,5,thetastart,optresults = useropt,control = default_control(negate=TRUE))
quadruntime <- difftime(Sys.time(),tm,units = 'secs')
cat("Run time for mass model quadrature:",quadruntime,"seconds.\n")
# Total runtime
aghqtime <- optruntime + quadruntime
## MCMC ----
tm <- Sys.time()
stanmod <- tmbstan(
ff,
chains = 4,
cores = 4,
iter = 1e04,
warmup = 1e03,
init = thetastart,
seed = 48645,
algorithm = "NUTS"
)
stantime <- difftime(Sys.time(),tm,units = 'secs')
# Save the traceplot
# pdf(file = file.path(plotpath,"stan-trace.pdf"),width = 7,height = 7)
# traceplot(stanmod)
# dev.off()
## TMB ----
tm <- Sys.time()
tmbsd <- TMB::sdreport(ff)
tmbtime <- difftime(Sys.time(),tm,units = "secs")
tmbsddat <- data.frame(var = paste0('theta',1:4),est = tmbsd$par.fixed,sd = sqrt(diag(tmbsd$cov.fixed)))
rownames(tmbsddat) <- NULL
# Times
# AGHQ
as.numeric(aghqtime) * stanmod@sim$iter / as.numeric(stantime)
# TMB
as.numeric(optruntime) * stanmod@sim$iter / as.numeric(stantime)
# Redefine parameters functions for plotting
get_psi0 <- function(theta)
(parambounds$Psi0[2] - parambounds$Psi0[1]) *
exp(-exp(theta)) + parambounds$Psi0[1]
get_theta1 <- function(Psi0)
log(-log( (Psi0 - parambounds$Psi0[1]) /
(parambounds$Psi0[2] - parambounds$Psi0[1]) ))
get_gamma <- function(theta)
(parambounds$gamma[2] - parambounds$gamma[1]) *
exp(-exp(theta)) + parambounds$gamma[1]
# Add a little buffer, for stability
get_theta2 <- function(gamma)
log(-log( (gamma - parambounds$gamma[1] + 1e-03) /
(parambounds$gamma[2] - parambounds$gamma[1] + 1e-03) ))
get_alpha <- function(theta)
exp(theta) + parambounds$alpha[1]
# Add a little buffer, for stability
get_theta3 <- function(alpha)
log(alpha - parambounds$alpha[1] + 1e-03)