speed-orthogonal.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from scipy.stats import linregress\n",
    "from scipy.odr import Model, Data, ODR\n",
    "from sklearn.decomposition import TruncatedSVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu, sigma = 0, 1 # mean and standard deviation\n",
    "n = 100000\n",
    "s1 = np.random.normal(mu, sigma, n)\n",
    "s2 = np.random.normal(mu, sigma, n)\n",
    "s3 = np.random.normal(mu, sigma, n)\n",
    "\n",
    "a = 0.8\n",
    "x1 = np.cumsum(s1)+s2\n",
    "x2 = a*np.cumsum(s1)+s3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 작은 예제 1\n",
    "x2 = np.array([0, 0, 1, 1])\n",
    "x1 = np.array([0, 1, 1, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# # 작은 예제 2\n",
    "# x2 = np.array([0, 1])\n",
    "# x1 = np.array([0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(x, y) = (x2, x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x11223b978>"
      ]
     },
     "execution_count": 187,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEPCAYAAABsj5JaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE2tJREFUeJzt3X+s5Hdd7/Hnq2w3rFCqbdcKKz2rFFxtWJoGlr2X3uus\nWDkkJkv4Q0tjkQa0BkETc8O2RtPjzTXQ/5TbVLLctRETtlVq0kWtFHBHU6F4ENoF3e0u4Fn6Q+qB\nC72Bu9671Ld/zHQ9np6fnz1nvmd2n49kkvnO9zPf85rpfPva7/c73++kqpAkabUu6DqAJGk8WSCS\npCYWiCSpiQUiSWpigUiSmlggkqQmnRdIkgNJnkpyZJH5NyR5ZHh7MMkrR51RkvRcnRcIcBfwhiXm\nfwX4r1X1KuB/AB8cSSpJ0pI2dR2gqh5MMrHE/IfmTD4EbFv/VJKk5WyELZDVeAdwf9chJEkbYAtk\npZLsAW4Cru06iyRpTAokyU5gPzBZVd9cYpwX9pKkVaqqtDxvo+zCyvD23BnJFcC9wI1V9eXlFlRV\nY3m77bbbOs9g/u5zmH88b+Oc/2x0vgWS5MNAD7g0yVeB24DNQFXVfuA3gUuAO5MEOF1Vu7rKK0ka\n6LxAquqGZeb/AvALI4ojSVqhjbIL67zX6/W6jnBWzN8t83dr3PO3ytnuA9tIktS59Hokab0locb8\nILokacxYIJKkJhaIJKmJBSJJamKBSJKaWCCSpCYWiCSpiQUiSWpigUiSmlggkqQmFogkqYkFIklq\nYoFIkppYIJKkJhaIJKmJBSJJamKBSJKaWCCSpCYWiCSpiQUiSWpigUiSmnReIEkOJHkqyZElxrw/\nyYkkDye5epT5JEkL67xAgLuANyw2M8kbgZdV1cuBm4EPjCrYqMzOzjI9Pc3s7GzXUSSt0vm8/nZe\nIFX1IPDNJYbsBT40HPsZ4OIkl48i2ygcPHgPExM7uO66X2JiYgcHD97TdSRJK3S+r7+pqq4zkGQC\n+GhV7Vxg3keB91bVp4bTnwDeU1WfW2BsbYTXs1Kzs7NMTOzg1KnDwE7gCFu27OHkyWNs3bq163iS\nlnCurL9JqKq0PHfTWofp2tTU1Jn7vV6PXq/XWZblzMzMsHnzdk6derY3d3LhhRPMzMyM1QdQOh+N\n6/rb7/fp9/trsqxx2AL5AHC4qu4ZTh8DfryqnlpgrFsgkkbiXFl/z2YLpPNjIEMZ3hZyCHgrQJLd\nwLcWKo9xtHXrVg4cuJMtW/bwohddw5Ytezhw4M6x+vBJ5yvX3w2wBZLkw0APuBR4CrgN2AxUVe0f\njrkDmAS+A9y00PGP4bix2gJ51uzsLDMzM2zfvv28+vBJ54JxX3/PZguk8wJZS+NaIJLUlXNhF5Yk\nacxYIJKkJhaIJKmJBSJJamKBSJKaWCCSpCYWiCSpiQUiSWpigUiSmlggkqQmFogkqYkFIklqYoFI\nkppYIJKkJhaIJKmJBSJJamKBSJKaWCCSpCYWiCSpiQUiSWpigUiSmlggkqQmFogkqUnnBZJkMsmx\nJMeT7Ftg/qVJ7k/ycJIvJHlbBzElSfOkqrr748kFwHHg9cCTwDRwfVUdmzPmNuD5VXVrksuAR4HL\nq+q7Cyyvunw9kjRuklBVaXlu11sgu4ATVXWyqk4DdwN75435GnDR8P5FwDcWKg9J0mht6vjvbwMe\nmzP9OINSmeuDwCeTPAm8EPjZEWWTJC2h6wJZiVuBR6pqT5KXAR9PsrOqvr3Q4KmpqTP3e70evV5v\nJCElaRz0+336/f6aLKvrYyC7gamqmhxO3wJUVd0+Z8yfA79dVX8znP4ksK+qPrvA8jwGIkmrMM7H\nQKaBK5NMJNkMXA8cmjfmKPCTAEkuB14BfGWkKSVJz9HpLqyqeibJu4AHGJTZgao6muTmwezaD7wX\nuCvJI0CA91TV/+4utSQJOt6FtdbchSVJqzPOu7AkSWPKApEkNbFAJElNLBBJUhMLRJLUxAKRJDWx\nQCRJTSwQSVITC0SS1MQCkSQ1sUAkSU0sEElSEwtEktTEApEkNbFAJElNLBBJUhMLRJLUxAKRJDWx\nQCRJTSwQSVITC0SS1MQCkSQ1sUAkSU06L5Akk0mOJTmeZN8iY3pJPp/ki0kOjzqjJOm5UlXd/fHk\nAuA48HrgSWAauL6qjs0ZczHwKeCnquqJJJdV1dcXWV51+XokadwkoarS8tyut0B2ASeq6mRVnQbu\nBvbOG3MDcG9VPQGwWHlIkkar6wLZBjw2Z/rx4WNzvQK4JMnhJNNJbhxZOknSojZ1HWAFNgHXAD8B\nvAD4dJJPV9WXFho8NTV15n6v16PX640goiSNh36/T7/fX5NldX0MZDcwVVWTw+lbgKqq2+eM2Qc8\nv6p+azj9v4D7q+reBZbnMRBJWoVxPgYyDVyZZCLJZuB64NC8MfcB1yZ5XpLvAV4LHB1xTknSPJ3u\nwqqqZ5K8C3iAQZkdqKqjSW4ezK79VXUsyceAI8AzwP6q+ocOY0uS6HgX1lpzF5Ykrc4478KSJI0p\nC0SS1MQCkSQ1sUAkSU0sEElSEwtEktTEApEkNbFAJElNmgokyXVrHUSSNF6azkRP8tWqumId8pwV\nz0SXpNU5mzPRF70WVpL5FzU8Mwu4tOWPSZLOHUtdTPG/AD8HfHve42HwS4KSpPPYUgXyEPB/q+qv\n5s9I8uj6RZIkjYNlj4Ek+bH5l09P0quq/noGa+ExEElanfW+Gu8fJdmXgS1J/ifw3pY/Jkk6d6yk\nQF4LvBT4FINfEHwSeN16hpIkbXwrKZDTwClgC/B84B+r6l/XNZUkacNbSYFMMyiQ1zD4ZtZbkvzx\nuqaSJG14KzmI/uqq+uy8x26sqj9c12QNPIguSatzNgfR/U10STqP+ZvokqSRs0AkSU0sEElSk84L\nJMlkkmNJjifZt8S41yQ5neTNo8wnSVpYpwWS5ALgDuANwFUMviK8Y5Fx7wM+NtqEkqTFdL0Fsgs4\nUVUnq+o0cDewd4Fx7wY+AvzzKMNJkhbXdYFsAx6bM/348LEzkrwEeFNV/R6DS8lLkjaApS7nvlH8\nDjD32MiSJTI1NXXmfq/Xo9frrUsoSRpH/X6ffr+/Jsvq9ETCJLuBqaqaHE7fAlRV3T5nzFeevQtc\nBnwH+MWqes4vJnoioSStztieiZ7kecCjwOuBfwL+FnhLVR1dZPxdwEer6k8WmW+BSNIqrMtvoo9C\nVT2T5F3AAwyOxxyoqqNJbh7Mrv3znzLykJKkBXktLEk6j3ktLEnSyFkgkqQmFogkqYkFIklqYoFI\nkppYIJKkJhaIJKmJBSJJamKBSJKaWCCSpCYWiCSpiQUiSWpigUiSmlggkqQmFogkqYkFIklqYoFI\nkppYIJKkJhaIJKmJBSJJamKBSJKaWCCSpCYWiCSpSecFkmQyybEkx5PsW2D+DUkeGd4eTPLKLnJK\nkv6jVFV3fzy5ADgOvB54EpgGrq+qY3PG7AaOVtXTSSaBqaravcjyqsvXI0njJglVlZbndr0Fsgs4\nUVUnq+o0cDewd+6Aqnqoqp4eTj4EbBtxRknSAroukG3AY3OmH2fpgngHcP+6JpIkrcimrgOsVJI9\nwE3AtUuNm5qaOnO/1+vR6/XWNZckjZN+v0+/31+TZXV9DGQ3g2Mak8PpW4CqqtvnjdsJ3AtMVtWX\nl1iex0AkaRXG+RjINHBlkokkm4HrgUNzByS5gkF53LhUeUiSRqvTXVhV9UySdwEPMCizA1V1NMnN\ng9m1H/hN4BLgziQBTlfVru5SS5Kg411Ya81dWJK0OuO8C0uSNKYsEElSEwtEktTEApEkNbFAJElN\nLBBJUhMLRJLUxAKRJDWxQCRJTSwQSVITC0SS1MQCkSQ1sUAkSU0sEElSEwtEktTEApEkNbFAJElN\nLBBJUhMLRJLUxAKRJDWxQCRJTSwQSVITC0SS1KTzAkkymeRYkuNJ9i0y5v1JTiR5OMnVo84oSXqu\nTgskyQXAHcAbgKuAtyTZMW/MG4GXVdXLgZuBD4w86DqbnZ1lenqa2dnZrqNIWqXzef3tegtkF3Ci\nqk5W1WngbmDvvDF7gQ8BVNVngIuTXD7amOvn4MF7mJjYwXXX/RITEzs4ePCeriNJWqHzff3tukC2\nAY/NmX58+NhSY55YYMxYmp2d5e1vfyenTh3m6af/jlOnDvP2t7/zvPyXjDRuXH9hU9cB1trU1NSZ\n+71ej16v11mW5czMzLB583ZOndo5fGQnF144wczMDFu3bu00m6Sljev62+/36ff7a7KsVNWaLKjp\njye7gamqmhxO3wJUVd0+Z8wHgMNVdc9w+hjw41X11ALLqy5fz2rNzs4yMbGDU6cOAzuBI2zZsoeT\nJ49t6A+gpHNn/U1CVaXluV3vwpoGrkwykWQzcD1waN6YQ8Bb4UzhfGuh8hhHW7du5cCBO9myZQ8v\netE1bNmyhwMH7hyrD590vnL97XgLBAZf4wV+l0GZHaiq9yW5mcGWyP7hmDuASeA7wE1V9blFljVW\nWyDPmp2dZWZmhu3bt59XHz7pXDDu6+/ZbIF0XiBraVwLRJK6Ms67sCRJY8oCkSQ1sUAkSU0sEElS\nEwtEktTEApEkNbFAJElNLBBJUhMLRJLUxAKRJDWxQCRJTSwQSVITC0SS1MQCkSQ1sUAkSU0sEElS\nEwtEktTEApEkNbFAJElNLBBJUhMLRJLUxAKRJDXprECSfF+SB5I8muRjSS5eYMwPJvnLJH+f5AtJ\nfqWLrJKk5+pyC+QW4BNV9SPAXwK3LjDmu8CvVdVVwH8CfjnJjhFmHJl+v991hLNi/m6Zv1vjnr9V\nlwWyF/iD4f0/AN40f0BVfa2qHh7e/zZwFNg2soQjNO4fQPN3y/zdGvf8rboskO+vqqdgUBTA9y81\nOMl24GrgM+ueTJK0rE3rufAkHwcun/sQUMBvLDC8lljOC4GPAL863BKRJHUsVYv+f3t9/3ByFOhV\n1VNJfgA4XFU/usC4TcCfAvdX1e8us8xuXowkjbGqSsvz1nULZBmHgLcBtwM/D9y3yLjfB/5hufKA\n9jdBkrR6XW6BXAL8EfBS4CTwM1X1rSQvBj5YVT+d5HXAXwNfYLCLq4Bfr6q/6CS0JOmMzgpEkjTe\nxvZM9HE9ETHJZJJjSY4n2bfImPcnOZHk4SRXjzrjUpbLn+SGJI8Mbw8meWUXORezkvd/OO41SU4n\nefMo8y1nhZ+fXpLPJ/liksOjzriYFXx2Lk1y//Bz/4Ukb+sg5qKSHEjyVJIjS4zZkOvuctmb19uq\nGssbg2Mn7xne3we8b4ExPwBcPbz/QuBRYEeHmS8AvgRMABcCD8/PA7wR+LPh/dcCD3X9Xq8y/27g\n4uH9yXHLP2fcJxl8eePNXede5ft/MfD3wLbh9GVd515F9tuA9z6bG/gGsKnr7HPyXcvgVIIji8zf\nyOvuctmb1tux3QJhPE9E3AWcqKqTVXUauJvB65hrL/AhgKr6DHBxksvZGJbNX1UPVdXTw8mH2Fgn\nfq7k/Qd4N4Ovjf/zKMOtwEry3wDcW1VPAFTV10eccTEryf414KLh/YuAb1TVd0eYcUlV9SDwzSWG\nbNh1d7nsrevtOBfIOJ6IuA14bM704zz3P9T8MU8sMKYrK8k/1zuA+9c10eosmz/JS4A3VdXvMThv\naSNZyfv/CuCSJIeTTCe5cWTplraS7B8ErkryJPAI8KsjyrZWNvK6uxorXm+7/BrvsjwRcXwl2QPc\nxGDTeZz8DoNdos/aaCWynE3ANcBPAC8APp3k01X1pW5jrcitwCNVtSfJy4CPJ9npOjs6q11vN3SB\nVNV1i80bHhC6vP79RMQFdzcMT0T8CPCHVbXYuSaj8gRwxZzpHxw+Nn/MS5cZ05WV5CfJTmA/MFlV\nS23yj9pK8r8auDtJGOyHf2OS01V1aEQZl7KS/I8DX6+qfwH+JclfA69icPyhSyvJ/jrgtwGq6stJ\n/hHYAXx2JAnP3kZed5fVst6O8y6sZ09EhDU6EXEEpoErk0wk2Qxcz+B1zHUIeCtAkt3At57dVbcB\nLJs/yRXAvcCNVfXlDjIuZdn8VfXDw9sPMfiHxzs3SHnAyj4/9wHXJnleku9hcDD36IhzLmQl2Y8C\nPwkwPHbwCuArI025vLD4VulGXndhiezN623X3w44i28VXAJ8gsE3qx4Avnf4+IuBPx3efx3wDINv\nfHwe+ByDdu0y9+Qw8wngluFjNwO/OGfMHQz+xfgIcE3X7/Vq8jPYj/2N4Xv9eeBvu8682vd/ztjf\nZwN9C2sVn5//xuCbWEeAd3edeRWfncuAjw4/90eAt3SdeV7+DwNPAv8P+CqDXT1jse4ul711vfVE\nQklSk3HehSVJ6pAFIklqYoFIkppYIJKkJhaIJKmJBSJJamKBSOssyauSfGp4ifKHk/xM15mkteB5\nINI6S/Jy4F9rcHmOFwN/x+BS5v+n42jSWXELRFpDSV49/FGezUlekOSLwIU1vDxEVf0Tg+u2be00\nqLQG3AKR1liS/w5sGd4eq6rb58zbBdxVVVd1lU9aKxaItMaSXMjg4oGngP9cw5VsuPvqMIML1k13\nGFFaE+7CktbeZQx+Qvki4PkASS5i8BO5t1oeOle4BSKtsST3AQeBHwJeAvwa8BfAfVX1/i6zSWtp\nQ/+glDRuhj8h+/+r6u4kFwB/w+C3L64Fvi/JTQx+PfNtVXWkw6jSWXMLRJLUxGMgkqQmFogkqYkF\nIklqYoFIkppYIJKkJhaIJKmJBSJJamKBSJKa/BsGeBu6j/xaIQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11223b828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pd.DataFrame({'x1': x1, 'x2': x2}).plot(kind='scatter', y='x1', x='x2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 기존 방법"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def orthoregress(x, y):\n",
    "    \"\"\"Perform an Orthogonal Distance Regression on the given data,\n",
    "    using the same interface as the standard scipy.stats.linregress function.\n",
    "    Arguments:\n",
    "    x: x data\n",
    "    y: y data\n",
    "    Returns:\n",
    "    [m, c, nan, nan, nan]\n",
    "    Uses standard ordinary least squares to estimate the starting parameters\n",
    "    then uses the scipy.odr interface to the ODRPACK Fortran code to do the\n",
    "    orthogonal distance calculations.\n",
    "    \"\"\"\n",
    "    linreg = linregress(x, y)\n",
    "    mod = Model(f)\n",
    "    dat = Data(x, y)\n",
    "    od = ODR(dat, mod, beta0=linreg[0:2])\n",
    "    out = od.run()\n",
    "\n",
    "    \n",
    "    return list(out.beta)\n",
    "\n",
    "\n",
    "def f(p, x):\n",
    "    \"\"\"Basic linear regression 'model' for use with ODR\"\"\"\n",
    "    return (p[0] * x) + p[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000 loops, best of 3: 529 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit orthoregress(x2, x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.0, 0.5]"
      ]
     },
     "execution_count": 190,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res1 = orthoregress(x2, x1)\n",
    "res1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 새로운 시도"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def ortho_tsvd(x2, x1):\n",
    "    X = np.hstack((x2[:, np.newaxis], x1[:, np.newaxis]))\n",
    "    tsvd = TruncatedSVD(n_components=1)\n",
    "    tsvd.fit(X)\n",
    "    beta = tsvd.components_[0][1] / tsvd.components_[0][0]\n",
    "    center = X.mean(axis=0)\n",
    "    intercept = center[1] - center[0] * beta\n",
    "    return tsvd.components_[0][1] / tsvd.components_[0][0], intercept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000 loops, best of 3: 429 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit ortho_tsvd(x2, x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.0000000000000007, -3.3306690738754696e-16)"
      ]
     },
     "execution_count": 193,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res2 = ortho_tsvd(x2, x1)\n",
    "res2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 새로운 시도 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def ortho_npsvd(x2, x1):\n",
    "    X = np.hstack((x2[:, np.newaxis], x1[:, np.newaxis]))\n",
    "    U, s, V = np.linalg.svd(X, full_matrices=False)\n",
    "    beta = V[0][1] / V[0][0]  # V 의 row vector 가 basis\n",
    "    center = X.mean(axis=0)\n",
    "    intercept = center[1] - center[0] * beta\n",
    "    return beta, intercept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 7.27 times longer than the fastest. This could mean that an intermediate result is being cached \n",
      "10000 loops, best of 3: 59.1 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit ortho_npsvd(x2, x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.99999999999999967, 1.6653345369377348e-16)"
      ]
     },
     "execution_count": 196,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res3 = ortho_npsvd(x2, x1)\n",
    "res3 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 검증"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def dist_point_line(coef, point):\n",
    "    return np.abs(point[1] - coef[0] * point[0] - coef[1]) / np.sqrt(coef[0] ** 2 + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.0, 0.5] 2.0\n",
      "(1.0000000000000007, -3.3306690738754696e-16) 1.41421356237\n",
      "(0.99999999999999967, 1.6653345369377348e-16) 1.41421356237\n"
     ]
    }
   ],
   "source": [
    "X = np.hstack((x2[:, np.newaxis], x1[:, np.newaxis]))\n",
    "for res in [res1, res2, res3]:\n",
    "    s = 0\n",
    "    for x in X:\n",
    "        s += dist_point_line(res, x)\n",
    "    print(res, s)"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}