Conditional_Expectation_Gaussian.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Introduction\n",
      "-----------------\n",
      "\n",
      "By this point, we have developed many tools to deal with computing the conditional expectation. In this section, we discuss a bizarre and amazing coincidence regarding Gaussian random variables and linear projection, a coincidence that is the basis for most of statistical signal processing."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Conditional Expectation by Optimization"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, let's consider the important case of the zero-mean bivariate Gaussian and try to find a  function $h$ that minimizes the mean squared error (MSE). Again,  trying to solve for the conditional expectation by minimizing the error over all possible functions $h$ is generally very, very hard. One alternative is to use parameters for the $h$ function and then just optimize over those. For example, we could assume that $h(Y)= \\alpha Y$ and then use calculus to find the $\\alpha$ parameter.\n",
      "\n",
      "Let's try this with the zero-mean bivariate Gaussian density,\n",
      "\n",
      "$$\\mathbb{E}((X-\\alpha Y )^2) = \\mathbb{E}(\\alpha^2 Y^2 - 2 \\alpha X Y + X^2 )$$\n",
      "\n",
      "and then differentiate this with respect to $\\alpha$ to obtain\n",
      "\n",
      "$$\\mathbb{E}(2 \\alpha Y^2 - 2 X Y  ) = 2 \\alpha \\sigma_y^2 - 2 \\mathbb{E}(XY) = 0$$\n",
      "\n",
      "Then, solving for $\\alpha$ gives\n",
      "\n",
      "$$ \\alpha = \\frac{ \\mathbb{E}(X Y)}{ \\sigma_y^2 } $$\n",
      "\n",
      "which means we that\n",
      "\n",
      "\\begin{equation}\n",
      "\\mathbb{ E}(X|Y) \\approx \\alpha Y =   \\frac{ \\mathbb{E}(X Y )}{ \\sigma_Y^2 } Y =\\frac{\\sigma_{X Y}}{ \\sigma_Y^2 } Y \n",
      "\\end{equation}\n",
      "\n",
      "where that last equality is just notation. Remember here we assumed a special linear form for $h=\\alpha Y$, but we did that for convenience. We still don't know whether or not this is the one true $h_{opt}$ that minimizes the MSE for all such functions.\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Conditional Expectation Using Direct Integration"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, let's try this again by computing  $ \\mathbb{E}(X|Y)$ in the case of the bivariate Gaussian distribution straight from the definition.\n",
      "\n",
      "\\begin{equation}\n",
      "\\mathbb{E}(X|Y)  = \\int_{\\mathbb{ R}} x \\frac{f_{X,Y}(x,y)}{f_Y(y)} dx\n",
      "\\end{equation}\n",
      "\n",
      "where \n",
      "\n",
      "$$ f_{X,Y}(x,y) = \\frac{1}{2\\pi |\\mathbf{R}|^{\\frac{1}{2}}} e^{-\\frac{1}{2} \\mathbf{v}^T \\mathbf{R}^{-1} \\mathbf{v} } $$ \n",
      "\n",
      "and where\n",
      "\n",
      "$$ \\mathbf{v}= \\left[ x,y \\right]^T$$ \n",
      "\n",
      "$$ \\mathbf{R} = \\left[ \\begin{array}{cc}\n",
      "\\sigma_{x}^2 & \\sigma_{xy}  \\\\\\\\\n",
      "\\sigma_{xy}  & \\sigma_{y}^2 \\\\\\\\\n",
      "\\end{array} \\right] $$ \n",
      "\n",
      "and with\n",
      "\n",
      "\\begin{eqnarray}\n",
      " \\sigma_{xy} &=& \\mathbb{E}(xy)   \\nonumber    \\\\\\\\\n",
      " \\sigma_{x}^2 &=& \\mathbb{E}(x^2) \\nonumber    \\\\\\\\ \n",
      " \\sigma_{y}^2 &=& \\mathbb{E}(y^2) \\nonumber      \n",
      "\\end{eqnarray}\n",
      "\n",
      "This conditional expectation (Eq. 4 above) is a tough integral to evaluate, so we'll do it with `sympy`.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sympy import Matrix, Symbol, exp, pi, simplify, integrate \n",
      "from sympy import stats, sqrt, oo, use\n",
      "from sympy.abc import y,x\n",
      "\n",
      "sigma_x = Symbol('sigma_x',positive=True)\n",
      "sigma_y = Symbol('sigma_y',positive=True)\n",
      "sigma_xy = Symbol('sigma_xy',real=True)\n",
      "fyy = stats.density(stats.Normal('y',0,sigma_y))(y)\n",
      " \n",
      "R = Matrix([[sigma_x**2, sigma_xy],\n",
      "            [sigma_xy,sigma_y**2]])\n",
      "fxy = 1/(2*pi)/sqrt(R.det()) * exp((-Matrix([[x,y]])*R.inv()* Matrix([[x],[y]]))[0,0]/2 )\n",
      "\n",
      "fcond = simplify(fxy/fyy)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Unfortunately, `sympy` cannot immediately integrate this without some hints. So, we need to define a positive variable ($u$) and substitute it into the integration"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "u=Symbol('u',positive=True) # define positive variable\n",
      "\n",
      "fcond2=fcond.subs(sigma_x**2*sigma_y**2-sigma_xy**2,u) # substitute as hint to integrate\n",
      "g=simplify(integrate(fcond2*x,(x,-oo,oo))) # evaluate integral\n",
      "gg=g.subs( u,sigma_x**2 *sigma_y**2 - sigma_xy**2 ) # substitute back in\n",
      "use( gg, simplify,level=2) # simplify exponent term"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "sigma_xy*y/sigma_y**2"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Thus, by direct integration using `sympy`, we found\n",
      "\n",
      "$$ \\mathbb{ E}(X|Y) = Y \\frac{\\sigma_{xy}}{\\sigma_{y}^{2}} $$ \n",
      "\n",
      "and this matches the prior result we obtained by direct minimization by assuming that $\\mathbb{E}(X|Y) = \\alpha Y$ and then solving for the optimal $\\alpha$!\n",
      "\n",
      "The importance of this result cannot be understated: the one true and optimal $h_{opt}$ *is a linear function* of $Y$. \n",
      "\n",
      "In other words, assuming a linear function, which made the direct search for an optimal $h(Y)$ merely convenient yields the optimal result! This is  a general result that extends for *all* Gaussian problems. The link between linear functions and optimal estimation of Gaussian random variables is the most fundamental result in statistical signal processing! This fact is exploited in everything from optimal filter design  to adaptive signal processing.\n",
      "\n",
      "We can easily extend this result to non-zero mean problems by inserting the means in the right places as follows:\n",
      "\n",
      "$$ \\mathbb{ E}(X|Y) = \\bar{X} + (Y-\\bar{Y}) \\frac{\\sigma_{xy}}{\\sigma_{y}^{2}}  $$\n",
      "\n",
      "where $\\bar{X}$ is the mean of $X$ (same for $Y$)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Summary\n",
      "-------------\n",
      "\n",
      "In this section, we showed that the conditional expectation for Gaussian random variables is a linear function, which, by a bizarre coincidence, is also the easiest one to work with. This result is fundamental to all optimal linear filtering problems (e.g. Kalman filter) and is the basis of most of the theory of stochastic processes used in signal processing. Up to this point, we have worked hard to illustrate all of the concepts we will need to unify our understanding of this entire field and figured out multiple approaches to these kinds of problems, most of which are far more difficult to compute. Thus, it is indeed just plain lucky that the most powerful distribution is the easiest to compute as a conditional expectation because it is a linear function. We will come back to this same result again and again as we work our way through these greater concepts."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### References \n",
      "\n",
      "This post was created using the [nbconvert](https://github.com/ipython/nbconvert) utility from the source [IPython Notebook](www.ipython.org) which is available for [download](https://github.com/unpingco/Python-for-Signal-Processing/blob/master/Conditional_Expectation_Gaussian.ipynb) from the main github [site](https://github.com/unpingco/Python-for-Signal-Processing) for this blog. "
     ]
    }
   ],
   "metadata": {}
  }
 ]
}