{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Bayesian signal reconstruction with Gaussian Random Fields (Wiener Filtering)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Florent Leclercq,
\n",
"Imperial Centre for Inference and Cosmology, Imperial College London,
\n",
"florent.leclercq@polytechnique.org"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.linalg\n",
"from matplotlib import pyplot as plt\n",
"np.random.seed(123456)\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"N=1000\n",
"epsilon=0.000000000001"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup signal covariance"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"t=np.array([100./(f*f) for f in range(1,N//2)])\n",
"signalcovar=np.concatenate([np.zeros(1),t,t[::-1],np.zeros(1)])"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAExNJREFUeJzt3X+MXWWdx/H3tzNQbFlsCwOpLVhYG5TVuOCEBdlsXPG3rpAsbiDu2rgk/YddUdko7G6Cm/1jNXEFTTZEAmp3Y0BEIoQ1uqRijFm361T8AVS2FaFUKh2hUCmlv+a7f8yZ9rbMdGbuuXPuPIf3K5nce8597jzPmdN+7nOf5/yIzESS1F4L+t0ASdLcMuglqeUMeklqOYNeklrOoJekljPoJanlDHpJajmDXpJazqCXpJYb7HcDAE455ZRctWpVv5shSUXZuHHjbzNzaLpy8yLoV61axcjISL+bIUlFiYjHZ1LOoRtJajmDXpJazqCXpJYz6CWp5Qx6SWq5aYM+Ir4UETsi4sGOdcsi4r6I2Fw9Lq3WR0R8ISK2RMTPIuK8uWy8JGl6M+nRfwV411HrrgXWZ+ZqYH21DPBuYHX1sxa4qTfNlCR1a9qgz8zvA88ctfoSYF31fB1wacf6f89x/wMsiYjlvWrsJG3jGxu3sWffwbmqQpLmxK4X93PPT59spK5ux+hPy8ztANXjqdX6FcATHeW2VeteIiLWRsRIRIyMjo521YgfPvo013z9p/zzfz7c1fslqV+uueOnfOS2B9j81O/mvK5eT8bGJOsmvft4Zt6cmcOZOTw0NO0ZvJN6/sUDAOzYtber90tSv2x/bg8AL+4fm/O6ug36pyaGZKrHHdX6bcDpHeVWAs18N5EkTarboL8HWFM9XwPc3bH+Q9XRNxcAz00M8UiS+mPai5pFxG3AW4BTImIbcD3waeCOiLgS2Ap8oCr+LeA9wBbgBeDDc9DmSUw6OiRJ81Y2GFvTBn1mXjHFSxdPUjaBq+o2aqYiJpsSkKRyNBFjnhkrSS1n0EtSyxn0ktRyrQj6Jic1JKkXmsytVgS9JGlqBr0ktZxBL0ktZ9BLUssVHfSeLiVJ0ys66CVJ0zPoJanlWhH0HkYvqTRN5lYrgl6SNDWDXpJazqCXpJYz6CWp5YoOeu87Iql03nhEklRbK4I+vU6xpMI0mVutCHpJ0tQMeklquaKD3slYSaWLBi7PWHTQS5KmZ9BLUssZ9JLUcga9JLVcK4Leo+glaWqtCHpJ0tQMeklqOYNeklquVtBHxMci4qGIeDAibouIEyLizIjYEBGbI+JrEXF8rxorSZq9roM+IlYAHwGGM/P1wABwOfAZ4IbMXA3sBK7sRUMnbUMDZ5RJ0lwq4TLFg8ArImIQWARsB94K3Fm9vg64tGYdkqQaug76zPw18FlgK+MB/xywEXg2Mw9UxbYBK+o2UpLUvTpDN0uBS4AzgVcBi4F3T1J00sPcI2JtRIxExMjo6Gi3zZAkTaPO0M3bgF9l5mhm7gfuAt4MLKmGcgBWAk9O9ubMvDkzhzNzeGhoqEYzJEnHUifotwIXRMSiiAjgYuBh4H7gsqrMGuDuek2cnjeYklSaJnOrzhj9BsYnXX8M/Lz6XTcDnwQ+HhFbgJOBW3vQzsl50I2kwjVx1M3g9EWmlpnXA9cftfpR4Pw6v1eS1DutODPWkRtJpckGk6sVQS9JmppBL0ktZ9BLUssZ9JLUckUHvUdXStL0ig56SSpdEydOGfSS1HIGvST1kT16SVJtrQj69KpmkgpTxEXNJEn1NXEpBINeklrOoJekPnIyVpJUW9FBH01csV+SCld00EuSpmfQS1IfNHlQuEEvSX3kZKwkqbaig96pWEmaXtFBL0ml88xYSVJtBr0k9ZGTsZKk2loR9F6lWFJpmry8eiuCXpJK1UTcG/SS1HJFB73XNJOk6RUd9JJUuibG6g16SWq5WkEfEUsi4s6I+EVEbIqICyNiWUTcFxGbq8elvWrs0TzaRlLpSpiM/Tzw7cx8LfBGYBNwLbA+M1cD66tlSVKfdB30EXES8CfArQCZuS8znwUuAdZVxdYBl9Zt5HSauFaEJPVSKdejPwsYBb4cEQ9ExC0RsRg4LTO3A1SPp0725ohYGxEjETEyOjraVQOMd0mlm++XQBgEzgNuysxzgd3MYpgmM2/OzOHMHB4aGqrRDEnSsdQJ+m3AtszcUC3fyXjwPxURywGqxx31mji1Jk8hlqS5MY8Pr8zM3wBPRMTZ1aqLgYeBe4A11bo1wN21WihJqmWw5vv/FvhqRBwPPAp8mPEPjzsi4kpgK/CBmnVMy469pOI0mFu1gj4zfwIMT/LSxXV+74zrb6ISSZpD830yVpJUgLKD3i69pMKVcGasJGmeM+glqeWKDnovfSCpdE7GSpJqKzroPX5eUum88YgkqbZWBL09e0mlKeUyxX1nwEsqncfRS5JqKzro7dBLKp2HV0qSajPoJanlig567zAlqXRNnOFfdNBP8FIIkkrTZEe16KA33iUVz8lYSVJdBr0ktVzRQe9crKTSeWasJKm2VgS9PXtJpfGiZjNmwksqm5dAkCTVVnTQO2QjqXSeGStJqs2gl6SWKzroHbmRVDonYyVJtRUd9E7GSipdEWfGRsRARDwQEfdWy2dGxIaI2BwRX4uI4+s3U5LUrV706K8GNnUsfwa4ITNXAzuBK3tQxzHZsZdUmiZHJGoFfUSsBN4L3FItB/BW4M6qyDrg0jp1HIs3HJFUuiZuQFK3R38j8AlgrFo+GXg2Mw9Uy9uAFTXrkCTV0HXQR8T7gB2ZubFz9SRFJ/24ioi1ETESESOjo6NdtcHJWEmlm++TsRcB74+Ix4DbGR+yuRFYEhGDVZmVwJOTvTkzb87M4cwcHhoaqtEMSdKxdB30mXldZq7MzFXA5cB3M/ODwP3AZVWxNcDdtVs5bWPmvAZJ6qkm5xjn4jj6TwIfj4gtjI/Z3zoHdQDmu6QWaCDIBqcvMr3M/B7wver5o8D5vfi9kqT6Cj8z1j69pLJ5mWJJUm0GvSS1nEEvSX3kZYolSbUVHfTOxUoqnT16SVJtrQh6r2IpqTTFXKa43wx4SaWb7xc1kyQVwKCXpJYrOug96kZS6Uq4w5QkaZ4rOujt0UsqnZOxM2TgSyqNh1dKknqm6KC3Iy+pdF4CQZJUW9FB7x2mJJXPwyslSTW1Iujt10vS1IoOegNeUumcjJUk1VZ20Null1Q4z4ydIY++kaSptSLox8x5SYU52GBwFR30E3eYskcvqTRjOZFfc19X0UE/wR69pNI0mVtFB/3EJ+GYPXpJhZkYiWji3tdFB/2EJse6JKkXDjbYQW1F0Nujl1SaIiZjI+L0iLg/IjZFxEMRcXW1fllE3BcRm6vHpb1r7pEm/kx26CWVZqJ/Ot8nYw8A12Tm64ALgKsi4hzgWmB9Zq4G1lfLc2rMpJdUmCJ69Jm5PTN/XD3/HbAJWAFcAqyriq0DLq3byKnbMP7o0I2k0hw6vLKBunoyRh8Rq4BzgQ3AaZm5HcY/DIBTp3jP2ogYiYiR0dHRWvXboZdUmqLuGRsRJwLfAD6ambtm+r7MvDkzhzNzeGhoqFYbPOpGUmmKOeomIo5jPOS/mpl3Vaufiojl1evLgR31mjg1z4yVVKrDZ8bO4+PoIyKAW4FNmfm5jpfuAdZUz9cAd3ffvJlp8pNRknqhydgarPHei4C/An4eET+p1v098Gngjoi4EtgKfKBeE6d3cGyua5CkcnUd9Jn5AyCmePnibn/v7NpwqC1NVCdJRfLMWElquaKD3jNjJZVuvp8ZO294ZqwkTa0VQe9RN5I0tbKDvgp4x+gllcrr0c/QmIdXStKUig76w5Ox9ugllcnJ2Bky6CVpai0J+n63QJLmr6KD3o68pNI5dCNJqq3ooPcaN5JK1HmSZzF3mJIkzdz+ho8Jb03QexkESaXYe8Cgn7HOaN/nReklFWLv/sN5Na/vMDXfvLj/YL+bIEkz0nReFR30nR+ETX8VkqRudeaVk7Gz0PlVSJLms70H7NF3pek/nCR1y8nYWej8yuPQjaRSHDEC4ZmxM2ePXlIpHLqZhc7Dkvbss0cvqQydR91445FZ2PXi/n43QZJmZNeeA43W15qgf26PQS+pDE3nVWuC/tkXDHpJZXh2z75Dz71M8Sx0/uEkaT6zR9+Fk04YZHTX3n43Q5Jm5Klde/m9EwYbq6/ooJ/4ynP6skVs27mnv42RpBnatnMPpy9dBHgJhBlbdfJitow+741IJM17B8eSrU/v5oxlixqrs+ignzj+dHjVUp7ZvY/NO57vc4sk6di27Hie3fsO8oaVrwQKnoyNiHdFxCMRsSUirp2LOjpd+PsnA/C9R3bMdVWSVMu3H/wNcDi3mtDzoI+IAeDfgHcD5wBXRMQ5va6n0+lLF3H+mcu4+fuP8uSzjtVLmp8ef3o3t/zgUf707CFWLHlFY/XOxbTv+cCWzHwUICJuBy4BHu51RZ1feT71Z3/AX3zxh7zzhu/z529ayZtevZTTly3i5MXHs3jhIIsXDrBwcKDXTZCkQzKTsYQDY2O8sPcgz+7Zz84X9rH16RfY+PhOvvnArxkYCP7xfYf7vk1cAmEugn4F8ETH8jbgj+agnkMi4JxXncQ3r7qIf/2vR7j9R1v5yn8/9pJyAwuCgQgWLKB6DAYXBAMLggUx/hNR/c4jfn8c8Xs6Fw+Xj0nWvfT98ZInDWlonrqp6fCmJt6b256G6mlsixrcpsbqSQ5mcnAMDo6NcXDscKiPjVG9NnVjFg4u4O3nnMbfveNsVp2ymB27Xmym4cxN0E8WYS/Z+ohYC6wFOOOMM7qq6KyhE3nvG5azoArS15x6Ijf95Zt4cf9BHnt6N088s4edL+zjhb0H2L3vIHv2HeRgJmNjyYGxrHbU+OPEczjyH87E04l1R/xHySPLjJfLSd935Lo8Yl1TmX/0B9ac1dNILUd+4M5pPc1U07r902Rl0VBFgwuO7CAOdHQUJ14biGBgASxeOMgrX3EcSxYdx4olizhraDHHDRweLV84OMB737C8kaNvotc9o4i4EPhUZr6zWr4OIDP/Zar3DA8P58jISE/bIUltFxEbM3N4unJzcdTNj4DVEXFmRBwPXA7cMwf1SJJmoOdDN5l5ICL+BvgOMAB8KTMf6nU9kqSZmZOLLWTmt4BvzcXvliTNTtFnxkqSpmfQS1LLGfSS1HIGvSS1nEEvSS3X8xOmumpExCjweJdvPwX4bQ+bUwK3+eXBbX55qLPNr87MoekKzYugryMiRmZyZlibuM0vD27zy0MT2+zQjSS1nEEvSS3XhqC/ud8N6AO3+eXBbX55mPNtLn6MXpJ0bG3o0UuSjqHooG/6JuRNiYjTI+L+iNgUEQ9FxNXV+mURcV9EbK4el1brIyK+UP0dfhYR5/V3C7oTEQMR8UBE3FstnxkRG6rt/Vp12WsiYmG1vKV6fVU/292tiFgSEXdGxC+qfX3hy2Aff6z6N/1gRNwWESe0cT9HxJciYkdEPNixbtb7NiLWVOU3R8SabttTbND34ybkDToAXJOZrwMuAK6qtu1aYH1mrgbWV8sw/jdYXf2sBW5qvsk9cTWwqWP5M8AN1fbuBK6s1l8J7MzM1wA3VOVK9Hng25n5WuCNjG97a/dxRKwAPgIMZ+brGb+M+eW0cz9/BXjXUetmtW8jYhlwPeO3Yj0fuH7iw2HWMrPIH+BC4Dsdy9cB1/W7XXO0rXcDbwceAZZX65YDj1TPvwhc0VH+ULlSfoCV1T/+twL3Mn4Tut8Cg0fvb8bvdXBh9XywKhf93oZZbu9JwK+ObnfL9/HE/aSXVfvtXuCdbd3PwCrgwW73LXAF8MWO9UeUm81PsT16Jr8J+Yo+tWXOVF9XzwU2AKdl5naA6vHUqlgb/hY3Ap8Axqrlk4FnM/NAtdy5TYe2t3r9uap8Sc4CRoEvV8NVt0TEYlq8jzPz18Bnga3Adsb320bavZ87zXbf9myflxz0M7oJecki4kTgG8BHM3PXsYpOsq6Yv0VEvA/YkZkbO1dPUjRn8FopBoHzgJsy81xgN4e/yk+m+G2uhh0uAc4EXgUsZnzY4mht2s8zMdV29mz7Sw76bcDpHcsrgSf71Jaei4jjGA/5r2bmXdXqpyJiefX6cmBHtb70v8VFwPsj4jHgdsaHb24ElkTExF3QOrfp0PZWr78SeKbJBvfANmBbZm6olu9kPPjbuo8B3gb8KjNHM3M/cBfwZtq9nzvNdt/2bJ+XHPStvQl5RARwK7ApMz/X8dI9wMTM+xrGx+4n1n+omr2/AHhu4itiCTLzusxcmZmrGN+P383MDwL3A5dVxY7e3om/w2VV+aJ6epn5G+CJiDi7WnUx8DAt3ceVrcAFEbGo+jc+sc2t3c9Hme2+/Q7wjohYWn0beke1bvb6PWFRc7LjPcD/Ab8E/qHf7enhdv0x41/Rfgb8pPp5D+Pjk+uBzdXjsqp8MH4E0i+BnzN+VEPft6PLbX8LcG/1/Czgf4EtwNeBhdX6E6rlLdXrZ/W73V1u6x8CI9V+/iawtO37GPgn4BfAg8B/AAvbuJ+B2xifh9jPeM/8ym72LfDX1fZvAT7cbXs8M1aSWq7koRtJ0gwY9JLUcga9JLWcQS9JLWfQS1LLGfSS1HIGvSS1nEEvSS33/5+XUQNFYYnyAAAAAElFTkSuQmCC\n",
"text/plain": [
"