mlvsmax.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Created by @Mehmet Ozan UNAL\n",
    "for BYM 509 @Istanbul Technical Univercity\n",
    "# Maximum Likelihood  Estimation vs Maximum a Posteriori Notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Create Simulation Enviroment\n",
    "\n",
    "In this example $Y=AX+N$ simulation enviroment will be created.\n",
    "In this equation:\n",
    "- Y is the observation\n",
    "- A is system geometry parameters\n",
    "- X is varible (the thing which is observed )\n",
    "- N is randomly distrubuted zero mean gaussian noise\n",
    "\n",
    "Size of this vector and matrices are:\n",
    "$$ Y_{mx1}=A_{mxn}X_{nx1} + N_{mx1} $$\n",
    "\n",
    "\n",
    "now contiune with code.\n",
    "\n",
    "#### Import libraries\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import misc\n",
    "\n",
    "#plt.rcParams['figure.figsize'] = [10, 5]\n",
    "%matplotlib notebook "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert to vector\n",
    "def toVector(a):\n",
    "    return a.reshape(-1,1) \n",
    "    \n",
    "def toImage(a):\n",
    "    return a.reshape(K,L)\n",
    "\n",
    "def convLog(arr):\n",
    "    sample10 = np.mean( arr[0:10] )\n",
    "    sample100 = np.mean( arr[0:100] )\n",
    "    sample1_000 = np.mean( arr[0:1000] )\n",
    "    sample10_000 = np.mean( arr[0:10000] )\n",
    "    return [sample10, sample100, sample1_000, sample10_000]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementation of ML and MAP at 100x100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = 100\n",
    "L = 100\n",
    "m = int(1.2*K*L)\n",
    "n = K*L\n",
    "\n",
    "# Create X read image\n",
    "X_image = misc.imread('test/test1.jpg', mode=\"L\") # grayscale mode\n",
    "X_image = misc.imresize(X_image, (K,L))/ 255 # normalize image\n",
    "X = toVector(X_image)\n",
    "A = np.random.rand(m,n)\n",
    "\n",
    "\n",
    "def calculateMaxlikelihood(A, X, NoisePower):\n",
    "    # every time new noise will be generated randomly\n",
    "    N =  (np.random.rand(m,1)-0.5) * NoisePower\n",
    "    Y = A@X + N\n",
    "    # return Xml\n",
    "    return  np.linalg.inv(A.transpose() @ A ) @ A.transpose() @ Y\n",
    "\n",
    "def calculateMaximumAPosteriori(A, X, NoisePower, B):\n",
    "    # every time new noise will be generated randomly\n",
    "    N =  (np.random.rand(m,1)-0.5) * NoisePower\n",
    "    Y = A@X + N\n",
    "    I1 = np.ones((n,n))\n",
    "    I2 = np.ones((n,1))\n",
    "    # return Xmap\n",
    "    return  np.linalg.inv(A.transpose() @ A + B * I1 ) @ ( A.transpose() @ Y + B*X.mean()*I2)\n",
    "\n",
    "print(\"Started.\")\n",
    "Xml =  calculateMaxlikelihood(A,X,1)\n",
    "print(\"ML done.\")\n",
    "Xmap01 =  calculateMaximumAPosteriori(A,X,1,0.1)\n",
    "print(\"MAP01 done\")\n",
    "Xmap03 =  calculateMaximumAPosteriori(A,X,1,0.3)\n",
    "print(\"MAP03 done\")\n",
    "Xmap05 =  calculateMaximumAPosteriori(A,X,1,0.5)\n",
    "print(\"MAP05 done\")\n",
    "Xmap10 =  calculateMaximumAPosteriori(A,X,1,1.0)\n",
    "print(\"MAP10 done\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "\n",
    "ax1 = fig.add_subplot(321)\n",
    "ax1.title.set_text('Original Image')\n",
    "ax1.imshow(toImage(X))\n",
    "\n",
    "ax2 = fig.add_subplot(322)\n",
    "ax2.title.set_text('ML Estimated Image')\n",
    "ax2.imshow(toImage(Xml))\n",
    "\n",
    "ax3 = fig.add_subplot(323)\n",
    "ax3.title.set_text('Map Estimated Image B = 0.1')\n",
    "ax3.imshow(toImage(Xmap01))\n",
    "\n",
    "\n",
    "ax4 = fig.add_subplot(324)\n",
    "ax4.title.set_text('Map Estimated Image B = 0.3')\n",
    "ax4.imshow(toImage(Xmap03))\n",
    "\n",
    "ax5 = fig.add_subplot(325)\n",
    "ax5.title.set_text('Map Estimated Image B = 0.5')\n",
    "ax5.imshow(toImage(Xmap05))\n",
    "\n",
    "ax6 = fig.add_subplot(326)\n",
    "ax6.title.set_text('Map Estimated Image B = 1.0')\n",
    "ax6.imshow(toImage(Xmap10))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Decrease sizes for fast calculations to 10x10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = 10\n",
    "L = 10\n",
    "m = int(1.2*K*L)\n",
    "n = K*L\n",
    "\n",
    "# Create X read image\n",
    "X_image = misc.imread('test/test1.jpg', mode=\"L\") # grayscale mode\n",
    "X_image = misc.imresize(X_image, (K,L))/ 255 # normalize image\n",
    "X = toVector(X_image)\n",
    "A = np.random.rand(m,n)\n",
    "\n",
    "\n",
    "def calculateMaxlikelihood(A, X, NoisePower):\n",
    "    # every time new noise will be generated randomly\n",
    "    N =  (np.random.rand(m,1)-0.5) * NoisePower\n",
    "    Y = A@X + N\n",
    "    # return Xml\n",
    "    return  np.linalg.inv(A.transpose() @ A ) @ A.transpose() @ Y\n",
    "\n",
    "def calculateMaximumAPosteriori(A, X, NoisePower, B):\n",
    "    # every time new noise will be generated randomly\n",
    "    N =  (np.random.rand(m,1)-0.5) * NoisePower\n",
    "    Y = A@X + N\n",
    "    I1 = np.ones((n,n))\n",
    "    I2 = np.ones((n,1))\n",
    "    # return Xmap\n",
    "    return  np.linalg.inv(A.transpose() @ A + B * I1 ) @ ( A.transpose() @ Y + B*X.mean()*I2)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Bias Calculation\n",
    "Bias is calculated as:\n",
    "\n",
    "Theoritical bias of maximum likelihood is zero. It can be calculated as:\n",
    "\n",
    "Therefore the real bias should be decrease with increase of tries\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculcateBias(x, xml):\n",
    "    B = x-xml\n",
    "    return B.mean()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Varience Calculation\n",
    "Varience of estimations calculated as:\n",
    "$$\n",
    "varience = \\frac1{n} \\sum\\limits_{i=1}^n(x_i - \\overline{x})^2 \n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculateVarience(x, xml):\n",
    "    meanxml = xml.mean()\n",
    "    return np.square(xml-meanxml).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculateMSE(x, xml):\n",
    "    return np.square(xml-x).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create Simulation Class\n",
    "- To easily run simulations this class is created. This is using as container to save simulation results to and object and run this object with function pointer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "class sim:\n",
    "    def __init__(self, name, func, noisePower=1, B=None):\n",
    "        self.biasArr = []\n",
    "        self.varArr = []\n",
    "        self.mseArr = []\n",
    "        self.name = name\n",
    "        self.calculate = func\n",
    "        self.noisePower = noisePower\n",
    "        self.B = B\n",
    "    \n",
    "    def run(self):\n",
    "        if self.B == None:\n",
    "            X_ = self.calculate(A,X, self.noisePower)\n",
    "        else:\n",
    "            X_ = self.calculate(A,X, self.noisePower, self.B)\n",
    "        self.biasArr.append(calculcateBias(X, X_))\n",
    "        self.varArr.append(calculateVarience(X, X_))\n",
    "        self.mseArr.append(calculateMSE(X,X_))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run Simulation 1000 times for ML and Map with b=0.2 b=0.4 b=0.8 b=1.0 b=5.0\n",
    "Run simulation 1000 times and save all bias, varience and mean square error to arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NOISE_POW = 1\n",
    "sim_maxLikelihood = sim(\"Maximum Likelihood\", calculateMaxlikelihood, noisePower=NOISE_POW)\n",
    "sim_map02 = sim(\"MAP B=0.2\", calculateMaximumAPosteriori, noisePower=NOISE_POW, B=0.2)\n",
    "sim_map04 = sim(\"MAP B=0.4\", calculateMaximumAPosteriori, noisePower=NOISE_POW, B=0.4)\n",
    "sim_map08 = sim(\"MAP B=0.8\", calculateMaximumAPosteriori, noisePower=NOISE_POW, B=0.8)\n",
    "sim_map10 = sim(\"MAP B=1.0\", calculateMaximumAPosteriori, noisePower=NOISE_POW, B=1.0)\n",
    "sim_map50 = sim(\"MAP B=5.0\", calculateMaximumAPosteriori, noisePower=NOISE_POW, B=5.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "for i in range(1000):\n",
    "    sim_maxLikelihood.run()\n",
    "    sim_map02.run()\n",
    "    sim_map04.run()\n",
    "    sim_map08.run()\n",
    "    sim_map10.run()\n",
    "    sim_map50.run()\n",
    "    if i % 100 == 0:\n",
    "        print(i,)\n",
    "print(\"Done.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Result for different B values\n",
    "#### Results Bias"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " \n",
    "plt.figure()\n",
    "biasAll = [ np.mean(sim_map02.biasArr), np.mean(sim_map04.biasArr),\n",
    "           np.mean(sim_map08.biasArr),np.mean(sim_map10.biasArr), np.mean(sim_map50.biasArr)  ]\n",
    "plt.plot([\"0.2\",\"0.4\",\"0.8\",\"1.0\",\"5.0\"], biasAll)\n",
    "plt.title('Bias & B')\n",
    "plt.xlabel('B values')\n",
    "plt.ylabel('Bias')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Result Varience"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "varAll = [ np.mean(sim_map02.varArr), np.mean(sim_map04.varArr),\n",
    "           np.mean(sim_map08.varArr),np.mean(sim_map10.varArr), np.mean(sim_map50.varArr)  ]\n",
    "plt.plot([\"0.2\",\"0.4\",\"0.8\",\"1.0\",\"5.0\"], varAll)\n",
    "plt.title('Varience & B')\n",
    "plt.xlabel('B values')\n",
    "plt.ylabel('Varience')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Result MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "mseAll = [ np.mean(sim_map02.mseArr), np.mean(sim_map04.mseArr),\n",
    "           np.mean(sim_map08.mseArr),np.mean(sim_map10.mseArr), np.mean(sim_map50.mseArr)  ]\n",
    "plt.plot([\"0.2\",\"0.4\",\"0.8\",\"1.0\",\"5.0\"], mseAll)\n",
    "plt.title('MSE & B')\n",
    "plt.xlabel('B values')\n",
    "plt.ylabel('MSE')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run Simulation 10,000 times for different power level noise\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_map_noise01 = sim(\"MAP N=0.1\", calculateMaximumAPosteriori, noisePower=0.1, B=0.3)\n",
    "sim_map_noise02 = sim(\"MAP N=0.2\", calculateMaximumAPosteriori, noisePower=0.2, B=0.3)\n",
    "sim_map_noise04 = sim(\"MAP N=0.4\", calculateMaximumAPosteriori, noisePower=0.4, B=0.3)\n",
    "sim_map_noise08 = sim(\"MAP N=0.8\", calculateMaximumAPosteriori, noisePower=0.8, B=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(10000):\n",
    "    sim_map_noise01.run()\n",
    "    sim_map_noise02.run()\n",
    "    sim_map_noise04.run()\n",
    "    sim_map_noise08.run()\n",
    "    if i % 1000 == 0:\n",
    "        print(i,)\n",
    "print(\"Done.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot([\"10\",\"100\",\"1000\",\"10000\"],convLog(sim_map_noise01.biasArr), label=\"map noise=0.1\")\n",
    "plt.plot([\"10\",\"100\",\"1000\",\"10000\"],convLog(sim_map_noise02.biasArr), label=\"map noise=0.2\")\n",
    "plt.plot([\"10\",\"100\",\"1000\",\"10000\"],convLog(sim_map_noise04.biasArr), label=\"map noise=0.4\")\n",
    "plt.plot([\"10\",\"100\",\"1000\",\"10000\"],convLog(sim_map_noise08.biasArr), label=\"map noise=0.8\")\n",
    "plt.legend(loc='lower right')\n",
    "plt.title('Bias & Number of Sim')\n",
    "plt.xlabel('Number of Sim')\n",
    "plt.ylabel('Bias')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot([\"10\",\"100\",\"1000\",\"10000\"],convLog(sim_map_noise01.varArr),label=\"map noise=0.1\")\n",
    "plt.plot([\"10\",\"100\",\"1000\",\"10000\"],convLog(sim_map_noise02.varArr),label=\"map noise=0.2\")\n",
    "plt.plot([\"10\",\"100\",\"1000\",\"10000\"],convLog(sim_map_noise04.varArr),label=\"map noise=0.4\")\n",
    "plt.plot([\"10\",\"100\",\"1000\",\"10000\"],convLog(sim_map_noise08.varArr),label=\"map noise=0.8\")\n",
    "plt.legend(loc='lower right')\n",
    "plt.title('Varience & Number of Sim')\n",
    "plt.xlabel('Number of Sim')\n",
    "plt.ylabel('Varience')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot([\"10\",\"100\",\"1000\",\"10000\"],convLog(sim_map_noise01.mseArr),label=\"map noise=0.1\")\n",
    "plt.plot([\"10\",\"100\",\"1000\",\"10000\"],convLog(sim_map_noise02.mseArr),label=\"map noise=0.2\")\n",
    "plt.plot([\"10\",\"100\",\"1000\",\"10000\"],convLog(sim_map_noise04.mseArr),label=\"map noise=0.4\")\n",
    "plt.plot([\"10\",\"100\",\"1000\",\"10000\"],convLog(sim_map_noise08.mseArr),label=\"map noise=0.8\")\n",
    "plt.legend(loc='lower right')\n",
    "plt.title('MSE & Number of Sim')\n",
    "plt.xlabel('Number of Sim')\n",
    "plt.ylabel('MSE')\n",
    "plt.show()"
   ]
  },
  {
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}