Newer
Older
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
{
"cells": [
{
"cell_type": "code",
"execution_count": 65,
"id": "ef2a766b-a3a1-431f-8dac-6c57c1157312",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.stats as stt\n",
"import statsmodels.api as sm\n",
"import statsmodels.formula.api as smf\n",
"from patsy import dmatrices\n",
"import pandas as pd\n",
"\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"import seaborn as sb\n",
"from statsmodels.graphics.factorplots import interaction_plot\n",
"\n",
"sb.set_style('whitegrid')\n",
"sb.set(font_scale=1.5)"
]
},
{
"cell_type": "markdown",
"id": "102b1593-da1a-485c-a8a3-0129619a69ac",
"metadata": {},
"source": [
"## Synthetic example\n",
"\n",
"We start with an example where we control the dependency from 2 predictor variables ($x1$ and $x2$) to a response variable ($y$).\n",
"\n",
"For further reference on the models and , see the docs about [ordinary least square](https://www.statsmodels.org/dev/examples/notebooks/generated/ols.html) estimation, as well as [patsy](https://patsy.readthedocs.io/en/latest/).\n",
"\n",
"### Exercise\n",
"- Generate $y$ such that it has a linear dependency with each of the variable, test a linear model on it.\n",
"- Add a non-linear dependence to $y$ with the product "
]
},
{
"cell_type": "code",
"execution_count": 57,
"id": "3576be4a-1c48-4f7d-abbc-f18064c930e6",
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>x1</th>\n",
" <th>x2</th>\n",
" <th>y</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1.951076</td>\n",
" <td>1.654086</td>\n",
" <td>0.061812</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>-0.024425</td>\n",
" <td>1.342077</td>\n",
" <td>0.091460</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>0.638307</td>\n",
" <td>-0.854780</td>\n",
" <td>1.663423</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>0.958588</td>\n",
" <td>1.135066</td>\n",
" <td>0.393107</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>-0.944212</td>\n",
" <td>-0.223619</td>\n",
" <td>-0.959402</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" x1 x2 y\n",
"0 1.951076 1.654086 0.061812\n",
"1 -0.024425 1.342077 0.091460\n",
"2 0.638307 -0.854780 1.663423\n",
"3 0.958588 1.135066 0.393107\n",
"4 -0.944212 -0.223619 -0.959402"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# number of samples\n",
"n = 50\n",
"\n",
"# exogeneous variables\n",
"x1 = stt.norm.rvs(size=n)\n",
"x2 = stt.norm.rvs(size=n)\n",
"\n",
"# error\n",
"e = stt.norm.rvs(size=n)\n",
"\n",
"# endogenous variable scaled by factor a\n",
"a1 = -0.5\n",
"a2 = 0.7\n",
"y = a1 * x1 + a2 * x2 + e\n",
"\n",
"# build dataframe\n",
"df = pd.DataFrame(np.column_stack((x1,x2,y)), columns=['x1','x2','y'])\n",
"df.head()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "6345745a-3046-48fe-b3b8-54e09c507dc9",
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAHMCAYAAAAarpbgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABGgklEQVR4nO3de3RU1f3//9dMboQwIQlXiRFQQ2IIfCopVy9YpYK2VVCLgFX5VMUC0t8SRWsV+6k31FqpUoVaq3IRxQvIEixWLloVQRovEJAASiAYIGISCAm5zZzfH3yTGpJMZiaTOefMPB9rdVVmdmbezk7klX32eW+HYRiGAAAA4DOn2QUAAADYDQEKAADATwQoAAAAPxGgAAAA/ESAAgAA8BMBCgAAwE8EKAAAAD8RoAAAAPxEgAIAAPATAQoAAMBP0WYXYJb169frww8/1Pbt23Xo0CGVlpYqOjpaqampGj58uCZPnqzU1FSzywQAABbkiNSz8K6//np9+umniomJUbdu3ZSSkqLS0lIdPHhQHo9H8fHx+utf/6rzzz/f7FIBAIDFRGyAeuutt9SjRw/l5OQoNja24fH9+/fr97//vbZs2aLk5GStX79eHTt2NLFSAABgNREboLw5cuSIzjvvPEnSc889p5EjR5pcEQAAsJKI3QPlTdeuXZWUlKSysjJVVVUF/DqGYcjjIZ/6yul08HmZjDkwH3NgPubAfGbOgdPpkMPhaHUcAaoZX3/9tcrKyuR0OpWVlRXw63g8hkpKKoJYWfiKjnYqOTlBx45Vqq7OY3Y5EYk5MB9zYD7mwHxmz0FKSoKioghQPjMMQyUlJcrNzdUTTzwhSfr1r3+ttLQ0kysDAABWE/EBauXKlbrrrrsaPXbmmWfqiSee0C9+8Ys2v350NK22fBEV5Wz0/wg95sB8zIH5mAPz2WUOIj5AdenSRYMGDZJhGDp06JAOHz6sgoICvf322xo8eLB69uwZ8Gs7nQ4lJycEsdrwl5gYb3YJEY85MB9zYD7mwHxWnwPuwjtFYWGhHn30Ua1du1Y9evTQ6tWr5XK5Anott9ujY8dOBLnC8BQV5VRiYryOHTsht5t9B2ZgDszHHJiPOTCf2XOQmBjv0+pXxK9AnSotLU1PP/20rrzySu3evVtLlizR1KlTA349NiH6x+328JmZjDkwH3NgPubAfFafA2tfYDRJVFSULrjgAklSXl6eydUAAACrIUC1oK6uTpLk8Vg3/QIAAHMQoJpRU1Oj999/X5La1AcKAACEp4gMUNu2bdNf/vIXFRQUNHlu7969mjp1qvbv36+OHTtq/PjxoS8QAABYWkRuIq+srNT8+fM1f/58paSk6LTTTlN0dLS+++47FRUVSZKSkpL0l7/8RT169DC5WsA8Ho+hXYVlKquoVlJCnPqlJcnpbL1DLwCEu4gMUJmZmbrvvvv06aefateuXdq3b5+qqqrUqVMn5eTk6IILLtC1116rlJQUs0sFTJObX6yla3ertLy64bFkV5wmjUpXTkZ3EysDAPPRB6odud0ezsLzUf3ZR6WlFZa+bTWc/XAONm8/pGdWtHwH6vRx2YSodsDPgfmYA/OZPQcnz8JrfYdTRO6BAtAyj8fQ0rW7vY55Ze1uTqsHENEIUAAayd9f2uiyXXNKyqu1q7AsNAUBgAURoAA0Una8xrdxFd5DFgCEMwIUgEaSOsX6Ni4hrp0rAQDrIkABaCTjjGQlu7yHoxTXyZYGABCpCFAAGnE6HZo0Kt3rmImj0ukHBSCiEaAANJGT0V3Tx2U3WYlKccXRwgAAFKGNNAG0Lieju85N70YncgBoBgEKQIucTocyeyebXQYAWA6X8AAAAPxEgAIAAPATAQoAAMBPBCgAAAA/EaAAAAD8RIACAADwEwEKAADATwQoAAAAP9FIE5bk8Rh0wAYAWBYBCpaTm1+spWt3q7S8uuGxZFecJo1K5ww2AIAlcAkPlrJlZ7GeWZHXKDxJUml5tZ5Zkafc/GKTKgMA4L8IULAMt8fQy+/mex3zytrd8niMEFUEAEDzCFCwjB3ffK+SU1aeTlVSXq1dhWWhKQgAgBYQoGAZJceqfBpXVuE9ZAEA0N4IULCMlMQOPo1LSohr50oAAPCOAAXLyDqzi1Jc3sNRiutkSwMAAMxEgIJlRDkdum50htcxE0el0w8KAGA6AhQsZXBmd00fl63kU1aiUlxxmj4umz5QAABLoJEmLCcno7vOTe9GJ3IAgGURoGBJTqdDmb2TzS4DAIBmEaAAcfYeAMA/BChEPM7eAwD4i03kiGi5+Zy9BwDwHwEKEcvjMbR07W6vYzh7DwDQHAIUItauwrImK0+n4uw9AEBzCFCIWL6eqcfZewCAUxGgELF8PVOPs/cAAKciQCFi9UtLatLx/FScvQcAaA4BChHL6XRo0qh0r2M4ew8A0BwCFCJaTgZn7wEA/EcjTUQ8zt4DAPiLAAWIs/cAAP7hEh4AAICfCFAAAAB+IkABAAD4iQAFAADgJwIUAACAnwhQAAAAfiJAAQAA+IkABQAA4CcCFAAAgJ8irhO5YRj6/PPPtX79euXm5uqbb77R8ePH5XK5lJWVpbFjx+oXv/iFHA6O8QAAAM2LuAC1adMmTZ48ueHPaWlpSk1N1bfffquPP/5YH3/8sVavXq158+YpNjbWvEIBAIBlRVyAMgxDp59+um688Ub97Gc/U5cuXRqee+uttzR79my9//77evrpp3XnnXeaWCkAALCqiNsDNXDgQK1Zs0Y33HBDo/AkSWPHjtX06dMlSa+//ro8Ho8ZJQIAAIuLuADVqVMnxcTEtPj8hRdeKEkqKytTSUlJqMoCAAA2EnEBqjXV1dUN/9yhQwcTKwEAAFZFgDrF6tWrJUmZmZnq1KmTydUAAAArirhN5N5s375dr776qiRpypQpQXnN6Ggyqi+iopyN/h+hxxyYjzkwH3NgPrvMgcMwDMPsIqzgyJEj+uUvf6mioiL99Kc/1V//+tc2v6ZhGPSTAgAgDBGgJJWXl+uGG27Qjh071L9/fy1atCgol+/cbo+OHTsRhArDX1SUU4mJ8Tp27ITcbu5+NANzYD7mwHzMgfnMnoPExHifVr8i/hJeRUWFbr75Zu3YsUPp6en6xz/+EdS9T3V1/AD6w+328JmZjDkwH3NgPubAfFafA2tfYGxnJ06c0K233qovvvhCffr00Ysvvqjk5GSzywIAABYXsQGqurpa06ZN05YtW5SamqqFCxeqW7duZpcFAABsICIDVG1trWbMmKGNGzeqZ8+eWrhwoXr27Gl2WQAAwCYiLkC53W7deeed+uCDD9StWzctXLhQaWlpZpcFAABsJOI2kf/zn//UmjVrJEmxsbG65557Whw7e/ZsZWVlhao0AABgExEXoGpqahr++dtvv9W3337b4tjy8vJQlAQAAGwm4gLUVVddpauuusrsMgAAgI1F3B4oAACAtiJAAQAA+IkABQAA4CcCFAAAgJ8ibhM5AOvyeAztKixTWUW1khLi1C8tSU6nw+yyAKAJAhQAS8jNL9bStbtVWl7d8FiyK06TRqUrJ6O7iZUBQFNcwgNgutz8Yj2zIq9ReJKk0vJqPbMiT7n5xSZVBgDNI0ABMJXHY2jp2t1ex7yydrc8HiNEFQFA6whQAEy1q7CsycrTqUrKq7WrsCw0BQGADwhQAExVVuE9PPk7DgBCgQAFwFRJCXFBHQcAoUCAAmCqfmlJSnZ5D0cprpMtDQDAKghQAEzldDo0aVS61zETR6XTDwqApRCgAJguJ6O7po/LbrISleKK0/Rx2fSBAmA5NNIEYAk5Gd11bno3OpEDsAUCFADLcDodyuydbHYZANAqAhTaBWeaAQDCGQEKQceZZgCAcMcmcgQVZ5oBACIBAQpBw5lmAIBIQYBC0HCmGdB+PB5DO/eVatOOQ9q5r5RfRACTsQcKQcOZZkD7YF8hYD2sQCFoONMMCD72FQLWRIBC0HCmGRBc7CsErIsAhaDhTDMguNhXCFgXAQpBxZlmQPCwrxCwLjaRI+g40wwIDvYVAtZFgEK74EwzoO3q9xV6u4zHvkLAHFzCAwCLYl8hYF0EKACwMPYVAtbEJTwAsDj2FQLWQ4ACABtgXyFgLVzCAwAA8BMBCgAAwE8EKAAAAD+xBwpAxPJ4DDZmAwgIAQpARMrNL9bStbsbNalMdsVp0qh0WgMAaBWX8ABEnNz8Yj2zIq9Jh+/S8mo9syJPufnFJlUGwC4IUAAiisdjaOna3V7HvLJ2tzweI0QVAbAjAhSAiLKrsMzr2XKSVFJerV2FZaEpCIAtEaAARJSyCu/hyd9xACITAQpARElKiGt9kB/jAEQmAhSAiNIvLanJwbynSnGdbGkAAC0hQAGIKE6nQ5NGpXsdM3FUOv2gAHhFgAIQcXIyumv6uOwmK1EprjhNH5dNHygAraKRJoCIlJPRXeemd6MTOYCAEKAARCyn06HM3slmlwHAhriEBwAA4CcCFAAAgJ8IUAAAAH6KyD1Q3333nTZu3Kht27YpLy9PX331laqqqtS/f38tX77c7PIAAIDFRWSAWr16tebMmWN2GQAAwKYiMkB16tRJI0aMUHZ2trKzs1VQUKAnn3zS7LIAAIBNBBygtmzZosGDBwezlpC55pprdM011zT8mct2AADAHwFvIr/++ut1+eWX66WXXlJZWVkQSwIAALC2Nt2F98033+ixxx7ThRdeqFmzZmnLli3BqgsAAMCyAg5Q//rXv3TzzTerS5cuqqmp0apVq3TDDTewKtWOPB5DO/eVatOOQ9q5r1Qej2F2SQAARCSHYRht+lu4rq5O69ev17Jly/TJJ5/I4/HI4XAoNjZWl156qcaPH2/5vVLLly/XPffcE/Q2Bm63R8eOnQjKa23ZWayX381XSXl1w2MprjhdNzpDgzPtf/BpVJRTiYnxOnbshNxuj9nlRCTmwHzMgfmYA/OZPQeJifGKimp9fanNAeqHioqK9Nprr2nFihU6fPjwyTdwONS3b1+NHz9eY8eOVVJSUrDeLmjaK0AZhiGHo+0Hk27cWqQ5C1u+PHrPjYM1YmCvNr8PrMXtMbTjm+9VcqxKKYkdlHVmF0Vx0C0AWEJQA1Q9j8ej999/X6+//rr+/e9/y+12N1qVmjBhgnJycoL9tgGz8gqUx2No5ryPGq08nSolMU5P3na+rU+RN/s3Dn95PIby95eq7HiNkjrFKuOM5KB+/masONptDsIRc2A+5sB8Zs+BrytQ7dIHyul06uKLL1Z0dLTKysr0xRdfyDAMVVdX6+2339aqVav0P//zP/r973+vgQMHtkcJllFX17bJ37mv1Gt4kqSSY9XasbckLE6Vd7s9bf7M2ltufrGWrt2t0h/MS7IrTpNGpSsno+3hJje/WM+syGvyeEl5tea9sVXTx2UH5X1aYoc5CHfMgfmYA/NZfQ6CfhZecXGxnn32WV188cW69dZb9fnnn8swDOXk5Ojee+/VRRddJIfDoS+++EKTJk3S5s2bg11CWCmr8B6e/B2HtqkPN6WnhNrS8mo9syJPufnFbXp9j8fQ0rW7vY55Ze1ubiAAAJMFZQXKMAx98MEHeu211xou2RmGoU6dOunKK6/UhAkTlJ6eLulk/6jCwkL93//9nz7++GM99dRTWrp0aTDKCEtJCXFBHYfA+Rpuzk3vFvDlvF2FZU3C2alKyqu1q7AsLFYcAcCu2hSgDh06pDfeeENvvvmmDh06pPrtVFlZWZo4caJ+/vOfKz4+vsnXpaWl6amnntLw4cOVn5/flhLCXr+0JCW74rz+pZriilO/tKTQFRWhQhFuWHEEAHsIOEDdeuut+uijj+TxeGQYhuLj43XZZZdpwoQJPu1r6tSpk7p166aDBw8GWkJEcDodmjQqvdk9MfUmjkq39QZyuwhFuGHFEQDsIeAA9cEHH0iSzjrrLF177bUaN26cXC6XX68xevRoUxpuHjx4UGPHjm34c01NjSQpPz9fQ4cObXj85ptv1i233BLq8prIyeiu6eOym2xcTnHFaWKQNi6jdaEIN6w4AoA9BBygLr/8ck2YMEFDhgwJ+M3vvvvugL+2Ldxud7PBra6urtHjVVVVoSuqFTkZ3XVuejftKixTWUW1khJO/iXKylPohCLcsOIIAPbQLn2gcJLb7VFJSUVQX9PjMcIyREVHO5WcnKDS0gpL37baUouBesFqMdBcq4T2XnG0yxyEM+bAfMyB+cyeg5SUBPP6QKF9tHf/IbQuVJdTWXEEAGsjQNlESysf9f2H2ru5Iv4rVOHG6XTQqgAALIoAZQOh6D8E/xBuAHuoq/No/WcHVFx2Qt2T4nXxoNMVHR30HtKIQAQoG6C5IgD477X1u/XulkL9cKfvsg17NHpwmsZfnG5eYQgLBCgboLkiAPjntfW7tebTwiaPG4YaHidEoS1Yx7QBmisCLfN4DO3cV6pNOw5p575SzgmE6uo8endL0/D0Q+9uKeQuO7QJK1A2QHNFoHncmYrmrP/sgFpr0GMYJ8ddOuSM0BSFsMMKlA3UN1f0huaKiDT1d6ae+otF/Z2pufnFJlUGsxWXnQjqOKA5BCibqO8/lOxqfJkuxRVHCwNEHF/vTOVyXmTqntT0EPu2jAOaE/AlvKKiIr/Gx8XFyeVyKTY2NtC3jHg0V4xM4dp9vi24MxXeXDzodC3bsMfrZTyH4+Q4IFABB6hLLrkkoK9LS0vThRdeqOuvv169e/cO9O0jFv2HIgt7fJrHnanwJjraqdGD05q9C6/e6MFp9INCmwT83WMYRkD/279/v15++WVdeeWV+uc//xnMfxcgrLDHp2Xcmfpf3IXYvPEXp2vMkDQ5TlmsdTikMUPoA4W2C3gFat26ddq2bZv+8Ic/yOl0auLEiRoyZIh69OghwzBUXFysTz/9VK+++qrcbrcefPBBpaWladu2bVq0aJH27Nmju+++W1lZWaxEAaeg+7x3dr8zNViXZVmh9G78xem66sKz6ESOduEwjNZu9mze/v37dfXVV+v000/XCy+8oOTk5i8rlZaW6te//rWKioq0fPlypaamqqamRjfccIO+/PJLTZo0SbNnz27Tv4RVud0elZRUmF2GLZh9+rbV7NxXqsdf+bzVcXdNPDdol3TtNgctnQ9Zz6o3V3gLPUP79/R5Duz67291dvs5CEdmz0FKSoKioloP2QHH8Pnz5+v48eN68MEHWwxPkpScnKwHHnhAR48e1YIFCyRJsbGxuuOOO2QYhjZt2hRoCUDYYo9P6+x4Z2prl2W37PTtsix3IQLmC/gS3saNG9WxY0dlZ2e3OnbAgAHq2LGjPvroo4bHBg0apJiYGB08eDDQEoCwxR4f39jpzlRfQs/L/8rXJUP7tPpa3IUImC/gAFVSUqLoaN+/3DAMff/99w1/joqKUseOHVVVVRVoCUDYsvsen1Cyy52pPoWeY9Xa8c33Or2L9/5ErFAC5gv4El5KSoqqqqq0ZcuWVsdu2bJFJ06caHSpr7a2VseOHfN6+Q+IVHSfDz++hpmSY63/UhnpK5TceQgrCHgF6rzzztPy5ct133336fnnn1daWlqz4w4cOKD77rtPDodD559/fsPjBQUFMgxDvXr1CrQEIKzV7/E5dcNxiitOE7nLynZ8DTMpiR1aHRPJK5TceQirCDhA3XbbbXr33Xe1f/9+/eIXv9Bll12mIUOGqHv37nI4HCouLtbmzZu1Zs0anThxQgkJCZo2bVrD169evVqSNGTIkLb/WwBhyk57fOCdT6EnMU5ZZ3bRsaOVXl+rfoXS21144bhC2dKdh/Wb8K168wDCU8BtDCTp888/14wZM3TkyBE5Tu1W9v8YhqGuXbvqqaeeUk5OTsPj77zzjr777jv95Cc/0RlnhOdp2LQx8J3Zt62COQiF1loPzLhmoC4d3tfnOWhuNSZcVyg9HkOz5m9sddXt8akj2hQc+Tkwn9lz4Gsbg4BXoCTp3HPP1T//+U8tXrxY7777rvbs2SO32y3p5Cbxs88+W6NHj9avfvUrJSYmNvrayy+/vC1vDQC209pl2cGZ/oWeSFqh5M5DWE2bApQkuVwuTZs2TdOmTVNtba2OHj0qwzCUlJSkmJiYYNQIAGEj2KHHanchttfh19x5CKtpc4D6oZiYGHXt2jWYLwkAYcdqoSdY2nODd6TfeQjrCeqBQG63WyUlJSopKWm4lAcACH/tffh1/SZ8b8L1zkNYU5tXoE6cOKFXX31Vq1atUn5+fqM9UJmZmfr5z3+ua6+9VvHx3hvDAQDsKRSHX0fqnYewrjYFqG+++UZTp07V/v37derNfHV1dcrLy9P27dv1yiuvaMGCBerbt2+bigVgjvba1xIqdq/f6kK1wZveaLCSgAPU8ePHddNNN+ngwYOKjo7WT3/6U40YMUI9e/aUJB06dEiffPKJ/vWvf2nfvn266aab9PbbbyshISFoxQNof3ZvXGj3+u0glBu8I+nOQ1hbwAFq4cKFOnjwoLp3766//e1vOuecc5qM+eUvf6mdO3dqypQpOnjwoBYtWqSpU6e2qWAAoWP3xoV2r98uQr3BO1w34cNeAt5Evm7dOjkcDj3wwAPNhqd6mZmZevDBB2UYht57771A3w5AiPm6r8Wq55DZvX47YYM3IlHAAWrfvn2KjY3VRRdd1OrYCy+8UHFxcdq3b1+gbwcgxPzZ12JFdq/fTjj8GpEo4ABVV1fnc6NMh8OhmJgY1dXVBfp2AELM7o0LrV6/x2No575SbdpxSDv3ldp+Jax+g/epK1EprjgulSIsBbwHqmfPntq/f792796t9HTvv3ns2rVLx48fV+/evQN9OwAhZvfGhVauv6WN7b8anaFLh9v3bmU2eCOSBLwCNWzYMBmGoT/+8Y+qrm75N7jq6mr98Y9/lMPh0PDhwwN9OwAhZvd9LVat31vDyXlvbNXGrUUhrSfY6jd4D8vqqczeyYQnhK2AA9TNN9+s2NhY5ebm6oorrtDrr7+uAwcOqLa2VrW1tSosLNTrr7+uK664Qrm5uYqJidFNN90UzNoBtCO772uxYv2+bGz/+8o821/OAyKBwzi1A6Yf3nnnHd11112qq6uTw9H8f4QMw1B0dLQef/xxXX755QEXakdut0clJRVml2EL0dFOJScnqLS0QnV1HrPLiUgtzUFzl5t+2LjQ6k0qW6s/lHbuK9Xjr3ze6rh7fjVI6acntX9BaIL/FpnP7DlISUlQVFTr60tt6kR++eWXq3fv3po7d64+/vjjJt3InU6nLrjgAt1+++1eWx0AsC5v+1rs0KTSSvtyfN7YfrymnSsB0FZtPguvf//+ev7551VeXq7t27erpKREkpSSkqL+/fvL5XK1uUgA5mqucaGdmlRapfGizxvbO8W2cyUA2qrNAaqey+XSsGHDgvVyACwsFIfHhqP6je3e+lN1TYpXxhnJ7IMCLC7gTeQAIhdNKgPjy8b2W67MJnQCNkCAAuA3qzeptDJvDSdnXDNQIwb2MqkyAP7w6RJesDaAOxwO7dixIyivBcA8Vm5SaQctbWyPjY0yuzQAPvIpQLWh0wGAMOTLXh4rN9m0AqtsbAcQGJ8C1KJFi5p9/MCBA3r00UdVWVmp0aNHa9iwYerZs6ck6fDhw9q0aZPeffdddezYUb/73e+UmpoavMoBmKZ+L09zd+HVs3KTTQBoq4AbaX7//fcaO3asoqOj9dxzz7V4Ht6ePXs0ZcoUud1urVixQikpKW0q2E5opOk7sxunIbA5sFKTSjMFq5koPwfmYw7MZ/YctHsjzfnz5+vIkSP6+9//7vUw4bPPPlsPPPCAbr75Zs2fP1/33ntvoG8JwGKs1KTSLHZoJgog+AK+C+/9999XXFyczj///FbHnn/++erQoYM2bNgQ6NsBsKhIPjzW28HAz6zIU25+sUmVAWhvAQeo4uJiRUX5fsdIVFSUvvvuu0DfDgAsxddmojTEBMJTwAEqMTFRlZWVystreRNpvby8PFVUVHCsC4CwQTNRILIFHKCGDh0qwzA0e/ZslZaWtjiurKxMs2fPlsPh0NChQwN9OwCwFJqJApEt4E3k06dP13vvvaedO3fq8ssv18SJEzV06FD16NFDDodDhw4d0ubNm/Xqq6+qpKREcXFxmjZtWjBrb7NNmzbpxRdf1JdffqnKykr16tVLY8aM0ZQpU9SxY0ezywNgYTQTBSJbwG0MJOnf//63Zs6cqePHj8vhaH7jqGEYSkhI0JNPPqmRI0cGXGiwLV68WA8//LAMw1DPnj2VkpKiPXv2qKamRmeddZaWLl2qpKSkNr0HbQx8Z/Ztq2AO/OXxGJo1f2OrzUQfnzrC5431zIH5mAPzmT0HvrYxaNNZeBdeeKFWr16tCRMmKDExUYZhNPpfYmKiJkyYoFWrVlkqPOXl5emRRx6RJD3wwAN6//33tWLFCq1du1b9+/fX119/rdmzZ5tcJQAr8+VgYJqJAuGrTStQpyosLFRJSYkkKSUlRWlpacF66aCaNm2a1q1bp7Fjx+qxxx5r9FxBQYEuu+wyeTwerVy5UpmZmQG/DytQvjP7Nw4wB4EKZjNR5sB8zIH5zJ6Ddm+k2Zy0tDTLhqZ6FRUV+vDDDyVJ48ePb/J8nz59NGzYMG3cuFFr1qxpU4ACEP7CuZlosDqsA+EoqAHKDr766ivV1NQoNjZWAwcObHZMTk6ONm7cqC+//DLE1QGwo3A8GJgO64B3QQlQHo9HBQUFOnr0qOrq6ryOHTx4cDDeMmB79+6VJPXq1UsxMTHNjjnjjDMajQWASFLfYf1U9R3Wp4/LDjhEWX1Vy+MxtG3PERUePCpXfIzl6oN1tClAFRcX68knn9S7776rqqqqVsc7HA7t2LGjLW/ZZkePHpUkde7cucUx9c/Vj22L6Og27dOPGPXXm3257oz2wRyYzwpz4PEYeqW1DuvrdmvwOT38DhZbdhbr5XfzVXLKfrHrRmdocKb5q1pbdhbr5X/lq+SYNeuLFFb4OfBFwAHq8OHDGj9+vIqLi+XrPvQg7lcPWHX1yR+MllafJCk2NrbR2EA5nQ4lJye06TUiTWJivNklRDzmwHxmzsG2PUcaBZzmlByrVlFplQac3dXn1924tUjz3tja9LXKqzXvja2658bBGjGwl9/1BovV64tEVv9vUcAB6q9//asOHz6shIQE3X777brkkkvUvXt3v87HM0Nc3MmmdrW1tS2OqampaTQ2UB6PoWPHKtv0GpEiKsqpxMR4HTt2Qm43d76YgTkwnxXmoPCgbyvvhQeP6vQuvv0F5/EY+tvypuHkh/62YqsyUhNNuVxm9foijdk/B4mJ8e17F96///1vORwOPfzwwxozZkygLxNyvlye8+Uyn6+4DdY/breHz8xkzIH5zJwDV3zLq/OnjvO1xp37Sn1a1dqxt8SUzfhWry9SWf2/RQFfYCwpKVFUVJRGjRoVzHraXZ8+fSRJRUVFLa5C7d+/v9FYAIgU/dKSlOzyvvqe4jq5+dtXVj830Or1wZoCDlBdunRRhw4dFB1tr04IWVlZiomJUU1NjbZubX7JNjc3V5L0ox/9KISVAYD52qPDutXPDbR6fbCmgAPU8OHDVVFRoYKCgiCW0/4SEhJ0/vnnS5Jee+21Js8XFBRo06ZNkmSrS5MAECw5Gd01fVx2k5WoFFdcQC0M2mNVK5isXh+sKeAA9Zvf/Ebx8fF64okngllPSEybNk0Oh0MrV67UsmXLGu4OLC4u1syZM+XxeDRq1Ci6kAOIWDkZ3fWnqSN018RzNeWKLN018Vw9PnVEQP2frH5uoNXrgzW16Sy8zZs367e//a2ysrJ06623auDAgerYsWMw62s3L730kh599FEZhqHTTjtNycnJ2rNnj2pqatS3b18tXbpUKSkpbXoPzsLzndlnH4E5sIJwn4NgnhvYHnLzi/XK2t1N+lRZpb5IYfbPga9n4QUcoM455xy/v8YKjTR/6JNPPtELL7ygrVu3qrKyUr169dKYMWM0ZcoUJSS0vX8TAcp3Zv/AgDmwgkiYA6t3Inc6HSoqraITuYnM/jlo98OErdAUs62GDx+u4cOHm10GAEQMq58b6HQ6NODsrjq9S3zYhlgER8ABatGiRcGsAwAAwDYCDlBDhgwJZh0AAAC2Ya8mTgBszer7XwDAV0ELUIZhqLS0VFVVVerViwMXATTW3B1Yya44TeIOJwA21OYAtX37ds2fP18bN27UiRMnmtxpd/ToUf35z3+WJN13332KjY1t61sCsJnc/GI9syKvyeOl5dV6ZkVeQM0Z0b5YLQw9PnN7aVOAeuutt3Tfffeprq6uxTGdO3fWgQMH9Mknn+jiiy/WRRdd1Ja3BGAzHo+hpWt3ex3zytrdOje9G39ZWASrhaHHZ24/AXci//rrrzV79mzV1dXp+uuv15tvvqnk5OZvTb3yyitlGIbWrVsXcKEA7GlXYVmjvxSaU1JerV2FZaEpCF7VrxaeOmf1q4W5+cUmVRa++MztKeAA9eKLL6q2tlbXXXed7r33XvXv319RUVHNjh02bJgk6Ysvvgj07QDYFCfd24evq4Uej/37AFoFn7l9BRygNm3aJIfDoVtuuaXVsT169FB8fLyKiooCfTsANsVJ9/bBamHo8ZnbV8ABqri4WPHx8erZs6dP4+Pi4lRdzW+YQKThpHv7YLUw9PjM7SvgABUbG6va2lqfjnSpqqpSeXm5OnXqFOjbAbApTrq3D1YLQ4/P3L4CDlCpqamqq6tTQUFBq2M/+OADud1unX322YG+HQAby8norunjspusRKW44mhhYCGsFoYen7l9BdzG4IILLlB+fr4WLVqkP/zhDy2OKy0t1Z/+9Cc5HA6NHDky0LcDWkTvFHvIyeiuc9O7MVcWVr9a2FzPrnqsFgYXn7l9BRygJk+erKVLl+rVV19Vly5dNHny5EbPV1VV6b333tPcuXNVVFSk5ORkTZw4sa31Ao3QO8VenE6HMns33+4E1lC/Wnjqz1WKK04Tg/RzxS89jYXiM0fwOQxfNjG1YMOGDfrtb3+ruro6RUdHyzAMud1unXnmmSosLGzYIxUbG6u//e1vGj58eDBrtzy326OSkgqzy7CF6GinkpMTVFpaobo6j09f01J363pcGvJPIHOA4LLSHLRXyLH6Lz1mzgHB8iSzfw5SUhIUFdX6DqeA90BJ0k9+8hO9/PLL6t+/v2pra1VXVyfDMPT111+rpqZGhmEoKytLS5YsibjwhPZF7xSgfdWvFg7L6qnM3slBC080jGxZe3zmaD9tPgtv4MCBeuONN7Rz507l5uaquLhYHo9HXbt21aBBgzRgwIBg1Ak04k/vFC4ZAebjSB+EmzYHqHqZmZnKzMwM1ssBXtE7BbAXfulBuGnTJTzALPROAeyFX3oQbghQsCV6pwD2wi89CDcEKNgS3a0Be+GXHoQbAhRsKxy6W3s8hnbuK9WmHYe0c18pdw0ibPFLD8JN0DaRA2awc3drq/fDAYKNhpEIJwQo2J4du1u31AS0vh+OXVbQAH/Z+Zce4IcIUECI0Q8Hkc6Ov/QAp2IPFBBi/vTDAQBYEwEKCDH64QCA/RGggBCjHw4A2B97oIAQq++H4+0yHv1wAKB5Ho9hiZsQCFBAiNX3w2nuLrx69MMBgKas1P6FS3iACcKhCSgAhFJ9+5dTV+/r27/k5heHtB5WoACT0A8HAHxjxfYvBCjARPTDAQJXV+fR+s8OqLjshLonxeviQacrOpoLK+HIn/YvofpvKgEKAGA7r63frXe3FMr4wfGRyzbs0ejBaRp/sfcz92A/Vmz/QoACANjKa+t3a82nhU0eNww1PE6ICi9WbP/CWicAwDbq6jx6d0vT8PRD724pVF2dJ0QVIRTq2794E+r2LwQoAIBtrP/sQKPLds0xjJPjED7q2794E+r2LwQoAIBtFJedCOo42IfV2r+wBwoAYBvdk+KDOg72YqX2LwQoAIBtXDzodC3bsMfrZTyH4+Q4hCertH/hEh4AwDaio50aPTjN65jRg9PoB4V2xwoUAMBW6lsUnNoHyuEQfaAQMgQoAIDtjL84XVddeBadyGEaAhQAwJaio526dMgZZpeBCEWAAgAgBDwewxJ3jyE4CFAAALSz3PxiLV27u9GBuMmuOE0alR7y/kUIDi4WAwDQjnLzi/XMirxG4UmSSsur9cyKPOXmF5tUGdqCAAUAQDvxeAwtXbvb65hX1u6Wx9PK+TSwHAIUAADtZFdhWZOVp1OVlFdrV2FZaApC0BCgAABoJ2UV3sOTv+NgHRG3ibyqqkofffSRtm3bpry8POXl5amsrEyS9NlnnykhIcHcAgEAYSMpIa71QX6Mg3VEXIDau3evpk+fbnYZAHzErd+ws35pSUp2xXm9jJfiOvl9DXuJuAAVHR2tgQMHasCAAcrOzlaXLl00ZcoUs8sC0Axu/YbdOZ0OTRqVrmdW5LU4ZuKodH4psKGIC1Dp6el6/fXXG/584MABE6sB0JL6W79PVX/r9/Rx2YQo2EJORndNH5fd5JeBFFecJvLLgG1FXIACYH2+3vp9bno3fnOHLeRkdNe56d24HB1GCFAALMefW78zeyeHqCqgbZxOB9+vYYQ2BgAsh1u/AVgdK1DtLDqajOqLqChno/9H6FlpDrokdvB5XDj9jFlpDiIVc2A+u8wBAaodOZ0OJSfTV8ofiYnxZpcQ8awwB0M7d1SXt3fo+6NVLY7pmhSvof9zuqLCcA+JFeYg0jEH5rP6HNgmQN1///1atmyZ3183ZMgQLV68uB0qap3HY+jYsUpT3ttuoqKcSkyM17FjJ+R2e8wuJyJZbQ4m/bSf5r2xtcXnJ45K17Gj4fXzZbU5iETMgfnMnoPExHifVr9sE6BcLpe6du3q99d17ty5HarxXV0dP4D+cLs9fGYms8ocnHt2V6+3fp97dldL1NkerDIHkYw5MJ/V58A2AWrWrFmaNWuW2WUACCFu/QZgVbYJUEA44XgS33Hrt3/43gJCgwAFhBjHk6A9eDyGVm3cq/f+c0AVVXUNj/O9BbQPa98jCISZ+uNJTm0SWX88SW5+sUmVwc5y84v1/z39od76qKBReJL43gLaS0SuQI0bN05FRUWSJMMwGh6/+OKLG/550KBBmj9/fshrQ/jieBK0h5bODDwV31tAcEVkgDp69KjKysqaPP7Dx44fPx66ghAROJ4EweZLKK/H9xYQXBEZoNavX292CYhAHE+CYPMllP8Q31tA8LAHCgiRpIS4oI4D/A1EfG8BwUOAAkKkX1qSkl3e/wJLcZ287RzwhT+BiO8tILgIUECIOJ0OTRqV7nXMxFHpbPKFz3wJ5fX43gKCiwAFhFBORndNH5fd5C+9FFecpo/LplcP/OJLKO/UIZrvLaAdROQmcsBMHE+CYKoP5ac2Z03oEK2f/jhNPx/Rh+8toB0QoAATcDwJgolQDoQeAQoAwoBdQzln98GuCFAAAFNwLiTsjE3kAICQ41xI2B0BCgAQUr6eC+nxGF7HAGYiQAEAQsqfcyEBqyJAAQBCinMhEQ7YRA4AJmju7rNIwbmQCAcEKAAIsZbuPvvV6AxdOryviZWFRv0RNN4u43F2H6yOS3gAIpbHY2jnvlJt2nFIO/eVhmTTsre7z+a9sVUbtxa1ew1m41xIhANWoABEJDN6EPly99nfV+bpiWkj2uX9raSlI2hSXHGaSB8o2AABCkDEqV8FOlV9D6L2OnzXl7vPjpSdUP7+UqWfnhT09w+11rqMcwQN7IwABSCi+NqD6Nz0bkH/i9znu8+O1wT1fc3g6wqf0+lQv7SkhhC1q7CMEAVbIEABiCj+9CAK9tlyPt991ik2qO8bav6s8HGcC+yKTeQAIoqZPYjq7z7zpmtSvDLOsN+hwPX86TLOcS6wMwIUgIhiZg8iX+4+u+XKbFtfvvJ1hW/n/lKOc4GtEaAARBRfVoHaswdR/d1np9aQ4orTjGsGasTAXu3yvqHi68rdzn2lHOcCW2MPFICIUr8K1NwenXrt3YOopbvPYmOj2u09Q8XnlTsfP16Oc4FVsQIFIOJ4WwVqrxYGp3I6HcrsnaxhWT2V2TvZ1pftfsjXFb7MNN/2eXGcC6yKFSgAEYkeRO3D1xW+zN7JHOcCW2MFCkDECtdVILP5ssLHcS6wO1agAABB58sKH8e5wM4IUACAdlG/wucNl1JhVwQoAICpfAlagNWwBwoAAMBPBCgAAAA/EaAAAAD8xB4oAIBleTwGG8xhSQQoAIAl5eYXN2lxkOyK0yRaHMACuIQHALCc3PxiPbMir0mn8tLyaj2zIk+5+cUmVQacRIACAFiKx2No6drdXse8sna3PB4jRBUBTRGgAACWsquwzOsZeZJUUl6tXYVloSkIaAZ7oAAAllJW4T08+TvOytgkb18EKACApSQlxLU+yI9xVsUmeXvjEh4AwFL6pSUp2eU9HKW4Tq7W2BWb5O2PAAUAsBSn06FJo9K9jpk4Kt22l7rYJB8eCFAAAMvJyeiu6eOym6xEpbjiNH1ctq0vcbFJPjywBwoAYEk5Gd11bnq3sNtkHUmb5MMZAQoAYFlOp0OZvZPNLiOoImWTfLjjEh4AACEUCZvkIwEBCgCAEAr3TfKRggAFAECIhfMm+UjBHigAAEwQrpvkIwUBCgAAk4TjJvlIwSU8AAAAP0XcClRBQYHee+89bd68Wfn5+SotLVVcXJz69u2rSy+9VNddd50SEhLMLhMAAFhYRAUot9ut0aNHN/y5W7duysjI0JEjR7Rt2zZt27ZNr7/+ul566SWlpqaaWCkAALCyiApQhmGoU6dOmjBhgq666iqdddZZDc998cUXuuOOO7R//37dfvvteu2110ysFAAAWFlEBaioqCitW7dOSUlJTZ770Y9+pD/96U+aOHGivvzyS3311Vc655xzQl8kAACwvIjaRO5wOJoNT/UGDRokl8slSdq7d2+IqgIAAHYTUQGqNW63W3V1dZKkDh06mFwNAACwKgLUD6xbt04nTpxQdHS0fvSjH5ldDgAAsKiI2gPlzfHjx/XYY49Jkq6++mqlpKQE5XWjo8movoiKcjb6f4Qec2A+5sB8dp0Dj8dQ/v5SlR2vUVKnWGWckWzbjuZ2mQOHYRiG2UWYze12a+rUqfrggw+UmpqqlStXNuyFagvDMORw2PMbGABgDxu3Fum5t7bp+6NVDY916dxBU8YO0IiBvUysLLzZJkDdf//9WrZsmd9fN2TIEC1evLjF5w3D0L333qs333xTnTt31pIlS9SvX7+2lNrA7fbo2LETQXmtcBcV5VRiYryOHTsht9tjdjkRiTkwH3NgPrvNwZadxZr3xtYWn59xzUANzrTXwcRmz0FiYrxPq1+2uYTncrnUtWtXv7+uc+fOXp9/6KGH9OabbyohIUHPP/980MJTvbo66/8AWonb7eEzMxlzYD7mwHx2mAOPx9CSd/O9jnn53Xz9z5ldbHk5z+pzYJsANWvWLM2aNSuor/nYY49pyZIlio+P13PPPaeBAwcG9fUBwM48HkO7CstUVlGtpIQ49UtLsuVfxOFqV2GZSsurvY4pKa/WrsIyDixuB7YJUME2d+5cvfDCC4qLi9P8+fP14x//2OySAMAycvOLtXTt7kZ/QSe74jRpVLpyMux1SShclVV4D0/+joN/rL3FvZ0sWLBACxYsUExMjObNm6fhw4ebXRIAWEZufrGeWZHXZHWjtLxaz6zIU25+sUmV4YeSEuKCOg7+ibgAtWjRIs2dO1fR0dGaO3euRo4caXZJAGAZHo+hpWt3ex3zytrd8nhscf9RWOuXlqRkl/dwlOI6eekVwRdRl/AOHz6sRx55RJKUkJCgF154QS+88EKzY6+++mpdc801oSwPAEzHvhr7cDodmjQqXc+syGtxzMRR6exbaycRFaBqa2tV37Xh6NGj+uyzz1ocO2LEiFCVBQCWwb4ae8nJ6K7p47Kb7FdLccVpIvvV2lVEBajTTz9d+fneb/kEgEjGvhr7ycnornPTu3HHZIhFVIACAHhXv6/G22U89tVYj9Pp4JJqiEXcJnIAQMvq99V4w74agAAFADhF/b6aU+/wSnHFafq4bPbVAOISHgCgGeyrAbwjQAEAmsW+GqBlXMIDAADwEwEKAADATwQoAAAAPxGgAAAA/ESAAgAA8BMBCgAAwE8EKAAAAD8RoAAAAPxEgAIAAPCTwzAMw+wiwpVhGPJ4+Hh9FRXllNvtMbuMiMYcmI85MB9zYD4z58DpdMjhaP3IIgIUAACAn7iEBwAA4CcCFAAAgJ8IUAAAAH4iQAEAAPiJAAUAAOAnAhQAAICfCFAAAAB+IkABAAD4iQAFAADgJwIUAACAnwhQAAAAfiJAAQAA+IkABQAA4CcCFAAAgJ+izS4AaElBQYHee+89bd68Wfn5+SotLVVcXJz69u2rSy+9VNddd50SEhLMLjNsVVVV6aOPPtK2bduUl5envLw8lZWVSZI+++wzPvsg2rRpk1588UV9+eWXqqysVK9evTRmzBhNmTJFHTt2NLu8sPbdd99p48aNDd/nX331laqqqtS/f38tX77c7PLCnmEY+vzzz7V+/Xrl5ubqm2++0fHjx+VyuZSVlaWxY8fqF7/4hRwOh9mlNuEwDMMwuwjgVG63W1lZWQ1/7tatm3r06KEjR47o0KFDkqQzzjhDL730klJTU80qM6x99dVXGjt2bLPPEaCCZ/HixXr44YdlGIZ69uyplJQU7dmzRzU1NTrrrLO0dOlSJSUlmV1m2HrppZc0Z86cJo8ToELjk08+0eTJkxv+nJaWpsTERH377bcNv7BddNFFmjdvnmJjY80psgWsQMGSDMNQp06dNGHCBF111VU666yzGp774osvdMcdd2j//v26/fbb9dprr5lYafiKjo7WwIEDNWDAAGVnZ6tLly6aMmWK2WWFlby8PD3yyCOSpAceeEDjx4+Xw+HQ4cOHNXXqVG3fvl2zZ8/WvHnzTK40fHXq1EkjRoxQdna2srOzVVBQoCeffNLssiKGYRg6/fTTdeONN+pnP/uZunTp0vDcW2+9pdmzZ+v999/X008/rTvvvNPESptiBQqWZBiGjh492uJv3p999pkmTpwo6eQP2TnnnBPC6iLTgQMHdMkll0hiBSpYpk2bpnXr1mns2LF67LHHGj1XUFCgyy67TB6PRytXrlRmZqZJVUaW5cuX65577mEFKkSOHz+uuLg4xcTENPv8ggULNHfuXCUlJemTTz6R02mdrdvWqQT4AYfD4fWyxaBBg+RyuSRJe/fuDVFVQPBUVFToww8/lCSNHz++yfN9+vTRsGHDJElr1qwJaW1AqHTq1KnF8CRJF154oSSprKxMJSUloSrLJwQo2JLb7VZdXZ0kqUOHDiZXA/jvq6++Uk1NjWJjYzVw4MBmx+Tk5EiSvvzyy1CWBlhGdXV1wz9b7b/1BCjY0rp163TixAlFR0frRz/6kdnlAH6rXznt1atXi7+Bn3HGGY3GApFm9erVkqTMzEx16tTJ5GoaI0DBdo4fP96wX+Tqq69WSkqKyRUB/jt69KgkqXPnzi2OqX+ufiwQSbZv365XX31Vkix5AwsBCrbidrs1c+ZMHThwQKmpqZo1a5bZJQEBqb804W3/R/1t2z+8jAFEgiNHjui2225TbW2tfvrTn+pnP/uZ2SU1QRsDBN3999+vZcuW+f11Q4YM0eLFi1t83jAMzZ49Wx988IE6d+6sBQsWNGwkx3+11+eP4IqLi5Mk1dbWtjimpqam0VggEpSXl+uWW25RUVGR+vfvr0cffdTskppFgELQuVwude3a1e+v83YpQ5Ieeughvfnmm0pISNDzzz+vfv36BVpiWGuvzx/B5cvlOV8u8wHhpKKiQjfffLN27Nih9PR0/eMf/7Dc3qd6BCgE3axZs4J+ae2xxx7TkiVLFB8fr+eee67Fu5bQPp8/gq9Pnz6SpKKiItXW1jZ7KW///v2NxgLh7MSJE7r11lv1xRdfqE+fPnrxxReVnJxsdlktYg8ULG/u3Ll64YUXFBcXp/nz5+vHP/6x2SUBbZaVlaWYmBjV1NRo69atzY7Jzc2VJO40Rdirrq7WtGnTtGXLFqWmpmrhwoXq1q2b2WV5RYCCpS1YsEALFixQTEyM5s2bp+HDh5tdEhAUCQkJOv/88yWp2eOICgoKtGnTJknSmDFjQlobEEq1tbWaMWOGNm7cqJ49e2rhwoXq2bOn2WW1igAFy1q0aJHmzp2r6OhozZ07VyNHjjS7JCCopk2bJofDoZUrV2rZsmWqP1mruLhYM2fOlMfj0ahRozjGBWHL7Xbrzjvv1AcffKBu3bpp4cKFSktLM7ssn3AWHizp8OHDGjlypAzDUOfOnRsdJnyqq6++Wtdcc00Iq4sc48aNU1FRkaT/nk8oqdExO4MGDdL8+fPNKC8svPTSS3r00UdlGIZOO+00JScna8+ePaqpqVHfvn21dOlSep21o4MHD2rs2LENf66pqVFlZaWio6MbbV6++eabdcstt5hQYXhbtWqV7rjjDklSamqqevTo0eLY2bNnKysrK1SltYpN5LCk2traht/Gjx49qs8++6zFsSNGjAhVWRHn6NGjKisra/L4Dx87fvx46AoKQ5MnT1ZGRoZeeOEFbd26Vd9//7169eqlMWPGaMqUKRza3M7cbnez3+N1dXWNHq+qqgpdURGkvlWHJH377bf69ttvWxxbXl4eipJ8xgoUAACAn9gDBQAA4CcCFAAAgJ8IUAAAAH4iQAEAAPiJAAUAAOAnAhQAAICfCFAAAAB+IkABAAD4iQAFAADgJwIUAACAnzgLDwBs7tixY/r000+1fft27dixQ9u3b9d3330nSZozZ46uuuoqkysEwg8BCgBsbu3atbrnnnvMLgOIKAQoAAgD3bp10znnnKOsrCz1799fM2bMMLskIKwRoADA5q644gou0wEhxiZyAKb67rvv9MQTT+iKK65QTk6OBgwYoEsuuUT33nuv9uzZ02T8c889p4yMDGVnZ2vr1q3NvuYHH3ygzMxMZWRk6O23327y/ObNm/Xb3/5WF1xwgbKzszV06FDdeOONevPNN+V2u73Wu2XLFv3mN7/R0KFDNXDgQI0ePVpz585VRUWFli9froyMDF188cV+fQbvvPOOMjIylJGRoffee6/ZMV999ZUGDBigjIwMLViwoNFz0dH8LgyEmsMwDMPsIgBEpg0bNmjmzJmqrKyUJMXExCgmJqbRnx966CGNHTu24WsMw9Cvf/1rbdy4UWlpaXrrrbfUqVOnhueLi4t15ZVXqqSkRGPHjtVjjz3W6D3nzJmjl156SZLkcDjkcrlUUVHREJyGDRumZ555ptFr1lu8eLEefvhh1f9n0+VyqaqqSrW1tTrrrLM0fvx4zZkzR6mpqVq/fr1fn8Xvf/97vfnmm0pKStJbb72l0047reG5yspKXXXVVdq7d6+GDBmihQsXyun0/vtvRkZGw78vq1NA8LECBcAUW7du1YwZM1RZWalrr71W77zzjr788kt9/vnn2rBhgyZNmqTa2lrde++92rZtW8PXORwOPf744+rSpYsKCwv1hz/8oeE5wzB09913q6SkRL1799b999/f6D2XLFnSEJ6uvfZaffjhh9qyZYv+85//6J577lF0dLQ2bdqk2bNnN6n3s88+0yOPPCLDMHTeeedpzZo1+s9//qMvvvhCTz31lL7//ns9++yzAX8es2fP1plnnqmysjLNmjWr0UrYgw8+qL179yopKUlPPPFEq+EJQPvjpxCAKR544AHV1tZq2rRpeuCBB3TWWWcpKipKktSrVy/94Q9/0PXXX6+6ujrNnz+/0dd269ZNc+bMkcPh0KpVq7RixQpJ0t///ndt3LhRMTEx+vOf/6yEhISGr6mqqtK8efMkST//+c/1wAMPqFu3bpKkjh07avLkyfrd734n6eQltR+GNkl6+umn5fF4dPbZZ2vBggXq27evpJOXz8aMGaOnnnpKR48eDfjziI+P19y5cxUbG6stW7Y0/Du/8847Wr58uSTp4YcfVo8ePQJ+DwDBQ4ACEHI7d+7Utm3bFBMTo1//+tctjqu/dPfJJ5802Zs0cuRITZ48WdLJMLZy5Uo9/fTTkqTbb79dAwYMaDT+448/VllZmSTptttua/b9Jk2a1BCqVq9e3fB4WVmZNm3aJEm66aabFBsb2+Rrhw0bph//+Mct/rv4IjMzU3fddZck6dlnn9Vbb73VsIp23XXXadSoUW16fQDBw85DACGXm5srSfJ4PBozZkyL4+pDU2VlpcrKytSlS5dGz8+cObOhgWR98Dj//PObDWV5eXmSpNNOO61h9ehUUVFRGjZsmN5+++2G8dLJDdz1+54GDx7cYr1DhgzRf/7znxaf98X111+vjz/+WBs2bNDdd98tSerXr1/DPwOwBlagAIRccXGxpJMB6ciRIy3+r7S0tOFrTpw40eR1YmNj9eijjzb82eVy6dFHH5XD4Wgy9vvvv5ekVi+B9ezZs9F4SSopKWn4Z29fH6zLa4888oji4uIknQx1Tz75ZMOfAVgDK1AAQq5+ZenMM8/UP//5zza91muvvdbwz8ePH9dXX33VcBmuOc2Fq9bG+XqzcrBual65cqWqq6slnfyscnNzlZ6eHpTXBhAcrEABCLmuXbtKkgoLCxtaFgRiw4YNWrx4saSTt+0bhqHf/e53OnLkSJOx9Zf/Dh486PU1Dx06JElKSUlp8rXSf1fPmuPtOV9t375df/7znyU1bkXQXE8sAOYhQAEIuUGDBkmSamtrW2wc2Zri4uKG89+uuuoqvfzyy0pNTdX333+vu+++u8lqUHZ2tqSTAWnv3r3Nvqbb7dbmzZslqdEm9HPOOadhRerTTz9tsSZvz/misrJSM2fOVG1trYYNG6Y33nhDAwcOVFVVlWbOnKmampo2vT6A4CFAAQi5AQMGKCsrS5L0l7/8pdEeo+bU3z1Xz+Px6K677lJpaan69Omj2bNny+Vy6c9//rOio6P10Ucf6cUXX2z0Needd56SkpIkSX/961+bfZ9XX321YRXpZz/7WcPjSUlJGjp0qCTpxRdfbDbI1PeTaosHH3xQBQUFSkpK0uOPP67Y2NiGdgz5+flNmoICMA8BCkDIORwO/fGPf1RsbKyKior0y1/+UmvWrGm0Ufzw4cNauXKl/vd//1dPPPFEo69//vnn9cknnzT0e+rYsaMk6dxzz9X06dMlSU8++aS2b9/e8DUdOnRoOGB31apVuv/++xsu9Z04cUKLFy/WnDlzJEmXX355w4pVvRkzZsjhcGjXrl2aOnWqCgoKJEl1dXX617/+pRkzZqhz584BfyarV69u6Pf0yCOPNGxIP+OMMxqahS5ZskQbNmxo9utLSkoa/a9eZWVlo8eb24wPwH8c5QLANB9//LFmzpzZsMIUFRUll8ul6urqRn/R//KXv9RDDz0k6WQH8/ou5XfddZduuummRq/p8Xh044036tNPP1WfPn20YsWKhoAlNT3KJTExURUVFaqrq5MkDR06VM8++2yzR7m89NJLDSFLkhITE1VVVaWamhr169dPV199tebMmaO+fftqzZo1Pn8OBw4c0NixY1VeXq7rrruuSQd1Sbrrrru0cuVKJScna+XKlU3u+KvfL9Wa2267rSFIAggcK1AATHPeeefpvffe0x133KGcnBy5XC6Vl5fL4XDo7LPP1jXXXKP58+c3HK1y/Pjxhj1C5513XrP9npxOp/70pz8pKSlJBQUFevDBBxs9f88992jhwoUaPXq0unbtqsrKSiUkJGjo0KF65JFH9OKLLzYbniRp8uTJWrx4sUaOHKnOnTururpaqampmjp1qpYtW9YwLjEx0efPoK6uTnfccYfKy8u99nu6//771bt3b5WWluruu++Wx+Px+T0ABB8rUAAQJHfccYdWrVqlq6++Wo888ojZ5QBoR6xAAUAQ7N27t+GOwgsuuMDkagC0NwIUAPjoqaee0pIlS1RUVNRwCa2yslLvvPOObrjhBlVXV+vMM8/kzDogAnAJDwB8NG3aNK1bt06SFBMTo4SEBB07dqwhTPXo0UPPP/+8+vXrZ2aZAEKAo1wAwEeTJ09W9+7d9fnnn+u7777T0aNHlZCQoD59+uiiiy7Sr371q4ZeUwDCGytQAAAAfmIPFAAAgJ8IUAAAAH4iQAEAAPiJAAUAAOAnAhQAAICfCFAAAAB+IkABAAD4iQAFAADgJwIUAACAn/5/eLVo0lcXI1kAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAHMCAYAAAAarpbgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABDzElEQVR4nO3deXjU1d3//9dkJYQJybAKoqAmQQgopLKJS5Uq0lqDuADW5XbBG6i9v2LBWsW22rpVRYsKWosIFJcqy6W4VBCtyiINIoQlgBIIBgiQTEiAbDOf3x/8JhKSTGbfPs/HdXldMnNm5ihnZl5zzvm8j8UwDEMAAADwWFy4OwAAABBtCFAAAABeIkABAAB4iQAFAADgJQIUAACAlwhQAAAAXiJAAQAAeIkABQAA4CUCFAAAgJcIUAAAAF5KCHcHwuXTTz/VF198oc2bN2v//v0qLy9XQkKCunfvrqFDh+q2225T9+7dw91NAAAQgSxmPQvv5ptv1tdff63ExER16tRJNptN5eXl2rdvn5xOp1JSUvTCCy9o+PDh4e4qAACIMKYNUEuWLFGXLl2Um5urpKSkhtv37Nmj3//+91q3bp0yMjL06aefqm3btmHsKQAAiDSmDVDuHDp0SBdeeKEk6ZVXXtEll1wS5h4BAIBIYto9UO507NhR6enpstvtqq6u9vl5DMOQ00k+NbO4OAtjwOQYA2AMRJe4OIssFkur7QhQzfjuu+9kt9sVFxenPn36+Pw8FotFR44cU329M4C9Q7RISIhTRkYqY8DEGANgDEQfmy1V8fEEKI8ZhqGysjLl5+fr6aefliTdfvvt6tGjR5h7BgAAIo3pA9TSpUs1bdq0RredddZZevrpp3X11Vf7/fzx8ZTaMivX3z1jwLwYA2AMxC7TB6gOHTpo4MCBMgxD+/fv14EDB1RUVKT33ntPF1xwgbp27erX86elpQSop4hWjAEwBsAYiD1chXeK4uJiPfHEE1q+fLm6dOmiZcuWyWq1+vx8R44cl8PBurcZxcfHKS0thTFgYowBMAaiT1paikczhqafgTpVjx499Le//U3XXHONduzYoQULFmjixIk+P5/D4WTjoMkxBsAYAGMg9rAo24z4+HhddNFFkqSCgoIw9wYAAEQaAlQL6uvrJUlOJ78YAABAYwSoZtTW1uqzzz6TJL/qQAEAgNhkygC1adMmPffccyoqKmpy365duzRx4kTt2bNHbdu21Q033BD6DgIAgIhmyk3kx44d06xZszRr1izZbDaddtppSkhI0MGDB1VSUiJJSk9P13PPPacuXbqEubcAgGjndBratrtc9qM1Sk9NVlaPdMXFtV7tGpHLlAGqd+/eeuihh/T1119r+/bt2r17t6qrq9WuXTvl5ubqoosu0o033iibzRburgIAotyqjSV6edFGlVXWNNyWYU3W+BGZys3uHMaewR/UgQqy8vKjXLpqUq4zsBgD5sUYwDc7D2nmOxtbvH/y6BxCVIQ5cRZe6zucTLkHCgCAYHM6Df3z40K3bd5YvkNOJ/MY0YgABQBAEGwvtjdatmtOWWWNthfbQ9MhBBQBCgCAILAfdR+evG2HyEKAAgAgCNJTkwPaDpGFAAUAQBBk9UiXzeo+HNmsJ0oaIPoQoAAACIK4OItuujLbbZtxIzKpBxWlCFAAAATJBb0764FbL2gyE2WzJlPCIMqZspAmAAChMqx/N2V3T9OWXWVUIo8hBCgAAIIsLs6i3mdmhLsbCCCW8AAAALxEgAIAAPASAQoAAMBLBCgAAAAvEaAAAAC8RIACAADwEgEKAADASwQoAAAAL1FIEwAABJ3TaWh7sT1mqrEToAAAQFDlF5Zq4fIdKq+sabgtw5qs8SMyo/Y8QJbwAABA0OQXlurFxQWNwpMklVfW6MXFBcovLA1Tz/xDgAIAAEHhdBpauHyH2zZvLN8hp9MIUY8ChwAFAACCYnuxvcnM06nKKmu0vdgemg4FEAEKAAAEhf2o+/DkbbtIQoACAABBkZ6aHNB2kYQABQAAgiKrR7oyrO7Dkc16oqRBtCFAAQCAoIiLs2j8iEy3bcaNyIzKelAEKAAAEDS52Z01eXROk5komzVZk0fnRG0dKAppAgCAoMrN7qwBmZ2oRA4AAOCNuDiLep+ZEe5uBAwBCgDQINbOKwOChQAFAJAUm+eVAcHCJnIAQMyeVwYECwEKAEwuls8rA4KFAAUAJhfL55UBwUKAAgCTi+XzyoBgIUABgMnF8nllQLAQoADA5GL5vDIgWAhQAGBysXxeGRAsBCgAQMyeVwYEC4U0AQCSYvO8MiBYCFAAgAaxdl4ZECws4QEAAHiJAAUAAOAlAhQAAICXCFAAAABeIkABAAB4iQAFAADgJQIUAACAlwhQAAAAXiJAAQAAeMl0lcgNw9A333yjTz/9VPn5+fr+++9VVVUlq9WqPn36KC8vT1dffbUsFo4uAAAAzTNdgFqzZo1uu+22hj/36NFD3bt31w8//KCvvvpKX331lZYtW6aZM2cqKSkpfB0FAAARy3QByjAMnX766br11lv185//XB06dGi4b8mSJZo+fbo+++wz/e1vf9Nvf/vbMPYUAABEKothGEa4OxFKVVVVSk5OVmJiYrP3z549WzNmzFB6erpWr16tuDj/tomVlx9Vfb3Tr+dAdEpIiFNGRipjwMQYA2AMRB+bLVXx8a1/95tuE3m7du1aDE+SdPHFF0uS7Ha7ysrKQtUtAAAQRUwXoFpTU1PT8O9t2rQJY08AAECkIkCdYtmyZZKk3r17q127dmHuDQAAiESm20TuzubNm/Xmm29KkiZMmBCQ5/RkHRWxyfV3zxgwL8YAGAOxy3SbyFty6NAhXX/99SopKdHPfvYzvfDCC+HuEgAAiFAEKEmVlZW65ZZbtGXLFvXt21fz5s0L2PLdkSPH5XBw5YUZxcfHKS0thTFgYowBMAaiT1paikczhqZfwjt69KjuvPNObdmyRZmZmfrHP/4R0L1PDoeTS1dNjjEAxgAYA7HH1Iuyx48f1913360NGzaoZ8+eeu2115SRkRHubgEAgAhn2gBVU1OjSZMmad26derevbtef/11derUKdzdAgAAUcCUAaqurk733HOPVq1apa5du+r1119X165dw90tAAAQJUwXoBwOh37729/q888/V6dOnfT666+rR48e4e4WAACIIqbbRP7hhx/qo48+kiQlJSXpgQceaLHt9OnT1adPn1B1DQAARAnTBaja2tqGf//hhx/0ww8/tNi2srIyFF0CAABRhjpQQcYJ3ObFKexgDIAxEH1stlSP6kCZbg8UAACAvwhQAAAAXiJAAQAAeIkABQAA4CXTXYUHALHO6TS0vdgu+9EapacmK6tHuuLiLOHuFhBTCFAAEEPyC0u1cPkOlVfWNNyWYU3W+BGZys3uHMaeAbGFJTwAiBH5haV6cXFBo/AkSeWVNXpxcYHyC0vD1DMg9hCgACAGOJ2GFi7f4bbNG8t3yOmk9B8QCAQoAIgB24vtTWaeTlVWWaPtxfbQdAiIcQQoAIgB9qPuw5O37QC4R4ACgBiQnpoc0HYA3CNAAUAMyOqRrgyr+3Bks54oaQDAfwQoAIgBcXEWjR+R6bbNuBGZ1IMCAoQABQAxIje7syaPzmkyE2WzJmvy6BzqQAEBRCFNAIghudmdNSCzE5XIgSAjQAFAjImLs6j3mRnh7gYQ0whQAADEMM5GDA4CFAAAMYqzEYOHTeQAAMQgzkYMLgIUAAAxhrMRg48ABSConE5D23aXa82W/dq2u5wPbCAEOBsx+NgDBSBo2H8BhAdnIwYfM1AAgoL9F0D4cDZi8BGgAAQc+y+A8OJsxOAjQAEIOPZfAOHF2YjBR4ACEHDsvwDCj7MRg4tN5AACjv0XQGTgbMTgIUABCDjX/gt3y3jsvwBCg7MRg4MlPAABx/4LALGOAAUgKNh/ASCWsYQHIGjYfwEgVhGgAAQV+y8AxCKW8AAAALxEgAIAAPASAQoAAMBL7IECACAGOJ0GF2yEEAEKAIAol19YqoXLdzQqXpthTdb4EZmUDAkSlvAAAIhi+YWlenFxQZPK/+WVNXpxcYHyC0vD1LPYRoACACBKOZ2GFi7f4bbNG8t3yOk0QtQj8yBAAQAQpbYX292eOSlJZZU12l5sD02HTIQABQBAlLIfdR+evG0HzxGgAACIUumpya038qIdPEeAAgAgSmX1SG9yYPepbNYTJQ0QWAQoAACiVFycReNHZLptM25EJvWggoAABQBAFMvN7qzJo3OazETZrMmaPDqHOlBBQiFNAACiXG52Zw3I7EQl8hAiQAEAEAPi4izqfWZGuLthGizhAQAAeIkABQAA4CUCFAAAgJdMuQfq4MGDWrVqlTZt2qSCggJt3bpV1dXV6tu3rxYtWhTu7gEAgAhnygC1bNkyPf744+HuBgAAiFKmDFDt2rXTsGHDlJOTo5ycHBUVFenZZ58Nd7cAAECU8DlArVu3ThdccEEg+xIy1113na677rqGP7NsBwAAvOHzJvKbb75Zo0aN0ty5c2W32wPYJQAAgMjm11V433//vZ588kldfPHFmjp1qtatWxeofgEAAEQsnwPUv//9b915553q0KGDamtr9f777+uWW25hVgoAAASN02lo2+5yrdmyX9t2l8vpNMLSD4thGH69cn19vT799FO99dZbWr16tZxOpywWi5KSknTFFVfohhtuiPi9UosWLdIDDzwQlDIGR44cl8PhDOhzIjrEx8cpLS2FMWBijAEwBgJr3bZS/fPjQpVV1jTcZrMm66Yrs3VB78AcmpyWlqL4+Nbnl/wOUCcrKSnR22+/rcWLF+vAgQMnXsBiUa9evXTDDTcoLy9P6enpgXq5gAlmgAIQexxOQ1u+P6yyI9WypbVRn7M6KJ5DW4GgWrWxRI+/3vJWoQduvUDD+ncLWX8CGqBcnE6nPvvsM/3rX//Sf/7zHzkcjkazUmPHjlVubm6gX9ZnzEAh0JxOQzt+qFBNvaHkBIsyu7fnVPQY4c0vYGYfwBgIDKfT0JSZXzZ6353KlpasZ3893O/PWk9noIJSByouLk6XXXaZEhISZLfbtWHDBhmGoZqaGr333nt6//33dd555+n3v/+9+vfvH4wuRAyHw6n6et40ZpJfWKqFy3eo/KQ3eoY1WeNHZCo3OzBTzAiP/MJSvbi4oMntZZU1mvnORk0endPs3zGfA2AM+Gfb7nK34UmSyo7UaMuuMvU+MyMkfQr4WXilpaV66aWXdNlll+nuu+/WN998I8MwlJubqwcffFCXXnqpLBaLNmzYoPHjx2vt2rWB7gIQNq4v2PJT3ujllTV6cXGB8gtLw9Qz+MvpNLRw+Q63bd5YviNsG1qBWGY/6j48edsuEAIyA2UYhj7//HO9/fbbDUt2hmGoXbt2uuaaazR27FhlZmZKOlE/qri4WH/84x/11Vdf6fnnn9fChQsD0Q0grDz9gh2Q2YnlvCi0vdjeJBifqqyyRtuL7SH7BQyYRXpqckDbBYJfAWr//v1655139O6772r//v1ybafq06ePxo0bp1/84hdKSUlp8rgePXro+eef19ChQ1VYWOhPF4CIwRdsbIvEX8CAWWT1SFeGNdntZ6zNmqysHukh65PPAeruu+/Wl19+KafTKcMwlJKSoquuukpjx471aF9Tu3bt1KlTJ+3bt8/XLgARhS/Y2BaJv4ABs4iLs2j8iMxm9yC6jBuRGdLZfZ8D1Oeffy5JOvvss3XjjTdq9OjRslqtXj3HlVdeGZaCm/v27VNeXl7Dn2trayVJhYWFGjx4cMPtd955p+66665Qdw9Rii/Y2BaJv4ABM8nN7qzJo3OaXKRjsyZrXBgu0vE5QI0aNUpjx47VoEGDfH7x+++/3+fH+sPhcDQb3Orr6xvdXl1dHbpOIerxBRvbIvEXMGA2udmdNSCzk7YX22U/WqP01BOfqeF43wWlDhR+VF5+lEtXTaSly9xdWrrMHdGjuTIVLf0CTkiIU0ZGaqPPAafTiIgPf4RGc2MAkc1mSw19JXI0xZvGfLz5gkV08jQEnfrlSY0w8yFARR8CVITgTWNOTqeh70oqVGdYlGgxdHY3KpGb0clfnms372d20oQIUNHH0wAVlErkgNnFxVl0bk8bH5yQRI2wSOfPsipLsuZFgAKAICvcU06NsAjlz7IqS7LmFvCjXAAAjdmraj1rR42wkPLn6CWObQIBCgCCLL1dkmftqBEWMv6cbRjN5yI6nYa27S7Xmi37tW13eUT2MVqwhAcAQZZ9RgY1wiKMP0cvefPYnLM7+N3XQGHJMbCYgQKAIHMV4XSHIpyh5c/RS9F4bBNLjoFHgAKAEHAdQ5FhbbxMZ7MmU8IgDPw5einajm2K5iXHSObzEl5JSYlX7ZOTk2W1WpWU5NleAACINZF0DIXZ+XP0UrQd2+TPcmUwRXsJCJ8D1OWXX+7T43r06KGLL75YN998s84880xfXx4AolJcnIVSBRHAn7MNo+1cxEhccoyF/Vg+L+EZhuHTP3v27NE///lPXXPNNfrwww8D+d8CAAHFFUuxzZ9l1Whako20JcdY2Y/l8wzUihUrtGnTJv3hD39QXFycxo0bp0GDBqlLly4yDEOlpaX6+uuv9eabb8rhcOjRRx9Vjx49tGnTJs2bN087d+7U/fffrz59+jATBSDixMIvZCn6l0mCzZ9l1WhZko2kJcdYqsrv81l4e/bs0ZgxY3T66adrzpw5yshofkq6vLxct99+u0pKSrRo0SJ1795dtbW1uuWWW/Ttt99q/Pjxmj59ul//EZGMYzzMizOwopfrF3JLPJ1hCPcYiJUQGM3CPQZcAjWm/bVtd7meeuObVttNGzcgbEvdnp6F5/MS3qxZs1RVVaVHH320xfAkSRkZGXrkkUdUUVGh2bNnS5KSkpJ03333yTAMrVmzxtcuAEDAxcoVS7GyTILAiJQlx0jcj+Urn5fwVq1apbZt2yonJ6fVtv369VPbtm315ZdfNtw2cOBAJSYmat++fb52AQACLlKvWPJGLC2TIHAiYckx0vZj+cPnAFVWVqaEBM8fbhiGDh8+3PDn+Ph4tW3bVtXV1b52AQACLhZ+IcdCCERwhPsq0Ejaj+Uvn5fwbDabqqurtW7dulbbrlu3TsePH2+01FdXV6cjR464Xf4DgFALxS/kYF/dFwshELEplqry+zwDdeGFF2rRokV66KGH9Oqrr6pHjx7Nttu7d68eeughWSwWDR8+vOH2oqIiGYahbt26+doFAAi4YP9CDsXG7lhaJkHsce3HOvV9YLMma1wUXeDgc4D69a9/rY8//lh79uzR1VdfrauuukqDBg1S586dZbFYVFpaqrVr1+qjjz7S8ePHlZqaqkmTJjU8ftmyZZKkQYMG+f9fAQABEswiiS1dCeXa2B2ozbyxtEyC2BQJ+7H85XMZA0n65ptvdM899+jQoUOyWJr/jzYMQx07dtTzzz+v3Nzchts/+OADHTx4UD/96U91xhln+NqFiBfuS1cRPpFy+TJ809xMkbe/kE8eA7W1Dk2dtarVUPPUxGEB+RKJlMvWzY7PgejjaRkDvwKUJFVWVmr+/Pn6+OOPtXPnTjkcDkknNomfc845uvLKK/WrX/1KaWlp/rxM1OJNY158cEY/f4tQnjwGCr47HPL6N4EIgfAPnwPRx9MA5fMSnovVatWkSZM0adIk1dXVqaKiQoZhKD09XYmJif4+PQCETSCvWArHxu5YWCaJJFR1x8n8DlAnS0xMVMeOHQP5lAAQE8K1sTvcl63HCqq641Q+lzFojsPhUFlZmcrKyhqW8gAAP27sdoeN3ZGJqu5ojt8zUMePH9ebb76p999/X4WFhY32QPXu3Vu/+MUvdOONNyolJcXvzgJAtArm1X0IHqq6oyV+Bajvv/9eEydO1J49e3TqXvT6+noVFBRo8+bNeuONNzR79mz16tXLr84CiB2u/STllTWqPFardqmJsrVrE9P7SmKl/o2ZUNUdLfE5QFVVVemOO+7Qvn37lJCQoJ/97GcaNmyYunbtKknav3+/Vq9erX//+9/avXu37rjjDr333ntKTU0NWOcBRKfm9pO4xPq+EjZ2RxequqMlPgeo119/Xfv27VPnzp318ssv69xzz23S5vrrr9e2bds0YcIE7du3T/PmzdPEiRP96jCA6NZafaJAF5WMRGzsjh5UdUdLfN5EvmLFClksFj3yyCPNhieX3r1769FHH5VhGPrkk098fTkAMcCT/SQubyzfEfAz4gBvsfkfLfE5QO3evVtJSUm69NJLW2178cUXKzk5Wbt37/b15QDEAE/2k7i49pUA4RRLh98isHwOUPX19R4XyrRYLEpMTFR9fb2vLwcgBni7T4R9JYgErs3/p85E2azJMb3UDPd83gPVtWtX7dmzRzt27FBmpvt0vn37dlVVVenMM8/09eUAxABv94mUlh0PUk8A77D5H6fyeQZqyJAhMgxDf/rTn1RT0/KvxJqaGv3pT3+SxWLR0KFDfX05ADHAk/0kJ1vy5S6KFCJiuDb/D+nTVb3PzCA8mZzPAerOO+9UUlKS8vPz9ctf/lL/+te/tHfvXtXV1amurk7FxcX617/+pV/+8pfKz89XYmKi7rjjjkD2HUCU8WQ/yanYTA4gElmMUytgeuGDDz7QtGnTVF9fL4ul+SRuGIYSEhL01FNPadSoUT53NFpxArd5cQp7y9zVgWrOtHEDovKyf8ZA6EXagb+Mgehjs6UqPr71+SW/KpGPGjVKZ555pmbMmKGvvvqqSTXyuLg4XXTRRbr33nvdljoAYC6u/SRLvvhe769u/epcNpPDExz4i1Dy+yy8vn376tVXX1VlZaU2b96ssrIySZLNZlPfvn1ltVr97iSA2BMXZ1GfnjaPAhRFCtGalgq0mqEwK8LD7wDlYrVaNWTIkEA9HQATcG0qd7eUR5FCtIYDfxEOPm8iBwB/UaQQgeDNgb9AoBCgAIQVRQrhLw78RTh4tIQXqA3gFotFW7ZsCchzAYgdFCmEPzjwF+HgUYDyo9IBAHjEVaQQ8BZ76RAOHgWoefPmNXv73r179cQTT+jYsWO68sorNWTIEHXt2lWSdODAAa1Zs0Yff/yx2rZtq9/97nfq3r174HoOAIB+3EvX3FV4LuylQ6D5XEjz8OHDysvLU0JCgl555ZUWz8PbuXOnJkyYIIfDocWLF8tms/nV4WhD8TTzCmUBvUgrHogTKKIYWs3VgbJZkzUujHWgGAPRJ+iFNGfNmqVDhw7p73//u9vDhM855xw98sgjuvPOOzVr1iw9+OCDvr4kgGZQPBA4gb10CCWfr8L77LPPlJycrOHDh7fadvjw4WrTpo1Wrlzp68sBaIareOCpez9cxQM5iBdmw4G/CBWfA1Rpaani4+M9bh8fH6+DBw/6+nIATuFp8UAO4gWAwPM5QKWlpenYsWMqKGh5055LQUGBjh49yrEuQABRPBAAwsfnADV48GAZhqHp06ervLy8xXZ2u13Tp0+XxWLR4MGDfX05AKegeCAAhI/Pm8gnT56sTz75RNu2bdOoUaM0btw4DR48WF26dJHFYtH+/fu1du1avfnmmyorK1NycrImTZoUyL77bc2aNXrttdf07bff6tixY+rWrZtGjhypCRMmqG3btuHuHuAWxQMBIHx8LmMgSf/5z380ZcoUVVVVyWJpfqOeYRhKTU3Vs88+q0suucTnjgba/Pnz9Ze//EWGYahr166y2WzauXOnamtrdfbZZ2vhwoVKT0/3+3W4dNW8gn35stNpaOqsVa0WD3xq4jA20oYJl7CDMRB9PC1j4NdZeBdffLGWLVumsWPHKi0tTYZhNPonLS1NY8eO1fvvvx9R4amgoECPPfaYJOmRRx7RZ599psWLF2v58uXq27evvvvuO02fPj3MvQTc4yBeAAgfv2agTlVcXKyysjJJks1mU48ePQL11AE1adIkrVixQnl5eXryyScb3VdUVKSrrrpKTqdTS5cuVe/evf16LX51mFeofnlGYvFAnMDsAxgD0SfohTSb06NHj4gNTS5Hjx7VF198IUm64YYbmtzfs2dPDRkyRKtWrdJHH33kd4ACgo3igYg2VM5HLAhogIoGW7duVW1trZKSktS/f/9m2+Tm5mrVqlX69ttvQ9w7wDccxItoQeV8xIqABCin06mioiJVVFSovr7ebdsLLrggEC/ps127dkmSunXrpsTExGbbnHHGGY3aIjj4FQqYi6ty/qlclfMnj86JyBDFZxWa41eAKi0t1bPPPquPP/5Y1dXVrba3WCzasmWLPy/pt4qKCklS+/btW2zjus/V1h+erKOa0bptpfrnx4UqO2Xfzk1XZuuC3pH3AeoL1989Y8C8GAM/cjoNvdFa5fwVO3TBuV0iKpz4+1nFGIhdPgeoAwcO6IYbblBpaak83YcewP3qPqupOfEmaGn2SZKSkpIatfVHWlqK388Ra1ZtLNHMdzY2ub2sskYz39moB269QMP6dwtDz4KDMQDGgLRp56FGIaQ5ZUdqVFJerX7ndAxRr9wL5GcVYyD2+BygXnjhBR04cECpqam69957dfnll6tz585enY8XDsnJJ4oK1tXVtdimtra2UVt/HDlyXA4HV164OJ2GXl7U9APpZC8v3qjs7mkR9SvUF/HxcUpLS2EMmBhj4EfF+zyb0S/eV6HTO4Q/bATqs4oxEH3S0lKCexXef/7zH1ksFv3lL3/RyJEjfX2akPNkec6TZT5PORxOLl09ybbd5R79Ct2yqyxmNkUzBsAYkKwpLc/6n9ouEv5fBfqzijEQe3xelC0rK1N8fLxGjBgRyP4EXc+ePSVJJSUlLc5C7dmzp1FbBA7ntwHmlNUjXRlW97P6NuuJDdqRgM8qtMbnANWhQwe1adNGCQnRVQmhT58+SkxMVG1trTZubH56Nj8/X5J0/vnnh7Bn5sD5bYA5RVvlfD6r0BqfA9TQoUN19OhRFRUVBbA7wZeamqrhw4dLkt5+++0m9xcVFWnNmjWSFFVLk9Ei2n6FAgic3OzOmjw6p8lngM2aHHElDPisQmt8DlD/+7//q5SUFD399NOB7E9ITJo0SRaLRUuXLtVbb73VcHVgaWmppkyZIqfTqREjRlCFPAii7VcogMDKze6sv04cpmnjBmjCL/to2rgBemrisIgKTxKfVWidX2fhrV27Vr/5zW/Up08f3X333erfv7/atm0byP4Fzdy5c/XEE0/IMAyddtppysjI0M6dO1VbW6tevXpp4cKFstlsfr8O5x81zwznt3EGFhgD0c/fzyrGQPTx9Cw8nwPUueee6/VjIqGQ5slWr16tOXPmaOPGjTp27Ji6deumkSNHasKECUpNTQ3Ia/CmaVmsV/flgxOMgdjgz2cVYyD6BP0w4UgoiumvoUOHaujQoeHuhmlxfhuAaMBnFZrjc4CaN29eIPsBAAAQNXwOUIMGDQpkPwAAAKJGdBVxAiJArO/dAgC0LmAByjAMlZeXq7q6Wt26xc5BsMDJmrsiJ8OarPExdPUgAKB1fgeozZs3a9asWVq1apWOHz/e5Eq7iooKPfPMM5Kkhx56SElJSf6+JBAW+YWlenFxQZPbyytr9OLigogrBGgGzAbCE4wTBINfAWrJkiV66KGHVF9f32Kb9u3ba+/evVq9erUuu+wyXXrppf68JBAWTqehhct3uG3zxvIdGpDZiQ/mEGE2EJ5gnCBYfK5E/t1332n69Omqr6/XzTffrHfffVcZGc1f5nnNNdfIMAytWLHC544C4bS92N7oA7g5ZZU12l5sD02HTM41G3jq34lrNjC/sDRMPUMkYZwgmHwOUK+99prq6up000036cEHH1Tfvn0VHx/fbNshQ4ZIkjZs2ODrywFhxcnskcPT2UCnM/pr1cF3jBMEm88Bas2aNbJYLLrrrrtabdulSxelpKSopKTE15cDwoqT2SMHs4HwBOMEweZzgCotLVVKSoq6du3qUfvk5GTV1PDrHNGJk9kjB7OB8ATjBMHmc4BKSkpSXV2dR0e6VFdXq7KyUu3atfP15YCw4mT2yMFsIDzBOEGw+Rygunfvrvr6ehUVFbXa9vPPP5fD4dA555zj68sBYZeb3VmTR+c0mYmyWZMpYRBCzAbCE4wTBJvPZQwuuugiFRYWat68efrDH/7QYrvy8nL99a9/lcVi0SWXXOLrywERITe7swZkdqKmTBi5ZgObq8nl4m42kJpA5uDvOAFa43OAuu2227Rw4UK9+eab6tChg2677bZG91dXV+uTTz7RjBkzVFJSooyMDI0bN87f/gJhx8nsoXdq6BmQ2UmTR+c0qe9jsyZrnJv6PtQEMhfXrLG34wTwhMXwZBNTC1auXKnf/OY3qq+vV0JCggzDkMPh0FlnnaXi4uKGPVJJSUl6+eWXNXTo0ED2PSqUlx9Vfb0z3N1AGCQkxCkjI5Ux4Cd3oceb2cCWKsm7BGMZljEQGcI568gYiD42W6ri41vf4eRXgJKkjRs36pFHHlFBQfMfTH369NEf//hH9e/f35+XiVq8acyLD07/BSr0OJ2Gps5a5faydps1WU9NHBbQL1bGABgD0cfTAOX3WXj9+/fXO++8o23btik/P1+lpaVyOp3q2LGjBg4cqH79+vn7EgBMKJDH53hTE4jlWQCe8DtAufTu3Vu9e/cO1NMBMLlAhh5qAgEINJ/LGABAMAUy9FATCECgEaAARKRAhh5qAgEINAIUgIgUyNBDJXkAgUaAAhCRAh16qCTvGafT0Lbd5VqzZb+27S6X0+nXhdpAzPK7jAHc49JV8+Ly5cBorg6UP4UQQ1kTKNrGAIVGAy/axgBCWAcK7vGmMS8+OAMnWo9fiaYxEI5Co2YQTWMAJ4SsDhQABBvH5wRXIGtuAWbBHigAMDlvam4BOIEABQAmR6FRwHsEKAAwOQqNAt5jDxQA+CBaN7Y3x1Vzq7XDlik0CvyIAAUAXoq1y/1dNbfcXYVHoVGgMZbwAMALrsv9T52tKa+s0YuLC5RfWBqmnvmHQqOAd5iBAgAPxfrl/rnZnTUgs1PMLE0CwUSAijGxtC8DiDTeXO4frXWrIr3mFp9xiBQEqBgSa/sygEjD5f7hxWccIgl7oGJErO7LACIJl/uHD59xiDQEqBjg6b4MTlUH/OO63N8dLvcPPD7jEIkIUDGAYxiA0HBd7u8Ol/sHHp9xiEQEqBjAvgwgdLjcP/T4jEMkYhN5DGBfBhBaXO4fWnzGIRIRoGIAxzAAoRfpl/vHEj7jEIlYwosB7MsAEMv4jEMkIkDFCPZlAIhlfMYh0rCEF0PYlwEglvEZh0hCgIox7MsAEMv4jEOkIEABIcIZXgAQOwhQQAhwhhcAxBY2kQNBtm4bZ3gBQKwhQAFB5HAa+ufHhW7bcIYXAEQfAhQQRFu+P6wyzvACgJhDgAKCqOxItUftOMMLAKKL6TaRV1dX68svv9SmTZtUUFCggoIC2e12SdL69euVmpoa3g4iptjS2njUjjO8ACC6mC5A7dq1S5MnTw53NxAl/C090OesDrJZk90u43GGFwBEH9MFqISEBPXv31/9+vVTTk6OOnTooAkTJoS7W4hAgSg9EB9n0U1XZmvmOxtbbMMZXgAQfSyGYZj68p+9e/fq8ssvlxScJbzy8qOqr3cG9DkRfPmFJ0oPtMSTs7cSEuKUkZGq8vKjWrt5f5MwZrMmaxx1oGLayWOAzwFzYgxEH5stVfHxrW8RN90MFNAap9PQwuU73LZ5Y/kODcjs5PHMEWd4AUBsIUABp9hebG9S9PJUrtID3pzJxRleABA7KGMAnMLTkgKUHgAA82IGKsg8WUdFZOngYemBDmltlJDQ8t+v6++eMWBejAEwBmIXASrI0tJSwt0FeGlw+7bq8N4WHa5ouQhmx/QUDT7vdMV7sIeJMQDGABgDsSdqAtTDDz+st956y+vHDRo0SPPnzw9Cjzxz5MhxORxceRFtxv8sq9XSA0cqjrl9jvj4OKWlpTAGTIwxAMZA9ElLS4mtq/CsVqs6duzo9ePat28fhN54zuFwculqFBpwTkdNHp3TYumBAed09PjvlTEAxgAYA7GHOlDUgYIb/lQip/4LGANgDEQf6kABARCrpQf8PaIGQHjxHg4/AhRgMoE4ogY/4osMocZ7ODIQoAATaemImvLKGr24uMCjI2rwI77IEGq8hyOHKQPU6NGjVVJSIkk6eQvYZZdd1vDvAwcO1KxZs0LeNyBYgnFEjZl58kU2uG/XMPQMsYr3cGQxZWWviooK2e122e12VVRUNNzuus1ut6uqqiqMPQQCz5sjauCep19kTqepr9FBgPEejiymnIH69NNPw90FIOQ4oiZwPP0iK9xTrmEd2oWoV4h1vIcjiylnoAAzSk9NDmg7M/P4i6yqNsg9gZnwHo4sBCjAJLJ6pCvD6v6D1WY9cRUZ3PP4i6xdUpB7AjPhPRxZCFCAScTFWTR+RKbbNuNGZLL51AOefpFlnxF7NcQQPryHIwsBCjCR3OzOmjw6p8mXv82azOXPXuCLDOHCezhymP4ol2CjfL95RfIRDhR/DIzm6kC5zkvMze4c0WMAoRGsMcB7OHg4ygVAi2L1iJpQy83urAGZnfgiQ8jxHg4/AhQQJE6noU07D6l4X4WsKYl8scYovsgAcyJAAUGQX1iqN5bvUBlHfABATGITORBgriM+yk4ptOg64iO/sDRMPQMABAoBCgggjvgAAHMgQAEBxFlVAGAOBCgggDirCgDMgU3kQAAF46wq6r0AQOQhQAEB5Driw90ynjdnVTVXqJGr+QAg/FjCMzmn09C23eVas2W/tu0uZ3OznwJ5xIfrar5TwxhX8wFA+DEDZWLMbgSH66yqU+tAnXzER2s8vZpvQGYnlvMAIAwIUCblmt04lWt2g0Mp/ZOb3VkXnNtFJeXVPlUi9+ZqPn+qYLO/CgB8Q4AyIWY3QiMuzqJ+53TU6R1SvD5ENBRX8zEDGVqEVSC2EKBMKFSzG/BdMK7mOxkzkKG1blupFnxcSFgFYgibyE2IWkWRz3U1nzveXM13Mqqlh9aqjSWa+c5GLgYAYgwByoSCPbsB/wXyar5TUS09dJxOQ68s2eS2DWEViE4EKBMK5uwGAsd1Nd+pf1c2a7JfS2zMQIZO4Z5yHa6odtuGsApEJ/ZAmZBrdqO5PTAuvs5uILBysztrQGangG4+ZgYydOxVtZ61I6wCUYcAZVKu2Y1Tr8LyplYRQiMuzhLQzfyBrpaOlqW3S/KsHWEViDoEKBMLxuwGIh8zkKGTfUaGOrRv43YZj7AKRCf2QJmca3ZjSJ+u6n1mBl+aJhGs/VVoLC7Oogl5/dy2IawC0cliGAaXfwRReflRr4soIjYkJMQpIyM1oscAxR2DyzUG/r16V5M6UCyXm0M0fA6gMZstVfHxrc8vsYQHmFig91eheRf07qzzzupAWAViCAEKAEKAsArEFvZAAQAAeIkABQAA4CUCFAAAgJfYAxUDuJIKAIDQIkBFufzC0ibVxDOsyRrP5dEAAAQNS3hRLL+wVC8uLmhyJEd5ZY1eXFyg/MLSMPUMAIDYRoCKUk6noYXLd7ht88byHXI6qZMKAECgEaCi1PZiu9vDYCWprLJG24vtoekQAAAmwh6oKGU/6j48edsumrGJHgAQagSoKJWemtx6Iy/aRSs20QMAwoElvCiV1SNdGVb34chmPTEbE6vYRA8ACBcCVJSKi7No/IhMt23GjciM2aUsNtEDAMKJABXFcrM7a/LonCYzUTZrsiaPzonpJSw20QMAwok9UFEuN7uzBmR2Mt0majbRAwDCiQAVA+LiLOp9Zka4uxFSbKIHAIQTS3iISmyiBwCEEwEKUcnsm+gBAOFFgELUMvMmegBAeLEHClHNrJvoAQDhRYBC1DPjJnoAQHixhAcAAOAl081AFRUV6ZNPPtHatWtVWFio8vJyJScnq1evXrriiit00003KTU1NdzdBAAAEcxiGIZpzrpwOBzq06dPw587deqkLl266NChQ9q/f78k6YwzztDcuXPVvXv3gLxmeflR1dc7A/JciC4JCXHKyEhlDJgYYwCMgehjs6UqPr71BTpTzUAZhqF27dpp7Nixuvbaa3X22Wc33Ldhwwbdd9992rNnj+699169/fbbYewpAACIZKaagTIMQxUVFUpPT2/2/vXr12vcuHGSpCVLlujcc8/1+zX51WFe/PIEYwCMgejj6QyUqTaRWyyWFsOTJA0cOFBWq1WStGvXrhD1CgAARBtTBajWOBwO1dfXS5LatGkT5t4AAIBIRYA6yYoVK3T8+HElJCTo/PPPD3d3AABAhDLVJnJ3qqqq9OSTT0qSxowZI5vNFpDn9WQdFbHJ9XfPGAgOp9NQ4Z5y2atqld4uSdlnZERcBXrGABgDsctUm8hb4nA4NHHiRH3++efq3r27li5d2rAXCkDkWbWxRK8s2aTDFdUNt3Vo30YT8vppWP9uYewZALOImgD18MMP66233vL6cYMGDdL8+fNbvN8wDD344IN699131b59ey1YsEBZWVn+dLWRI0eOy+Hgygszio+PU1paCmMgwNZtK9XMdza2eP891/XXBb0j4yBpxgAYA9EnLS0ltupAWa1WdezY0evHtW/f3u39f/7zn/Xuu+8qNTVVr776akDDkyQ5HE4uXTU5xkDgOJ2GFnxc6LbNPz8u1HlndYio5TzGABgDsSdqAtTUqVM1derUgD7nk08+qQULFiglJUWvvPKK+vfvH9DnDzSn09D2YrvsR2uUnpqsrB7pEfUlAQTb9mK7yitr3LYpq6zR9mI7B0wDCKqoCVCBNmPGDM2ZM0fJycmaNWuWfvKTn4S7S27lF5Zq4fIdjb48MqzJGj8iU7nZkbFcAQSb/aj78ORtOwDwlSkvC5g9e7Zmz56txMREzZw5U0OHDg13l9zKLyzVi4sLmvzyLq+s0YuLC5RfWBqmngGhlZ6aHNB2AOAr0wWoefPmacaMGUpISNCMGTN0ySWXhLtLbjmdhhYu3+G2zRvLd8jpjIprAQC/ZPVIV4bVfTiyWU8sbwNAMJlqCe/AgQN67LHHJEmpqamaM2eO5syZ02zbMWPG6Lrrrgtl95rFng/gR3FxFo0fkakXFxe02GbciEz2BgIIOlMFqLq6OrmqNlRUVGj9+vUtth02bFiouuUWez6AxnKzO2vy6JwmewJt1mSNY08ggBAxVYA6/fTTVVjo/hLoSMOeD6Cp3OzOGpDZiatSAYSNqQJUNHLt+XC3jMeeD5hRXJyFZWsAYWO6TeTRxrXnwx32fAAAEFoEqCjg2vNx6tVHNmuyJo/OYc8HAAAhxhJelGDPBwAAkYMAFUXY8wEAQGRgCQ8AAMBLBCgAAAAvEaAAAAC8RIACAADwEgEKAADASwQoAAAALxGgAAAAvESAAgAA8BIBCgAAwEsWwzCMcHciljkcznB3AWEUHx/HGDA5xgAYA9ElLs4ii6X1Y9IIUAAAAF5iCQ8AAMBLBCgAAAAvEaAAAAC8RIACAADwEgEKAADASwQoAAAALxGgAAAAvESAAgAA8BIBCgAAwEsEKAAAAC8RoAAAALxEgAIAAPASAQoAAMBLBCgAAAAvJYS7A2ZRVFSkTz75RGvXrlVhYaHKy8uVnJysXr166YorrtBNN92k1NTUcHcTQVZdXa0vv/xSmzZtUkFBgQoKCmS32yVJ69evZwzEiDVr1ui1117Tt99+q2PHjqlbt24aOXKkJkyYoLZt24a7ewiygwcPatWqVQ3v861bt6q6ulp9+/bVokWLwt09BIjFMAwj3J2IdQ6HQ3369Gn4c6dOndSlSxcdOnRI+/fvlySdccYZmjt3rrp37x6ubiIEtm7dqry8vGbvI0DFhvnz5+svf/mLDMNQ165dZbPZtHPnTtXW1urss8/WwoULlZ6eHu5uIojmzp2rxx9/vMntBKjYwgxUCBiGoXbt2mns2LG69tprdfbZZzfct2HDBt13333as2eP7r33Xr399tth7CmCLSEhQf3791e/fv2Uk5OjDh06aMKECeHuFgKkoKBAjz32mCTpkUce0Q033CCLxaIDBw5o4sSJ2rx5s6ZPn66ZM2eGuacIpnbt2mnYsGHKyclRTk6OioqK9Oyzz4a7WwgwZqBCwDAMVVRUtPirc/369Ro3bpwkacmSJTr33HND2DuE0969e3X55ZdLYgYqFkyaNEkrVqxQXl6ennzyyUb3FRUV6aqrrpLT6dTSpUvVu3fvMPUSobZo0SI98MADzEDFGDaRh4DFYnE7ZT9w4EBZrVZJ0q5du0LUKwCBdPToUX3xxReSpBtuuKHJ/T179tSQIUMkSR999FFI+wYg8AhQEcDhcKi+vl6S1KZNmzD3BoAvtm7dqtraWiUlJal///7NtsnNzZUkffvtt6HsGoAgIEBFgBUrVuj48eNKSEjQ+eefH+7uAPCBa/a4W7duSkxMbLbNGWec0agtgOhFgAqzqqqqhr0SY8aMkc1mC3OPAPiioqJCktS+ffsW27juc7UFEL0IUGHkcDg0ZcoU7d27V927d9fUqVPD3SUAPqqpqZGkFmefJCkpKalRWwDRizIGrXj44Yf11ltvef24QYMGaf78+S3ebxiGpk+frs8//1zt27fX7NmzGzaSI/IEaxwgdiQnJ0uS6urqWmxTW1vbqC2A6EWAaoXValXHjh29fpy7aXxJ+vOf/6x3331XqampevXVV5WVleVrFxECwRoHiB2eLM95sswHIDoQoFoxderUgC+tPfnkk1qwYIFSUlL0yiuvtHjFDiJHMMYBYkvPnj0lSSUlJaqrq2t2KW/Pnj2N2gKIXuyBCrEZM2Zozpw5Sk5O1qxZs/STn/wk3F0CEAB9+vRRYmKiamtrtXHjxmbb5OfnSxJX2wIxgAAVQrNnz9bs2bOVmJiomTNnaujQoeHuEoAASU1N1fDhwyWp2SOZioqKtGbNGknSyJEjQ9o3AIFHgAqRefPmacaMGUpISNCMGTN0ySWXhLtLAAJs0qRJslgsWrp0qd566y25TsoqLS3VlClT5HQ6NWLECI5xAWIAZ+GFwIEDB3TJJZfIMAy1b9++0WHCpxozZoyuu+66EPYOoTZ69GiVlJRI+vGcREmNjvsZOHCgZs2aFY7uwU9z587VE088IcMwdNpppykjI0M7d+5UbW2tevXqpYULF1LvLcbt27dPeXl5DX+ura3VsWPHlJCQoHbt2jXcfuedd+quu+4KQw8RCGwiD4G6urqGX6IVFRVav359i22HDRsWqm4hTCoqKmS325vcfvJtVVVVoesQAuq2225Tdna25syZo40bN+rw4cPq1q2bRo4cqQkTJnBgtAk4HI5m3+P19fWNbq+urg5dpxBwzEABAAB4iT1QAAAAXiJAAQAAeIkABQAA4CUCFAAAgJcIUAAAAF4iQAEAAHiJAAUAAOAlAhQAAICXCFAAAABeIkABAAB4ibPwACDKHThwQMuXL9fatWu1detWHThwQJLUsWNHnX/++br++us1dOjQMPcSiC2chQcAUWzfvn366U9/qpM/ylNSUmQYRqPDaseMGaNHH31U8fHx4egmEHOYgQKAKOZwOGQYhoYOHaq8vDwNHTpUXbp0kdPp1Pfff69nn31WK1as0LvvvqvOnTvr//2//xfuLgMxgRkoAIhilZWV2rNnj/r27dvs/YZh6K677tIXX3yhtm3bas2aNUpOTg5xL4HYwyZyAGF18OBBPf300/rlL3+p3Nxc9evXT5dffrkefPBB7dy5s0n7V155RdnZ2crJydHGjRubfc7PP/9cvXv3VnZ2tt57770m969du1a/+c1vdNFFFyknJ0eDBw/WrbfeqnfffVcOh8Ntf9etW6f//d//1eDBg9W/f39deeWVmjFjho4ePapFixYpOztbl112mVf/Dz744ANlZ2crOztbn3zySbNttm7dqn79+ik7O1uzZ89uuN1qtbYYniTJYrFozJgxkqRjx47pu+++86pvAJpHgAIQNitXrtQVV1yhv//97yosLFRNTY0SEhK0d+9evfPOO8rLy9OSJUsaPeauu+7SsGHDVFdXpylTpqiqqqrR/aWlpfrd734nwzCUl5enq6++utH9jz/+uG655RZ9/PHHOnjwoFJSUlRZWak1a9bo97//vW6//fYmz+kyf/583XzzzVq5cqXsdruSkpL0ww8/aPbs2br++ut15MgRn/4/jBo1qiHkPPTQQ9q3b1+j+48dO6Z7771XtbW1GjRokCZMmODV858849RaQATgGQIUgLDYuHGj7rnnHh07dkw33nijPvjgA3377bf65ptvtHLlSo0fP151dXV68MEHtWnTpobHWSwWPfXUU+rQoYOKi4v1hz/8oeE+wzB0//33q6ysTGeeeaYefvjhRq+5YMECzZ07V5J044036osvvtC6dev03//+Vw888IASEhK0Zs0aTZ8+vUl/169fr8cee0yGYejCCy/URx99pP/+97/asGGDnn/+eR0+fFgvvfSSz/8/pk+frrPOOkt2u11Tp05tFHQeffRR7dq1S+np6Xr66acVF+fdR/fXX38tSUpMTFSvXr187iOAkxgAEAZjxowxsrKyjOeee67FNo8++qiRlZVlTJw4scl9n332mZGdnW1kZWUZixYtMgzDMF5++WUjKyvL6Nu3r7Fx48ZG7Y8fP24MGjTIyMrKMqZMmdLs682bN8/IysoysrKymjz+1ltvNbKysoxRo0YZNTU1TR67evXqhsf+9Kc/bfW/vzlbt241cnJyjKysLGPmzJmGYRjGsmXLGp73k08+8fo59+zZY5x33nlGVlaWMXXqVJ/6BaApZqAAhNy2bdu0adMmJSYm6vbbb2+xXV5eniRp9erVTZaeLrnkEt12222SpEceeURLly7V3/72N0nSvffeq379+jVq/9VXX8lut0uSfv3rXzf7euPHj1enTp0kScuWLWu43W63a82aNZKkO+64Q0lJSU0eO2TIEP3kJz9p8b/FE71799a0adMkSS+99JKWLFnSMIt20003acSIEV49X3V1tf7v//5Px48fV3p6uu677z6/+gfgR5QxABBy+fn5kiSn06mRI0e22M4Vmo4dOya73a4OHTo0un/KlCn6+uuvtXnz5obgMXz48GZDWUFBgSTptNNOa3EZKz4+XkOGDNF7773X0F46sYHb+P8vWL7gggta7O+gQYP03//+t8X7PXHzzTfrq6++0sqVK3X//fdLkrKyshr+3VP19fW67777tHnzZiUmJuqZZ55Rly5d/OobgB8xAwUg5EpLSyWdCEiHDh1q8Z/y8vKGxxw/frzJ8yQlJemJJ55o+LPVatUTTzwhi8XSpO3hw4clqdUQ0bVr10btJamsrKzh3909PlAB5bHHHmvY+B0fH69nn33Wq9IDDodDU6dO1fLly5WQkKCnn35aw4cPD0jfAJzADBSAkHPNLJ111ln68MMP/Xqut99+u+Hfq6qqtHXr1oZluOY0F65aa2d4WC7P03atWbp0qWpqaiSd+H+Vn5+vzMxMjx7rCk8ffPCB4uPj9de//tXtLB8A3zADBSDkOnbsKEkqLi7WsWPHfH6elStXav78+ZKk7OxsGYah3/3udzp06FCTtq7lv1NLBJxq//79kiSbzdbksdKPs2fNcXefpzZv3qxnnnlG0on/JulE6YXmamKdyuFw6Le//a2WLVvWEJ5GjRrld58ANEWAAhByAwcOlCTV1dW1WDiyNaWlpXrggQckSddee63++c9/qnv37jp8+LDuv//+JrNBOTk5kk4EpF27djX7nA6HQ2vXrpWkRpvQzz333IYZKVdJgOa4u88Tx44d05QpU1RXV6chQ4bonXfeUf/+/VVdXa0pU6aotra2xce6wtPJM08///nP/eoPgJYRoACEXL9+/dSnTx9J0nPPPddoj1FzXFfPuTidTk2bNk3l5eXq2bOnpk+fLqvVqmeeeUYJCQn68ssv9dprrzV6zIUXXqj09HRJ0gsvvNDs67z55psNs0gnh4/09HQNHjxYkvTaa681G2Rc9aT88eijj6qoqEjp6el66qmnlJSUpGeeeUapqakqLCzUk08+2ezjHA6H7rvvPn3wwQcNe54IT0BwEaAAhJzFYtGf/vQnJSUlqaSkRNdff70++uijRhvFDxw4oKVLl+p//ud/9PTTTzd6/KuvvqrVq1c3XF3Wtm1bSdKAAQM0efJkSdKzzz6rzZs3NzymTZs2uueeeyRJ77//vh5++OGGpb7jx49r/vz5evzxxyWdqAzumrFyueeee2SxWLR9+3ZNnDhRRUVFkk5c7fbvf/9b99xzj9q3b+/z/5Nly5Zp0aJFkk5sIndtSD/jjDMaioUuWLBAK1eubPQ4h8OhadOm6cMPP2wITyzbAcHHYcIAwuarr77SlClTGmaY4uPjZbVaVVNT0yhMXX/99frzn/8s6UQFc1eV8mnTpumOO+5o9JxOp1O33nqrvv76a/Xs2VOLFy9uCFjSif1ErmrkFotFaWlpOnr0qOrr6yVJgwcP1ksvvaR27do16e/cuXMbQpYkpaWlqbq6WrW1tcrKytKYMWP0+OOPq1evXvroo488/v+wd+9e5eXlqbKyUjfddFOTCuqSNG3aNC1dulQZGRlaunRpQ8Bat26dfvWrX0k6UWm8tRD34IMPErCAACBAAQirI0eO6M0339Rnn32m7777TpWVlUpOTla3bt10/vnn6/LLL9eFF16o5ORkVVVVKS8vT8XFxbrwwgv1j3/8o9mr6vbv369rrrlGdrtd1157baPQI0lr1qzRwoULtX79etntdrVt21a9e/fWNddco7y8PMXHx7fY36+//lqvvvqqNmzYoOrqanXr1k0jR47UhAkT9Pbbb+vxxx/Xeeed1+jqQHfq6+t10003acOGDcrKytI777zTbMmCqqoqXXvttdq9e7eGDh2qOXPmKC4uTmvXrtUtt9zi0WtJJwLktdde63F7AM0jQAFAgNx33316//33NWbMGD322GPh7g6AIGIPFAAEwK5duxquKLzooovC3BsAwUaAAgAPPf/881qwYIFKSkrkdDolnSg98MEHH+iWW25RTU2NzjrrLK/PrAMQfVjCAwAPTZo0SStWrJB0YsN2amqqjhw50hCmunTpoldffVVZWVnh7CaAEOAoFwDw0G233abOnTvrm2++0cGDB1VRUaHU1FT17NlTl156qX71q1811JoCENuYgQIAAPASe6AAAAC8RIACAADwEgEKAADASwQoAAAALxGgAAAAvESAAgAA8BIBCgAAwEsEKAAAAC8RoAAAALz0/wFluKBfoH8kOgAAAABJRU5ErkJggg==",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure()\n",
"plt.scatter(x1, y)\n",
"plt.xlabel('exog x1')\n",
"plt.ylabel('endog y')\n",
"\n",
"plt.figure()\n",
"plt.scatter(x2, y)\n",
"plt.xlabel('exog x2')\n",
"plt.ylabel('endog y')\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "d30de122-7250-4c31-9e6d-f2e22370b242",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y R-squared: 0.484\n",
"Model: OLS Adj. R-squared: 0.462\n",
"Method: Least Squares F-statistic: 22.01\n",
"Date: Thu, 06 Jul 2023 Prob (F-statistic): 1.80e-07\n",
"Time: 00:36:28 Log-Likelihood: -62.306\n",
"No. Observations: 50 AIC: 130.6\n",
"Df Residuals: 47 BIC: 136.3\n",
"Df Model: 2 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept -0.0413 0.125 -0.332 0.741 -0.292 0.209\n",
"x1 -0.4097 0.118 -3.468 0.001 -0.647 -0.172\n",
"x2 0.5906 0.136 4.335 0.000 0.316 0.865\n",
"==============================================================================\n",
"Omnibus: 1.032 Durbin-Watson: 2.003\n",
"Prob(Omnibus): 0.597 Jarque-Bera (JB): 0.986\n",
"Skew: 0.167 Prob(JB): 0.611\n",
"Kurtosis: 2.398 Cond. No. 1.47\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"source": [
"# define linear model\n",
"lm = smf.ols('y ~ x1 + x2', df)\n",
"\n",
"# fit model to data\n",
"lmf = lm.fit()\n",
"\n",
"# summary\n",
"print(lmf.summary())"
]
},
{
"cell_type": "markdown",
"id": "776bee7d-176d-443d-90ff-256ab36bb5e2",
"metadata": {},
"source": [
"The construction of the model relies on the design matrix created by `patsy.dmatrix`."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "9237dd64-3e28-4a6e-8a26-fb8fa63d9842",
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Intercept</th>\n",
" <th>x1</th>\n",
" <th>x2</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1.0</td>\n",
" <td>-0.598903</td>\n",
" <td>-1.107940</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1.0</td>\n",
" <td>1.484895</td>\n",
" <td>-0.308404</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1.0</td>\n",
" <td>-1.279121</td>\n",
" <td>-0.071924</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>1.0</td>\n",
" <td>0.796557</td>\n",
" <td>-1.757192</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1.0</td>\n",
" <td>0.772527</td>\n",
" <td>0.049753</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Intercept x1 x2\n",
"0 1.0 -0.598903 -1.107940\n",
"1 1.0 1.484895 -0.308404\n",
"2 1.0 -1.279121 -0.071924\n",
"3 1.0 0.796557 -1.757192\n",
"4 1.0 0.772527 0.049753"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"y_aff, X_aff = dmatrices('y ~ x1 + x2', data=df, return_type='dataframe')\n",
"X_aff.head()"
]
},
{
"cell_type": "markdown",
"id": "b9596e22-c8b4-43ea-88ca-894c1cb4bbff",
"metadata": {},
"source": [
"## Interactions between predictors\n",
"\n",
"Now we change the model to incorporate a non-linear interaction between $x1$ and $x2$."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "0ae0753d-9007-489a-a102-5a1a1f1b97c9",
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>x1</th>\n",
" <th>x2</th>\n",
" <th>y</th>\n",
" <th>y2</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>-0.598903</td>\n",
" <td>-1.107940</td>\n",
" <td>-0.983083</td>\n",
" <td>-0.717664</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1.484895</td>\n",
" <td>-0.308404</td>\n",
" <td>-2.069647</td>\n",
" <td>-2.252826</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>-1.279121</td>\n",
" <td>-0.071924</td>\n",
" <td>0.927597</td>\n",
" <td>0.964397</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>0.796557</td>\n",
" <td>-1.757192</td>\n",
" <td>-0.702349</td>\n",
" <td>-1.262231</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>0.772527</td>\n",
" <td>0.049753</td>\n",
" <td>0.529516</td>\n",
" <td>0.544891</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" x1 x2 y y2\n",
"0 -0.598903 -1.107940 -0.983083 -0.717664\n",
"1 1.484895 -0.308404 -2.069647 -2.252826\n",
"2 -1.279121 -0.071924 0.927597 0.964397\n",
"3 0.796557 -1.757192 -0.702349 -1.262231\n",
"4 0.772527 0.049753 0.529516 0.544891"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# coefficient for the interaction, which is simply a second-order polynomial in x1, x2\n",
"a12 = 0.4\n",
"y2 = a1 * x1 + a2 * x2 + a12 * x1 * x2 + e\n",
"\n",
"# build dataframe\n",
"df['y2'] = y2\n",
"df.head()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "fc12371e-4da7-4497-aebb-a55a7aae0daf",
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Intercept</th>\n",
" <th>x1</th>\n",
" <th>x2</th>\n",
" <th>x1:x2</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1.0</td>\n",
" <td>-0.598903</td>\n",
" <td>-1.107940</td>\n",
" <td>0.663548</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1.0</td>\n",
" <td>1.484895</td>\n",
" <td>-0.308404</td>\n",
" <td>-0.457947</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1.0</td>\n",
" <td>-1.279121</td>\n",
" <td>-0.071924</td>\n",
" <td>0.092000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>1.0</td>\n",
" <td>0.796557</td>\n",
" <td>-1.757192</td>\n",
" <td>-1.399704</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1.0</td>\n",
" <td>0.772527</td>\n",
" <td>0.049753</td>\n",
" <td>0.038436</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Intercept x1 x2 x1:x2\n",
"0 1.0 -0.598903 -1.107940 0.663548\n",
"1 1.0 1.484895 -0.308404 -0.457947\n",
"2 1.0 -1.279121 -0.071924 0.092000\n",
"3 1.0 0.796557 -1.757192 -1.399704\n",
"4 1.0 0.772527 0.049753 0.038436"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"y_aff, X_aff = dmatrices('y2 ~ x1 * x2', data=df, return_type='dataframe')\n",
"X_aff.head()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "8c0ae1bb-4493-4a1f-bb2a-d467d216fc62",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y2 R-squared: 0.477\n",
"Model: OLS Adj. R-squared: 0.443\n",
"Method: Least Squares F-statistic: 14.01\n",
"Date: Thu, 06 Jul 2023 Prob (F-statistic): 1.28e-06\n",
"Time: 00:36:41 Log-Likelihood: -62.280\n",
"No. Observations: 50 AIC: 132.6\n",
"Df Residuals: 46 BIC: 140.2\n",
"Df Model: 3 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept -0.0335 0.131 -0.256 0.799 -0.297 0.230\n",
"x1 -0.4150 0.122 -3.405 0.001 -0.660 -0.170\n",
"x2 0.5960 0.140 4.259 0.000 0.314 0.878\n",
"x1:x2 0.4273 0.126 3.389 0.001 0.174 0.681\n",
"==============================================================================\n",
"Omnibus: 0.909 Durbin-Watson: 1.989\n",
"Prob(Omnibus): 0.635 Jarque-Bera (JB): 0.917\n",
"Skew: 0.164 Prob(JB): 0.632\n",
"Kurtosis: 2.423 Cond. No. 1.71\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"source": [
"# define linear model\n",
"lm2 = smf.ols('y2 ~ x1 * x2', df)\n",
"\n",
"# fit model to data\n",
"lmf2 = lm2.fit()\n",
"\n",
"# summary\n",
"print(lmf2.summary())"
]
},
{
"cell_type": "markdown",
"id": "3761c74f-1dc0-47a0-aeb3-ccb928875dad",
"metadata": {},
"source": [
"This model accurately captures the coefficients used to generate $y2$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "dfc36822-9af7-4eab-9e25-baaa0b9b6eb9",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y2 R-squared: 0.347\n",
"Model: OLS Adj. R-squared: 0.319\n",
"Method: Least Squares F-statistic: 12.48\n",
"Date: Thu, 06 Jul 2023 Prob (F-statistic): 4.48e-05\n",
"Time: 00:36:52 Log-Likelihood: -67.852\n",
"No. Observations: 50 AIC: 141.7\n",
"Df Residuals: 47 BIC: 147.4\n",
"Df Model: 2 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept -0.1561 0.139 -1.121 0.268 -0.436 0.124\n",
"x1 -0.3311 0.132 -2.509 0.016 -0.597 -0.066\n",
"x2 0.5106 0.152 3.354 0.002 0.204 0.817\n",
"==============================================================================\n",
"Omnibus: 0.529 Durbin-Watson: 2.231\n",
"Prob(Omnibus): 0.768 Jarque-Bera (JB): 0.634\n",
"Skew: -0.016 Prob(JB): 0.728\n",
"Kurtosis: 2.449 Cond. No. 1.47\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"source": [
"# model without interaction\n",
"lm = smf.ols('y2 ~ x1 + x2', df)\n",
"\n",
"# fit model to data\n",
"lmf = lm.fit()\n",
"\n",
"# summary\n",
"print(lmf.summary())"
]
},
{
"cell_type": "markdown",
"id": "edf8ac33-7145-4bed-b9d9-c0362104debc",
"metadata": {},
"source": [
"Note that correctly estimating the interaction requires including the linear dependencies."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "db858a08-c017-4c7b-9f2a-24d9fbe5dfbd",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y2 R-squared: 0.026\n",
"Model: OLS Adj. R-squared: 0.005\n",
"Method: Least Squares F-statistic: 1.257\n",
"Date: Thu, 06 Jul 2023 Prob (F-statistic): 0.268\n",
"Time: 00:37:28 Log-Likelihood: -77.858\n",
"No. Observations: 50 AIC: 159.7\n",
"Df Residuals: 48 BIC: 163.5\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept -0.0803 0.173 -0.465 0.644 -0.428 0.267\n",
"x1:x2 0.1794 0.160 1.121 0.268 -0.142 0.501\n",
"==============================================================================\n",
"Omnibus: 0.414 Durbin-Watson: 2.034\n",
"Prob(Omnibus): 0.813 Jarque-Bera (JB): 0.567\n",
"Skew: 0.064 Prob(JB): 0.753\n",
"Kurtosis: 2.494 Cond. No. 1.35\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"source": [
"# model with only the interaction\n",
"lm3 = smf.ols('y2 ~ x1 : x2', df)\n",
"\n",
"# fit model to data\n",
"lmf3 = lm3.fit()\n",
"\n",
"# summary\n",
"print(lmf3.summary())"
]
},
{
"cell_type": "markdown",
"id": "793b95f7-1df7-4b58-858d-309f2ab93717",
"metadata": {},
"source": [
"### EXERCISE\n",
"- Redo the fit of $y$ by the model with interactions (`lm2`)\n",
"- with the below data, check relationships and interactions between the cognitive scores.\n",
"- Adapt the synthetic model to the data and compare the perf with the model fitted on data.\n",
"\n",
"Dataset taken from R dataset repository ([link](https://github.com/vincentarelbundock/Rdatasets/)), use e.g.\n",
"`df = sm.datasets.get_rdataset(\"NeuroCog\", \"heplots\").data`"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "e1625945-157d-4716-a6fd-5f287a4b17d4",
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Unnamed: 0</th>\n",
" <th>Dx</th>\n",
" <th>Speed</th>\n",
" <th>Attention</th>\n",
" <th>Memory</th>\n",
" <th>Verbal</th>\n",
" <th>Visual</th>\n",
" <th>ProbSolv</th>\n",
" <th>SocialCog</th>\n",
" <th>Age</th>\n",
" <th>Sex</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>14</td>\n",
" <td>Schizophrenia</td>\n",
" <td>19</td>\n",
" <td>9</td>\n",
" <td>19</td>\n",
" <td>33</td>\n",
" <td>24</td>\n",
" <td>39</td>\n",
" <td>28</td>\n",
" <td>44</td>\n",
" <td>Female</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>15</td>\n",
" <td>Schizophrenia</td>\n",
" <td>8</td>\n",
" <td>25</td>\n",
" <td>15</td>\n",
" <td>28</td>\n",
" <td>24</td>\n",
" <td>40</td>\n",
" <td>37</td>\n",
" <td>26</td>\n",
" <td>Male</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>16</td>\n",
" <td>Schizophrenia</td>\n",
" <td>14</td>\n",
" <td>23</td>\n",
" <td>15</td>\n",
" <td>20</td>\n",
" <td>13</td>\n",
" <td>32</td>\n",
" <td>24</td>\n",
" <td>55</td>\n",
" <td>Female</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>17</td>\n",
" <td>Schizophrenia</td>\n",
" <td>7</td>\n",
" <td>18</td>\n",
" <td>14</td>\n",
" <td>34</td>\n",
" <td>16</td>\n",
" <td>31</td>\n",
" <td>36</td>\n",
" <td>53</td>\n",
" <td>Male</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>18</td>\n",
" <td>Schizophrenia</td>\n",
" <td>21</td>\n",
" <td>9</td>\n",
" <td>35</td>\n",
" <td>28</td>\n",
" <td>29</td>\n",
" <td>45</td>\n",
" <td>28</td>\n",
" <td>51</td>\n",
" <td>Male</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Unnamed: 0 Dx Speed Attention Memory Verbal Visual \n",
"0 14 Schizophrenia 19 9 19 33 24 \\\n",
"1 15 Schizophrenia 8 25 15 28 24 \n",
"2 16 Schizophrenia 14 23 15 20 13 \n",
"3 17 Schizophrenia 7 18 14 34 16 \n",
"4 18 Schizophrenia 21 9 35 28 29 \n",
"\n",
" ProbSolv SocialCog Age Sex \n",
"0 39 28 44 Female \n",
"1 40 37 26 Male \n",
"2 32 24 55 Female \n",
"3 31 36 53 Male \n",
"4 45 28 51 Male "
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# load NeuroCog data \n",
"df = pd.read_csv('NeuroCog_dataset.csv', sep=',')\n",
"\n",
"df.head()"
]
},
{
"cell_type": "markdown",
"id": "04a9875a-c70f-4d36-8046-d5651092109a",
"metadata": {},
"source": [
"## "
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}