gaussian discriminant analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian Discriminant Analysis\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "GDA is an interesting classification method. Unlike those optimization-based models, it is a generative learning algorithm based on probabilities. It requires a stronger assumption than Logistic regression. To our surprise, Logistic does not require any distribution assumption, not even exponential family distribution. That's why we should always try Logistic first. \n",
    "\n",
    "The foundation of GDA is Bayes' theorem. Let's breakdown what each element means.\n",
    "\n",
    "$$P (y \\mid x)=\\frac {P (x \\mid y) \\cdot P (y)} {P (x)}$$\n",
    "\n",
    "In a binary GDA, the assumption is that $x$ is always independent of $y$ which yields $P (x)=1$. $P (y)$ denotes the probability density function of Bernoulli distribution which is $P (y)=\\phi^{y}+(1-\\phi)^{1-y}$ where $\\phi$ specifies the portion of each label. $P (x \\mid y)$ denotes the probability density function of multivariate Gaussian distribution (because we assume it follows Gaussian distribution to use GDA). To make prediction, we merely compute the probability for each label and pick the label associated with the higher probability. It's literally a walk in the park.\n",
    "\n",
    "Details of classification and discriminative algorithm such as Logistic\n",
    "\n",
    "https://github.com/je-suis-tm/machine-learning/blob/master/newton%20method%20for%20logistic%20regression.ipynb\n",
    "\n",
    "Reference to some of the mathematics\n",
    "\n",
    "https://en.wikipedia.org/wiki/Bayes%27_theorem\n",
    "https://en.wikipedia.org/wiki/Multivariate_normal_distribution\n",
    "https://en.wikipedia.org/wiki/Covariance_matrix\n",
    "\n",
    "Reference to the lecture material\n",
    "\n",
    "https://see.stanford.edu/materials/aimlcs229/cs229-notes2.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as lda\n",
    "import os\n",
    "os.chdir('d:/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plz note all calculations involved are linear algebra\n",
    "#the sum of x**2 in matrix form is x'x where prime stands for transposed\n",
    "def gaussian_discriminant_analysis(y,x):\n",
    "    \n",
    "    #split test and train as usual\n",
    "    x_train,x_test,y_train,y_test=train_test_split(x,y,test_size=0.3)\n",
    "    \n",
    "    #denote phi as the probability of y==1\n",
    "    #we just calculate the number of samples where y==1\n",
    "    #divided by the total number samples\n",
    "    phi=len(x_train[y_train==1])/len(x_train)\n",
    "    \n",
    "    #calculate the 'mean vector' for two scenarios\n",
    "    #y==0 and y==1\n",
    "    #note that we use pandas dataframe\n",
    "    #so the results are not a vector\n",
    "    #and we dont need to transpose it\n",
    "    #cuz later when we need to input a transposed vector\n",
    "    #we can directly input mean0/mean1\n",
    "    mean0=np.mean(x_train[y_train==0],0)\n",
    "    mean1=np.mean(x_train[y_train==1],0)\n",
    "    \n",
    "    #calculate the difference between x and mean for two scenarios\n",
    "    #concatenate x for two scenarios together\n",
    "    dif=pd.concat([x_train[y_train==0]-mean0,x_train[y_train==1]-mean1])\n",
    "    \n",
    "    #calculate the covariance matrix\n",
    "    #we use a list to append all covariance/variance\n",
    "    #later we reshape it into a 4 by 4 matrix\n",
    "    #our x is four-dimensional so we should have 4 by 4\n",
    "    #this can be done via np.cov\n",
    "    arr=list(dif.columns)\n",
    "    cov=[]    \n",
    "    for i in range(len(dif.columns)):\n",
    "        for j in range(len(dif.columns)):\n",
    "            cov.append( \n",
    "                       ( \n",
    "                        np.mat(dif[arr[i]])*np.mat(\n",
    "                dif[arr[j]]).T).item()/len(arr) \n",
    "                      )\n",
    "    cov=np.mat(cov).reshape(4,4)    \n",
    "    \n",
    "    #now we have mean for two scenarios,covariance matrix and phi\n",
    "    #we use bayesian conditional probability formula\n",
    "    #to calculate the probability of y==0 and y==1 for each x\n",
    "    #and we use the larger probability of the two to forecast y\n",
    "    p0=[]\n",
    "    p1=[]   \n",
    "        \n",
    "    #we calculate the probability for each element in x matrix\n",
    "    #if we use the whole matrix to do algebra\n",
    "    #lets forget anything before np.e as we can see determinant there\n",
    "    #which indicates everything before np.e becomes scalar\n",
    "    #for (x-mean).T*cov.I*(x-mean)\n",
    "    #(x-mean).T is n by 4 matrix\n",
    "    #cov.I is 4 by 4 matrix\n",
    "    #we get n by 4 matrix\n",
    "    #then multiplied by (x-mean) which is 4 by n matrix\n",
    "    #we would end up with n by n matrix in the end\n",
    "    #oh,god,thats not what we want\n",
    "    #when (x-mean).T is 1 by 4 matrix\n",
    "    #after multiplication by cov.I which is 4 by 4 matrix\n",
    "    #now we get 1 by 4 matrix\n",
    "    #eventually times (x-mean).T which is 4 by 1 matrix\n",
    "    #we get scalar,a 1 by 1 matrix!!!\n",
    "    for k in range(len(x_test)):        \n",
    "        probability0=phi/( \n",
    "                          np.linalg.det(2*np.pi*cov)**(0.5) \n",
    "                         ) \ \n",
    "        *np.exp( \n",
    "                -0.5*(np.mat(x_test.iloc[k]-mean0))*cov.I*(\n",
    "                    np.mat(x_test.iloc[k]-mean0)).T \n",
    "               )       \n",
    "        probability1=phi/( \n",
    "                          np.linalg.det(2*np.pi*cov)**(0.5) \n",
    "                         ) \ \n",
    "        *np.exp( \n",
    "                -0.5*(np.mat(x_test.iloc[k]-mean1))*cov.I*(\n",
    "                    np.mat(x_test.iloc[k]-mean1)).T \n",
    "               )       \n",
    "        p0.append(probability0.item())\n",
    "        p1.append(probability1.item())\n",
    "    \n",
    "    #here we use numpy sign\n",
    "    #numpy sign treats positive number as 1\n",
    "    #it treats zero as 0\n",
    "    #however,it treats negative number as -1\n",
    "    #we use a map function to convert -1 to 0\n",
    "    forecast=np.sign(np.subtract(p1,p0))\n",
    "    forecast=list(map(lambda x: 0 if x<0 else int(x),forecast))\n",
    "    \n",
    "    #accuracy\n",
    "    print('test accuracy: {}%'.format(len(\n",
    "        y_test[forecast==y_test])/len(y_test)*100))\n",
    "    \n",
    "    #just too lazy to write codes for plotting\n",
    "\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#linear discriminant analysis from sklearn\n",
    "def LDA(y,x):\n",
    "    \n",
    "    x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3)    \n",
    "    m=lda().fit(x_train,y_train)    \n",
    "    forecast=m.predict(x_test)    \n",
    "    print('test accuracy: {}%'.format(len(\n",
    "        y_test[forecast==y_test])/len(y_test)*100))\n",
    "    \n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#iris data manipulation as usual\n",
    "#simplify data into a binary classification problem\n",
    "df=pd.read_csv('iris.csv')\n",
    "df=df[df['type']!='Iris-versicolor']\n",
    "df['y']=np.select([df['type']=='Iris-setosa',df['type']=='Iris-virginica'],[0,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#prepare x and y\n",
    "x=pd.concat([df[i] for i in df.columns[:4]],axis=1)\n",
    "y=df['y']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test accuracy: 100.0%\n"
     ]
    }
   ],
   "source": [
    "gaussian_discriminant_analysis(y,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test accuracy: 100.0%\n"
     ]
    }
   ],
   "source": [
    "LDA(y,x)"
   ]
  },
  {
   "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.6.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}