CCA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CCA\n",
    "\n",
    "canonical correlation analysis\n",
    "\n",
    "Let $X^{(1)}$ be a p by 1 vector and $X^{(2)}$ be a q by 1 vector, $p<=q$\n",
    "\n",
    "Let $E(X^{1}) = \\mu^1_{(p*1)}; Cov(X^1) = \\Sigma_{11(p*p)}$\n",
    "$E(X^{2}) = \\mu^2_{(q*1)}; Cov(X^2) = \\Sigma_{22(q*q)}$\n",
    "\n",
    "and\n",
    "\n",
    "$Cov(X^1,X^2) = \\Sigma_{12} = \\Sigma_{21}^T$,\n",
    "$E[[X^1 - \\mu^1][X^2 - \\mu^2]^T]$ is a p by q matrix.\n",
    "\n",
    "So:\n",
    "\n",
    "<img src=\"ex.png\" alt=\"Drawing\" style=\"width: 300px;\"/>\n",
    "<img src=\"cov.png\" alt=\"Drawing\" style=\"width: 300px;\"/>\n",
    "\n",
    "Set linear combinations:\n",
    "$U=a^T X^1, V = b^TX^2$\n",
    "\n",
    "We can obtain:\n",
    "\\begin{align*}\n",
    "Var(U) &= a^T Cov(X^1)a = a^T \\Sigma_{11}a\\\\\n",
    "Var(V) &= b^T Cov(X^2)b = b^T \\Sigma_{22}b\\\\\n",
    "Cov(U,V) &= a^T Cov(X^1,X^2)b = a^T \\Sigma_{12}b\n",
    "\\end{align*}\n",
    "\n",
    "We seek coefficient vectors a and b such that\n",
    "\n",
    "$Corr(U,V) = \\frac{a^T\\Sigma_{12}b}{\\sqrt{a^T\\Sigma_{11}a}\\sqrt{b^T\\Sigma_{22}b}}$\n",
    "is as large as possible.\n",
    "\n",
    "Variance of U and V = 1.\n",
    "\n",
    "\n",
    "Define:\n",
    "\n",
    "kth canonical variate pair = the pair of linear combination $U_k$ and $V_k$ having\n",
    "unit variances, which maximize the correlation among all ** choices uncorrelated with\n",
    "the previous k - 1 canonical variate pairs.** Note: uncorrelated does not equal to independent!\n",
    "\n",
    "<img src=\"eig.png\" alt=\"Drawing\" style=\"width: 300px;\"/>\n",
    "\n",
    "Result:\n",
    "\n",
    "$Max_{a,b} Corr(U,V) = \\rho^*_1$\n",
    "\n",
    "attained by the linear combinations\n",
    "\n",
    "$U_1 = e_1^T \\Sigma_{11}^{-1/2}X^1$ and $V_1 = f_1^T \\Sigma_{22}^{-1/2} X^2$\n",
    "\n",
    "The canonical variates have the properties\n",
    "$Var(U_k) = Var(V_k) = 1$\n",
    "\n",
    "$Cov(U_k,U_l) = Cov(V_k,V_l) = Cov(U_k,V_l) = 0$ when $k != l$\n",
    "\n",
    "### Canonical variates from standardized variables\n",
    "\n",
    "$Z^1 = [Z_1^1,Z_2^1,...,Z_p^1]^T$ and\n",
    "\n",
    "$Z^2 = [Z_1^2, Z_2^2,...,Z_q^2]^T$\n",
    "\n",
    "Here, $Cov(Z^1) = \\rho_{11}, Cov(Z^2) = \\rho_{22}$, and\n",
    "\n",
    "$Cov(Z^1,Z^2) = \\rho_{12}=\\rho_{21}^T$.\n",
    "\n",
    "Then use\n",
    "$\\rho_{11}^{-1/2}\\rho_{12}\\rho_{22}^{-1}\\rho_{21}\\rho_{11}^{-1/2}$\n",
    "find eigenvalues and eigenvectors. Eigenvalues are the same as before but the eigenvectors are different.\n",
    "\n",
    "<img src=\"z.png\" alt=\"Drawing\" style=\"width: 300px;\"/>\n",
    "\n",
    "$a^T_k(X^1-\\mu^1) = a_k^T V_{11}^{1/2}Z^1$\n",
    "where $V_{11}$ is the diagonal matrix with ith element $Var(X_i^1)$.\n",
    "\n",
    "### Interpretation:\n",
    "\n",
    "Look at the linear combinatin: $U_1 = a_{11}X_1^1 + a_{12}X_2^1+..$\n",
    "1. $Corr(U_1, X_1^1)$... Here the correlation != the coefficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import dot\n",
    "import scipy.linalg as la\n",
    "from sklearn.cross_decomposition import CCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def var_single(X):\n",
    "    \"\"\"\n",
    "    Find variance matrix\n",
    "    \"\"\"\n",
    "    v = np.cov(X)\n",
    "    r = np.matrix('v,1;1,v')\n",
    "    return r\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "http://users.stat.umn.edu/~helwig/notes/cancor-Notes.pdf\n",
    "reference of R code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cca(X,Y):\n",
    "    \"\"\"\n",
    "    Canonical Correlatio Analysis\n",
    "    \n",
    "    Input:\n",
    "    X: observation matrix X, every column is one data point\n",
    "    Y: observation matrix Y, every column is one data point\n",
    "    \n",
    "    Output:\n",
    "    basis in X space, basis in Y space, correlation\n",
    "    \"\"\"\n",
    "    # find variance and covariance matrix\n",
    "    if len(X) == 1:\n",
    "        cov_xx = var_single(X)\n",
    "    else:\n",
    "        cov_xx = np.cov(X)\n",
    "    if len(Y) == 1:\n",
    "        cov_yy = var_single(Y)\n",
    "    else:\n",
    "        cov_yy = np.cov(Y)\n",
    "    n = len(X)\n",
    "\n",
    "    cov_xy = np.cov(X, Y)[:n,n:]    \n",
    "    cov_yx = np.transpose(cov_xy)\n",
    "    # eigen\n",
    "    cov_xx_evalue,cov_xx_evector = la.eig(cov_xx)\n",
    "    cov_xx_isqrt = dot(dot(cov_xx_evector,np.diag(1/np.sqrt(cov_xx_evalue))),np.transpose(cov_xx_evector))\n",
    "    \n",
    "    cov_yy_evalue, cov_yy_evector = la.eig(cov_yy)\n",
    "    cov_yy_isqrt = dot(dot(cov_yy_evector,np.diag(1/np.sqrt(cov_yy_evalue))), np.transpose(cov_yy_evector))\n",
    "    a = la.inv(cov_yy)\n",
    "    # Xmat and Ymat\n",
    "    Xmat = dot(dot(dot(dot(cov_xx_isqrt,cov_xy),la.inv(cov_yy)),cov_yx),cov_xx_isqrt)\n",
    "    ymat = dot(dot(dot(dot(cov_yy_isqrt,cov_yx),la.inv(cov_xx)),cov_xy),cov_yy_isqrt)\n",
    "    \n",
    "    r1=la.eig(Xmat)\n",
    "    r2=la.eig(Ymat)\n",
    "    \n",
    "    return r1,r2\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a=[-11.25,7.43, 15.48, 2.27, -48.90, -15.13, 49.28, 4.7, 61.32, -268.95, 8488]\n",
    "b=[-10.87, 7.45, 14.97, 1.97, -47.71, -14.46, 44.36, 5.1, 61.76, -273.02, 8399]\n",
    "c=[-11.18, 7.44, 14.20, 1.97, -48.29, -14.81, 43.66, 5.2, 64.16, -263.20, 8328]\n",
    "d=[-10.62, 7.38, 15.02, 2.03, -49.06, -14.72, 44.80, 4.9, 64.04, -285.11, 8306]\n",
    "e=[-11.02, 7.43, 12.92, 1.97, -47.44, -14.40, 41.20, 5.2, 57.46, -256.64, 8286]\n",
    "f=[-10.83, 7.72, 13.58, 2.12, -48.34, -14.18, 43.06, 4.9, 52.18, -274.07, 8272]\n",
    "g=[-11.18, 7.05, 14.12, 2.06, -49.34, -14.39, 41.68, 5.7, 61.60, -291.20, 8216]\n",
    "h=[-11.05, 6.95, 15.34, 2.00, -48.21, -14.36, 41.32, 4.8, 63.00, -265.86, 8189]\n",
    "i=[-11.15, 7.12, 14.52, 2.03, -49.15, -14.66, 42.36,4.9, 66.46, -269.62, 8180]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = np.vstack((a,b,c,d,e,f,g,h,i))\n",
    "run100=data[:,0]\n",
    "long_jump = data[:,1]\n",
    "shot = data[:,2]\n",
    "high_jump = data[:,3]\n",
    "run400 = data[:,4]\n",
    "hurdle = data[:,5]\n",
    "discus = data[:,6]\n",
    "pole_vault = data[:,7]\n",
    "javelin = data[:,8]\n",
    "run1500 = data[:,9]\n",
    "score = data[:,10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "X: shot, discus, javelin, pole_vault\n",
    "\n",
    "Y: run100,run400,run1500,hurdle,long_jump,high_jump"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = np.vstack((shot, discus, javelin,pole_vault))\n",
    "Y = np.vstack((run100, run400,run1500,hurdle,long_jump,high_jump))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "r1,r2=cca(X,Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.60595938+0.j,  0.94686012+0.j,  1.00000000+0.j,  1.00000000+0.j])"
      ]
     },
     "execution_count": 206,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r1[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.77843392+0.j,  0.97306738+0.j,  1.00000000+0.j,  1.00000000+0.j])"
      ]
     },
     "execution_count": 207,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sqrt(r1[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  7.78433925e-01 +0.00000000e+00j,\n",
       "         9.73067377e-01 +0.00000000e+00j,\n",
       "         0.00000000e+00 +7.82039080e-09j,\n",
       "         1.42187762e-08 +0.00000000e+00j,\n",
       "         1.00000000e+00 +0.00000000e+00j,   1.00000000e+00 +0.00000000e+00j])"
      ]
     },
     "execution_count": 208,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sqrt(r2[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}