data-science/mvn-conditional.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Intro\n",
    "\n",
    "In this notebook we will learn about the [conditional multivariate normal (MVN) distribution](https://en.wikipedia.org/wiki/Multivariate_normal_distribution). In particular, we want to estimate the expected value (or the mean) of some subset of variables given that another subset has been conditioned on. Though the notation is quasi-dense, it is not terribly difficult to produce a conditional MVN from a marginal MVN distribution. \n",
    "\n",
    "# Case 1\n",
    "\n",
    "* $X_0 \\rightarrow X_1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.random import normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10000, 2)\n",
      "(2,)\n",
      "(2, 2)\n",
      "mean [-0.01217851  0.97154671]\n",
      "cov [[0.99624904 1.98586371]\n",
      " [1.98586371 4.94146349]]\n"
     ]
    }
   ],
   "source": [
    "N = 10000\n",
    "x0 = normal(0, 1, N)\n",
    "x1 = normal(1 + 2 * x0, 1, N)\n",
    "\n",
    "X = np.hstack([x0.reshape(-1, 1), x1.reshape(-1, 1)])\n",
    "M = np.mean(X, axis=0)\n",
    "S = np.cov(X.T)\n",
    "\n",
    "print(X.shape)\n",
    "print(M.shape)\n",
    "print(S.shape)\n",
    "print('mean', M)\n",
    "print('cov', S)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.20168259547735531"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M[0] + S[0,1] / S[1,1] * (0.5 - M[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.9924929599471097"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M[1] + S[1,0] / S[0,0] * (0.5 - M[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.19817481229753187"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S[0,0] - S[0,1] / S[1,1] * S[1,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9829606417418284"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S[1,1] - S[1,0] / S[0,0] * S[1,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Case 2\n",
    "\n",
    "* $X_0 \\rightarrow X_1 \\rightarrow X_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import namedtuple\n",
    "from numpy.linalg import inv\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "COV = namedtuple('COV', 'C11 C12 C21 C22 C22I')\n",
    "\n",
    "def to_row_indices(indices):\n",
    "    return [[i] for i in indices]\n",
    "\n",
    "def to_col_indices(indices):\n",
    "    return indices\n",
    "\n",
    "def get_covariances(i1, i2, S):\n",
    "    r = to_row_indices(i1)\n",
    "    c = to_col_indices(i1)\n",
    "    C11 = S[r,c]\n",
    "    \n",
    "    r = to_row_indices(i1)\n",
    "    c = to_col_indices(i2)\n",
    "    C12 = S[r,c]\n",
    "    \n",
    "    r = to_row_indices(i2)\n",
    "    c = to_col_indices(i1)\n",
    "    C21 = S[r,c]\n",
    "    \n",
    "    r = to_row_indices(i2)\n",
    "    c = to_col_indices(i2)\n",
    "    C22 = S[r,c]\n",
    "    \n",
    "    C22I = inv(C22)\n",
    "    \n",
    "    return COV(C11, C12, C21, C22, C22I)\n",
    "\n",
    "def compute_means(a, M, C, i1, i2):\n",
    "    a = np.array([2.0])\n",
    "    return M[i1] + C.C12.dot(C.C22I).dot(a - M[i2])\n",
    "\n",
    "def compute_covs(C):\n",
    "    return C.C11 - C.C12.dot(C.C22I).dot(C.C21)\n",
    "\n",
    "def update_mean(m, a, M, i1, i2):\n",
    "    v = np.copy(M)\n",
    "    for i, mu in zip(i1, m):\n",
    "        v[i] = mu\n",
    "    for i, mu in zip(i2, a):\n",
    "        v[i] = mu\n",
    "    return v\n",
    "\n",
    "def update_cov(c, S, i1, i2):\n",
    "    m = np.copy(S)\n",
    "    rows, cols = c.shape\n",
    "    for row in range(rows):\n",
    "        for col in range(cols):\n",
    "            m[i1[row],i1[col]] = c[row,col]\n",
    "    for i in i2:\n",
    "        m[i,i] = 0.01\n",
    "    return m\n",
    "\n",
    "def update_mean_cov(v, iv, M, S):\n",
    "    if v is None or iv is None or len(v) == 0 or len(iv) == 0:\n",
    "        return np.copy(M), np.copy(S)\n",
    "    i2 = iv.copy()\n",
    "    i1 = [i for i in range(S.shape[0]) if i not in i2]\n",
    "    \n",
    "    C = get_covariances(i1, i2, S)\n",
    "    m = compute_means(v, M, C, i1, i2)\n",
    "    c = compute_covs(C)\n",
    "    M_u = update_mean(m, v, M, i1, i2)\n",
    "    S_u = update_cov(c, S, i1, i2)\n",
    "    return M_u, S_u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean [-0.01313569  0.96903015  2.93044396]\n",
      ">\n",
      "cov [[ 0.99331094  1.96882311  3.92097634]\n",
      " [ 1.96882311  4.87809341  9.73351352]\n",
      " [ 3.92097634  9.73351352 20.43568728]]\n",
      ">\n",
      "corr [[1.         0.89441492 0.87027596]\n",
      " [0.89441492 1.         0.9748773 ]\n",
      " [0.87027596 0.9748773  1.        ]]\n"
     ]
    }
   ],
   "source": [
    "N = 10000\n",
    "x0 = normal(0, 1, N)\n",
    "x1 = normal(1 + 2 * x0, 1, N)\n",
    "x2 = normal(1 + 2 * x1, 1, N)\n",
    "\n",
    "X = np.hstack([x0.reshape(-1, 1), x1.reshape(-1, 1), x2.reshape(-1, 1)])\n",
    "M = np.mean(X, axis=0)\n",
    "S = np.cov(X.T)\n",
    "\n",
    "print('mean', M)\n",
    "print('>')\n",
    "print('cov', S)\n",
    "print('>')\n",
    "print('corr', np.corrcoef(X.T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean [0.40296894 2.         4.98759174]\n",
      ">\n",
      "cov [[ 1.98683996e-01  1.96882311e+00 -7.51882725e-03]\n",
      " [ 1.96882311e+00  1.00000000e-02  9.73351352e+00]\n",
      " [-7.51882725e-03  9.73351352e+00  1.01390146e+00]]\n",
      ">\n",
      "corr [[ 1.         -0.02477834  0.77462271]\n",
      " [-0.02477834  1.          0.05241713]\n",
      " [ 0.77462271  0.05241713  1.        ]]\n"
     ]
    }
   ],
   "source": [
    "M_u, S_u = update_mean_cov(np.array([2.0]), [1], M, S)\n",
    "\n",
    "print('mean', M_u)\n",
    "print('>')\n",
    "print('cov', S_u)\n",
    "print('>')\n",
    "print('corr', np.corrcoef(np.random.multivariate_normal(M_u, S_u, N*10).T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Case 2\n",
    "\n",
    "* $X_0 \\leftarrow X_1 \\rightarrow X_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean [ 0.96073372 -0.01011371  0.98507138]\n",
      ">\n",
      "cov [[16.71912352  3.92724763  7.87120405]\n",
      " [ 3.92724763  0.9809372   1.96828134]\n",
      " [ 7.87120405  1.96828134  4.9446487 ]]\n",
      ">\n",
      "corr [[1.         0.96975255 0.86569856]\n",
      " [0.96975255 1.         0.8937146 ]\n",
      " [0.86569856 0.8937146  1.        ]]\n"
     ]
    }
   ],
   "source": [
    "N = 10000\n",
    "\n",
    "x1 = normal(0, 1, N)\n",
    "x0 = normal(1 + 4.0 * x1, 1, N)\n",
    "x2 = normal(1 + 2.0 * x1, 1, N)\n",
    "\n",
    "X = np.hstack([x0.reshape(-1, 1), x1.reshape(-1, 1), x2.reshape(-1, 1)])\n",
    "M = np.mean(X, axis=0)\n",
    "S = np.cov(X.T)\n",
    "\n",
    "print('mean', M)\n",
    "print('>')\n",
    "print('cov', S)\n",
    "print('>')\n",
    "print('corr', np.corrcoef(X.T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean [9.00835823 2.         5.01842773]\n",
      ">\n",
      "cov [[ 0.99612527  3.92724763 -0.00894181]\n",
      " [ 3.92724763  0.01        1.96828134]\n",
      " [-0.00894181  1.96828134  0.99523029]]\n",
      ">\n",
      "corr [[1.         0.11353019 0.5455407 ]\n",
      " [0.11353019 1.         0.08049836]\n",
      " [0.5455407  0.08049836 1.        ]]\n"
     ]
    }
   ],
   "source": [
    "M_u, S_u = update_mean_cov(np.array([2.0]), [1], M, S)\n",
    "\n",
    "print('mean', M_u)\n",
    "print('>')\n",
    "print('cov', S_u)\n",
    "print('>')\n",
    "print('corr', np.corrcoef(np.random.multivariate_normal(M_u, S_u, N*10).T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Case 2\n",
    "\n",
    "* $X_0 \\rightarrow X_1 \\leftarrow X_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean [3.11610448e-03 1.00024244e+00 4.90651340e-05]\n",
      ">\n",
      "cov [[ 9.88652679e-01  1.94567309e+00 -7.66954037e-03]\n",
      " [ 1.94567309e+00  1.36958157e+01  2.92099886e+00]\n",
      " [-7.66954037e-03  2.92099886e+00  9.73747382e-01]]\n",
      ">\n",
      "corr [[ 1.          0.52875446 -0.00781672]\n",
      " [ 0.52875446  1.          0.79986056]\n",
      " [-0.00781672  0.79986056  1.        ]]\n"
     ]
    }
   ],
   "source": [
    "N = 10000\n",
    "\n",
    "x0 = normal(0, 1, N)\n",
    "x2 = normal(0, 1, N)\n",
    "x1 = normal(1 + 2 * x0 + 3 * x2, 1, N)\n",
    "\n",
    "X = np.hstack([x0.reshape(-1, 1), x1.reshape(-1, 1), x2.reshape(-1, 1)])\n",
    "M = np.mean(X, axis=0)\n",
    "S = np.cov(X.T)\n",
    "\n",
    "print('mean', M)\n",
    "print('>')\n",
    "print('cov', S)\n",
    "print('>')\n",
    "print('corr', np.corrcoef(X.T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean [0.14514499 2.         0.21327409]\n",
      ">\n",
      "cov [[ 0.71224389  1.94567309 -0.42263635]\n",
      " [ 1.94567309  0.01        2.92099886]\n",
      " [-0.42263635  2.92099886  0.35076628]]\n",
      ">\n",
      "corr [[1.         0.01068818 0.53011777]\n",
      " [0.01068818 1.         0.00570366]\n",
      " [0.53011777 0.00570366 1.        ]]\n"
     ]
    }
   ],
   "source": [
    "M_u, S_u = update_mean_cov(np.array([2.0]), [1], M, S)\n",
    "\n",
    "print('mean', M_u)\n",
    "print('>')\n",
    "print('cov', S_u)\n",
    "print('>')\n",
    "print('corr', np.corrcoef(np.random.multivariate_normal(M_u, S_u, N*10).T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Links\n",
    "\n",
    "* https://stackoverflow.com/questions/22927181/selecting-specific-rows-and-columns-from-numpy-array\n",
    "* https://docs.scipy.org/doc/numpy/reference/generated/numpy.delete.html\n",
    "* https://stackoverflow.com/questions/38713746/python-numpy-conditional-simulation-from-a-multivatiate-distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00311610448033705\n",
      "1.0002424415942128\n",
      "4.906513396201674e-05\n"
     ]
    }
   ],
   "source": [
    "for m in M_u:\n",
    "    print(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}