Example_1_Exploration.ipynb 46.4 KB
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
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploration\n",
    "The first step in the `CORBASS` algorithm is an exploration of the parameter space. This way the region of interest, i.e. the region where the probability mass is concentrated, is found. Here we give an example how to run the `corbass.exploration` module and how to interpret the output. First we import the function `explore` from the module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import sys\n",
    "import os\n",
    "# relative import\n",
    "sys.path.append(os.path.abspath('') + '/../')\n",
    "from corbass.exploration import explore"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After the import we load the example parameter file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from corbass.utils import load\n",
    "pars = load('./Example_Parfile.py')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
46
    "Next we use the `explore` function to explore the parameter space. This is done via grid evaluation. The model parameter bounds and the number of gridpoints are specified in the parameter file. Then for each point in the grid, i.e. each choice of model parameters, the posterior probability is calculated. `explore` returns an instance of `ExplorationResult`, containing the posterior on the grid and a scale. We use the `write` member function to write the results to the disk, to use them in [the next example](https://gitext.gfz-potsdam.de/arthus/corbass/blob/master/examples/Example_1_Exploration.ipynb)."
47
48
49
50
51
52
53
54
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
Maximilian Schanner's avatar
Typo    
Maximilian Schanner committed
55
   "outputs": [
56
57
58
59
60
61
62
63
64
65
66
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exploring 1450, this may take a while...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
Maximilian Schanner's avatar
Maximilian Schanner committed
67
      "100% (15625 of 15625) |##################| Elapsed Time: 0:43:00 Time:  0:43:00\n"
Maximilian Schanner's avatar
typo    
Maximilian Schanner committed
68
     ]
69
70
71
72
73
74
75
76
77
78
79
80
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exploring 1550, this may take a while...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
Maximilian Schanner's avatar
Maximilian Schanner committed
81
      "100% (15625 of 15625) |##################| Elapsed Time: 0:58:53 Time:  0:58:53\n"
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
     ]
    }
   ],
   "source": [
    "for bin_name, data in pars.data:\n",
    "    print(f\"Exploring {bin_name}, this may take a while...\")\n",
    "    rslt = explore(pars.bin_fname(bin_name), data, pars.n_ex, pars.bounds,\n",
    "                   r_ref=pars.r_ref)\n",
    "    rslt.write()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can load the results for one of the bins and analyse them. The standard output name can be accessed from the `ExplorationResult` class, the output location is specified in the parameter file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from corbass.exploration import ExplorationResult\n",
    "import numpy as np\n",
    "\n",
    "with np.load(f'{pars.bin_fname(1450)}{ExplorationResult.suffix}') as fh:\n",
    "    posterior = fh['posterior']\n",
    "    scale = fh['scale']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
118
    "To demonstrate how this result is used, we focus on the marginal posterior for the error level and the residual scaling. First we get the corresponding parameter space arrays and build a meshgrid:"
119
120
121
122
123
124
125
126
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
127
    "from corbass.utils import bounds2grid\n",
128
129
    "lams, errs, ress = bounds2grid(pars.bounds, pars.n_ex)\n",
    "ER, RE = np.meshgrid(errs[0], ress[0])"
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we integrate out the non-dipole scaling, by means of a simple Riemann sum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
145
    "marginal = posterior.sum(axis=0)*lams[1]"
146
147
148
149
150
151
152
153
154
155
156
157
158
159
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot the marginal distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
160
161
162
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEGCAYAAACgt3iRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAabklEQVR4nO3df9Rd1V3n8feHECAFUoj8MBJsaY11oCqUDEVbK4gWrFqwIzWzqsRazZRBrc5oC85MrePSqWPtUuwCZVVLaHEwxVYyrNKaRIrTGQqGlpZfRWKDEEFif0FaMZDwnT/uDr2GJ8+9zz25eZ775P1a66y7z757n7PPcyDfe87Ze59UFZIkdXHQbDdAkjT5DCaSpM4MJpKkzgwmkqTODCaSpM4Onu0GjMshObQO4/DZboYkTZTtfPkLVXXsTOvN22ByGIfz8pwz282QpImyoa7/+1HqeZtLktSZwUSS1JnBRJLUmcFEktSZwUSS1JnBRJLUmcFEktSZwUSS1JnBRJLUmcFEktSZwUSS1JnBRJLUmcFEktSZwUSS1JnBRJLU2ViDSZKjklyf5HNJ7kvyXUmWJFmf5IH2eXRf+cuSbE5yf5Jz+/JPT3JX++7yJBlnuyVJMzPuK5PfBz5aVd8GfCdwH3ApsLGqlgMb2zpJTgZWAqcA5wFXJFnQtnMlsBpY3pbzxtxuSdIMjC2YJFkMvAr4Y4CqeqqqvgKcD6xpxdYAF7T0+cB1VbWjqrYAm4EzkiwFFlfVrVVVwDV9dSRJc8A4r0xeBPwT8L4kn07y3iSHA8dX1aMA7fO4Vv4E4OG++ltb3gktvWf+cyRZnWRTkk1Ps2PfHo0kaa/GGUwOBl4GXFlVpwFfo93S2oupnoPUNPnPzay6qqpWVNWKhRw60/ZKkkY0zmCyFdhaVbe19evpBZfH2q0r2ue2vvIn9tVfBjzS8pdNkS9JmiPGFkyq6h+Bh5O8pGWdA9wLrANWtbxVwA0tvQ5YmeTQJCfRe9B+e7sVtj3Jma0X10V9dSRJc8DBY97+zwPXJjkE+DzwRnoBbG2SNwEPARcCVNU9SdbSCzg7gUuqalfbzsXA1cAi4Ka2SJLmiPQ6SM0/i7OkXp5zZrsZkjRRNtT1d1TVipnWcwS8JKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKkzg4kkqTODiSSpM4OJJKmzsQaTJA8muSvJnUk2tbwlSdYneaB9Ht1X/rIkm5Pcn+TcvvzT23Y2J7k8ScbZbknSzOyPK5Ozq+rUqlrR1i8FNlbVcmBjWyfJycBK4BTgPOCKJAtanSuB1cDytpy3H9otSRrSbNzmOh9Y09JrgAv68q+rqh1VtQXYDJyRZCmwuKpuraoCrumrI0maA8YdTAr4yyR3JFnd8o6vqkcB2udxLf8E4OG+ultb3gktvWf+cyRZnWRTkk1Ps2MfHoYkaToHj3n7r6iqR5IcB6xP8rlpyk71HKSmyX9uZtVVwFUAi7NkyjKSpH1vrFcmVfVI+9wGfBg4A3is3bqifW5rxbcCJ/ZVXwY80vKXTZEvSZojxhZMkhye5MjdaeDVwN3AOmBVK7YKuKGl1wErkxya5CR6D9pvb7fCtic5s/XiuqivjiRpDpj2NleS1w2xjX+pqo9MkX888OHWi/dg4E+r6qNJ/gZYm+RNwEPAhQBVdU+StcC9wE7gkqra1bZ1MXA1sAi4qS2SpDkivQ5Se/ky+SK9q4DpxnW8qqpevK8b1tXiLKmX55zZboYkTZQNdf0dfUM5hjboAfxNVfXT0xVI8oGZ7lSSNL8MembynkEbqKqf2EdtkSRNqEHB5Ir90gpJ0kRzokdJUmeDnpm8KMm6vX1ZVa/dx+2RJE2gQcHkn4Df3R8NkSRNrkHBZHtV3bJfWiJJmliDgsmD+6MRmp8OOuSQ0SoePOKUcbueGanaM089Ndr+arT9SfPRtP/XVtWzI+CTfDfwwv46VXXN2FomSZoYQ/0ETPJ+4MXAncDuKU52v1tEknSAG/Z+wgrg5Jpu7hVJ0gFr2HEmdwPfOM6GSJIm17BXJscA9ya5Hb7+CkPHmUiSYPhg8o5xNkKSNNmGCiaONZEkTWfaZyZJbhy0gWHKSJLmt0FXJq+cbm4uei/NOnkftkfzyYiDD/OCZSPV23XUopHqHfzY4yPVe+bRx0ar9+STI9WT5rJB/7efP8Q2Rhw+LEmaLwaNgPdZiSRpIN9nIknqzGAiSerMYCJJ6myoYJLkh5N8OsmXkjyRZHuSJ8bdOEnSZBi27+bvAa8D7nKyR0nSnoa9zfUwcLeBRJI0lWGvTN4KfCTJLfzriR7fPZZWaV6oHTsGF5rCQSP+Znn41YePVO/JE0d7I+Q33rx0pHpL/mrLSPV2bfunkerVrl2DC0kdDRtMfhP4KnAYMOK7WCVJ89WwwWRJVb16rC2RJE2sYZ+ZbEhiMJEkTWnYYHIJ8NEkT9o1WJK0p4HBJEmAU6rqoKpaVFWLq+rIqlo8zA6SLGhjVG5s60uSrE/yQPs8uq/sZUk2J7k/ybl9+acnuat9d3lrkyRpjhgYTFp34A932MdbgPv61i8FNlbVcmBjWyfJycBK4BTgPOCKJAtanSuB1cDytpzXoT2SpH1s2Ntcn0zyb2e68STLgB8C3tuXfT6wpqXXABf05V9XVTuqaguwGTgjyVJgcVXd2gLbNX11JElzwLC9uc4G3pzkQeBr9F6KVVX1HQPq/R69MSpH9uUdX1WP0tvAo0mOa/knAJ/sK7e15T3d0nvmP0eS1fSuYDiM5w0+KknSPjFsMPnBmW44yQ8D26rqjiRnDVNliryaJv+5mVVXAVcBLM4SR+tL0n4yVDCpqr9P8kpgeVW9L8mxwBEDqr0CeG2S19Ab7Lg4yQeAx5IsbVclS4FtrfxW4MS++suAR1r+sinyNceNOvL6mS0Pj1Tv2E8fPbjQFF73+g0j1XvLax4cqd63XPfmkeq95KojBxeawq7ND864Tu18eqR96cA17KzBvwa8DbisZS0EPjBdnaq6rKqWVdUL6T1Y/6uq+glgHbCqFVsF3NDS64CVSQ5NchK9B+23t1ti25Oc2XpxXdRXR5I0Bwx7m+tHgdOATwFU1SNJRvuZBO8E1iZ5E/AQcGHb5j1J1gL3AjuBS6pq90/bi4GrgUXATW2RJM0RwwaTp6qqkhRAkhnNqFdVHwc+3tJfBM7ZS7nfpDcP2J75m4CXzmSfkqT9Z9iuwWuT/BFwVJKfBTbwr7v7SpIOYMM+gH9Xkh8AngBeAry9qtaPtWWSpIkxVDBJ8ttV9TZg/RR5kqQD3LC3uX5girwZjz2RJM1P016ZJLkY+I/Ai5J8tu+rI4H/O86GSZImx6DbXH9Krxvu/6BNyNhsr6ovja1VOqA9s+NfRqr3vI99dnChKXxoxFf1bP+Nvx6p3i0X/s5I9c7e8daR6i3/k2dmXGfXlodG2peDHQ9c0waTqnoceBz49/unOZKkSTTsMxNJkvbKYCJJ6sxgIknqbFBvru1MPd377veZDPXqXknS/DboAfyokzlKkg4gw070CEB7K+Jhu9erarT+g5KkeWXY95m8NskDwBbgFuBBnAZektQMe2XyG8CZwIaqOi3J2Tj2RHPMqIMdF924aaR6t37ulJHqffAHzx6p3kHfMFI1/vlblsy4zvO+9JWR9rXzS18eqR4184GVmluG7c31dHsPyUFJDqqqm4FTx9guSdIEGfbK5CtJjgD+Grg2yTZ6b0OUJGnoK5PzgSeBXwI+Cvwd8CPjapQkabIM+3Ksr/WtrhlTWyRJE2rYl2P1D148BFgIfM1Bi5IkGP7K5F8NXkxyAXDGWFokSZo4I83NVVV/AXzfPm6LJGlCDXub63V9qwcBK5h6zi5J0gFo2K7B/T23dtIbAX/+Pm+NNAtq166R6u28f/NI9b7pwa0j1TvomJkPPgTg0ENGqzeCLFgwUr3a6aDFSTfsM5M3jrshkqTJNWgK+j9gmttZVfUL+7xFkqSJM+gB/CbgDnozBb8MeKAtpwKj3RuQJM07g95nsgYgyU8BZ1fV0239D4G/HHvrJEkTYdiuwd8E9I81OaLlSZI0dG+udwKfTnJzW/9e4B1jaZEkaeIM25vrfUluAl7esi6tqn8cX7MkSZNk2ttcSb6tfb6M3m2th9vyTS1vurqHJbk9yWeS3JPk11v+kiTrkzzQPo/uq3NZks1J7k9ybl/+6Unuat9dniSjH7IkaV8bdGXyn4DVwO9O8V0x/ZQqO4Dvq6qvJlkIfKJd3bwO2FhV70xyKXAp8LYkJwMrgVPoBa4NSb61qnYBV7Z2fBL4CHAevjZYE2rUN0I+8w+PjLbDzHzWpBw02u+1esaJMQ5Ug3pzrW6fM37PaFUV8NW2urAtRW/k/Fktfw3wceBtLf+6qtoBbEmyGTgjyYPA4qq6FSDJNcAFGEwkac4Y6idLkguTHNnS/zXJh5KcNkS9BUnuBLYB66vqNuD4qnoUoH0e14qfQO8W2m5bW94JLb1n/lT7W51kU5JNT7NjmEOTJO0Dw17//req2p7klcC59K4o/nBQparaVVWnAsvoXWW8dJriU11X1zT5U+3vqqpaUVUrFnLooOZJkvaRYYPJ7tHuPwRcWVU30HtJ1lCq6iv0bmedBzyWZClA+9zWim0FTuyrtgx4pOUvmyJfkjRHDBtM/iHJHwGvBz6S5NBBdZMcm+Soll4EfD/wOWAdsKoVWwXc0NLrgJVJDk1yErAcuL3dCtue5MzWi+uivjqSpDlg2EGLr6d3VfGuqvpKu6L4lQF1lgJrkiygF3jWVtWNSW4F1iZ5E/AQcCFAVd2TZC1wL71p7i9pPbkALgauBhbRe/Duw3dJmkPS63Q1RMHe85LlbQDjscARVbVlrK3rYHGW1Mtzzmw3Q5p9k9A1uHyfyVyxoa6/o6pWzLTesL25fo1e993LWtZC4AMz3ZkkaX4a9jbXjwKnAZ8CqKpHdncVljTHjfCrv3zBhGZo2Ovfp9ogxAJIcvj4miRJmjTDBpO1rTfXUUl+FtgAvHd8zZIkTZJhZw1+V5IfAJ4AXgK8varWj7VlkqSJMewzE1rwWA/PTpPyhqq6dmwtkyRNjEEDDxe3aeHfk+TV6fk54PP0xp5IkjTwyuT9wJeBW4GfoTdQ8RDg/Kq6c8xtkyRNiEHB5EVV9e0ASd4LfAH45qraPvaWSZImxqDeXE/vTrSpTbYYSCRJexp0ZfKdSZ5o6QCL2nrovf9q8VhbJ0maCIPetLhgfzVEkjS5Zj4DnCRJezCYSJI6M5hIkjozmEiSOjOYSJI6M5hIkjozmEiSOjOYSJI6M5hIkjozmEiSOjOYSJI6M5hIkjozmEiSOjOYSJI6M5hIkjozmEiSOjOYSJI6M5hIkjobWzBJcmKSm5Pcl+SeJG9p+UuSrE/yQPs8uq/OZUk2J7k/ybl9+acnuat9d3mSjKvdkqSZG+eVyU7gP1fVvwHOBC5JcjJwKbCxqpYDG9s67buVwCnAecAVSXa/g/5KYDWwvC3njbHdkqQZGlswqapHq+pTLb0duA84ATgfWNOKrQEuaOnzgeuqakdVbQE2A2ckWQosrqpbq6qAa/rqSJLmgP3yzCTJC4HTgNuA46vqUegFHOC4VuwE4OG+altb3gktvWf+VPtZnWRTkk1Ps2NfHoIkaRpjDyZJjgD+HPjFqnpiuqJT5NU0+c/NrLqqqlZU1YqFHDrzxkqSRjLWYJJkIb1Acm1VfahlP9ZuXdE+t7X8rcCJfdWXAY+0/GVT5EuS5ohx9uYK8MfAfVX17r6v1gGrWnoVcENf/sokhyY5id6D9tvbrbDtSc5s27yor44kaQ44eIzbfgXwk8BdSe5seb8KvBNYm+RNwEPAhQBVdU+StcC99HqCXVJVu1q9i4GrgUXATW2RJM0R6XWQmn8WZ0m9POfMdjMkaaJsqOvvqKoVM63nCHhJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmcGE0lSZwYTSVJnBhNJUmdjCyZJ/iTJtiR39+UtSbI+yQPt8+i+7y5LsjnJ/UnO7cs/Pcld7bvLk2RcbZYkjWacVyZXA+ftkXcpsLGqlgMb2zpJTgZWAqe0OlckWdDqXAmsBpa3Zc9tSpJm2diCSVX9NfClPbLPB9a09Brggr7866pqR1VtATYDZyRZCiyuqlurqoBr+upIkuaIg/fz/o6vqkcBqurRJMe1/BOAT/aV29rynm7pPfOnlGQ1vasYgK9uqOvv31cN34+OAb4w243YDw6E4zwQjhE8zvnmJaNU2t/BZG+meg5S0+RPqaquAq7aV42aDUk2VdWK2W7HuB0Ix3kgHCN4nPNNkk2j1Nvfvbkea7euaJ/bWv5W4MS+csuAR1r+sinyJUlzyP4OJuuAVS29CrihL39lkkOTnETvQfvt7ZbY9iRntl5cF/XVkSTNEWO7zZXkfwFnAcck2Qr8GvBOYG2SNwEPARcCVNU9SdYC9wI7gUuqalfb1MX0eoYtAm5qy3w20bfpZuBAOM4D4RjB45xvRjrO9DpJSZI0OkfAS5I6M5hIkjozmMySJOe1qWM2J7l0iu/PSvJ4kjvb8vbZaGcXU02ps8f3aVPkbE7y2SQv299t7GqIY5z48wiQ5MQkNye5L8k9Sd4yRZn5cD6HOc6JPqdJDktye5LPtGP89SnKzPxcVpXLfl6ABcDfAS8CDgE+A5y8R5mzgBtnu60dj/NVwMuAu/fy/WvodagIcCZw22y3eQzHOPHnsR3HUuBlLX0k8LdT/Dc7H87nMMc50ee0nZ8jWnohcBtwZtdz6ZXJ7DgD2FxVn6+qp4Dr6E0pM6/U1FPq9DsfuKZ6PgkctXsc0qQY4hjnhap6tKo+1dLbgft47mwU8+F8DnOcE62dn6+21YVt2bMn1ozPpcFkdpwAPNy3vrdpYr6rXYrelOSU/dO0/WrYv8Okm1fnMckLgdPo/aLtN6/O5zTHCRN+TpMsSHInvYHj66uq87mcK9OpHGiGmSbmU8ALquqrSV4D/AW9wZzzyYymy5lQ8+o8JjkC+HPgF6vqiT2/nqLKRJ7PAcc58ee0euP4Tk1yFPDhJC+tqv7nfjM+l16ZzI69TR/zrKp6YvelaFV9BFiY5Jj918T9YuDfYdLNp/OYZCG9f2CvraoPTVFkXpzPQcc5n85pVX0F+DjPfbXHjM+lwWR2/A2wPMlJSQ6h9y6Xdf0Fknzj7heBJTmD3rn64n5v6XitAy5qPUfOBB6vNqv0fDFfzmM7hj8G7quqd++l2MSfz2GOc9LPaZJj2xUJSRYB3w98bo9iMz6X3uaaBVW1M8nPAR+j17PrT6o3pcyb2/d/CPwYcHGSncCTwMpq3SwmRaaeUmchPHuMH6HXa2Qz8M/AG2enpaMb4hgn/jw2rwB+Erir3WsH+FXgm2H+nE+GO85JP6dLgTXpvYDwIGBtVd24x78/Mz6XTqciSerM21ySpM4MJpKkzgwmkqTODCaSpM4MJpKkzgwmmnVJdvXNwHpnpphFeR/u66wkN871bfZt+6eSvKel35zkon203e9pM8be2cYazKTuBUlO3hft0PzhOBPNBU9W1anTFUiyoL7+KmeSHFxVOwdteNhyk6D1/99X3gC8q6reN0LdC4Ab6b1meyjz6Txoal6ZaM5K8mCStyf5BHBhko8n+a0ktwBvSfKCJBvb+xY2JvnmVu/qJO9OcjPw29Ns//D03kfyN0k+neT8ln9b/+R9bb+n7638NNs/Jb33RtzZ2ri85V/U1j+T5P0t70fafj+dZEOS46fY3juS/HJfm367bf9vk3xPy39ekrVt+3/Wtrlij+38DPB64O1Jrk1yRPv7fSrJXf3HtWdbk3w38Frgd9pxvTjJqUk+2cp9OMnRfW189nxN97fSPDDbc+u7uAC7gDv7lh9v+Q8Cb+0r93Hgir71/w2saumfBv6ipa+m98t5wRT7Oov2Lgrgt4CfaOmj6L274nDgl4Bfb/lLgb8dUP7Zbe6xrz8A3tDShwCLgFOA+4FjWv6S9nk0Xx9E/DPA77b0TwHvael3AL/c97fYXeY1wIaW/mXgj1r6pcBOYMUUbbsa+LGWPhhY3NLH0Bv1nGna+mzdtv5Z4Htb+r8DvzfV+XKZ34u3uTQXTHeb68+mWf8u4HUt/X7gf/Z998Hquy22F68GXrv71z5wGL1pM9YC6+lNjfJ64IMDyu/NrcB/SbIM+FBVPZDk+4Drq+oLAFW1+10oy4A/S++dEYcAWwa0HWD3JIR3AC9s6VcCv9+2fXeSzw6xnQC/leRVwDP0pho/HthbW79eMXk+cFRV3dKy1vD1vxc89/xpnjKYaK772oD1fv1zA01XbrcA/66q7n/OF8kXk3wH8OPAf5iu/FS3pACq6k+T3Ab8EPCxdnspTD2V9x8A766qdUnOoncVMsiO9rmLr/+/PNXU4YO8ATgWOL2qnk7yIL1Aube2zsQw50HzgM9MNMn+H70Zl6H3D+InZlj/Y8DPJ8/OAHta33fXAW8Fnl9Vdw1R/jmSvAj4fFVdTm8W1u8ANgKvT/INrcySVvz5wD+09KoZHke/T9C7mqL1uPr2Ieo8H9jWAsnZwAta/t7aup3eK22pqseBL+9+ZkNvksRb0AHHYKK5YNEeXYPfOWS9XwDe2G7l/CQzf8j7G/Rm+P1skrvb+m7X0wtUa4csP5UfB+5Ob/bZb6P3GtR7gN8EbknyGWD3NOfvAD6Y5P8AX5jhcfS7Aji2/U3eRu95xuMD6lwLrEiyiV5Q/hzANG29DviV1lngxfSC3++0fZ5K77mJDjDOGizNI+lNK76wqv6l/UO/EfjWqnpqlpumec5nJtL88jzg5vTeFhjgYgOJ9gevTCRJnfnMRJLUmcFEktSZwUSS1JnBRJLUmcFEktTZ/wfWw/4XCJdnIAAAAABJRU5ErkJggg==\n",
163
164
165
166
167
168
169
170
171
172
173
174
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
175
176
177
    "plt.pcolormesh(ER, RE, marginal.T)\n",
    "plt.xlabel('Error level scaling factor')\n",
    "plt.ylabel('Residual term [nT]')"
178
179
180
181
182
183
184
185
186
187
188
189
190
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By calculating the empirical mean and covariance of the individual components, we find the region where the probability mass is concentrated. Therefore the function `get_integ_bounds` from `corbass.integration` is available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
Maximilian Schanner's avatar
Maximilian Schanner committed
191
   "outputs": [],
192
193
194
195
196
197
198
199
200
201
   "source": [
    "from corbass.integration import get_integ_bounds\n",
    "integ_bounds, moments = get_integ_bounds(pars.bounds, posterior,\n",
    "                                         moments_out=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
Maximilian Schanner's avatar
Maximilian Schanner committed
202
    "For illustration purposes, we show the marginal distribution for the error scaling and the residual (blue), together with the Gaussian from the empirical mean and variance (orange). The black lines show the interval which will be integrated in the next step of `CORBASS`."
203
204
205
206
   ]
  },
  {
   "cell_type": "code",
Maximilian Schanner's avatar
Maximilian Schanner committed
207
   "execution_count": 9,
208
   "metadata": {},
Maximilian Schanner's avatar
Maximilian Schanner committed
209
   "outputs": [
Maximilian Schanner's avatar
Maximilian Schanner committed
210
211
212
213
214
215
216
217
218
219
220
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWsAAAEGCAYAAACjLLT8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd5gkVbn48e/pnKcnb2LzsoQlgyTFRVAvkkQRTCh48d6rYkbB8FPEe1GMBDEAKoIgAgKCCEpwCRKWvEtcYGFzmtw51fn9UT3LsDuhe6aqa7rn/TzPPNPTVX3q7e7qd6qrznmP0lojhBBicnM5HYAQQoixSbIWQog6IMlaCCHqgCRrIYSoA5KshRCiDnjsaLStrU3PnTvXjqYd8fLLLwOwePFihyN502SMqVaefPLJLq11uxPbbrR9ezS13sem8j49aLR925ZkPXfuXJ544gk7mnbE0qVLAVi2bJmjcQw1GWOqFaXUGqe23Wj79mhqvY9N5X160Gj7tpwGEUKIOiDJWggh6oAkayGEqAOSrIUQog5IshZCiDogyVoIIeqAJGshhKgDkqxrpGRorl++lnzRcDoUIUQdkmRdI4+9spH1fz2fu599zelQhLBfPg2/fju8dIfTkTQMSda18vytnO29kdblP3U6EiHst+pO2LwS/FGnI2kYkqxrZGCgH4D78ns6HIkQNbDiBvP3wCZn42ggkqxrRA2so6Dd/Kl7ASVDplITDSzVhX71HvPmXec5G0sDkWRdI4HUBvqIsF/xGdasX+90OELYRj93M8oo8nBpD1S2z+lwGoYk6xq5Kn8UD0fezdW+C9m64m6nwxHCNn/bHOe3xWN4gt0J6TQYJadDagiSrGtgIFvgX5mFbNz/bDLaB2secTokIWxxx4pNfP6RME/t8XU6OjrNO7P9zgbVICRZ18C6rn6OcD3LoliJ1/y70d77lNMhCWG5J9f0cN0N13HszAw//dA+EIibC+RUiCUkWddA9/rXuNp3IYsTj9LTegBzC69RTMvRhmgca7pTfPoPT3Ch9wouCv2OgNfNuvYjeXfpEmjaxenwGoIk6xpIbDYHwjTPWIhn3mG4lWbjcw84HJUQ1uhN5Tnj94+zp36FWXoT3v0+AkAg2swrhTby2u1whI1BknUNFLrfACAybT7TlxzJsbkLeJQlzgYlhEV+8a9XWdeb5ie7vQxuP+xxAgDt3gz/476N1LpnHI6wMUiyrgFX/zpKuCA6gznT2ljrW8izG5JOhyWEJTb3Z5nb7KNzzR2w+BgINAEQ9xQ513s9hTXLHY6wMUiyroFgeiN9nnZwe3C5FCd0bOXQl38IxbzToQkxYQPZAnt71pm9PvY+dfv9wVgLAIVkj1OhNRRJ1jYzDM1Psifw94Xnbb/vgKYBjsv+jfyGpx2LSwirJLJFtkb3gLNXwcKjt98fiUTJazelVK+D0TUOSdY225rI8XJxGsw9fPt9sV3fAUDX88scikoI6yQyeaIBD4RawOPbfn9TyEc/YUoZ6bpnBUnWNlvf1ccp7n+xyLN1+32LFy7kdaOT0hsyOEbUv8Mz9/G1DV+EVNdb7o8FvQzoMEoGxVhCkrXNujas5kfeK5ifXrH9vlnNQZ517UFL91NgyGQEor7tU1hBR349hFrfcn8s4OXk/He5a9F3HYqssUiytllq62oA4jMWbL9PKUVXy/6kDS8ktzgVmhATVigZRHSSjK8FlHrLsoDXTdoTpycnacYK8irarNC9BgBf69y33N+/68kckruETKDDgaiEsEYqVyRKmqJ3+EkGTvQ/yYFrrqxxVI1JkrXNPAPrMHBBbOZb7t9rlxZKBrywacChyISYuES2SFSlMfyxYZcf7nqew7bdUOOoGpMka5uFMhvp87aD2/uW+/eeFecM953MuuUkhyITYuIGsgVWGvNJte077PKiL0rQSIKWCTcmSpK1jbKFEt9Mf5S/7XXZTss6Y36ifg+dfU9D/wYHohNi4hLZIt8snsnW/b847PKirwk3BuQSNY6s8UiyttGGvgx9RInM2n2nZUopkp0HmX+slS58oj4lskUAs5/1MAy/OfRcyqROnCRrG63r6ufLnpvYTa8ednl8/v4kdYD86w/XODIhrJFODvCY/7PMeG3489IqGMdAQVauzUyUJGsb9W58nS96bmZW9pVhl++1SytPGYsovP7vGkcmhDVyyV46VR8B9/DnpDdOO5JFuWswOvascWSNR5K1jVJbXwcgOm3BsMuXzGziL6V3sC6wu8xTJ+rS4CQa/mjzsMtjoSAl7SKRK9YyrIY0/IkmYYlij9nHWjXPHnZ5e9TP49GjMWItXOqSAu2i/hTS5rlob2j4ZN3qyXKB50pyq3ywz7trGVrDkSNrG3kT683zdbFZI66z96w4K9f1yrBzUZeMwSJNI/SzjgR8fNRzH6WNMgHBREmytonWmmBmMwlv21sqke3oyOg67k6dTPblu2sYnRDW6CpF+Jf7MIhOG3Z5JBqnpBXFlNS0nihJ1jbpzxT4Su5M/nroTaOuF4234VUlUj2baxSZENZ5Xi3gZ/FvQfOcYZc3hXwMEH7zCFyMmyRrm6ztSQOKzs7hjzgGheKdAGT7t466nhCTUSJbHLGPNUAs6KFfhyEjExBMlCRrm6zvSnCh53J2y4x+ri7W3EpBuykmukZdT4jJ6KS+q7h804dGHE7eFPSylTiFkgw3nyhJ1jbp3vwGp3qW0VnaNOp6rRE/vUQppbbVKDIhrBMoDJiVUXcojzoo4vfw4cJ3uXX++bUNrAFJsrZJptzHOtA2d9T1WsI+rikezeuRA2oQlRDW8peS5N2REZcrpYgFvfRnCjWMqjFJsraJ0bfWvBEfvo/1oIjfw6/5IMuj76pBVEJYxzA0QSNFwTtysgY41fMAJ6w+rzZBNTBJ1jbxJdabN5pG7mMN5pFHW8hNpk/OWYv6ksybtayL3uH7WA9a5N7EvgPLpEzqBEmytkHJ0BSySfr8M8DjH3P9b7iu5pxXPlKDyISwTiJb5L7SfmyacfSo6xW9TXgpQCFTo8gakyRrG2zqz/CDwke48113VbR+wd9M2EhASeoniPqRyBb4Tel4Nu92+qjrGQEpk2oFSdY2MPtYw+zWcEXrlwLlWaEzMspL1I9EtoiPwqj9rAEYTNYyMGZCJFnbYEN3kj94f8iuPcsqWl+Fy8k6JeetRf1IZHK85D+dRS9cMup6RriddbodXcrXKLLGJMnaBj2b1/JO9wqaVWUF192RdgDyCRnFKOpHJtmPS2m8oaZR1xvoPJh35C4m175XjSJrTJKsbZDtMvtYe0aol7Aj3bqAnxQ+RL939KHpQkwmuaR5WsMXjo+6XixgThYtfa0nRpK1DXSvWceaeGXJOtS6C78oncQ27wwboxLCWoWUmaz9keFrWQ9q8Ra4xnsBxorRi5qJ0UmytoE/VZ6tfIw+1oNawj6m0U2iS2Y5F/WjlDaLM/nCoyfrSDjCO9zPYXQNP72dqIwka4sVSga9OcWm6BLwBip6TGvEx13+c2l/evQLNUJMJlt0M7/jRFTr/FHXi4UDDOggpZRU3psISdYW603nubx0PPccdm3Fj2kJ++nWMUh32xiZENZaozv5feh0aJ476npNQS8DhCEryXoiJFlbrDdlXkRpDo88O8yOmoJeeonizkqyFvWjkB5gujc35jDypqCXAR2GTH+NImtMkqwt1pvKcZfvHBZvvLXix7hdioSrCV9OjjxE/Ti892Zu6P8wFHOjrhcNeHjemEOfp71GkTUmSdYWS/R1s5trHVGqq4OQ8cYJFmSEl6gf3nyCgvKOeW3G63Zxnutz/HXWV2sUWWMaY5yoqFam35xEINDUUdXjHou+m9WlfTjLjqCEsIG3mCDriuCtYN0mqWk9YXJkbbF8wkzWoXh1X/m2thzIrfoIO0ISwhb+UoqcZ/Ra1oM+5vonZ71yps0RNTZJ1hYrlpO1L1Zdsp4eKDAj+RzkknaEJYSltNaEjCQFT7Si9VvcGebmV415fluMTJK1xboLXp5w7QWR6oaOLym9yNXGNyltft6myISwTjpf4obiO3lh9ocrWr/kHyyTKj1CxkuStcWeVEs4r/kH0DSzqsd5o+aReKp3sx1hCWGpRLbIHcYhbJ77/soeECjXD5EyqeMmydpiPak8zaHK+1gP8sfNC5KZPqm8Jya/ZK7APLWJFleqovVVsJysZQKCcZNkbbGP9f6K87qq76IUaTZPm+QGJFmLyW8gW+Rvvm+y5NXLK1rfiE7n36U9KapK+o6I4UiytlhrcRNhXf1cc/GmOBnt236BUojJLJHOElY53KHRy6MOKrbtzscK32KgeU+bI2tckqwtVCwZRI0B8v7KduChWiM+vlT4HC90HGtDZEJYK5swp6BzB0efeGBQU9A8oh6QvtbjJsnaQn2ZAi0kKAVaqn5sc8jHP4yDWO2aa31gQlgsl6ps4oFBTT7FMt+X8T7+KzvDamiSrC3Ul87TohIQaq36sT6Pi/0CG4ltfsSGyISw1uDEA4FoZQcmsUiQDtWH7pea7eMlydpCPck89xn7kZt+wLge/xn37Rz3xgUWRyWE9bboFr5Z+E/8s/ataP2moJd+whjSdW/cJFlbqDdT4KuFz1Dc85RxPT7rayZclJ1ZTH5bjRi3e9+Lq6WyqetiAS/9OiyDYiZAkrWFepPmUNpqalkPVQy0ENQZKGStDEsI6yU2c5BvLZQqu2A4eGTtzkmyHi9J1hbyb1rOC/4zaOt6clyPN4Llc93pLgujEsJ6u/Xcy+/yZ0MuUdH6Aa+LR/RerAktsTmyxiUlUi1USnYRUjkIx8b1eFe4DQCd2oaqcLJdIZzgzg+YN/yV7etKKf7oP5UtbdM4zMa4GpkcWVtIp8rTco2jNwhAqvNAPpz/NolwZecBhXCKp5AgqwLgrvx4Lxb0MpCVftbjJcnaQq7sxJJ1qHkajxp70FPwWxiVENbzFpNk3ZXVsh50mnE7P33lGCgVbYqqsUmytpA320tGBcec5mgkrSEXx7seJr3uGYsjE8JagWKSXJXJ2uv1EtBZyA3YFFVjk2RtoaeNhTwUP2Hcj28OB7jIexnBVbdbGJUQ1tJa85viMdw/94vVPW57mVSZGHo8JFlb6JbC23hgzhfG/fiWSIBeohgp6Q0iJq9c0eDx4iK2TatyGrqgTEAwEZKsLVIyNPlMgubQ+EtAtkZ89OgoSrruiUkskS1yiOsFZhXXVfU4d6gZAC2jGMdFkrVF+jMF7vKewwmvnz/uNkI+D30qhicrXxPF5JXIFrjUewl7r7+2qsfppln8qXgkGV/1hc6EJGvL9KbzNKskKjSxHTHpjuPP91gUlRDWS2SLxMigAtWVAnbHZ/ON4qfpjS22KbLGJoNiLNI3kGCBytBVHtgyXtfFzqQ55OXHFsUlhNWSqRR+VcAdqm7wl1nTWjOQTDMzHrQnuAYmR9YWSfWZM7z4ohNL1oXYbFblx9dPW4haGJx4wFPhLDGDmoJeXvB/iqZHfmhHWA1PjqwtMjjRrb+pfULt7ObZRLTvPsjtA/6oFaEJYalccnDigeaqHhcLekkRoJSWC4zjIUfWFtlWDHFx8SRCs/aeUDuL9eucVbwKpEi7mKS2uVr5WP4beBcurepxTUEvAzokM5yPkyRri6wzmrlMn0pw+sQunrij5pF5TibOFZNUb8HDv429CLdVV2wsFjDLpLokWY+LJGuLZAe6mR9KoybYjj/WAUCyZ/PEgxLCBu7+tZzkW467kKrqcdGAh34dfrNin6iKJGuLHLD5Bu4qfAqMiRWpCcbNZJ3t22JFWEJYblr3Y/zcdVHVw8ZdLsW97sN5uvm9NkXW2OQCo0U82R6SKkLEPf4RjADRFjNZ5wfkNIiYnNRgIaZA9XXblwXfTSrawjEWxzQVyJG1RQKFPtKepgm3E49GOTx7MSvmnGZBVEJYz51PYKDAV31vpZaAguRW0NqGyBqbJGuLhIr9ZL3V9TsdTmvYzwba2ZaRt0ZMTubEA0FwVb+PnlL6Oz9f9yEpkzoOkhEsYBiaqNFP3j/xmgexoIeTPQ8ye/X1FkQmhPV8xQRZ9/jGAJT85W+fUsypapKsLTCQLXB58VhWzz55wm0ppTjRu5wlm262IDIhrHcZp3Dt3P8b12NVsPztU7rvVU2StQV60wX+aryd5JyjLWkv7W0mWJDKe2JyWpVrJtEyvlnKXYND1OXIumqSrC3Qk0ixp3qDNm/OkvbyvmYipT65CCMmnULJ4N2lh9gtM76p5wZrWhfScjBSLUnWFsh2r+cO/zeZu/UeS9orBlrwUoRcwpL2hLBKIlvkbM8N7LPttnE93tU0ix8VTiERWWBxZI1vzH7WSqlLKmhnQGv9bQviqUvpchGnYHn04UTpwdnRMz3j6ssqKiP7dvUS2QIxlaY7ML5uqoGmNn5Zej8fCM1DpiCoTiWDYk4EvjPGOucCU3aHLibMabhCzdYk6/WzjmPhqt14MTabiQ2xEWOQfbtKiUyBmaTpGWeyjgW9zKCLdPd66NjN4ugaWyXJ+uda6z+MtoJSqrpaiQ2mWJ7gNhS3Jlk3R8MU8dCbztMRDVjSphiW7NtVSiYH8Cij6okHBjUFvdzq/w7FJ94Du19hcXSNbcxz1lrri6xYp6GluwFQIWsmDejw5vie5/dkXrnfkvbE8GTfrl4uaV4Y9IbG9z8sFvDSqyMgFxirNqELjEqp46wKpJ497t6PnwS/AFXOSTeS5rCfT3ruRq9/ypL2RPVk3x5et27iiNzPMXY/cVyPbwp66SOCSyaFrtpEe4McZEkUde6lwnQej79vXMNvhxNvbiGnPRSTUszJQbJvD2Mgr1mrOwnHxzcjUjzkpU9H8OQkWVdrQtlFa/1dqwKpZ22JF1jiXmdZey0RPz3EoHwuXNSe7NvDc/e8xn+7bydSHF+y9bpdpNwxfIV+iyNrfBWXSFVKBYDPAm8HNPAQ8Cutddam2OrGf6aupAUvcIol7cWDXl7WUQLpHkvaE6OTfbty0b4X+Ib3T5D/MtA5rjaWBY4i2Xwon7A2tIZXzZH11cCewKXAL4DdgWvsCKqeaK2JGgMUfNZ1GvC4XfS6mimWJjaRgaiY7NsV0oM1Pfzj7/+/LrYf/3S/06KIpo5qJh9YrLXeZ8jf/1JKPWt1QPVmIFukWSXoDlrbxf87kfPYvbOJyyxtVYxA9u0KqcFRtRMYrDXTn8fX9zzk9wJfyKLIGl81R9ZPK6UOGfxDKXUw8G/rQ6ovfakczSTAom57g1ojfrpT1tQaEWOSfbtC7vwAJVzgHX+SPYiV/CzxNeh+1cLIGl8lw81XYp7H8wKfUEqtLf89B3jB3vAmv/7eLuYoA3e0zdJ2365WcuDWv0DuL+AfX+1gMTrZt6vnLSRIuyJE1finhnaHzQMbne6Z8ATTU0klp0Gkv+koevJuTsufyzkL32dpuzN9SQ4vPgaJLZKs7SP7dpUu9ZzOfZ1n8KMJtOEtH9jkE134rQlrShgzWWut1yilXMAKrfX4itg2sO6s4kFjb77faW0VMU9sGgC5/k342xZa2rYwyb5dve6cGyIT+xYZjJl9tNP92yRZV6Gic9ZaawN4Vik12+Z46k6h+w3e61pOs9va88uB5hkAJLZtsLRd8Vayb1fn+OxtHJa8e0JthMoDanID3VaENGVU0xtkOvC8Umo5kBq8U2t9guVR1ZHYlkf5je8iDON0sLDoY7jNTNbp3o2WtSlGJPt2BUqG5gPci3dgAXD2uNtpjoX5XP4LfHLa0UyzLryGV02y/p5tUdSzchEn1wS/Gu6opbWTdUY7yZz0ta4B2bcrkMwWiak0qQn0sQZoDvm4wziEoz2zLIpsaqg4WWut31ICTil1OPBRYEqXhnNnesjjweeLWNpuRyzIQfmLOb9zT/awtGWxI9m3KzOQLdDExJN1S9jHErUaz4Y87CcJu1JV1QZRSu2rlPqRUuoN4H+BF22Jqo54870kXE0wga5Mw2kN+3C7FFsHpK91Lci+PbZkNk+ELCo4sWQdC3g513s9B7z0E4simxoq6We9K/Bh4CNAN/BnQGmtj7Q5troQyPeR9sSxdkgMuFyKc4K3cuALPfDemyxuXYDs29VKJ/pwKY07OLHSCi6XIu2O4S1YV/xsKqjkNMhLwIPA8VrrVwGUUl+2Nao68iM+wdJdQnzJhrbnePqYl3zShpZFmezbVegzgszP/pFb9jlk7JXHkPU0ESw8Z0FUU0clp0E+CGzGrJdwhVLqKJCBR2AWcXo+00K2zZ6zyvlAG7FSHxglW9oXsm9XI5EtYuAiEgpOuK2CP07ISIJhWBDZ1FDJtF63aK1PBXYDlgFfBjqVUr9SSr3H5vgmtWSuyIe4h0XFV2xp3wh34MYAKZVqC9m3q2NsfYn/8/yW1tz6CbdV8jfjwoDBKn5iTBVfYNRap7TW12qtjwNmAc9gzvw8ZfUlM1zg/S27DdhT88cVNXuhFgc22dK+MMm+XaGe1/mY515iOjX2umNY3f4uPu0+HyzuRdXIxkzWSqmdJgLUWvdorX+jtX7XSOtMBQM9WwFwW9zHepC7ZQ5PGwvpk+p7tpB9uzqFtHkU7Ao2Tbgtd/Mu3JdZhOHyTritqaKSC4y7K6VWjLJcARN/9+pQqm8LAL7Y+OajG4tn1n6clD+f24KLsOffwZQn+3YVXOnynKDhifd96vDlOJaHSG5ZTGy6tXV1GlUlyXq3CtaZklfAsv3mzhtssidZd8QCANLX2j6yb1fBm+2mgBdvID7htqZ7klziu4xtry0GSdYVqajqXi0CqUeFhDmhbbjZngoHHVE/N/i+R/Cpd8EeP7RlG1OZ7NvVKeaz9HnbabdgAFiofICTHZBJoSs1odnNp7rnggfxrvxPCc+0p+teW8RPO334+l6zpX0hqvH94ie4dM8bLGkrFm/F0IpCUirvVUqS9QRsy7noDczG7QvY0r7P46LP1Yw3s82W9oWoVKFk0J8p0BKxpgJ1cyRIP2FKKemWWilJ1hMwY+sDnOGdWG3fsSS9rQTzcvQhnNWbzvNT7694W9/fLWmvJeyjT4dRGUnWlZJkPQFLeu/l1MJfbd1Gxt9GtCjJWjirO5HjeNfDdOYnPiAGIORzc5ZxNnfN+Jwl7U0FkqwnwJ/vI+u1t2dXd9OeLFd7y5Bz4aj+3i58qoQ71mFJe0opukMLWFOYeM+SqUKS9ThprYkVu8n77e0BvW6XE/mv7Bcw5K0SDkr3mqNoA03W9Xx6h28Ve22+2bL2Gp1kgHEayBbppItiZLqt2+mI+ikamt503tbtCDGaXN9mAIIt1u3vR/EYJ/f8xrL2Gp0k63Ha0t1Li0pCfBdbtzOfdTzp/2/SK2+3dTtCjCaZzrBOtxNpnWlZm0V/MyGdgaIciFRCkvU4bUzB4uxVZPY709btNMXbaFUJMj0yca5wztPefXm/51e4p1k4pmBwEoNMr3VtNjBJ1uO0uT9LDh8dbVbPEfNWze3mkUyhXyrvCef0JPO0hH2WtukKtwBQTMooxkpIsh4n95qH+LbnGjp8BVu309EcoVtH0Ykttm5HiNG8c8sf+G7+Z5a26Y2aF+eT/ZKsKyHJepyatj3OmZ478fntGb04KOB106PieNJbbd2OEKOZnX2RecZaS9ssTD+IQ7OXsjW2l6XtNipJ1uPkS22kV8XBY83w29E86DuClV7ZoYVzosVeMr4WS9tsampiE630ZLWl7TYqSdbjFMltod/XWZNt3d32Cf7sPr4m2xJiR8WSQdzooxC0dkxBS9DNWe5bUK8/YGm7jUqS9Tg1F7eSDtQmWXfE/PQNJEDLEYiovd50gTY1gA5ZnKwjAT7vuYXI+vstbbdRSbIeh2SuiM/IU4xY1+d0NO9L38Y9mVPQ2f6abE+IoboHEjyv51Bo2dXSduNhH31EpOtehSRZj8Pm/izvyF/M6gO+VZPtectzPCa7pa+1qL2eDHwofx7pJR+3tF2/x80AUVxZSdaVkGQ9Dpv7swBMi4dqsj1/eSaa/m3WVDwTohrdKXOEYWvE2n7WACl3FG++z/J2G5Ek63HIvfEIl3kvYpaqTf/QULkeQ6p7Q022J8RQgTX/4h++r9NmUXnUoTKeJgIFOb1XiUomzBU7UFue51j3crKx2hxZN7XPAiBfLqYjRE31r2Oxaz2laMzypq/q/AabU2BvVfjGIEfW46AG1lPATaB5Rk2219Y+jcuLx/K6V2aBFrWnUua0cu5Iu+VthyNNdKWlVnslJFmPgz+1iV5XC7jcNdleNODlZ+o0VriX1GR7QgzlyWxjgAh4rD9nvY/xImdlfgn5tOVtNxpJ1uMQzm2h32vNjBmVUEqxS0SR6ZXeIKL2/LluBtzNtrS9i97AR9TdZAeknMJY5Jz1OPQW/WSbZtV0mxeUfkbbmi7gqZpuV4hVehfSkU7s2OM9EbNqZaJHCpWNRY6sq5QtlDg9+xUe3/PbNd1uLtBGrCQzQYvau6jwAf4178u2tO2PmufBU31yZD0WSdZV2t7HuilY0+2WQh3EjT6ZOFfUlAb6MgVawvYULAs1mQO+sv3dtrTfSCRZV6n/jae52fcdFhVW1XbD0Q7cSpPpk6+LonaKxRLP+M7ksG57JraNtLRT0G4y6QFb2m8kkqyrlNnyKvu7XqUlVtsja2+TOTCmd+u6mm5XTG2lQo4mlSYYsKdue6x1JotyV/Nsm1SVHIsk6yoVe80C7C3T59V0u+4Z+3J+4TS2Gk013a6Y2orlyWx9cXsqTMbDfpRS9KTtnXGpEUhvkCqpgY1k8RK2aecdSdOMBfyudAwHFJvYt6ZbFlOZLpjJOtQ83Zb23S7F//PfQOcbc21pv5FIsq6SL7WJblcbM5Wq6XY7ogHmq42ktrUA9nxwhNiRLpnJOtJq32jdt7tWkpNSCmOS0yBVWlNqZVVo/5pvtznk5Tbf/2PeK1fVfNti6srh5ZbS4TS12ZesM54YgbwUcxqLJOsq/bD4Uf4x79yab1cpRY8rjie9rebbFlNXUkX4nvdLeAIR27aR88YJliRZj0WSdRXyRYOuZI5pTfbOaD6ShKeFQE6StaidYqlEa9j6miBv2YY/TtSQrntjkWRdha4NqwHIJZIAABoRSURBVHnEfxYHZv7tyPYzvjYiBRk8IGpnRm41F+W+a+s2CqFpDOjadoWtR5Ksq9C/5XWmqx5i0agj288H24kbMgWSqB23LlLyhm3dxsu7fpp35C6mZMiE0KORZF2F1LY3AIh1znVk+2tmnsA5+TPJF2TIuagNjy6QC1g7q/mOmkPmaZZiSZL1aCRZV6HYY44ebJkx35kAZu7H341D2FaeE08IO2mt8VDCCFo/6cBQcwqvcZX3Qkq5pK3bqXeSrKugBjaQ0EGiMXtq+45leqDAYa7n6NkiczEK+5XKA2JUxN4j6yavwVL3s+hi1tbt1DtJ1lV4ldnc6zsSVeMBMYNmGpu4zncBxTWPOrJ9MbUUDM1G3Up+mr3jCiLN5kQeRrFo63bqnYxgrMKNHE2k8z94v0Pbj3eY5d9zfZscikBMJQXtZq3uwD3rAFu3E20xk7UuSX2Q0ciRdRW29aUc62MN0Nw2HUMr9IAMzRX2KxaLeDBoDXtt3U60qZWSVmBIsh6NJOsKacPg/vypHJuyp65vJTw+P12qBffAWsdiEFOHN7OVA10v0+rN2bod5XLzkmsBRS3paDTy6lSomM/iUQaBaKujcWwOLaQ9+ZKjMYgpopTHQBGP27/PfyX2c7rc9vY6qXeSrCtUzGcACLTNdjSOF3f7Ap/JfJaiDCAQNlOlAkU8eD1u27fVHPZSkH7Wo5JkXaFSwfwqGO2o7aQDO5q+28G8pGeTzsmVc2EvpQsUVW36IJyRuZqZhTU12Va9kmRdIV00k3Xr9LmOxrHX9Agnu++nkOpzNA7R+NxGkVKNknUHPYR0uibbqleSrCuUUQF+axxHvCnmaBzNkQDf8v4Jb1aq7wl7bVOtJDy1uUZjBJpxU0JOhIxMknWF+olyTfRMxwbEbKcUG0OL8RblKETYa6tuIuur0WjdUAtuDEpFqXszEknWFTIKOWbE7O1vWqls214EyFEsyY4t7GEUcniNDD5XbY51PRHzCH5wgl6xM0nWFVpYeJnPZK90OgwAwnMPQKHJpqRgu7DHwKZX2FutJlKjGVyiHXMZ0CEyUlFyRJKsK2EUcWNQjM50OhIAZux+KACFbMrhSESjSnWbxcKUx95ZYgbNOvgkXmIuA3lJSSOR2iAVKOSzeAF3fJbToQAQmzaf59y74dN+p0MRDSrdY5Y0cNcoWfs8LsJ+D6mcDDkfiSTrCpTKyTrQ6uyAmO2Uwu8PkMxKX2thj3y5/ozLU7sDgtl6I6qQIl808HnkCHtH8opUoFQw6+zGOp0dEDNUiyvN7NJauvsTTociGpCR2IYGvN7aHFkDeLw+wmR5ad3Wmm2znkiyrkBWhVinO2idPsfpULYLuTWtaoDXX3jC6VBEA1oZeyer9Qw87tqlCE8whkKz/oWHa7bNeiLJugIJw88m2miNhpwOZTtf2Jy0t3+1JGthvZfVPHpVnFoOK/CGmgDIv7G8dhutI5Ksx5LuwUj3EA+6cbkcHhAzhNsXooQb15YVTociGlCk6xnC7hpfE3F7KSgfsZ5na7vdOiHJegzbnrqN+cYaWn2Tr/9n3h2kI/Gi02GIBvTJjeczXdf+3HE60Mmd2SV0Je2toV2PJFmPIbHyTgq4CUebnA5lJyVvlFTJxdaBjNOhiEaiNbFSL9pV+xG7rqaZ3FhaytNrpVDZjiRZj8Yo0bHt3yRUlIDX/pq+1dLxOZyS/y7PbZSRjMJC+SQB8o4k67DfzUxXL6+8uqrm257sJFmPIrdmOREjQckfdzqUYYV8bpSCFetrMyRYTA1G92oAtLt23fYGudAs832Rmauurvm2JztJ1qPY+PQ/KGmFz+GpvEbidimuD/+MRSt/6nQoooHknr8DQytKPgdO/SkXW0MLmZZ4jmLJqP32JzFJ1qO42vNBji/9mGjIuRnNx9LuzTFr4Gm0lkrAwhrP7fIRzih8HY/PmXIGuc79WcJrrNok3xiHkmQ9ivtXddM+fx9cTtewHkW2bQmLjDfY0if1rYU1bnspxWPu/WgKOlMSOLbwUMIqx+svyhiCoSRZj2DbE7dyRt+lvHtB0OlQRhWcewAhlWP1S884HYpoAKVnb6RlxeUctbgDt0PjCloXHwZAavVjjmx/spJkPYLEk3/mGPdyDt9jrtOhjGr64kMA6F/9uMORiEaQevAyjio+wHH7zHAsBtW6gMvbv8GNA7s7FsNkJMl6OIZB+5Z/85RnP+a1R52OZlSB6bvzsOdgXuiXAopigvo3EOt6mvs4mCN363AuDqUo7HEyj3cH6EvLzDGDJFkPI7f2CaJGP/2zljodytjcHm5e/GOu614sFxnFhJRevB2A1IJjHB9X8La2HB9138vK1esdjWMykWQ9jE1P/g1DK6btf4zToVRk71lN5FN9bJSLjGICBp66mVXGTA468BCnQ2GJZwMXeH/L5hcfcTqUSUOS9TBW9+S4T+/Pgbvv6nQoFTki/xArA2fy+FNy9VyMk1Fia6rE3eow3rm43eloCM49CAC9Xq7FDJJkPYz/GziGP8z5AUHf5BtiPpzZi/cHYOODf2RrIutwNKIe5Q3Fyamv89oen8PvmQT7fbCZrf7ZtPevxDDk9B5Ist7J+s3beG1bkqWLHbzAUiXXtD1ILDqRM/VfuPLG25wOR9Shh194g0S2yPH7TI5JoQHS7fuyRK9i9TaZDQkkWe8kd9uXucP3TZZOgq+C1YiedBEFX5wT3/hf7l651ulwRD3J9nPYLYdwZuA+Dl/Y5nQ024XmH0y7GuDFl192OpRJQZL1UIZB++YH2eidzfy2sNPRVCfUgv+kS5jj6uKav97FQFZmiRaVyb94Jz6do3n+AZNqotq2Qz7O4fyOh7ucGfY+2Uyed2YSKLz4N2JGH/0z34maxEPMR+LZ4zhe//jDPJSayYV3vuR0OKJO9DzxFzbrZvY6+CinQ3kLVyjOgtm7sPz1Hkpy3lqS9SDj6Wtx3XQ6LxqzaT3gA06HM257LZzLGYfNJfn4dTz+6ianwxGTXT5Ny8b7ud91MIctnHyn/j7T/Din9l7OWdc9RbYw+WZrqiVJ1kB/usC99/6Dh4u7c92ev+GIveY7HdKEfG2PXi72/ZKX//z/pvwOLkaXe+kf+HSOxLxjajqTeaUOberjvzx3cNjLP+BTVzxAb2rqjmicfO9OLZWKvLLqBY7/xUOc1XsK69/3B84/5VDHCthYJbDg7Wye9wE+nP8LN94mvUPEyB7MzOP8wmnsceh/OB3K8I74Ghz2eU5z3813t3yer172J9b1TM3BX1M3WedTbL78A0SuPQ5PIcF1/3U4Hzl0YV2eqx7OtFN+TtLbyuHPnsslPz2P3z2wis390gdbAJleePRX9D7wa65ameP20Ps5eEGn01ENz+OD9/wvfPwvzA+muSx9Nv/9y9t5bsPUq3U9Jar/aMNgIFdi29ZNqGeuhe7XiG5ZTnt2Lb9t+hzXf/o9dEQn7wQD4xKMEzjlCjpv+k9OS/6e/f9+AOf//RU+N+0F9l44hwX7HkE0GiMc8BHyunHV+bcJMQatSb/+KH0P/Ib2NXfg1XmWlw7kocIufPGoRZP/2+TCo/Ge9Sg9T9xK/6PtnPqbR/jm4WGmtTTT0tpKR0ucjqgf7yQ8lWOVmibrC+96ifc8/mmade9b7n/Ssx+/DpwJwMWprxHSb/2a82/PoVwV+DgAv0l+HjeD0/2YV4jv8S7lT/5T8BpZrkx9ATclXBi4dQk3JX5rHM+l+ePooJflgR/SqyO8rqdxz6IfcfpHzpxU3ZWsFNj1SDj3VcKJjdyTj3PHik184N9nM+uJjVAemZ7THu40DuAc11cJ+91cWfwmcZ0ov7LmB/gRz0FcEfwUAJclv0KA3PZtaOABz+FcE/goaM0VqbN2iuOf3qO40f8BAjrDpamzd1p+m+9Ybve9j7jRy4/T3wbgZt+JHPPJc9h9eszKl8QWxZLBV356BV/JXrrTsosDn+U5z57sV3yGz2av2Gn5hcGv8Kp7AYcWHuNTuZ3nHTw/+A3WuWdxZOF+Ppq74S3LFPCN0HlsdbVzTP6fnJy/FQCXNvBSxEOBTwYuok/H+HThj5xe+guGDnCTPoIXZ36Q2Xscwj8WtbNrZ8SaF8JukQ46l/4XN++f5cdXXsOpD38djzJzQU57GCDI992f43HfweyjX+Cc/C8BNbgbo1H8PPA5nvfswQHFp/lM9sqdNvHD4Fd5zT2fwwqPckbump2Wnxf8JhvcM3lXYRkfyd240/JzQ+fT7Wrlffm7OCl/+07LvxS+kLa2Dq785EFVP/2aJuvOqJ9EZA66tMMEtIEZLGo2d5j+zXPJ6sxbFruD01kUN5f3GvNQGAy+AxqFLzyTxbEoLh1kw9a90MqNVm4M5UErFzOb9+Fbc3anM+bnqeDTtLdPY8+Yn/0nw7Bau7lc0DSLBcAXjloEhz/Gxmf+Se/a5ynl0xj5DB7fbE5p2oVUrsjA+sUUjCRobSZsrXEFZ7E4bpaK7dfzSOnBizxmSveEprOoKQJa07N53k4h+MPTWBSL4DU89GzZeXkw0sGiaIRwqUTPVnN5ONrheOW3Siml6Gxrpqdn5+fW3tzCokCEtmwrPb07L5/eGkf5IjSn2+jp33n5jLY4AW+EaKqDnoGdl89pj9PsiRJJdtKVWAiYn4mSy0dJedins42cJ0o+dTh3uRcTf9upfGDhLpNjSPk4dcYC/PgTS0m+dgnJgV7SiT5yyT4K6T7mRxfj9rYyPTONnq5dy/vwm3M5dra0UvRHacu00dM3n8F9eNC01jguX4TmdOuw78estiZC3gjRVPuw78fs9hgtngjhZAc9iZ2Xz+uI0Rwf3xgOZUdZzQMPPFA/8UTjFBVaunQpAMuWLXM0jqEmY0y1opR6Umt9oBPbbrR9ezS13sem8j49aLR9uzG//wshRIORZC2EEHVAkrUQQtQBSdZCCFEHJFkLIUQdkGQthBB1QJK1EELUAUnWQghRB2wZFKOU2gasAdqALss3UDsSv3NGi32O1tqR4stD9m2n1PN7OlSjPA+w9rmMuG/bkqy3N67UE06NNLOCxO+ceo7dTo3yujTK84DaPRc5DSKEEHVAkrUQQtQBu5P15Ta3bzeJ3zn1HLudGuV1aZTnATV6LraesxZCCGENOQ0ihBB1QJK1EELUAduStVLqP5RSLyulXlVKnWvXduyglPqdUmqrUuo5p2OpllJqF6XUv5RSLyqlnldKfdHpmKqhlAoopZYrpZ4tx/89p2Oy00jvl1KqRSl1t1LqlfLv5iGP+Ub5c/WyUuq9Q+4/QCm1srzsEuXA7M9KKbdS6mml1N/q/HnElVI3KaVeKr83hzr+XLTWlv8AbuA1YD7gA54F9rBjWzbFfwSwP/Cc07GMI/bpwP7l21FgVZ299gqIlG97gceAQ5yOq9bvF/Aj4Nzy/ecCF5Zv71H+PPmBeeXPmbu8bDlwaPk1vBM4xoHn8xXgOuBv5b/r9Xn8ATizfNsHxJ1+LnYdWb8NeFVrvVprnQeuB060aVuW01o/APQ4Hcd4aK03aa2fKt9OAC8CM52NqnLalCz/6S3/NOxV8FHerxMxEwbl3+8v3z4RuF5rndNavw68CrxNKTUdiGmtH9Fmlrh6yGNqQik1CzgWGDoTbT0+jxjmAdtvAbTWea11Hw4/F7uS9Uxg3ZC/11NHCaNRKKXmAvthHp3WjfJX6WeArcDdWuu6in+8dni/OrXWm8BM6EBHebWRPlszy7d3vL+WLgK+DhhD7qvH5zEf2Ab8vnxK50qlVBiHn4tdyXq48zINe3Q0GSmlIsBfgC9prQecjqcaWuuS1npfYBbmEcoSp2OyWxXv10ifLUc/c0qp44CtWusnK33IMPc5/jzKPJinQX+ltd4PSGGe9hhJTZ6LXcl6PbDLkL9nARtt2pbYgVLKi/nBv1ZrfbPT8YxX+avnMuA/HA7FViO8X1vKX6Mp/95avn+kz9b68u0d76+Vw4ETlFJvYJ72fJdS6o/U3/OgHMP6Id/obsJM3o4+F7uS9ePAIqXUPKWUD/gwcJtN2xJDlK82/xZ4UWv9M6fjqZZSql0pFS/fDgJHAy85G5V9Rnm/bgM+Wb79SeCvQ+7/sFLKr5SaBywClpe/lieUUoeU2/zEkMfYTmv9Da31LK31XMzP+31a64/X2/MoP5fNwDql1OLyXUcBL+D0c7Hxaur7MK9svwZ8q9ZXcycY+5+ATUAB87/jfzodUxWxvx3zq9YK4Jnyz/ucjquK+PcGni7H/xzwHadjcuL9AlqBe4FXyr9bhjzmW+XP1csM6V0AHFh+zV4DfkF5hLIDz2kpb/YGqcvnAewLPFF+X24Fmp1+LjLcXAgh6oCMYBRCiDogyVoIIeqAJGshhKgDkqyFEKIOSLIWQog6MCmTtVKqpJR6ZsiPbVX7lFJLByuETeY2h7R9ulLqF+Xb/6OU+oRF7b6jXPXtmXL/5moe+36l1B5WxCGEGN6kTNZARmu975CfH+64glLKvcPfnkoarnS9eqC1/rXW+mqLmvsY8JPy652p8rHvx6w8VrFGeh9qYcgBzHNKqdsHBw6No50rh/vHOvQgYJztJoe5L66U+ux426xy+0uVUv1Kqb+Psd4t5dfx1fL6gweEhymlrlVK9SilTq5FzNWarMl6WEqpN5RS31FKPQR8SCm1TCl1gVLqfuCLSqk5Sql7lVIryr9nlx93lVLqZ0qpfwEXjtJ+WJm1rB8vF3A5sXz/Y0qpPYest6xcp3bY9Udpf09l1mp+phzjovL9nyj//axS6pryfceXt/u0UuoepVTnMO2dp5Q6e0hMF5bbX6WUekf5/pBS6oZy+38ut3ngDu2cCZwCfKe8w0bKr99TyqzFe+KQdd8Sq1LqMOAE4Mfl57VAKbWvUurR8nq3qHLd3x3fr9FeK7GTwQOYJZgVIT83nka01mdqrV+wNrQRxYGqkrUyjTcvPai1ft9oK2itT9Jm3Zkzy+sPHhA+rLX+GJN5pLUTI5wqGD1U4s3RXM8Ap5bvfwP4+pD1lgG/HPL37cAny7c/Bdxavn0V8DfKNWZHGW11AfDx8u045gjMMPBl4Hvl+6cDq8ZYf3ubO2zrUuBj5ds+IAjsiTnqqa18f0v5dzNvzpF5JvDT8u3TgV+Ub58HnD3ktRhc533APeXbZwO/Kd9eAhSBA4eJ7Srg5PJtD2ZpR4A2zJKPapRYtz+2/PcK4J3l2+cDFw33fslPVZ+J5JDb/7PDfv81zBIPK4bsp2HgDsw6y8/x5mdo2eD7D5xR3mfvB64Ysl/t+H4my78jmCP3ngJWAicOF9+Q+64HMpif4R+PEutczNKwv8QcvToHSGIeWD0J3INZdnkZsBo4YZhtLWXIZ6789zLMuh4vAdcyZPTgjusP9zmYbD+T9atoRpv//Ybz51H+PhT4QPn2NZjFwgfdqLUujbHd92AWozm7/HcAmA3cANwNfBfzCPTGMdYfySPAt5RZ9/dmrfUrSql3ATdprbsAtNaDdbRnAX9WZsEYH/D6GLEDDBYBehLzAwDmcOaLy20/p5RaUUE7CrhAKXUEZrnLmUAnMFKsbz5QqSYgrrW+v3zXH3jz9YKd3z9RhfLpv6Mo11pWSr0HsxbF2zDft9vK71s7sFFrfWx5vaYd2pkOfA84AOgH/oWZKEeTBU7SWg8opdqAR5VSt+lylhvGucCSwc/yKLGuBRYDZ2itP1teNwws01qfo5S6Bfhf4N2Yp9v+QGVHwPthHmBsBP6NWWzqoQoeNylN1mQ9mtQYfw81dCcabb1BCvig1vrlnRYo1a2U2hs4Ffjv0dYf7pQFgNb6OqXUY5gF2v9RPv2gGL5s4qXAz7TWtymllmIeRY8lV/5d4s33djzTCH0M88N+gNa6oMxKaoFRYq1GJe+D2FlQmTW+52L+M767fP97yj+DiTaCmRAfBH6ilLoQ8wjywR3aOxgzGW4DUEr9Gdh1jBhG+ie+ucLnMFKsa4E1WutHh6ybB+4q314J5Mr74krePBAZy3Kt9XqAIa9d3SbrujpnXYGHMSt+gZlwqn1j/gF8XilznjSl1H5Dll2PWVi9SWu9soL1d6KUmg+s1lpfgnlksDfm18pTlFKt5XVayqs3ARvKtz+5Y1tVeAjz2wDKvLC0VwWPacKsTVxQSh2J+bWUUWJNYE5Jhda6H+gdPGcOnIb5NVtMzOC3zTmY37QGz1kr4Af6zXOvC7XWv9Var8I8al4J/EAp9Z1h2hzpH2+Rcm4o79u+8v1D/4nvC2zB/CdeqWFjLS/b8Z94YcgRu0H5QERrbVD5QWZuyO2hBzB1abIm66B6a9e9nXqDjOALwBnlr/qnUf1FrO9jTiO1QpmT5X5/yLKbMP8R3FDh+sM5FXiu/F9+N+BqrfXzwP8B9yulngUGy2SeB9yolHoQ6KryeQz1S6C9/Jqcg3musH+Mx1wLHKiUegLzA/oSwCixXg98rXwxdAHmP5cfl7e5L+Z5a2GB8j/DLwBnK7MO9j+ATylz8gKUUjOVUh1KqRlAWmv9R+AnmPWYh3oMWKqUai2386Ehy97ATPRgTlnlLd8e6Z/4SLb/Ey8bNtZKn/uU5/RJc/mx9wdz8uJA+fYCzA+iz+m45Kfq9zG5w9+3A6eVb38R8wh6JeZ1kQXAe3mz7OrjvHlRcRnDX2C8mDcvMHYCj2JO9voD3rzA2FZu/wnMeRZfBOYOF9+QOK/DvMD541FincsOk1Pz1guq51G+kD7Sthj+AuPQv38BnD7S8iH3X8UkvcAoJVIbnFIqinnxyIv5NfQcrfWdzkYlhLXK13XO1lofN8F2rsJM4jdZEZeV6vocjhibNmfMPnDMFYWob3lgiVLq73qMvtYjUUpdCxyGecpz0pEjayGEqAOT9QKjEEKIISRZCyFEHZBkLYQQdUCStRBC1IH/DxUPJnJyIRWPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
Maximilian Schanner's avatar
Maximilian Schanner committed
221
222
    }
   ],
223
224
225
226
   "source": [
    "from scipy.stats import norm\n",
    "fig, ax = plt.subplots(1, 2)\n",
    "\n",
227
    "ax[0].plot(errs[0], marginal.sum(axis=1)*ress[1])\n",
Maximilian Schanner's avatar
Maximilian Schanner committed
228
229
230
    "ax[0].plot(errs[0], norm.pdf(errs[0], loc=moments[1, 0], scale=np.sqrt(moments[1, 1])), ls='--')\n",
    "ax[0].axvline(integ_bounds[1, 0], color='black')\n",
    "ax[0].axvline(integ_bounds[1, 1], color='black')\n",
231
    "ax[0].set_xlabel('Error level scaling factor')\n",
232
233
234
    "ax[0].set_ylabel('[Arb.]')\n",
    "ax[0].set_yticks([])\n",
    "\n",
235
236
    "ax[1].plot(ress[0], marginal.sum(axis=0)*errs[1])\n",
    "ax[1].plot(ress[0], norm.pdf(ress[0], loc=moments[2, 0], scale=np.sqrt(moments[2, 1])), ls='--')\n",
Maximilian Schanner's avatar
Maximilian Schanner committed
237
238
    "ax[1].axvline(integ_bounds[2, 0], color='black')\n",
    "ax[1].axvline(integ_bounds[2, 1], color='black')\n",
239
    "ax[1].set_xlabel('Residual term [nT]')\n",
240
241
242
243
244
245
246
247
248
249
250
251
252
    "ax[1].set_ylabel('[Arb.]')\n",
    "ax[1].set_yticks([])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we draw the found area of interest on top of the marginal distribution, to show wich area would be integrated over in the following step of `CORBASS`:"
   ]
  },
  {
   "cell_type": "code",
Maximilian Schanner's avatar
Maximilian Schanner committed
253
   "execution_count": 10,
254
   "metadata": {},
Maximilian Schanner's avatar
Maximilian Schanner committed
255
   "outputs": [
256
257
258
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEGCAYAAACgt3iRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAayklEQVR4nO3df7RV5X3n8feH3/jjigxKCfgzQ5OBNPUHRVLTVGNUJmmF2mDoSgJJTWmoTWxn0gidqbHTlYxOE1dqXJiwkgY0ZsgN1ci4xASo2nEGJahERCTSYJRAJcZfaJSf3/6xn6s7l3PP2efse+6959zPa62zzj7PeZ69v/tuuN+7936eZysiMDMzK2NIfwdgZmatz8nEzMxKczIxM7PSnEzMzKw0JxMzMyttWH8H0CwjNDJGcXR/h2Fm1lL28sJzEXFCve3aNpmM4mjO0QX9HYaZWUtZGyt/2kg7X+YyM7PSnEzMzKw0JxMzMyvNycTMzEpzMjEzs9KcTMzMrDQnEzMzK83JxMzMSnMyMTOz0pxMzMysNCcTMzMrzcnEzMxKczIxM7PSnEzMzKw0JxMzMyutqclE0hhJKyU9IWmrpHdJGitpjaQn0/vxufqLJW2XtE3SxbnysyVtTt/dIEnNjNvMzOrT7DOTfwDujoi3A78JbAUWAesiYjKwLn1G0hRgLjAVmAkskTQ0recmYAEwOb1mNjluMzOrQ9OSiaQO4D3ANwAiYn9EvAjMApanasuB2Wl5FrAiIvZFxA5gOzBd0gSgIyLWR0QAN+famJnZANDMM5PTgZ8D35T0iKSvSzoaGB8RuwHS+4mp/kTgmVz7nalsYlruXn4ESQskbZS08QD7endvzMysR81MJsOAs4CbIuJM4FXSJa0eVLoPElXKjyyMWBoR0yJi2nBG1huvmZk1qJnJZCewMyIeTJ9XkiWXZ9OlK9L7nlz9k3LtJwG7UvmkCuVmZjZANC2ZRMS/Ac9IelsqugB4HFgFzE9l84E70vIqYK6kkZJOI7vRviFdCtsraUbqxTUv18bMzAaAYU1e/6eAWyWNAH4CfJwsgXVKuhx4GpgDEBFbJHWSJZyDwBURcSitZyGwDBgNrE4vMzMbIJR1kGo/HRob5+iC/g7DzKylrI2VD0XEtHrbeQS8mZmV5mRiZmalOZmYmVlpTiZmZlaak4mZmZXmZGJmZqU5mZiZWWlOJmZmVpqTiZmZleZkYmZmpTmZmJlZaU4mZmZWmpOJmZmV5mRiZmalOZmYmVlpTiZmZlaak4mZmZXmZGJmZqU5mZiZWWlOJmZmVpqTiZmZleZkYmZmpTmZmJlZaU4mZmZWmpOJmZmV5mRiZmalOZmYmVlpTiZmZlaak4mZmZXmZGJmZqU5mZiZWWlNTSaSnpK0WdImSRtT2VhJayQ9md6Pz9VfLGm7pG2SLs6Vn53Ws13SDZLUzLjNzKw+fXFmcn5EnBER09LnRcC6iJgMrEufkTQFmAtMBWYCSyQNTW1uAhYAk9NrZh/EbWZmBfXHZa5ZwPK0vByYnStfERH7ImIHsB2YLmkC0BER6yMigJtzbczMbABodjIJ4AeSHpK0IJWNj4jdAOn9xFQ+EXgm13ZnKpuYlruXH0HSAkkbJW08wL5e3A0zM6tmWJPXf25E7JJ0IrBG0hNV6la6DxJVyo8sjFgKLAXo0NiKdczMrPc19cwkInal9z3A7cB04Nl06Yr0vidV3wmclGs+CdiVyidVKDczswGiaclE0tGSju1aBi4CHgNWAfNTtfnAHWl5FTBX0khJp5HdaN+QLoXtlTQj9eKal2tjZmYDQNXLXJIuLbCO1yPirgrl44HbUy/eYcC3I+JuST8EOiVdDjwNzAGIiC2SOoHHgYPAFRFxKK1rIbAMGA2sTi8zMxsglHWQ6uFL6RdkZwHVxnW8JyLe2tuBldWhsXGOLujvMMzMWsraWPlQbihHYbVuwK+OiD+uVkHSt+rdqJmZtZda90xurLWCiPhIL8ViZmYtqlYyWdInUZiZWUvzRI9mZlZarXsmp0ta1dOXEXFJL8djZmYtqFYy+Tnwpb4IxMzMWletZLI3Iu7rk0jMzKxl1UomT/VFENaehowY0VjDYQ1OGXfocEPNDu/f39j2orHtmbWjqv9rI+KNEfCSfhs4Nd8mIm5uWmRmZtYyCv0JKOkW4K3AJqBripOuZ4uYmdkgV/R6wjRgSlSbe8XMzAatouNMHgN+rZmBmJlZ6yp6ZjIOeFzSBnjzEYYeZ2JmZlA8mVzTzCDMzKy1FUomHmtiZmbVVL1nIunOWisoUsfMzNpbrTOTd1ebm4vsoVlTejEeaycNDj7UKZMaandozOiG2g179qWG2h3e/Wxj7V57raF2ZgNZrf/tswqso8Hhw2Zm1i5qjYD3vRIzM6vJzzMxM7PSnEzMzKw0JxMzMyutUDKR9HuSHpH0vKSXJe2V9HKzgzMzs9ZQtO/ml4FLgc2e7NGa7e4tn+/vEAaFC4fM6e8QrI0Uvcz1DPCYE4mZmVVS9Mzks8Bdku7jVyd6vL4pUVlbiH37aleq4lMfvLGu+rsuOr6h7bw+/mBD7cZtHNpQu44f7mqo3eEXXmioXXR7AuWSjdc1tB6zaoomk88DrwCjgAafxWpmZu2qaDIZGxEXNTUSMzNrWUXvmayV5GRiZmYVFU0mVwB3S3rNXYPNzKy7mslEkoCpETEkIkZHREdEHBsRHUU2IGloGqNyZ/o8VtIaSU+m9+NzdRdL2i5pm6SLc+VnS9qcvrshxWRmZgNEzWSSugPfXmIbVwJbc58XAesiYjKwLn1G0hRgLjAVmAkskdTVXeYmYAEwOb1mlojHzMx6WdHLXA9I+q16Vy5pEvAB4Ou54lnA8rS8HJidK18REfsiYgewHZguaQLQERHrU2K7OdfGzMwGgKK9uc4HPinpKeBVsodiRUS8s0a7L5ONUTk2VzY+InaTrWC3pBNT+UTggVy9nansQFruXn4ESQvIzmAYxVG198rMzHpF0WTyn+tdsaTfA/ZExEOSzivSpEJZVCk/sjBiKbAUoENjPVrfzKyPFEomEfFTSe8GJkfENyWdABxTo9m5wCWS3k822LFD0reAZyVNSGclE4A9qf5O4KRc+0nArlQ+qUK5DXBx6FCp9od31fdY3DHbCvUJOcJ7L9zQULs/Orexf4azf/CHDbU75bbGzrYPP7O7YrmGDe+xTRw80NC2bPAqOmvw54CrgMWpaDjwrWptImJxREyKiFPJbqz/c0R8BFgFzE/V5gN3pOVVwFxJIyWdRnajfUO6JLZX0ozUi2tero2ZmQ0ARS9z/QFwJvAwQETsknRs9SY9uhbolHQ58DQwJ61zi6RO4HHgIHBFRHT9absQWAaMBlanl5mZDRBFk8n+iAhJASDp6Ho2EhH3Avem5V8AF/RQ7/Nk84B1L98IvKOebZqZWd8p2jW4U9LXgDGS/gRYy6929zUzs0Gs6A34L0q6EHgZeBtwdUSsaWpkZmbWMgolE0nXRcRVwJoKZWZmNsgVvcx1YYWyuseemJlZe6p6ZiJpIfBnwOmSHs19dSzw/5oZmJmZtY5al7m+TdYN93+SJmRM9kbE802LygyI/fvrqj9q/ZMNbeefmdFQu18ufLihdkvfd0tD7f50/7yG2p28qvJkEENPO7nHNod2PN3QtjzYcfCqmkwi4iXgJeCP+iYcMzNrRUXvmZiZmfXIycTMzEpzMjEzs9Jq9ebaS+Xp3rueZ9LYNK1mZtZWat2Ab3QyRzMzG0SKTvQIQHoq4qiuzxHRWP9BMzNrK0WfZ3KJpCeBHcB9wFN4GngzM0uKnpn8HTADWBsRZ0o6H489sQGm3kGOXUbd/3hD7R596vSG2q09d1pD7YYc11AzXj+p8q3NX/7HsT22Oer5Fxva1sHnX2ioHXG4sXY2YBTtzXUgPYdkiKQhEXEPcEYT4zIzsxZS9MzkRUnHAP8C3CppD9nTEM3MzAqfmcwCXgP+Ergb+Ffg95sVlJmZtZaiD8d6NfdxeZNiMTOzFlX04Vj5wYsjgOHAqx60aGZmUPzM5FcGL0qaDUxvSkRmZtZyGpqbKyK+B7y3l2MxM7MWVfQy16W5j0OAaVSes8vMzAahol2D8z23DpKNgJ/V69GY9YM41NiAuUM/3dlQu3G79zTUbsiYBm9RDh9esfiobT9vbH1VaOjQhtrFQQ9abHVF75l8vNmBmJlZ66o1Bf1XqHI5KyI+3esRmZlZy6l1A34j8BDZTMFnAU+m1xnAoeaGZmZmraLW80yWA0j6GHB+RBxIn78K/KDp0ZmZWUso2jX4LUB+rMkxqczMzKxwb65rgUck3ZM+/y5wTVMiMkuWbLyuv0Mws4KK9ub6pqTVwDmpaFFE/FvzwjIzs1ZSqzfX2yPiCUlnpaJn0vtbJL0lIh6u0nYU2ZT1I9N2VkbE5ySNBb4DnEo2XuWyiHghtVkMXE52c//TEfH9VH42sAwYDdwFXBkRHjTZpi4cMqe/Q2iqISNH1a5Uqd24nh9mVdXIEY21M6tDrTOT/wIsAL5U4bug+pQq+4D3RsQrkoYD96ezm0uBdRFxraRFwCLgKklTgLnAVLL7MWsl/XpEHAJuSnE8QJZMZuLHBluLOrzv9cba/WxXYxtU/bMmaYga2lQc9t94g1Wt3lwL0vv59a44nTm8kj4OT68gGzl/XipfDtwLXJXKV0TEPmCHpO3AdElPAR0RsR5A0s3AbJxMzMwGjEJ/skiaI+nYtPzfJd0m6cwC7YZK2gTsAdZExIPA+IjYDZDeT0zVJ/LmZTSAnalsYlruXl5pewskbZS08QD7iuyamZn1gqLnv38TEXslvRu4mOyM4qu1GkXEoYg4A5hEdpbxjirVK51XR5XySttbGhHTImLacEbWCs/MzHpJ0WTSNdr9A8BNEXEH2UOyComIF8kuZ80EnpU0ASC9d816txM4KddsErArlU+qUG5mZgNE0WTyM0lfAy4D7pI0slZbSSdIGpOWRwPvA54AVgHzU7X5wB1peRUwV9JISacBk4EN6VLYXkkzJAmYl2tjZmYDQNFBi5eRnVV8MSJeTGcUf1WjzQRguaShZImnMyLulLQe6JR0OfA0MAcgIrZI6gQeJ5vm/orUkwtgIW92DV6Nb76bmQ0oKjpcI90vmZwGMJ4AHBMRO5oaXQkdGhvn6IL+DsOs/7VC1+Dw80wGirWx8qGImFZvu6K9uT5H1n13cSoaDnyr3o2ZmVl7KnqZ6w+AM4GHASJiV1dXYTMb4Br4qz/8gAmrU9Hz3/1pEGIASDq6eSGZmVmrKZpMOlNvrjGS/gRYC3y9eWGZmVkrKTpr8BclXQi8DLwNuDoi1jQ1MjMzaxlF75mQkscaeGOalA9HxK1Ni8zMzFpGrYGHHZIWS7pR0kXK/DnwE7KxJ2ZmZjXPTG4BXgDWA58gG6g4ApgVEZuaHJuZmbWIWsnk9Ij4DQBJXweeA06OiL1Nj8zMzFpGrd5cB7oW0tQmO5xIzMysu1pnJr8p6eW0LGB0+iyy5191NDU6MzNrCbWetDi0rwIxM7PWVf8McGZmZt04mZiZWWlOJmZmVpqTiZmZleZkYmZmpTmZmJlZaU4mZmZWmpOJmZmV5mRiZmalOZmYmVlpTiZmZlaak4mZmZXmZGJmZqU5mZiZWWlOJmZmVpqTiZmZleZkYmZmpTmZmJlZaU1LJpJOknSPpK2Stki6MpWPlbRG0pPp/fhcm8WStkvaJuniXPnZkjan726QpGbFbWZm9WvmmclB4L9GxH8CZgBXSJoCLALWRcRkYF36TPpuLjAVmAkskdT1DPqbgAXA5PSa2cS4zcysTk1LJhGxOyIeTst7ga3ARGAWsDxVWw7MTsuzgBURsS8idgDbgemSJgAdEbE+IgK4OdfGzMwGgD65ZyLpVOBM4EFgfETshizhACemahOBZ3LNdqayiWm5e3ml7SyQtFHSxgPs681dMDOzKpqeTCQdA/wT8BcR8XK1qhXKokr5kYURSyNiWkRMG87I+oM1M7OGNDWZSBpOlkhujYjbUvGz6dIV6X1PKt8JnJRrPgnYlconVSg3M7MBopm9uQR8A9gaEdfnvloFzE/L84E7cuVzJY2UdBrZjfYN6VLYXkkz0jrn5dqYmdkAMKyJ6z4X+CiwWdKmVPbXwLVAp6TLgaeBOQARsUVSJ/A4WU+wKyLiUGq3EFgGjAZWp5eZmQ0QyjpItZ8OjY1zdEF/h2Fm1lLWxsqHImJave08At7MzEpzMjEzs9KcTMzMrDQnEzMzK83JxMzMSnMyMTOz0pxMzMysNCcTMzMrzcnEzMxKczIxM7PSnEzMzKw0JxMzMyvNycTMzEpzMjEzs9KcTMzMrDQnEzMzK83JxMzMSnMyMTOz0pxMzMysNCcTMzMrzcnEzMxKczIxM7PSnEzMzKw0JxMzMyvNycTMzEpzMjEzs9KcTMzMrDQnEzMzK83JxMzMSnMyMTOz0pxMzMystKYlE0n/KGmPpMdyZWMlrZH0ZHo/PvfdYknbJW2TdHGu/GxJm9N3N0hSs2I2M7PGNPPMZBkws1vZImBdREwG1qXPSJoCzAWmpjZLJA1NbW4CFgCT06v7Os3MrJ81LZlExL8Az3crngUsT8vLgdm58hURsS8idgDbgemSJgAdEbE+IgK4OdfGzMwGiGF9vL3xEbEbICJ2SzoxlU8EHsjV25nKDqTl7uUVSVpAdhYD8MraWLmttwLvQ+OA5/o7iD4wGPZzMOwjeD/bzdsaadTXyaQnle6DRJXyiiJiKbC0t4LqD5I2RsS0/o6j2QbDfg6GfQTvZ7uRtLGRdn3dm+vZdOmK9L4nle8ETsrVmwTsSuWTKpSbmdkA0tfJZBUwPy3PB+7Ilc+VNFLSaWQ32jekS2J7Jc1Ivbjm5dqYmdkA0bTLXJL+N3AeME7STuBzwLVAp6TLgaeBOQARsUVSJ/A4cBC4IiIOpVUtJOsZNhpYnV7trKUv09VhMOznYNhH8H62m4b2U1knKTMzs8Z5BLyZmZXmZGJmZqU5mfQTSTPT1DHbJS2q8P15kl6StCm9ru6POMuoNKVOt++VpsjZLulRSWf1dYxlFdjHlj+OAJJOknSPpK2Stki6skKddjieRfazpY+ppFGSNkj6UdrHv61Qp/5jGRF+9fELGAr8K3A6MAL4ETClW53zgDv7O9aS+/ke4CzgsR6+fz9ZhwoBM4AH+zvmJuxjyx/HtB8TgLPS8rHAjyv8m22H41lkP1v6mKbjc0xaHg48CMwoeyx9ZtI/pgPbI+InEbEfWEE2pUxbicpT6uTNAm6OzAPAmK5xSK2iwD62hYjYHREPp+W9wFaOnI2iHY5nkf1saen4vJI+Dk+v7j2x6j6WTib9YyLwTO5zT9PEvCudiq6WNLVvQutTRX8Ora6tjqOkU4Ezyf6izWur41llP6HFj6mkoZI2kQ0cXxMRpY/lQJlOZbApMk3Mw8ApEfGKpPcD3yMbzNlO6poup0W11XGUdAzwT8BfRMTL3b+u0KQlj2eN/Wz5YxrZOL4zJI0Bbpf0jojI3/er+1j6zKR/9DR9zBsi4uWuU9GIuAsYLmlc34XYJ2r+HFpdOx1HScPJfsHeGhG3VajSFsez1n620zGNiBeBezny0R51H0snk/7xQ2CypNMkjSB7lsuqfAVJv9b1IDBJ08mO1S/6PNLmWgXMSz1HZgAvRZpVul20y3FM+/ANYGtEXN9DtZY/nkX2s9WPqaQT0hkJkkYD7wOe6Fat7mPpy1z9ICIOSvpz4PtkPbv+MbIpZT6Zvv8q8EFgoaSDwGvA3EjdLFqFKk+pMxze2Me7yHqNbAd+CXy8fyJtXIF9bPnjmJwLfBTYnK61A/w1cDK0z/Gk2H62+jGdACxX9gDCIUBnRNzZ7fdP3cfS06mYmVlpvsxlZmalOZmYmVlpTiZmZlaak4mZmZXmZGJmZqU5mVi/k3QoNwPrJlWYRbkXt3WepDsH+jpz6/6YpBvT8iclzeul9f5OmjF2UxprUE/b2ZKm9EYc1j48zsQGgtci4oxqFSQNjTcf5YykYRFxsNaKi9ZrBan/f2/5MPDFiPhmA21nA3eSPWa7kHY6DlaZz0xswJL0lKSrJd0PzJF0r6QvSLoPuFLSKZLWpectrJN0cmq3TNL1ku4Brquy/qOVPY/kh5IekTQrlT+Yn7wvbffsnupXWf9UZc+N2JRinJzK56XPP5J0Syr7/bTdRyStlTS+wvqukfSZXEzXpfX/WNLvpPKjJHWm9X8nrXNat/V8ArgMuFrSrZKOST+/hyVtzu9X91gl/TZwCfD3ab/eKukMSQ+kerdLOj4X4xvHq9rPytpAf8+t75dfwCFgU+71oVT+FPDZXL17gSW5z/8HmJ+W/xj4XlpeRvaX89AK2zqP9CwK4AvAR9LyGLJnVxwN/CXwt6l8AvDjGvXfWGe3bX0F+HBaHgGMBqYC24BxqXxsej+eNwcRfwL4Ulr+GHBjWr4G+EzuZ9FV5/3A2rT8GeBrafkdwEFgWoXYlgEfTMvDgI60PI5s1LOqxPpG2/T5UeB30/L/AL5c6Xj51d4vX+aygaDaZa7vVPn8LuDStHwL8L9y3303cpfFenARcEnXX/vAKLJpMzqBNWRTo1wGfLdG/Z6sB/6bpEnAbRHxpKT3Aisj4jmAiOh6Fsok4DvKnhkxAthRI3aArkkIHwJOTcvvBv4hrfsxSY8WWI+AL0h6D3CYbKrx8UBPsb7ZUDoOGBMR96Wi5bz584Ijj5+1KScTG+herfE5Lz83ULV6XQT8YURsO+IL6ReS3gl8CPjTavUrXZICiIhvS3oQ+ADw/XR5SVSeyvsrwPURsUrSeWRnIbXsS++HePP/cqWpw2v5MHACcHZEHJD0FFmi7CnWehQ5DtYGfM/EWtn/J5txGbJfiPfX2f77wKekN2aAPTP33Qrgs8BxEbG5QP0jSDod+ElE3EA2C+s7gXXAZZL+Q6ozNlU/DvhZWp5f537k3U92NkXqcfUbBdocB+xJieR84JRU3lOse8keaUtEvAS80HXPhmySxPuwQcfJxAaC0d26Bl9bsN2ngY+nSzkfpf6bvH9HNsPvo5IeS5+7rCRLVJ0F61fyIeAxZbPPvp3sMahbgM8D90n6EdA1zfk1wHcl/V/guTr3I28JcEL6mVxFdj/jpRptbgWmSdpIlpSfAKgS6wrgr1JngbeSJb+/T9s8g+y+iQ0ynjXYrI0om1Z8eES8nn7RrwN+PSL293No1uZ8z8SsvRwF3KPsaYECFjqRWF/wmYmZmZXmeyZmZlaak4mZmZXmZGJmZqU5mZiZWWlOJmZmVtq/AzDvHHS3250/AAAAAElFTkSuQmCC\n",
Maximilian Schanner's avatar
Maximilian Schanner committed
259
260
261
262
263
264
265
266
267
268
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
269
270
271
272
273
274
275
276
277
278
   "source": [
    "from matplotlib.patches import Rectangle\n",
    "fig, ax = plt.subplots()\n",
    "ax.pcolormesh(ER, RE, marginal.T)\n",
    "\n",
    "rect = Rectangle((moments[1, 0] - 2*np.sqrt(moments[1, 1]),\n",
    "                  moments[2, 0] - 2*np.sqrt(moments[2, 1])),\n",
    "                 4*np.sqrt(moments[1, 1]), 4*np.sqrt(moments[2, 1]), \n",
    "                 fill=True, facecolor=(1, 1, 1, 0.2),\n",
    "                 lw=2, ec='white')\n",
279
280
281
    "ax.add_patch(rect)\n",
    "ax.set_xlabel('Error level scaling factor')\n",
    "ax.set_ylabel('Residual term [nT]')"
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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",
301
   "version": "3.7.3"
302
303
304
305
306
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}