Skip to content

Instantly share code, notes, and snippets.

@jamesmcm

jamesmcm/gist:fd78e66963136c7914dd

Created July 8, 2014 10:05
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save jamesmcm/fd78e66963136c7914dd to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save jamesmcm/fd78e66963136c7914dd to your computer and use it in GitHub Desktop.
Download ZIP
Raw
gistfile1.txt
{
"metadata": {
"name": "",
"signature": "sha256:73456cfce5c2f85c61f377c85ca744d41c8bca573e845d3214eb893e513723da"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "code",
"collapsed": false,
"input": [
"#Write linear model to generate gene expression data from SNPs\n",
"%pylab inline\n",
"\n",
"import numpy as np\n",
"import GPy\n",
"import pylab as pb\n",
"#Columns are samples\n",
"X = np.genfromtxt('ALLsnpspanama.csv', delimiter=',')[1:,1:]\n",
"#Y = WX + N , W is sparse\n",
"print X.shape\n",
"\n",
"SNPsconsider=np.array(range(66))*200\n",
"\n",
"ng=1\n",
"W=np.zeros((ng,13314))\n",
"\n",
"\n",
"W[:,[2000]]=5\n",
"\n",
"Y=np.dot(W,X) + np.random.normal(size=(ng,182))\n",
"\n",
"print np.dot(W,X).shape\n",
"print np.random.normal(size=(1,182)).shape\n",
"pb.matshow(Y)\n",
"pb.colorbar()\n",
"print Y.shape\n",
"pb.matshow(X[1000,:].reshape((1,182)))\n",
"pb.colorbar()\n",
"pb.matshow(X[453,:].reshape((1,182)))\n",
"pb.colorbar()\n",
"#fe=open(\"ALLexprpanama.csv\",\"r\")\n",
"\n",
"\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n",
"(13314, 182)"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"(1, 182)\n",
"(1, 182)\n",
"(1, 182)"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n"
]
},
{
"output_type": "stream",
"stream": "stderr",
"text": [
"WARNING: pylab import has clobbered these variables: ['step']\n",
"`%matplotlib` prevents importing * from pylab and numpy\n"
]
},
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 14,
"text": [
"<matplotlib.colorbar.Colorbar instance at 0x8257dd0>"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAyQAAACDCAYAAAB8+9/PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE3dJREFUeJzt3X1wVFWax/HfhYRXkYKaEEO6y1gkkAQZ0hLFNyQjJAxI\nGHAZCnAHhKCzuKj4NqP7JuigQZ1yULTWQQdx3A047GJwlnTxUumUMMtkpiC4axiT0cRNIolSCJKE\n0KTT+0ewO/d2J3S6kzTB76fqVPW599zznHT3DXk4995jeL1erwAAAAAgCgZEewAAAAAAvrtISAAA\nAABEDQkJAAAAgKghIQEAAAAQNSQkAAAAAKKGhAQAAABA1JCQAAAAAOiWTz75RA6Hw1dGjhypV155\nJay+DNYhAQAAABCutrY2JSYmqrS0VHa7vdvHM0MCAAAAIGz79+/XuHHjwkpGJBISAAAAABHYvn27\nli5dGvbxXLIFAAAAQJI01DDUEmT7qFGjdOrUqYDtbrdbiYmJKi8vV1xcXFgxY8I6CgAAAMAVp0XS\nLyR9Jqmqw/bir78O2r6oqEhTpkwJOxmRSEgAAAAAdHC1pIyL5VvFnbQtKCjQkiVLIorHJVsAAAAA\nJEmGYWhLkO33SbKmDU1NTbr22mtVVVWlESNGhB2TGRIAAAAAPkNDbDd8+HCdPHky4ngkJAAAAAB8\nru7jeCQkAAAAAHxCnSHpKSQkAAAAAHz6OiFhYUQAAAAAPkODlGBOnz6thQsXKi0tTenp6Tp8+HBY\n8ZghAQAAAOAT6vOyHn74Yc2ZM0c7d+5Ua2urmpqaworHY38BAAAASGp/7G9tkO02mR/7e+bMGTkc\nDn322WcRx+SSLQAAAAA+Q2MCi1VVVZXi4uK0YsUK3XDDDbrvvvvU3NwcVjwSEgAAAAA+I4YHFqvW\n1lYdOXJEDzzwgI4cOaLhw4crPz8/rHjcQwIAAADAJ3aw5HK3l87YbDbZbDbdeOONkqSFCxeSkAAA\nAADoAUOkrCFSVodN6//P3OSaa66R3W5XRUWFxo8fr/3792vixIlhheOmdgAAAACS2m9q96YE2V5p\nvqldko4dO6ZVq1bJ7XZr3Lhx2rp1q0aOHNn9mCQkAAAAAKSLCUlGkO1lgQlJT+GSLQAAAAB+Q/o2\nHAkJAAAAAL/BfRuOhAQAAACAX5DH/PYmEhIAAAAAfiFespWUlKSrr75aAwcOVGxsrEpLS8MKR0IC\nAAAAwC/ES7YMw5DL5dLo0aMjCkdCAgAAAMCvG5ds9cSTtwZE3AMAAACAK8eQICUIwzA0c+ZMZWZm\nasuWLWGHY4YEAAAAgF+Il2wdOnRICQkJ+uqrr5Sdna3U1FRNmzat2+FISAAAAAD4DZZcn7eXriQk\nJEiS4uLitGDBApWWlpKQAAAAAIjQVVLWxPbyrfUHzU2am5vl8Xg0YsQINTU1ae/evXr66afDCkdC\nAgAAAMAvhEu2GhoatGDBAklSa2ur7rnnHuXk5IQVzvD2xK3xAAAAAPo9wzDk3Rxk+5qeeaJWMMyQ\nAAAAAPBjpXYAAAAAURPiSu09hYQEAAAAgF+Ij/3tKSyMCAAAAMAvxIURJcnj8cjhcCg3NzfscMyQ\nAAAAAPDrxj0kmzZtUnp6us6ePRt2OGZIAAAAAPgNDlKCqK2t1Z49e7Rq1aqInsBFQgIAAADAL8RL\nth555BG9+OKLGjAgspSChAQAAACA3/AgxeL3v/+9xowZI4fDEfH6JNxDAgAAAMBvsOT67/bSmT/8\n4Q/avXu39uzZo5aWFn3zzTdatmyZ3nnnnW6HY6V2AAAAAJIurtT+ZZDtYzpfqb2kpEQvvfSSPvjg\ng7BiMkMCAAAAwOd8GCu1G4YRdrweuYfE6XQqNTVVKSkp2rhxY090CfQrSUlJ+v73vy+Hw6GbbrpJ\nknTq1CllZ2dr/PjxysnJ0enTp6M8SqB3rFy5UvHx8Zo0aZJvW1ff/+eff14pKSlKTU3V3r17ozFk\noNcEOx/WrVsnm80mh8Mhh8OhoqIi3z7OB1yOzg8eFFC6Mn36dO3evTvseBEnJB6PR2vWrJHT6VR5\nebkKCgp0/PjxSLsF+hXDMORyuXT06FGVlpZKkvLz85Wdna2KigrNmDFD+fn5UR4l0DtWrFghp9Np\n2tbZ97+8vFw7duxQeXm5nE6nHnjgAbW1tUVj2ECvCHY+GIahRx99VEePHtXRo0c1e/ZsSZwPuHy5\nBw4KKL0p4oSktLRUycnJSkpKUmxsrBYvXqzCwsKeGBvQr1ivq9y9e7eWL18uSVq+fLnef//9aAwL\n6HXTpk3TqFGjTNs6+/4XFhZqyZIlio2NVVJSkpKTk31JPHAlCHY+SMGvved8wOXqvAYHlN4UcUJS\nV1cnu93uq9tsNtXV1UXaLdCvGIahmTNnKjMzU1u2bJEkNTQ0KD4+XpIUHx+vhoaGaA4R6FOdff+/\n+OIL2Ww2Xzv+zcB3xauvvqrJkycrLy/Pdwkj5wMuV+c0NKD0pogTkkhuYAGuFIcOHdLRo0dVVFSk\n1157TR9++KFpv2EYnCv4zrrU959zA1e61atXq6qqSmVlZUpISNBjjz3WaVvOB1wOQpkhaWlp0dSp\nU5WRkaH09HQ99dRTYceLOCFJTExUTU2Nr15TU2PK9oHvgoSEBElSXFycFixYoNLSUsXHx6u+vl6S\ndOLECY0ZMyaaQwT6VGfff+u/GbW1tUpMTIzKGIG+MmbMGF9ivmrVKt9lWZwPuFy5NSigWA0ZMkTF\nxcUqKyvTRx99pOLiYh08eDCseBEnJJmZmaqsrFR1dbXcbrd27NihefPmRdot0G80Nzfr7NmzkqSm\npibt3btXkyZN0rx587Rt2zZJ0rZt2zR//vxoDhPoU519/+fNm6ft27fL7XarqqpKlZWVvifTAVeq\nEydO+F7v2rXL9wQuzgdcrpo1LKAEM2xY+3a32y2Px6PRo0eHFS/idUhiYmK0efNmzZo1Sx6PR3l5\neUpLS4u0W6DfaGho0IIFCyRJra2tuueee5STk6PMzEwtWrRIb731lpKSkvTee+9FeaRA71iyZIlK\nSkp08uRJ2e12PfPMM3ryySeDfv/T09O1aNEipaenKyYmRq+//jqXqOCKYj0f1q9fL5fLpbKyMhmG\noeuuu05vvPGGJM4HXL6CzYgE09bWphtuuEGffvqpVq9erfT09LDisVI7AAAAAEnt9zEVebP0ketr\nfeTyryH1b+s/73Sl9jNnzmjWrFnKz89XVlZW92OSkAAAAACQ2hOSQm9OwPYfGXs7TUgk6dlnn9XQ\noUP1+OOPdztmj6zUDgAAAODKEMo9JCdPnvQ9wvrcuXPat2+fHA5HWPEumZA4nU6lpqYqJSVFGzdu\nDCsIAAAAgP7hvAYFFKsTJ07ozjvvVEZGhqZOnarc3FzNmDEjrHhdXrLl8Xg0YcIE7d+/X4mJibrx\nxhtVUFDATesAAADAFcgwDP3a+5OA7fcbv+3ykq1IdDlDUlpaquTkZCUlJSk2NlaLFy9WYWFhrwwE\nAAAAQPQ1a2hA6U1dPva3rq5OdrvdV7fZbPrjH/9oamNMyJIqSnplcAAAAEC/NWK6vN+4oj2KbnMH\nWZm9N3WZkIT0LOyKEulNr/RClvQzl+pXmY+Jtyy9YPy7Zaonw1z98zrz8Zna63t9cqD5jv/vFZr7\n8v7FMt5dlrFmWupzLPU9Xe9fmvOWqV4we6W5wRrL8XMrzeN7abx5/13+l8Y/m3+Wwp3mn2Xehq7H\n9lnGNab6uNUnTPW2Deb+jN9aunv4P0z1IiPXVPf+2HztoPE3ls/x5gvm+vWx/mMXmmMv3Wp+H182\n8kz1WnNPynzTHKvpYXN/w5sqzGP9L/P7bPzJ8j251vI9MX9M0pvqkjGsQ39zzfu8d7b3nfVPkusX\nUtXd5v3X/YOls3stfX9meV+/Z+n/XyyfY7qlvWXtRW+Z5Wcda4n3iuV419u+l7Xex0y7Eg+dMh/7\nkvnYp983x1rvtPRt/W1js4z1by0/258tx1vap9UcMdXLHVPMx5cVm+qb9APf65HeRaZ9y3dZflFZ\n/o/F2GYZS6u56k2yvM+/tBxvmeLO/6G5/ZN/Z+n/fXO17IT5Oz35NfOX1jjc4fgWy9hmWt7XGEus\n/eZqxXZz+5RfmfdbP8fNf28+fx98w3wCHf9pkiRpWVa93nFdozTjUdP+fD1sqv/c8rvJ2OAf7wBX\nk2lf2/zh5rEfNo99/GLLuf+/lvfi+q7fd2+e5XMdY67qdstYZ1j6m2np735L/D0d2h+29P0X83/+\neX96s6l+wXJLZ9wQ8+/80wcSTPWDd5mqmvYry3vzrmVsO93m/QMs14/faxmv5eoOY5S/f++Dlr5v\nM8e+vuBPpvr/fGleGLA23ty3/YeW99m5zlxfaK57nzXHr7JcdZ5c32iqt5WZv1dWf+5w/mb+tc20\nr3qcebDXpn1lPvji3yVZyyTXO5LxA8vnMNfyXtks+4st+2db3ovrLYP9V0vd8neKd7XlO/6fHfp+\nxNx3scvcNrfxS1O98ao4U32Rd5upvuMn95rqRpLlZ9th7v+C+Wuh2J+b68Y0y/F3WN6b8eb9jacH\n+l4P/0fz52b9m8pYbj52Tc0LpvqrN5sHY/zY8jns7PDa8nupvwh2z0hv6jIhSUxMVE1Nja9eU1Mj\nm80W2LBwnXSyWipcp0OSbuvZMQIAAACXv7+6pBpXtEcRsXOdrMzeUU1NjZYtW6Yvv/xShmHo/vvv\n10MPPRRWvC4TkszMTFVWVqq6ulpjx47Vjh07VFBQENjwR+ukT1zSj9bptg/WhzUQAAAAoF9LzpLs\nWf56Xf/8uziUGZLY2Fi9/PLLysjIUGNjo6ZMmaLs7OywHn51yYURi4qKtHbtWnk8HuXl5empp54y\n7c/KylJJCfeQAAAAAB1Nnz5dLpcr2sPoFsMw9KD3hYDtrxo/6/IpW/Pnz9eDDz4Y1qN/WakdAAAA\ngKT2hOR+r/UGQunXxtpOE5Lq6mpNnz5dH3/8sa666qpux+zyki0AAAAA3y3nNEz1rk/U4Prkkm0b\nGxu1cOFCbdq0KaxkRGKGBAAAAMBFhmFokfftgO3vGfcGzJBcuHBBc+fO1ezZs7V27dqwYzJDAgAA\nAMAnlHVIvF6v8vLylJ6eHlEyIl1ipXYAAAAA3y3NGhZQrA4dOqR3331XxcXFcjgccjgccjqdYcVj\nhgQAAACAjzuEx/7efvvtamtru2S7UJCQAAAAAPA5H8IlWz2JhAQAAACATygzJD2Je0gAAAAA+IRy\nD8nKlSsVHx+vSZMmRRyPhAQAAACAj1uDAorVihUrwr6J3YpLtgAAAAD4hHIPybRp01RdXd0j8UhI\nAAAAAPicaxvap/FISAAAAAD4nG/hKVsAAAAAosTdMkhtBz9U26GDfRKPhAQAAACAj7txmJQxq718\n64X8XovHU7YAAAAA+LXEBBaLJUuW6NZbb1VFRYXsdru2bt0adjjD6/V6IxkvAAAAgCuDYRjSsSDp\nwWRDvZU2cMkWAAAAAL9zfRuOhAQAAACAX1PfhiMhAQAAAODX0rfhuKkdAAAAgF9LkBKE0+lUamqq\nUlJStHHjxrDDkZAAAAAA8GsMUiw8Ho/WrFkjp9Op8vJyFRQU6Pjx42GFIyEBAAAA4BfCDElpaamS\nk5OVlJSk2NhYLV68WIWFhWGFIyEBAAAA4BdCQlJXVye73e6r22w21dXVhRWOm9oBAAAA+IXw2F/D\nMHosHAkJAAAAAL8mSdUu6XNXp00SExNVU1Pjq9fU1Mhms4UVjpXaAQAAAEi6OPPxeJD04CXzSu2t\nra2aMGGCDhw4oLFjx+qmm25SQUGB0tLSuh2TGRIAAAAAfucv3SQmJkabN2/WrFmz5PF4lJeXF1Yy\nIjFDAgAAAOAiwzCkFUHSg63mGZKexAwJAAAAAL8+XqmdhAQAAACAHwkJAAAAgKg527fhWBgRAAAA\ngN/5IKUbfve732nixIkaOHCgjhw5csn2JCQAAAAA/M4FKd0wadIk7dq1S3fccUdI7blkCwAAAIBf\nN2dErFJTU7vVnoQEAAAAgF9j34YjIQEAAADgF8JTtrKzs1VfXx+w/bnnnlNubm63wpGQAAAAAPBr\nkXTBJbW6Om2yb9++HgtHQgIAAADA76wkZV0s31ofVlehrO7OU7YAAAAA+LUGKd2wa9cu2e12HT58\nWHfddZdmz57dZXvDG0raAgAAAOCKZxiGpGDpgRHSbEc4uGQLAAAAQAcX+jQaCQkAAACADs72aTQS\nEgAAAAAdNPdpNG5qBwAAANDBuSAldE888YTS0tI0efJk3X333Tpz5kyX7UlIAAAAAHRwNkgJXU5O\njj7++GMdO3ZM48eP1/PPP99lexISAAAAAB1ENkOSnZ2tAQPa04ypU6eqtra2y/YkJAAAAAA6iCwh\n6eg3v/mN5syZ02UbbmoHAAAA0ME3l2yRnZ2t+vr6gO3PPfeccnNzJUkbNmzQoEGDtHTp0i77IiEB\nAAAA0ME5SeWSjnfaYt++fV328Pbbb2vPnj06cODAJaORkAAAAADo4Jyk6y6Wb+0K+Win06kXX3xR\nJSUlGjJkyCXbG97eWgMeAAAAQL9iGIakV4LseUihpg0pKSlyu90aPXq0JOmWW27R66+/3ml7ZkgA\nAAAAdHDpe0i6UllZ2a32zJAAAAAAkPTtDEmgUaNG6dSpU70SkxkSAAAAAJIU8mVZPYl1SAAAAABE\nDQkJAAAAgKghIQEAAAAQNSQkAAAAAKKGhAQAAABA1JCQAAAAAIia/wcNZfndBgU8hwAAAABJRU5E\nrkJggg==\n",
"text": [
"<matplotlib.figure.Figure at 0x47b6850>"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAyYAAACICAYAAAASyf8xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFnhJREFUeJzt3X1wVNX9x/HPxV0wPPx4NpLdbRfcNbsMEFcT0KFUdKwB\nZtxC7TihnYqY4g5tftZpZ6y1nalYR41jO1piO9HBB9oa4kwti51kcUIbVBxYRVB/hWpwjG5WzRgh\n8hAhZs3vD3A3u3v3IYawGN6vmTNy7/3u+V7DXpgv555zjP7+/n4BAAAAQAGNKvQNAAAAAACFCQAA\nAICCozABAAAAUHAUJgAAAAAKjsIEAAAAQMFRmAAAAAAoOAoTAAAAAINy8803q7i4WHPnzs0Yc+ut\nt8rtdqusrEx79uzJ2SeFCQAAAIBBWb16tUKhUMbrTU1NOnDggNra2vToo49q7dq1OfukMAEAAAAw\nKIsWLdLkyZMzXt+yZYtWrVolSVqwYIG6u7vV2dmZtU8KEwAAAACnVTQalcPhiB/b7XZ1dHRk/QyF\nCQAAAIDTrr+/P+nYMIys8RQmAAAAACRJRYYhw6RNmDBhUP3YbDZFIpH4cUdHh2w2W9bPWL7SHQMA\nAAAYcY5Lusfk/G+OHh1UP36/X3V1daqqqtLOnTs1adIkFRcXZ/0MhQkAAACAuP/JI2blypXavn27\nurq65HA4tG7dOn3++eeSpEAgoGXLlqmpqUkul0vjxo3TE088kbNPoz/15S8AAAAA5yTDMPSYyfk1\nSp8zcroxYgIAAAAgrqhAeSlMAAAAAMTl8yrXcGBVLgAAAABxRSbNTCgUksfjkdvtVm1tbdr1Q4cO\nacWKFSorK9OCBQv0n//8J2teChMAAAAAcfkUJrFYTDU1NQqFQtq3b58aGhq0f//+pJh7771Xl156\nqV5//XVt3LhRP/vZz7LmpTABAAAAEJdPYRIOh+VyueR0OmW1WlVVVaVgMJgUs3//fl111VWSpNLS\nUrW3t+vjjz/OmJfCBAAAAEDcBJOWKhqNyuFwxI/tdrui0WhSTFlZmZ599llJJwuZ9957Tx0dHRnz\nUpgAAAAAiBtr0lIZhpGznzvuuEPd3d3y+Xyqq6uTz+fTeeedlzGeVbkAAAAAxBVZpJe+kHYM3LYk\nZQsTm82mSCQSP45EIrLb7UkxEyZM0OOPPx4/njlzpmbNmpUxLxssAgAAAJB0ciSkd2L6+dGfJm+w\n2NfXp9LSUm3btk0lJSWaP3++Ghoa5PV64zGffvqpioqKNHr0aD322GPasWOHnnzyyYy5GTEBAAAA\nEGcdkzvGYrGorq5OlZWVisViqq6ultfrVX19vSQpEAho3759uummm2QYhubMmaMNGzZk7ZMREwAA\nAACSTo6Y9H/D5Pz7ySMmw4EREwAAAAAJeYyYDAcKEwAAAAAJ4wqTlsIEAAAAQML5hUnLPiYAAAAA\nEsaYNBOhUEgej0dut1u1tbVp17u6urRkyRJdcsklmjNnTtYVuSQmvwMAAAA4xTAM9S8zOd+UPPk9\nFouptLRULS0tstlsqqioSFsu+K677tKJEyd03333qaurS6Wlpers7JTFYv7SFiMmAAAAABLON2kp\nwuGwXC6XnE6nrFarqqqqFAwGk2JmzJihw4cPS5IOHz6sqVOnZixKJOaYAAAAABgoj1W5otGoHA5H\n/Nhut2vXrl1JMWvWrNHVV1+tkpISHTlyRM8880zWPhkxAQAAAJAwzqSlMAwjZzf33nuvLrnkEn3w\nwQfau3evfvrTn+rIkSMZ4xkxAQAAAJBwvtTaIbVGM4fYbDZFIpH4cSQSkd1uT4p5+eWX9etf/1qS\ndNFFF2nmzJl66623VF5ebtonhQkAAACAhDHS4otOti+teyU5pLy8XG1tbWpvb1dJSYkaGxvV0NCQ\nFOPxeNTS0qKFCxeqs7NTb731lmbNmpUxLYUJAAAAgIQ85phYLBbV1dWpsrJSsVhM1dXV8nq9qq+v\nlyQFAgHdeeedWr16tcrKyvTFF1/ogQce0JQpUzL2yXLBAAAAACSdWi74XpPzdyYvFzwcGDEBAAAA\nkJDHiMlwoDABAAAAkFCgwoTlggEAAAAk5LFcsCSFQiF5PB653W7V1tamXX/wwQfl8/nk8/k0d+5c\nWSwWdXd3Z0zLHBMAAAAAkk7NMWkwOb8yeY5JLBZTaWmpWlpaZLPZVFFRoYaGBnm9XtN+//nPf+qh\nhx5SS0tLxtyMmAAAAABIGGPSUoTDYblcLjmdTlmtVlVVVSkYDGbs8umnn9bKlSuzpqUwAQAAAJBw\nvklLEY1G5XA44sd2u13RqPmOjD09Pdq6dauuv/76rGmZ/A4AAAAgIcOckoEMw8i7u+eee07f+ta3\nNGnSpKxxFCYAAAAAEsZIrbul1tcyh9hsNkUikfhxJBKR3W43jd20aVPO17gkJr8DAAAAOMUwDPXv\nNTl/SfLk976+PpWWlmrbtm0qKSnR/PnzTSe/f/rpp5o1a5Y6OjpUVFSUNTcjJgAAAAAS8niVy2Kx\nqK6uTpWVlYrFYqqurpbX61V9fb0kKRAISJI2b96sysrKnEWJxIgJAAAAgFMMw1D/+ybnv5E8YjIc\nGDEBAAAAkGCyCteZQGECAAAAIO5EHq9yDYfTso9Jru3ogZHO6XRq3rx58vl8mj9/viTp4MGD+s53\nvqOLL75Y1157rbq7uwt8l8DwuPnmm1VcXKy5c+fGz2X7/t93331yu93yeDx6/vnnC3HLwLAxex7u\nuusu2e12+Xw++Xw+NTc3x6/xPOBsdGLM6LRmJp8aoLW1VT6fT3PmzNHixYuz5h3yHJPBbkcPjEQz\nZ87U7t27NWXKlPi522+/XdOmTdPtt9+u2tpaHTp0SPfff38B7xIYHi+++KLGjx+vG2+8UW+++aak\nzN//ffv26Qc/+IFeeeUVRaNRXXPNNXr77bc1ahT7/WJkMHse1q1bpwkTJujnP/95UizPA85GhmHo\n4/7xaeenG0eT5pjkUwN0d3dr4cKF2rp1q+x2u7q6ujRt2rSMuYf8zR/sdvTASJVa42/ZskWrVq2S\nJK1atUqbN28uxG0Bw27RokWaPHly0rlM3/9gMKiVK1fKarXK6XTK5XIpHA6f8XsGhovZ8yCZTxrm\necDZ6oTGpLVU+dQATz/9tK6//vr4/ibZihLpNBQmg9mOHhipDMPQNddco/Lycj322GOSpM7OThUX\nF0uSiouL1dnZWchbBM6oTN//Dz74IGkDLv7OwLli/fr1KisrU3V1dfzVRp4HnK0+U1FaS5VPDdDW\n1qaDBw/qqquuUnl5uf7yl79kzTvkwmQw29EDI9WOHTu0Z88eNTc365FHHtGLL76YdN0wDJ4VnLNy\nff95NjDSrV27Vu+++6727t2rGTNm6Be/+EXGWJ4HnA3yGTHJ57v6+eef67XXXlNTU5O2bt2q3/3u\nd2pra8sYP+RVuQazHT0wUs2YMUOSNH36dK1YsULhcFjFxcX66KOPdOGFF+rDDz/UBRdcUOC7BM6c\nTN//1L8zOjo6ZLPZCnWbwBkx8M//H//4x7ruuusk8Tzg7NWr0Xq19Zhebe3JGJNPDeBwODRt2jQV\nFRWpqKhI3/72t/X666/L7Xab9jnkEZPy8nK1tbWpvb1dvb29amxslN/vH2q3wNdGT0+Pjhw5Ikk6\nduyYnn/+ec2dO1d+v19PPfWUJOmpp57S8uXLC3mbwBmV6fvv9/u1adMm9fb26t1331VbW1t8JTtg\npPrwww/jv/7HP/4RX7GL5wFnqx6N1ezF03XjXd+Mt1T51ADf/e539dJLLykWi6mnp0e7du3S7Nmz\nM+Yd8ohJpu3ogXNFZ2enVqxYIUnq6+vTD3/4Q1177bUqLy/XDTfcoA0bNsjpdOqZZ54p8J0Cw2Pl\nypXavn27urq65HA4dPfdd+uOO+4w/f7Pnj1bN9xwg2bPni2LxaI//elPvLqCESX1eVi3bp1aW1u1\nd+9eGYahmTNnqr6+XhLPA85evTJfHnigTDXAl9/vQCAgj8ejJUuWaN68eRo1apTWrFmTtTAZ8nLB\nAAAAAEYGwzDU3L847fxSo9V0dbnTiZ3fAQAAAMTlM2IyHChMAAAAAMT1aGxB8uac/J7PVvMAAAAA\nRoYTGp3WzOSqE1pbWzVx4kT5fD75fD7dc889WfNmHTGJxWKqqalJ2mre7/czuR0AAAAYoXpN9i1J\nlW+dcOWVV2rLli155c06YpLPVvMAAAAARo4eFaW1VPnWCYOZMJ91xMRsq/ldu3YlxTgNQ+/lnQ4A\nAAA4V3xT/f3thb6JQctnxCSfOsEwDL388ssqKyuTzWbTgw8++NX3MclnLe33JP1W0pOSbpK0Tr/N\n+ZmBfqt1ecfm6ju1r7MtPtXAzw8m9qvcy2Dvbaj5B17P9dkzbag/y3z+356U+fMwmJ9bPtdTne58\nQ8mdy+nuL1f/Q3neBtP36YgfiuH+uX7V/E/q9DwT2WIH63T/vgzleRvuZ3uw/Q13/sEYat+n+zt/\ner6TT0q66Yz/WZTr80PJfbo/P5z9Dfefg6nZvo4yzSkZKJ864dJLL1UkEtHYsWPV3Nys5cuX6+23\n384Yn7UwyWereUlqldR96r9SuyRnzhsFAAAARpb2U+3r7TON1Xut7XqvNfN7UfnUCRMmTIj/eunS\npfrJT36igwcPasqUKaZ9Zi1MBm41X1JSosbGRjU0NKTFLdbJ34LFkrZTlAAAAOCc5FTyP9BvL8xt\nDNEJjdaFiy/WhYsvjp97ad0LSTH51AmdnZ264IILZBiGwuGw+vv7MxYlUh47vzc3N+u2226LbzX/\nq1/9Kun64sWLtX371/OHDgAAAAyXK6+8Uq2trYW+jUExDEP/2/9A2vn1xu1pE9nN6oT6+npJUiAQ\n0COPPKI///nPslgsGjt2rP7whz/o8ssvz5w7V2ECAAAA4NxgGIZu6X8o7fyjxm2DWmHrq2DndwAA\nAABxnxVo53cKEwAAAABx+azKNRyybrAIAAAA4NzSqzFpzUwoFJLH45Hb7VZtbW3G/l555RVZLBY9\n++yzWfMyYgIAAAAgriePV7lisZhqamrU0tIim82miooK+f1+eb3etLhf/vKXWrJkSc45KoyYAAAA\nAIjr1ei0liocDsvlcsnpdMpqtaqqqkrBYDAtbv369fr+97+v6dOn58xLYQIAAAAg7oTGpLVU0WhU\nDocjfmy32xWNRtNigsGg1q5dKyn3bvG8ygUAAAAgzmyEJFWuIkOSbrvtNt1///0yDEP9/f05X+Wi\nMAEAAAAQ16OxOtb6qnpaX80YY7PZFIlE4seRSER2uz0pZvfu3aqqqpIkdXV1qbm5WVarVX6/37RP\nNlgEAAAAIOnkSMhF/f+Xdv4dY07SiEdfX59KS0u1bds2lZSUaP78+WpoaEib/P6l1atX67rrrtP3\nvve9jLkZMQEAAAAQZzanJJXFYlFdXZ0qKysVi8VUXV0tr9er+vp6SVIgEBh0XkZMAAAAAEg6OWIy\nNdaRdv6T8+w554gMFSMmAAAAAOJOHM89YjIcKEwAAAAAxPUez70q13BgHxMAAAAAcb1Hx6Y1M6FQ\nSB6PR263W7W1tWnXg8GgysrK5PP5dNlll+lf//pX1rzMMQEAAAAg6dT+JG+ZlAelRtIck1gsptLS\nUrW0tMhms6mioiJtVa5jx45p3LhxkqQ333xTK1as0IEDBzLmZsQEAAAAQMJxk5YiHA7L5XLJ6XTK\narWqqqpKwWAwKebLokSSjh49qmnTpmVNS2ECAAAAIOEzk5YiGo3K4XDEj+12u6LRaFrc5s2b5fV6\ntXTpUv3xj3/MmpbCBAAAAEDCMZOWwjCMvLpavny59u/fr+eee04/+tGPssayKhcAAACAhOOS3miV\n3mzNGGKz2RSJROLHkUhEdrs9Y/yiRYvU19enTz75RFOnTjWNoTABAAAAkHBc0sWLT7YvPb0uKaS8\nvFxtbW1qb29XSUmJGhsb1dDQkBTzzjvvaNasWTIMQ6+99pokZSxKJAoTAAAAAAMdzR1isVhUV1en\nyspKxWIxVVdXy+v1qr6+XpIUCAT097//XRs3bpTVatX48eO1adOmrH2yXDAAAAAASafmjtSblAeB\n5OWChwMjJgAAAAASTJYHPhMoTAAAAAAkmCwPfCawXDAAAACAhDyWC5akUCgkj8cjt9ut2tratOt/\n+9vfVFZWpnnz5mnhwoV64403sqZlxAQAAABAQh4jJrFYTDU1NWppaZHNZlNFRYX8fr+8Xm88Ztas\nWXrhhRc0ceJEhUIh3XLLLdq5c2fGPhkxAQAAAJBwwqSlCIfDcrlccjqdslqtqqqqUjAYTIq54oor\nNHHiREnSggUL1NHRkTUthQkAAACAhKMmLUU0GpXD4Ygf2+12RaPRjF1u2LBBy5Yty5qWV7kAAAAA\nJOSxKpdhGHl39+9//1uPP/64duzYkTWOwgQAAABAwnFJH7dKXa0ZQ2w2myKRSPw4EonIbrenxb3x\nxhtas2aNQqGQJk+enDUtGywCAAAAkHRqJGSpSXnQnLzBYl9fn0pLS7Vt2zaVlJRo/vz5amhoSJr8\n/v777+vqq6/WX//6V11++eU5czNiAgAAACDBZE5JKovForq6OlVWVioWi6m6ulper1f19fWSpEAg\noLvvvluHDh3S2rVrJUlWq1XhcDhjn4yYAAAAAJB0asSk3KQ8eDV5xGQ4MGICAAAAIMFkeeAzgcIE\nAAAAQEIer3INBwoTAAAAAAl5LBc8HNhgEQAAAEDCcZNmIhQKyePxyO12q7a2Nu36f//7X11xxRU6\n//zz9fvf/z5nWkZMAAAAACQcyR0Si8VUU1OjlpYW2Ww2VVRUyO/3Jy0XPHXqVK1fv16bN2/OKy0j\nJgAAAAAS+kxainA4LJfLJafTKavVqqqqKgWDwaSY6dOnq7y8XFarNa+0FCYAAAAABiUajcrhcMSP\n7Xa7otHokPrkVS4AAAAAA3yeM8IwjNOelcIEAAAAwABHJL0kaUfGCJvNpkgkEj+ORCKy2+1Dykph\nAgAAAGCAHkmXnmpfeiApory8XG1tbWpvb1dJSYkaGxvV0NBg2lu+O8ZTmAAAAAAY4LOcERaLRXV1\ndaqsrFQsFlN1dbW8Xq/q6+slSYFAQB999JEqKip0+PBhjRo1Sg8//LD27dun8ePHm/Zp9OdbwgAA\nAAAY0U7OHdltcuWyvEc+vipGTAAAAAAMkHvEZDhQmAAAAAAYoDCFCfuYAAAAABjgM5OWLhQKyePx\nyO12q7a21jTm1ltvldvtVllZmfbs2ZM1K4UJAAAAgAEOm7RksVhMNTU1CoVC2rdvnxoaGrR///6k\nmKamJh04cEBtbW169NFHtXbt2qxZKUwAAAAADJB7xCQcDsvlcsnpdMpqtaqqqkrBYDApZsuWLVq1\napUkacGCBeru7lZnZ2fGrBQmAAAAAAbIXZhEo1E5HI74sd1uVzQazRnT0dGRMSuT3wEAAAAMkP7q\nVqqTywrnlrrEcLbPUZgAAAAAGOA3aWdSN0W02WyKRCLx40gkIrvdnjWmo6NDNpstY1Ze5QIAAAAg\n6eQIh1k7cuRIUlx5ebna2trU3t6u3t5eNTY2yu/3J8X4/X5t3LhRkrRz505NmjRJxcXFGXMzYgIA\nAABgUCwWi+rq6lRZWalYLKbq6mp5vV7V19dLkgKBgJYtW6ampia5XC6NGzdOTzzxRNY+jf7h3lse\nAAAAAHLgVS4AAAAABUdhAgAAAKDgKEwAAAAAFByFCQAAAICCozABAAAAUHAUJgAAAAAKjsIEAAAA\nQMFRmAAAAAAouP8HE7/o5pQM9Z4AAAAASUVORK5CYII=\n",
"text": [
"<matplotlib.figure.Figure at 0x47c7e90>"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAyYAAACICAYAAAASyf8xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFo9JREFUeJzt3X1sW9X9x/HPLXZL+vDrM6GxvbnFJnbVNhiSFtR1FMRI\nWwmvHRNKN41SsmJ1y4+hTWKMTRplCAhiE1vDpoDKQ7eRBmmsLlPionRLgaLWUFrgt3aQIgKOgYjQ\nhj6ENsTk90eLHdvXDyFJ3YX3SzpS7/XX59wmvkm++t5zjtHf398vAAAAACigMYW+AAAAAAAgMQEA\nAABQcCQmAAAAAAqOxAQAAABAwZGYAAAAACg4EhMAAAAABUdiAgAAAGBQbr75ZhUXF2v+/PkZY269\n9Va53W6VlZVp3759OfskMQEAAAAwKGvXrlUoFMr4elNTkw4dOqS2tjY98sgjWr9+fc4+SUwAAAAA\nDMqSJUs0derUjK9v27ZNa9askSQtWrRI3d3d6uzszNoniQkAAACAYRWNRuVwOOLHdrtdHR0dWd9D\nYgIAAABg2PX39ycdG4aRNZ7EBAAAAIAkqcgwZJi0SZMmDaofm82mSCQSP+7o6JDNZsv6HsuXumIA\nAAAAo85JSfeYnP/V8eOD6sfv96uurk5VVVXavXu3pkyZouLi4qzvITEBAAAAEPc/ecSsXr1aO3fu\nVFdXlxwOhzZs2KDPPvtMkhQIBLRixQo1NTXJ5XJpwoQJevzxx3P2afSnPvwFAAAA4CvJMAw9anJ+\nndLnjAw3KiYAAAAA4ooKNC6JCQAAAIC4fB7lGgmsygUAAAAgrsikmQmFQvJ4PHK73aqtrU17/ciR\nI1q1apXKysq0aNEi/fvf/846LokJAAAAgLh8EpNYLKaamhqFQiEdOHBADQ0NOnjwYFLMvffeq0sv\nvVSvvfaaNm/erJ/85CdZxyUxAQAAABCXT2ISDoflcrnkdDpltVpVVVWlYDCYFHPw4EFdddVVkqTS\n0lK1t7fro48+yjguiQkAAACAuEkmLVU0GpXD4Ygf2+12RaPRpJiysjI988wzkk4nMu+++646Ojoy\njktiAgAAACBuvElLZRhGzn7uuOMOdXd3y+fzqa6uTj6fT+edd17GeFblAgAAABBXZJFe/FzaNXDb\nkpQtTGw2myKRSPw4EonIbrcnxUyaNEmPPfZY/Hj27NmaM2dOxnHZYBEAAACApNOVkN7J6efHfpK8\nwWJfX59KS0u1Y8cOlZSUaOHChWpoaJDX643HfPLJJyoqKtLYsWP16KOPateuXXriiScyjk3FBAAA\nAECcdVzuGIvForq6OlVWVioWi6m6ulper1f19fWSpEAgoAMHDuimm26SYRiaN2+eNm3alLVPKiYA\nAAAAJJ2umPR/zeT8e8kVk5FAxQQAAABAQh4Vk5FAYgIAAAAgYUJhhiUxAQAAAJBwfmGGZR8TAAAA\nAAnjTJqJUCgkj8cjt9ut2tratNe7urq0bNkyXXLJJZo3b17WFbkkJr8DAAAAOMMwDPWvMDnflDz5\nPRaLqbS0VC0tLbLZbKqoqEhbLviuu+7SqVOndN9996mrq0ulpaXq7OyUxWL+0BYVEwAAAAAJ55u0\nFOFwWC6XS06nU1arVVVVVQoGg0kxs2bN0tGjRyVJR48e1fTp0zMmJRJzTAAAAAAMlMeqXNFoVA6H\nI35st9u1Z8+epJh169bp6quvVklJiY4dO6ann346a59UTAAAAAAkTDBpKQzDyNnNvffeq0suuUTv\nv/++9u/frx//+Mc6duxYxngqJgAAAAASzpdaO6TWaOYQm82mSCQSP45EIrLb7UkxL730kn75y19K\nki666CLNnj1bb775psrLy037JDEBAAAAkDBOWnrR6faFDS8nh5SXl6utrU3t7e0qKSlRY2OjGhoa\nkmI8Ho9aWlq0ePFidXZ26s0339ScOXMyDktiAgAAACAhjzkmFotFdXV1qqysVCwWU3V1tbxer+rr\n6yVJgUBAd955p9auXauysjJ9/vnneuCBBzRt2rSMfbJcMAAAAABJZ5YLvtfk/J3JywWPBComAAAA\nABLyqJiMBBITAAAAAAkFSkxYLhgAAABAQh7LBUtSKBSSx+OR2+1WbW1t2usPPvigfD6ffD6f5s+f\nL4vFou7u7ozDMscEAAAAgKQzc0waTM6vTp5jEovFVFpaqpaWFtlsNlVUVKihoUFer9e033/84x96\n6KGH1NLSknFsKiYAAAAAEsaZtBThcFgul0tOp1NWq1VVVVUKBoMZu3zqqae0evXqrMOSmAAAAABI\nON+kpYhGo3I4HPFju92uaNR8R8aenh5t375d119/fdZhmfwOAAAAICHDnJKBDMPIu7tnn31W3/jG\nNzRlypSscSQmAAAAABLGSa17pdZXM4fYbDZFIpH4cSQSkd1uN43dsmVLzse4JCa/AwAAADjDMAz1\n7zc5f0ny5Pe+vj6VlpZqx44dKikp0cKFC00nv3/yySeaM2eOOjo6VFRUlHVsKiYAAAAAEvJ4lMti\nsaiurk6VlZWKxWKqrq6W1+tVfX29JCkQCEiStm7dqsrKypxJiUTFBAAAAMAZhmGo/z2T819LrpiM\nBComAAAAABJMVuE6G0hMAAAAAMSdyuNRrpEwLPuY5NqOHhjtnE6nFixYIJ/Pp4ULF0qSDh8+rG99\n61u6+OKLde2116q7u7vAVwmMjJtvvlnFxcWaP39+/Fy2z/99990nt9stj8ej5557rhCXDIwYs/vh\nrrvukt1ul8/nk8/nU3Nzc/w17geci06NG5vWzOSTA7S2tsrn82nevHlaunRp1nGHPMdksNvRA6PR\n7NmztXfvXk2bNi1+7vbbb9eMGTN0++23q7a2VkeOHNH9999fwKsERsYLL7ygiRMn6sYbb9Qbb7wh\nKfPn/8CBA/re976nl19+WdFoVNdcc43eeustjRnDfr8YHczuhw0bNmjSpEn66U9/mhTL/YBzkWEY\n+qh/Ytr5mcbxpDkm+eQA3d3dWrx4sbZv3y673a6uri7NmDEj49hD/uQPdjt6YLRKzfG3bdumNWvW\nSJLWrFmjrVu3FuKygBG3ZMkSTZ06Nelcps9/MBjU6tWrZbVa5XQ65XK5FA6Hz/o1AyPF7H6QzCcN\ncz/gXHVK49JaqnxygKeeekrXX399fH+TbEmJNAyJyWC2owdGK8MwdM0116i8vFyPPvqoJKmzs1PF\nxcWSpOLiYnV2dhbyEoGzKtPn//3330/agIvfGfiq2Lhxo8rKylRdXR1/tJH7AeeqT1WU1lLlkwO0\ntbXp8OHDuuqqq1ReXq4///nPWccdcmIymO3ogdFq165d2rdvn5qbm/Xwww/rhRdeSHrdMAzuFXxl\n5fr8c29gtFu/fr3eeecd7d+/X7NmzdLPfvazjLHcDzgX5FMxyeez+tlnn+nVV19VU1OTtm/frt/8\n5jdqa2vLGD/kVbkGsx09MFrNmjVLkjRz5kytWrVK4XBYxcXF+vDDD3XhhRfqgw8+0AUXXFDgqwTO\nnkyf/9TfGR0dHbLZbIW6TOCsGPjz/4c//KGuu+46SdwPOHf1aqxeaT2hV1p7MsbkkwM4HA7NmDFD\nRUVFKioq0je/+U299tprcrvdpn0OuWJSXl6utrY2tbe3q7e3V42NjfL7/UPtFviv0dPTo2PHjkmS\nTpw4oeeee07z58+X3+/Xk08+KUl68skntXLlykJeJnBWZfr8+/1+bdmyRb29vXrnnXfU1tYWX8kO\nGK0++OCD+L///ve/x1fs4n7AuapH4zV36UzdeNfX4y1VPjnAt7/9bb344ouKxWLq6enRnj17NHfu\n3IzjDrlikmk7euCrorOzU6tWrZIk9fX16fvf/76uvfZalZeX64YbbtCmTZvkdDr19NNPF/hKgZGx\nevVq7dy5U11dXXI4HLr77rt1xx13mH7+586dqxtuuEFz586VxWLRH//4Rx5dwaiSej9s2LBBra2t\n2r9/vwzD0OzZs1VfXy+J+wHnrl6ZLw88UKYc4IvPdyAQkMfj0bJly7RgwQKNGTNG69aty5qYDHm5\nYAAAAACjg2EYau5fmnZ+udFqurrccGLndwAAAABx+VRMRgKJCQAAAIC4Ho0vyLg5J7/ns9U8AAAA\ngNHhlMamNTO58oTW1lZNnjxZPp9PPp9P99xzT9Zxs1ZMYrGYampqkraa9/v9TG4HAAAARqlek31L\nUuWbJ1x55ZXatm1bXuNmrZjks9U8AAAAgNGjR0VpLVW+ecJgJsxnrZiYbTW/Z8+epBjDcEp6N+8B\nAQAAgK+Gr6u/v73QFzFo+VRM8ssTDL300ksqKyuTzWbTgw8++OX3MclvLe13Jf1a0hOSbtKvtSHp\n1Q36ddLxUF5PfW2oUsdONdz/l6Fc20hfy2C/Frlk+74N9tpzXctwX3uqbP1n6vsJSTcNsq98+h/O\nz9hg+xvu78tg3z+Uz2Cu9w/2WlMN9v4a6fhs7x3pz1BmTyif3xGpBsYP99dxuH8HDHb8kbyWkf6/\npcr2/uG+v4b6syXVUH6HfNmv2xM6/Tui0N+nwTjbP+OH+hkfzLWkGt6v8/D+7Xq2ZJpTMlA+ecKl\nl16qSCSi8ePHq7m5WStXrtRbb72VMT5rYpLPVvOntUrqltSqdknOnJcJAAAAjC7tkk7/Xfzf7VON\n17ut7Xq3NfNTUfnkCZMmTYr/e/ny5frRj36kw4cPa9q0aaZ9Zk1MBm41X1JSosbGRjU0NJhELtXp\nb8VSObUzW5cAAADAqOSUdPrv4i/8d/5dfEpjdeHSi3Xh0ovj517c8HxSTD55Qmdnpy644AIZhqFw\nOKz+/v6MSYmUx87vzc3Nuu222+Jbzf/iF79Ien3p0qXaufO/84sOAAAAjJQrr7xSra2thb6MQTEM\nQ//b/0Da+Y3G7WkT2c3yhPr6eklSIBDQww8/rD/96U+yWCwaP368fve73+nyyy/PPHauxAQAAADA\nV4NhGLql/6G0848Ytw1qha0vg53fAQAAAMR9WqCd30lMAAAAAMTlsyrXSMi6wSIAAACAr5ZejUtr\nZkKhkDwej9xut2prazP29/LLL8tiseiZZ57JOi4VEwAAAABxPXk8yhWLxVRTU6OWlhbZbDZVVFTI\n7/fL6/Wmxf385z/XsmXLcs5RoWICAAAAIK5XY9NaqnA4LJfLJafTKavVqqqqKgWDwbS4jRs36rvf\n/a5mzpyZc1wSEwAAAABxpzQuraWKRqNyOBzxY7vdrmg0mhYTDAa1fv16Sbl3i+dRLgAAAABxZhWS\nVLmSDEm67bbbdP/998swDPX39+d8lIvEBAAAAEBcj8brROsr6ml9JWOMzWZTJBKJH0ciEdnt9qSY\nvXv3qqqqSpLU1dWl5uZmWa1W+f1+0z7ZYBEAAACApNOVkIv6/y/t/NvGvKSKR19fn0pLS7Vjxw6V\nlJRo4cKFamhoSJv8/oW1a9fquuuu03e+852MY1MxAQAAABBnNqcklcViUV1dnSorKxWLxVRdXS2v\n16v6+npJUiAQGPS4VEwAAAAASDpdMZke60g7//F59pxzRIaKigkAAACAuFMnc1dMRgKJCQAAAIC4\n3pO5V+UaCexjAgAAACCu9/j4tGYmFArJ4/HI7XartrY27fVgMKiysjL5fD5ddtll+uc//5l1XOaY\nAAAAAJB0Zn+SN03Sg1IjaY5JLBZTaWmpWlpaZLPZVFFRkbYq14kTJzRhwgRJ0htvvKFVq1bp0KFD\nGcemYgIAAAAg4aRJSxEOh+VyueR0OmW1WlVVVaVgMJgU80VSIknHjx/XjBkzsg5LYgIAAAAg4VOT\nliIajcrhcMSP7Xa7otFoWtzWrVvl9Xq1fPly/eEPf8g6LIkJAAAAgIQTJi2FYRh5dbVy5UodPHhQ\nzz77rH7wgx9kjWVVLgAAAAAJJyW93iq90ZoxxGazKRKJxI8jkYjsdnvG+CVLlqivr08ff/yxpk+f\nbhpDYgIAAAAg4aSki5eebl94akNSSHl5udra2tTe3q6SkhI1NjaqoaEhKebtt9/WnDlzZBiGXn31\nVUnKmJRIJCYAAAAABjqeO8Risaiurk6VlZWKxWKqrq6W1+tVfX29JCkQCOhvf/ubNm/eLKvVqokT\nJ2rLli1Z+2S5YAAAAACSzswdqTdJDwLJywWPBComAAAAABJMlgc+G0hMAAAAACSYLA98NrBcMAAA\nAICEPJYLlqRQKCSPxyO3263a2tq01//617+qrKxMCxYs0OLFi/X6669nHZaKCQAAAICEPComsVhM\nNTU1amlpkc1mU0VFhfx+v7xebzxmzpw5ev755zV58mSFQiHdcsst2r17d8Y+qZgAAAAASDhl0lKE\nw2G5XC45nU5ZrVZVVVUpGAwmxVxxxRWaPHmyJGnRokXq6OjIOiyJCQAAAICE4yYtRTQalcPhiB/b\n7XZFo9GMXW7atEkrVqzIOiyPcgEAAABIyGNVLsMw8u7uX//6lx577DHt2rUraxyJCQAAAICEk5I+\napW6WjOG2Gw2RSKR+HEkEpHdbk+Le/3117Vu3TqFQiFNnTo167BssAgAAABA0plKyHKT9KA5eYPF\nvr4+lZaWaseOHSopKdHChQvV0NCQNPn9vffe09VXX62//OUvuvzyy3OOTcUEAAAAQILJnJJUFotF\ndXV1qqysVCwWU3V1tbxer+rr6yVJgUBAd999t44cOaL169dLkqxWq8LhcMY+qZgAAAAAkHSmYlJu\nkh68klwxGQlUTAAAAAAkmCwPfDaQmAAAAABIyONRrpFAYgIAAAAgIY/lgkcCGywCAAAASDhp0kyE\nQiF5PB653W7V1tamvf6f//xHV1xxhc4//3z99re/zTksFRMAAAAACcdyh8RiMdXU1KilpUU2m00V\nFRXy+/1JywVPnz5dGzdu1NatW/MalooJAAAAgIQ+k5YiHA7L5XLJ6XTKarWqqqpKwWAwKWbmzJkq\nLy+X1WrNa1gSEwAAAACDEo1G5XA44sd2u13RaHRIffIoFwAAAIABPssZYRjGsI9KYgIAAABggGOS\nXpS0K2OEzWZTJBKJH0ciEdnt9iGNSmICAAAAYIAeSZeeaV94ICmivLxcbW1tam9vV0lJiRobG9XQ\n0GDaW747xpOYAAAAABjg05wRFotFdXV1qqysVCwWU3V1tbxer+rr6yVJgUBAH374oSoqKnT06FGN\nGTNGv//973XgwAFNnDjRtE+jP98UBgAAAMCodnruyF6TVy7Lu/LxZVExAQAAADBA7orJSCAxAQAA\nADBAYRIT9jEBAAAAMMCnJi1dKBSSx+OR2+1WbW2tacytt94qt9utsrIy7du3L+uoJCYAAAAABjhq\n0pLFYjHV1NQoFArpwIEDamho0MGDB5NimpqadOjQIbW1temRRx7R+vXrs45KYgIAAABggNwVk3A4\nLJfLJafTKavVqqqqKgWDwaSYbdu2ac2aNZKkRYsWqbu7W52dnRlHJTEBAAAAMEDuxCQajcrhcMSP\n7Xa7otFozpiOjo6MozL5HQAAAMAA6Y9upTq9rHBuqUsMZ3sfiQkAAACAAX6VdiZ1U0SbzaZIJBI/\njkQistvtWWM6Ojpks9kyjsqjXAAAAAAkna5wmLVjx44lxZWXl6utrU3t7e3q7e1VY2Oj/H5/Uozf\n79fmzZslSbt379aUKVNUXFyccWwqJgAAAAAGxWKxqK6uTpWVlYrFYqqurpbX61V9fb0kKRAIaMWK\nFWpqapLL5dKECRP0+OOPZ+3T6B/pveUBAAAAIAce5QIAAABQcCQmAAAAAAqOxAQAAABAwZGYAAAA\nACg4EhMAAAAABUdiAgAAAKDgSEwAAAAAFByJCQAAAICC+39M7PPmfg6vrQAAAABJRU5ErkJggg==\n",
"text": [
"<matplotlib.figure.Figure at 0x47ca3d0>"
]
}
],
"prompt_number": 14
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#Do GPy ML based regression\n",
"#Start with no SNPs, add SNPs by greedy forward selection - add best SNP at each step, is there a maxima?\n",
"#Y ~ N(0, s^2 XX^T + e^2 I )\n",
"\n",
"X=X.transpose()\n",
"Y=Y.transpose()\n",
"#print X.shape\n",
"#Y=Y.reshape((182,1))\n",
"#print Y.shape\n",
"chosenSNPS=[]\n",
"X2=np.zeros((182,0))\n",
"#print X2.shape\n",
"#print X2.dtype\n",
"k = GPy.kern.linear(input_dim=0)\n",
"m = GPy.models.GPRegression(X2,Y,k)\n",
"m.optimize(messages=False)\n",
"#m['.*variance']=1\n",
"#m['linear_variance']=10000\n",
"bestll=m.log_likelihood()\n",
"print m\n",
"print bestll\n",
"#provides initial parameters for 0 SNPS\n",
"oldparams=m['.*variance']\n",
"\n",
"#print SNPsconsider\n",
"#print np.array(set(SNPsconsider).difference(set(chosenSNPS)))\n",
"for step in range(100):\n",
" hasbeenadded=False\n",
" #bestll=np.inf\n",
" bestSNP=np.nan\n",
" X2=np.hstack([X2, np.zeros((182, 1))])\n",
" tempbestll=bestll\n",
" #for SNP in np.array(list(set(range(X.shape[1])).difference(set(chosenSNPS)))):\n",
" for SNP in np.array(list(set(SNPsconsider).difference(set(chosenSNPS)))):\n",
" X2[:,-1]=X[:,SNP]\n",
" m.X=X2\n",
" m.kern.X=X2\n",
" m.input_dim=step+1\n",
" m.kern.input_dim=step+1\n",
" #m['.*variance']=1\n",
" #m['linear_variance']=10000\n",
" #m.optimize(messages=False)\n",
" #print m\n",
" m.update_likelihood_approximation()\n",
" #print m.log_likelihood()\n",
" if m.log_likelihood() > tempbestll:\n",
" hasbeenadded=True\n",
" tempbestll=m.log_likelihood()\n",
" bestSNP=SNP\n",
" #print SNP\n",
" if hasbeenadded==False:\n",
" print \"converged: \" + str(len(chosenSNPS)) + \"SNPs\" \n",
" X2=X2[:,:-1]\n",
" break\n",
" X2[:,-1]=X[:,bestSNP]\n",
"\n",
" #print \"X2\" + str(X2.shape)\n",
" k = GPy.kern.linear(input_dim=step+1)\n",
" tempm = GPy.models.GPRegression(X2,Y,k)\n",
" tempm.optimize(messages=False)\n",
" print tempm.log_likelihood()\n",
" if tempm.log_likelihood()>bestll:\n",
" bestll=tempm.log_likelihood()\n",
" m=tempm\n",
" oldparams=m['.*variance']\n",
" chosenSNPS.append(bestSNP)\n",
" else:\n",
" print \"converged: \" + str(len(chosenSNPS)) + \"SNPs\" \n",
" X2=X2[:,:-1]\n",
" break\n",
" #print X2.shape\n",
" \n",
"pb.matshow(X2[:,:].transpose())\n",
"pb.colorbar()\n",
"pb.show()\n",
"\n",
"#m.plot()\n",
"#print m\n",
"\n",
"#print m.predict(np.array([1]))\n",
"#print m.predict(np.array([2]))\n",
"\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"Log-likelihood: -4.993e+02\n",
"\n",
" Name | Value | Constraints | Ties | prior \n",
"------------------------------------------------------------------\n",
" linear_variance | 1.0000 | (+ve) | | \n",
" noise_variance | 14.1390 | (+ve) | | \n",
"\n",
"-499.300482466\n",
"-264.808464272"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"converged: 1SNPs"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAyYAAACICAYAAAASyf8xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFoFJREFUeJzt3X1wVNX9x/HPxV0wPPx4NpLdbRfcNbsMEFcT0KFUdKwB\nZtxC7TihnYqY4g5tftZpZ6y1nalYR41jO7bEdqKDD7Q1xJlaFjvJ4oQ2qDiwiqD+CtXgGN2smjFC\n5CFCzJrfH+BudvfuQ0zC0vh+zZwZ773fnHNJ9ga+fu85x+jv7+8XAAAAABTQmELfAAAAAACQmAAA\nAAAoOBITAAAAAAVHYgIAAACg4EhMAAAAABQciQkAAACAgiMxAQAAADAoN998s4qLizV//vyMMbfe\neqvcbrfKysq0b9++nH2SmAAAAAAYlLVr1yoUCmW83tTUpEOHDqmtrU2PPPKI1q9fn7NPEhMAAAAA\ng7JkyRJNnTo14/Vt27ZpzZo1kqRFixapu7tbnZ2dWfskMQEAAAAwrKLRqBwOR/zYbrero6Mj69eQ\nmAAAAAAYdv39/UnHhmFkjScxAQAAACBJKjIMGSZt0qRJg+rHZrMpEonEjzs6OmSz2bJ+jeVL3TEA\nAACAUeekpHtMzv/q+PFB9eP3+1VXV6eqqirt3r1bU6ZMUXFxcdavITEBAAAAEPc/ecSsXr1aO3fu\nVFdXlxwOhzZs2KDPPvtMkhQIBLRixQo1NTXJ5XJpwoQJevzxx3P2afSnvvwFAAAA4CvJMAw9anJ+\nndLnjAw3KiYAAAAA4ooKNC6JCQAAAIC4fF7lGgmsygUAAAAgrsikmQmFQvJ4PHK73aqtrU27fuTI\nEa1atUplZWVatGiR/v3vf2cdl8QEAAAAQFw+iUksFlNNTY1CoZAOHDighoYGHTx4MCnm3nvv1aWX\nXqrXXntNmzdv1k9+8pOs45KYAAAAAIjLJzEJh8NyuVxyOp2yWq2qqqpSMBhMijl48KCuuuoqSVJp\naana29v10UcfZRyXxAQAAABA3CSTlioajcrhcMSP7Xa7otFoUkxZWZmeeeYZSacTmXfffVcdHR0Z\nxyUxAQAAABA33qSlMgwjZz933HGHuru75fP5VFdXJ5/Pp/POOy9jPKtyAQAAAIgrskgvfi7tGrht\nScoWJjabTZFIJH4ciURkt9uTYiZNmqTHHnssfjx79mzNmTMn47hssAgAAABA0ulKSO/k9PNjP0ne\nYLGvr0+lpaXasWOHSkpKtHDhQjU0NMjr9cZjPvnkExUVFWns2LF69NFHtWvXLj3xxBMZx6ZiAgAA\nACDOOi53jMViUV1dnSorKxWLxVRdXS2v16v6+npJUiAQ0IEDB3TTTTfJMAzNmzdPmzZtytonFRMA\nAAAAkk5XTPq/ZnL+veSKyUigYgIAAAAgIY+KyUggMQEAAACQMKEww5KYAAAAAEg4vzDDso8JAAAA\ngIRxJs1EKBSSx+OR2+1WbW1t2vWuri4tW7ZMl1xyiebNm5d1RS6Jye8AAAAAzjAMQ/0rTM43JU9+\nj8ViKi0tVUtLi2w2myoqKtKWC77rrrt06tQp3Xffferq6lJpaak6OztlsZi/tEXFBAAAAEDC+SYt\nRTgclsvlktPplNVqVVVVlYLBYFLMrFmzdPToUUnS0aNHNX369IxJicQcEwAAAAAD5bEqVzQalcPh\niB/b7Xbt2bMnKWbdunW6+uqrVVJSomPHjunpp5/O2icVEwAAAAAJE0xaCsMwcnZz77336pJLLtH7\n77+v/fv368c//rGOHTuWMZ6KCQAAAICE86XWDqk1mjnEZrMpEonEjyORiOx2e1LMSy+9pF/+8peS\npIsuukizZ8/Wm2++qfLyctM+SUwAAAAAJIyTll50un1hw8vJIeXl5Wpra1N7e7tKSkrU2NiohoaG\npBiPx6OWlhYtXrxYnZ2devPNNzVnzpyMw5KYAAAAAEjIY46JxWJRXV2dKisrFYvFVF1dLa/Xq/r6\neklSIBDQnXfeqbVr16qsrEyff/65HnjgAU2bNi1jnywXDAAAAEDSmeWC7zU5f2fycsEjgYoJAAAA\ngIQ8KiYjgcQEAAAAQEKBEhOWCwYAAACQkMdywZIUCoXk8XjkdrtVW1ubdv3BBx+Uz+eTz+fT/Pnz\nZbFY1N3dnXFY5pgAAAAAkHRmjkmDyfnVyXNMYrGYSktL1dLSIpvNpoqKCjU0NMjr9Zr2+49//EMP\nPfSQWlpaMo5NxQQAAABAwjiTliIcDsvlcsnpdMpqtaqqqkrBYDBjl0899ZRWr16ddVgSEwAAAAAJ\n55u0FNFoVA6HI35st9sVjZrvyNjT06Pt27fr+uuvzzosk98BAAAAJGSYUzKQYRh5d/fss8/qG9/4\nhqZMmZI1jsQEAAAAQMI4qXWv1Ppq5hCbzaZIJBI/jkQistvtprFbtmzJ+RqXxOR3AAAAAGcYhqH+\n/SbnL0me/N7X16fS0lLt2LFDJSUlWrhwoenk908++URz5sxRR0eHioqKso5NxQQAAABAQh6vclks\nFtXV1amyslKxWEzV1dXyer2qr6+XJAUCAUnS1q1bVVlZmTMpkaiYAAAAADjDMAz1v2dy/mvJFZOR\nQMUEAAAAQILJKlxnA4kJAAAAgLhTebzKNRKGZR+TXNvRA6Od0+nUggUL5PP5tHDhQknS4cOH9a1v\nfUsXX3yxrr32WnV3dxf4LoGRcfPNN6u4uFjz58+Pn8v2+b/vvvvkdrvl8Xj03HPPFeKWgRFj9jzc\nddddstvt8vl88vl8am5ujl/jecC56NS4sWnNTD45QGtrq3w+n+bNm6elS5dmHXfIc0wGux09MBrN\nnj1be/fu1bRp0+Lnbr/9ds2YMUO33367amtrdeTIEd1///0FvEtgZLzwwguaOHGibrzxRr3xxhuS\nMn/+Dxw4oO9973t6+eWXFY1Gdc011+itt97SmDHs94vRwex52LBhgyZNmqSf/vSnSbE8DzgXGYah\nj/onpp2faRxPmmOSTw7Q3d2txYsXa/v27bLb7erq6tKMGTMyjj3kT/5gt6MHRqvUHH/btm1as2aN\nJGnNmjXaunVrIW4LGHFLlizR1KlTk85l+vwHg0GtXr1aVqtVTqdTLpdL4XD4rN8zMFLMngfJfNIw\nzwPOVac0Lq2lyicHeOqpp3T99dfH9zfJlpRIw5CYDGY7emC0MgxD11xzjcrLy/Xoo49Kkjo7O1Vc\nXCxJKi4uVmdnZyFvETirMn3+33///aQNuPg7A18VGzduVFlZmaqrq+OvNvI84Fz1qYrSWqp8coC2\ntjYdPnxYV111lcrLy/XnP/8567hDTkwGsx09MFrt2rVL+/btU3Nzsx5++GG98MILSdcNw+BZwVdW\nrs8/zwZGu/Xr1+udd97R/v37NWvWLP3sZz/LGMvzgHNBPhWTfD6rn332mV599VU1NTVp+/bt+s1v\nfqO2traM8UNelWsw29EDo9WsWbMkSTNnztSqVasUDodVXFysDz/8UBdeeKE++OADXXDBBQW+S+Ds\nyfT5T/07o6OjQzabrVC3CZwVA3////CHP9R1110niecB565ejdUrrSf0SmtPxph8cgCHw6EZM2ao\nqKhIRUVF+uY3v6nXXntNbrfbtM8hV0zKy8vV1tam9vZ29fb2qrGxUX6/f6jdAv81enp6dOzYMUnS\niRMn9Nxzz2n+/Pny+/168sknJUlPPvmkVq5cWcjbBM6qTJ9/v9+vLVu2qLe3V++8847a2triK9kB\no9UHH3wQ/++///3v8RW7eB5wrurReM1dOlM33vX1eEuVTw7w7W9/Wy+++KJisZh6enq0Z88ezZ07\nN+O4Q66YZNqOHviq6Ozs1KpVqyRJfX19+v73v69rr71W5eXluuGGG7Rp0yY5nU49/fTTBb5TYGSs\nXr1aO3fuVFdXlxwOh+6++27dcccdpp//uXPn6oYbbtDcuXNlsVj0xz/+kVdXMKqkPg8bNmxQa2ur\n9u/fL8MwNHv2bNXX10viecC5q1fmywMPlCkH+OLzHQgE5PF4tGzZMi1YsEBjxozRunXrsiYmQ14u\nGAAAAMDoYBiGmvuXpp1fbrSari43nNj5HQAAAEBcPhWTkUBiAgAAACCuR+MLMm7Oye/5bDUPAAAA\nYHQ4pbFpzUyuPKG1tVWTJ0+Wz+eTz+fTPffck3XcrBWTWCymmpqapK3m/X4/k9sBAACAUarXZN+S\nVPnmCVdeeaW2bduW17hZKyb5bDUPAAAAYPToUVFaS5VvnjCYCfNZKyZmW83v2bMnKcYwnJLezXtA\nAAAA4Kvh6+rvby/0TQxaPhWT/PIEQy+99JLKyspks9n04IMPfvl9TPJbS/tdSb+W9ISkm/Rrbcga\nvUG/zno99esHxme7ZnZ9uOW691yy3d9w/1kG299Qfi65jPSfZajXh3o/Zn0/IemmYe57YP9fNn4w\n4w32a4fyGTkbXz+Un/twf58L+XMc6liDGXtgf0/o9DMx3L+bsvU10t/3VIX8PTrU7+twPz+D+bs7\n1XB+Jka6vy/7fXpC+T0Pw/17NtVQ/l1S6N/RI/kZHt5neWT/bTpSMs0pGSifPOHSSy9VJBLR+PHj\n1dzcrJUrV+qtt97KGJ81Mclnq/nTWiV1S2pVuyRnztsEAAAARpv2M+2/26car3db2/Vua+a3ovLJ\nEyZNmhT/7+XLl+tHP/qRDh8+rGnTppn2mTUxGbjVfElJiRobG9XQ0GASuVSnfwhL5dTObF0CAAAA\no5RTyf+L/r/z38WnNFYXLr1YFy69OH7uxQ3PJ8Xkkyd0dnbqggsukGEYCofD6u/vz5iUSHns/N7c\n3KzbbrstvtX8L37xi6TrS5cu1c6d/53fdAAAAGCkXHnllWptbS30bQyKYRj63/4H0s5vNG5Pm8hu\nlifU19dLkgKBgB5++GH96U9/ksVi0fjx4/W73/1Ol19+eeaxcyUmAAAAAL4aDMPQLf0PpZ1/xLht\nUCtsfRns/A4AAAAg7tMC7fxOYgIAAAAgLp9VuUZC1g0WAQAAAHy19GpcWjMTCoXk8XjkdrtVW1ub\nsb+XX35ZFotFzzzzTNZxqZgAAAAAiOvJ41WuWCymmpoatbS0yGazqaKiQn6/X16vNy3u5z//uZYt\nW5ZzjgoVEwAAAABxvRqb1lKFw2G5XC45nU5ZrVZVVVUpGAymxW3cuFHf/e53NXPmzJzjkpgAAAAA\niDulcWktVTQalcPhiB/b7XZFo9G0mGAwqPXr10vKvVs8r3IBAAAAiDOrkKTKlWRI0m233ab7779f\nhmGov78/56tcJCYAAAAA4no0XidaX1FP6ysZY2w2myKRSPw4EonIbrcnxezdu1dVVVWSpK6uLjU3\nN8tqtcrv95v2yQaLAAAAACSdroRc1P9/aeffNuYlVTz6+vpUWlqqHTt2qKSkRAsXLlRDQ0Pa5Pcv\nrF27Vtddd52+853vZBybigkAAACAOLM5JaksFovq6upUWVmpWCym6upqeb1e1dfXS5ICgcCgx6Vi\nAgAAAEDS6YrJ9FhH2vmPz7PnnCMyVFRMAAAAAMSdOpm7YjISSEwAAAAAxPWezL0q10hgHxMAAAAA\ncb3Hx6c1M6FQSB6PR263W7W1tWnXg8GgysrK5PP5dNlll+mf//xn1nGZYwIAAABA0pn9Sd40SQ9K\njaQ5JrFYTKWlpWppaZHNZlNFRUXaqlwnTpzQhAkTJElvvPGGVq1apUOHDmUcm4oJAAAAgISTJi1F\nOByWy+WS0+mU1WpVVVWVgsFgUswXSYkkHT9+XDNmzMg6LIkJAAAAgIRPTVqKaDQqh8MRP7bb7YpG\no2lxW7duldfr1fLly/WHP/wh67AkJgAAAAASTpi0FIZh5NXVypUrdfDgQT377LP6wQ9+kDWWVbkA\nAAAAJJyU9Hqr9EZrxhCbzaZIJBI/jkQistvtGeOXLFmivr4+ffzxx5o+fbppDIkJAAAAgISTki5e\nerp94akNSSHl5eVqa2tTe3u7SkpK1NjYqIaGhqSYt99+W3PmzJFhGHr11VclKWNSIpGYAAAAABjo\neO4Qi8Wiuro6VVZWKhaLqbq6Wl6vV/X19ZKkQCCgv/3tb9q8ebOsVqsmTpyoLVu2ZO2T5YIBAAAA\nSDozd6TeJD0IJC8XPBKomAAAAABIMFke+GwgMQEAAACQYLI88NnAcsEAAAAAEvJYLliSQqGQPB6P\n3G63amtr067/9a9/VVlZmRYsWKDFixfr9ddfzzosFRMAAAAACXlUTGKxmGpqatTS0iKbzaaKigr5\n/X55vd54zJw5c/T8889r8uTJCoVCuuWWW7R79+6MfVIxAQAAAJBwyqSlCIfDcrlccjqdslqtqqqq\nUjAYTIq54oorNHnyZEnSokWL1NHRkXVYEhMAAAAACcdNWopoNCqHwxE/ttvtikajGbvctGmTVqxY\nkXVYXuUCAAAAkJDHqlyGYeTd3b/+9S899thj2rVrV9Y4EhMAAAAACSclfdQqdbVmDLHZbIpEIvHj\nSCQiu92eFvf6669r3bp1CoVCmjp1atZh2WARAAAAgKQzlZDlJulBc/IGi319fSotLdWOHTtUUlKi\nhQsXqqGhIWny+3vvvaerr75af/nLX3T55ZfnHJuKCQAAAIAEkzklqSwWi+rq6lRZWalYLKbq6mp5\nvV7V19dLkgKBgO6++24dOXJE69evlyRZrVaFw+GMfVIxAQAAACDpTMWk3CQ9eCW5YjISqJgAAAAA\nSDBZHvhsIDEBAAAAkJDHq1wjgcQEAAAAQEIeywWPBDZYBAAAAJBw0qSZCIVC8ng8crvdqq2tTbv+\nn//8R1dccYXOP/98/fa3v805LBUTAAAAAAnHcofEYjHV1NSopaVFNptNFRUV8vv9ScsFT58+XRs3\nbtTWrVvzGpaKCQAAAICEPpOWIhwOy+Vyyel0ymq1qqqqSsFgMClm5syZKi8vl9VqzWtYEhMAAAAA\ngxKNRuVwOOLHdrtd0Wh0SH3yKhcAAACAAT7LGWEYxrCPSmICAAAAYIBjkl6UtCtjhM1mUyQSiR9H\nIhHZ7fYhjUpiAgAAAGCAHkmXnmlfeCApory8XG1tbWpvb1dJSYkaGxvV0NBg2lu+O8aTmAAAAAAY\n4NOcERaLRXV1daqsrFQsFlN1dbW8Xq/q6+slSYFAQB9++KEqKip09OhRjRkzRr///e914MABTZw4\n0bRPoz/fFAYAAADAqHZ67shekyuX5V35+LKomAAAAAAYIHfFZCSQmAAAAAAYoDCJCfuYAAAAABjg\nU5OWLhQKyePxyO12q7a21jTm1ltvldvtVllZmfbt25d1VBITAAAAAAMcNWnJYrGYampqFAqFdODA\nATU0NOjgwYNJMU1NTTp06JDa2tr0yCOPaP369VlHJTEBAAAAMEDuikk4HJbL5ZLT6ZTValVVVZWC\nwWBSzLZt27RmzRpJ0qJFi9Td3a3Ozs6Mo5KYAAAAABggd2ISjUblcDjix3a7XdFoNGdMR0dHxlGZ\n/A4AAABggPRXt1KdXlY4t9QlhrN9HYkJAAAAgAF+lXYmdVNEm82mSCQSP45EIrLb7VljOjo6ZLPZ\nMo7Kq1wAAAAAJJ2ucJi1Y8eOJcWVl5erra1N7e3t6u3tVWNjo/x+f1KM3+/X5s2bJUm7d+/WlClT\nVFxcnHFsKiYAAAAABsVisaiurk6VlZWKxWKqrq6W1+tVfX29JCkQCGjFihVqamqSy+XShAkT9Pjj\nj2ft0+gf6b3lAQAAACAHXuUCAAAAUHAkJgAAAAAKjsQEAAAAQMGRmAAAAAAoOBITAAAAAAVHYgIA\nAACg4EhMAAAAABQciQkAAACAgvt/THvt5pvqBbMAAAAASUVORK5CYII=\n",
"text": [
"<matplotlib.figure.Figure at 0x4ecf950>"
]
}
],
"prompt_number": 15
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print chosenSNPS\n",
"print m\n",
"pb.matshow(X[:,chosenSNPS].transpose())\n",
"pb.colorbar()\n",
"pb.show()\n",
"pb.matshow(X[:,1000].reshape(1,182))\n",
"pb.colorbar()\n",
"pb.show()\n",
"#print m.plot()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[2000]\n",
"\n",
"Log-likelihood: -2.683e+02\n",
"\n",
" Name | Value | Constraints | Ties | prior \n",
"------------------------------------------------------------------\n",
" linear_variance | 24.1732 | (+ve) | | \n",
" noise_variance | 1.0300 | (+ve) | | \n",
"\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAyYAAACICAYAAAASyf8xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFoFJREFUeJzt3X1wVNX9x/HPxV0wPPx4NpLdbRfcNbsMEFcT0KFUdKwB\nZtxC7TihnYqY4g5tftZpZ6y1nalYR41jO7bEdqKDD7Q1xJlaFjvJ4oQ2qDiwiqD+CtXgGN2smjFC\n5CFCzJrfH+BudvfuQ0zC0vh+zZwZ773fnHNJ9ga+fu85x+jv7+8XAAAAABTQmELfAAAAAACQmAAA\nAAAoOBITAAAAAAVHYgIAAACg4EhMAAAAABQciQkAAACAgiMxAQAAADAoN998s4qLizV//vyMMbfe\neqvcbrfKysq0b9++nH2SmAAAAAAYlLVr1yoUCmW83tTUpEOHDqmtrU2PPPKI1q9fn7NPEhMAAAAA\ng7JkyRJNnTo14/Vt27ZpzZo1kqRFixapu7tbnZ2dWfskMQEAAAAwrKLRqBwOR/zYbrero6Mj69eQ\nmAAAAAAYdv39/UnHhmFkjScxAQAAACBJKjIMGSZt0qRJg+rHZrMpEonEjzs6OmSz2bJ+jeVL3TEA\nAACAUeekpHtMzv/q+PFB9eP3+1VXV6eqqirt3r1bU6ZMUXFxcdavITEBAAAAEPc/ecSsXr1aO3fu\nVFdXlxwOhzZs2KDPPvtMkhQIBLRixQo1NTXJ5XJpwoQJevzxx3P2afSnvvwFAAAA4CvJMAw9anJ+\nndLnjAw3KiYAAAAA4ooKNC6JCQAAAIC4fF7lGgmsygUAAAAgrsikmQmFQvJ4PHK73aqtrU27fuTI\nEa1atUplZWVatGiR/v3vf2cdl8QEAAAAQFw+iUksFlNNTY1CoZAOHDighoYGHTx4MCnm3nvv1aWX\nXqrXXntNmzdv1k9+8pOs45KYAAAAAIjLJzEJh8NyuVxyOp2yWq2qqqpSMBhMijl48KCuuuoqSVJp\naana29v10UcfZRyXxAQAAABA3CSTlioajcrhcMSP7Xa7otFoUkxZWZmeeeYZSacTmXfffVcdHR0Z\nxyUxAQAAABA33qSlMgwjZz933HGHuru75fP5VFdXJ5/Pp/POOy9jPKtyAQAAAIgrskgvfi7tGrht\nScoWJjabTZFIJH4ciURkt9uTYiZNmqTHHnssfjx79mzNmTMn47hssAgAAABA0ulKSO/k9PNjP0ne\nYLGvr0+lpaXasWOHSkpKtHDhQjU0NMjr9cZjPvnkExUVFWns2LF69NFHtWvXLj3xxBMZx6ZiAgAA\nACDOOi53jMViUV1dnSorKxWLxVRdXS2v16v6+npJUiAQ0IEDB3TTTTfJMAzNmzdPmzZtytonFRMA\nAAAAkk5XTPq/ZnL+veSKyUigYgIAAAAgIY+KyUggMQEAAACQMKEww5KYAAAAAEg4vzDDso8JAAAA\ngIRxJs1EKBSSx+OR2+1WbW1t2vWuri4tW7ZMl1xyiebNm5d1RS6Jye8AAAAAzjAMQ/0rTM43JU9+\nj8ViKi0tVUtLi2w2myoqKtKWC77rrrt06tQp3Xffferq6lJpaak6OztlsZi/tEXFBAAAAEDC+SYt\nRTgclsvlktPplNVqVVVVlYLBYFLMrFmzdPToUUnS0aNHNX369IxJicQcEwAAAAAD5bEqVzQalcPh\niB/b7Xbt2bMnKWbdunW6+uqrVVJSomPHjunpp5/O2icVEwAAAAAJE0xaCsMwcnZz77336pJLLtH7\n77+v/fv368c//rGOHTuWMZ6KCQAAAICE86XWDqk1mjnEZrMpEonEjyORiOx2e1LMSy+9pF/+8peS\npIsuukizZ8/Wm2++qfLyctM+SUwAAAAAJIyTll50un1hw8vJIeXl5Wpra1N7e7tKSkrU2NiohoaG\npBiPx6OWlhYtXrxYnZ2devPNNzVnzpyMw5KYAAAAAEjIY46JxWJRXV2dKisrFYvFVF1dLa/Xq/r6\neklSIBDQnXfeqbVr16qsrEyff/65HnjgAU2bNi1jnywXDAAAAEDSmeWC7zU5f2fycsEjgYoJAAAA\ngIQ8KiYjgcQEAAAAQEKBEhOWCwYAAACQkMdywZIUCoXk8XjkdrtVW1ubdv3BBx+Uz+eTz+fT/Pnz\nZbFY1N3dnXFY5pgAAAAAkHRmjkmDyfnVyXNMYrGYSktL1dLSIpvNpoqKCjU0NMjr9Zr2+49//EMP\nPfSQWlpaMo5NxQQAAABAwjiTliIcDsvlcsnpdMpqtaqqqkrBYDBjl0899ZRWr16ddVgSEwAAAAAJ\n55u0FNFoVA6HI35st9sVjZrvyNjT06Pt27fr+uuvzzosk98BAAAAJGSYUzKQYRh5d/fss8/qG9/4\nhqZMmZI1jsQEAAAAQMI4qXWv1Ppq5hCbzaZIJBI/jkQistvtprFbtmzJ+RqXxOR3AAAAAGcYhqH+\n/SbnL0me/N7X16fS0lLt2LFDJSUlWrhwoenk908++URz5sxRR0eHioqKso5NxQQAAABAQh6vclks\nFtXV1amyslKxWEzV1dXyer2qr6+XJAUCAUnS1q1bVVlZmTMpkaiYAAAAADjDMAz1v2dy/mvJFZOR\nQMUEAAAAQILJKlxnA4kJAAAAgLhTebzKNRKGZR+TXNvRA6Od0+nUggUL5PP5tHDhQknS4cOH9a1v\nfUsXX3yxrr32WnV3dxf4LoGRcfPNN6u4uFjz58+Pn8v2+b/vvvvkdrvl8Xj03HPPFeKWgRFj9jzc\nddddstvt8vl88vl8am5ujl/jecC56NS4sWnNTD45QGtrq3w+n+bNm6elS5dmHXfIc0wGux09MBrN\nnj1be/fu1bRp0+Lnbr/9ds2YMUO33367amtrdeTIEd1///0FvEtgZLzwwguaOHGibrzxRr3xxhuS\nMn/+Dxw4oO9973t6+eWXFY1Gdc011+itt97SmDHs94vRwex52LBhgyZNmqSf/vSnSbE8DzgXGYah\nj/onpp2faRxPmmOSTw7Q3d2txYsXa/v27bLb7erq6tKMGTMyjj3kT/5gt6MHRqvUHH/btm1as2aN\nJGnNmjXaunVrIW4LGHFLlizR1KlTk85l+vwHg0GtXr1aVqtVTqdTLpdL4XD4rN8zMFLMngfJfNIw\nzwPOVac0Lq2lyicHeOqpp3T99dfH9zfJlpRIw5CYDGY7emC0MgxD11xzjcrLy/Xoo49Kkjo7O1Vc\nXCxJKi4uVmdnZyFvETirMn3+33///aQNuPg7A18VGzduVFlZmaqrq+OvNvI84Fz1qYrSWqp8coC2\ntjYdPnxYV111lcrLy/XnP/8567hDTkwGsx09MFrt2rVL+/btU3Nzsx5++GG98MILSdcNw+BZwVdW\nrs8/zwZGu/Xr1+udd97R/v37NWvWLP3sZz/LGMvzgHNBPhWTfD6rn332mV599VU1NTVp+/bt+s1v\nfqO2traM8UNelWsw29EDo9WsWbMkSTNnztSqVasUDodVXFysDz/8UBdeeKE++OADXXDBBQW+S+Ds\nyfT5T/07o6OjQzabrVC3CZwVA3////CHP9R1110niecB565ejdUrrSf0SmtPxph8cgCHw6EZM2ao\nqKhIRUVF+uY3v6nXXntNbrfbtM8hV0zKy8vV1tam9vZ29fb2qrGxUX6/f6jdAv81enp6dOzYMUnS\niRMn9Nxzz2n+/Pny+/168sknJUlPPvmkVq5cWcjbBM6qTJ9/v9+vLVu2qLe3V++8847a2triK9kB\no9UHH3wQ/++///3v8RW7eB5wrurReM1dOlM33vX1eEuVTw7w7W9/Wy+++KJisZh6enq0Z88ezZ07\nN+O4Q66YZNqOHviq6Ozs1KpVqyRJfX19+v73v69rr71W5eXluuGGG7Rp0yY5nU49/fTTBb5TYGSs\nXr1aO3fuVFdXlxwOh+6++27dcccdpp//uXPn6oYbbtDcuXNlsVj0xz/+kVdXMKqkPg8bNmxQa2ur\n9u/fL8MwNHv2bNXX10viecC5q1fmywMPlCkH+OLzHQgE5PF4tGzZMi1YsEBjxozRunXrsiYmQ14u\nGAAAAMDoYBiGmvuXpp1fbrSari43nNj5HQAAAEBcPhWTkUBiAgAAACCuR+MLMm7Oye/5bDUPAAAA\nYHQ4pbFpzUyuPKG1tVWTJ0+Wz+eTz+fTPffck3XcrBWTWCymmpqapK3m/X4/k9sBAACAUarXZN+S\nVPnmCVdeeaW2bduW17hZKyb5bDUPAAAAYPToUVFaS5VvnjCYCfNZKyZmW83v2bMnKcYwnJLezXtA\nAAAA4Kvh6+rvby/0TQxaPhWT/PIEQy+99JLKyspks9n04IMPfvl9TPJbS/tdSb+W9ISkm/Rrbcga\nvUG/zno99esHxme7ZnZ9uOW691yy3d9w/1kG299Qfi65jPSfZajXh3o/Zn0/IemmYe57YP9fNn4w\n4w32a4fyGTkbXz+Un/twf58L+XMc6liDGXtgf0/o9DMx3L+bsvU10t/3VIX8PTrU7+twPz+D+bs7\n1XB+Jka6vy/7fXpC+T0Pw/17NtVQ/l1S6N/RI/kZHt5neWT/bTpSMs0pGSifPOHSSy9VJBLR+PHj\n1dzcrJUrV+qtt97KGJ81Mclnq/nTWiV1S2pVuyRnztsEAAAARpv2M+2/26car3db2/Vua+a3ovLJ\nEyZNmhT/7+XLl+tHP/qRDh8+rGnTppn2mTUxGbjVfElJiRobG9XQ0GASuVSnfwhL5dTObF0CAAAA\no5RTyf+L/r/z38WnNFYXLr1YFy69OH7uxQ3PJ8Xkkyd0dnbqggsukGEYCofD6u/vz5iUSHns/N7c\n3KzbbrstvtX8L37xi6TrS5cu1c6d/53fdAAAAGCkXHnllWptbS30bQyKYRj63/4H0s5vNG5Pm8hu\nlifU19dLkgKBgB5++GH96U9/ksVi0fjx4/W73/1Ol19+eeaxcyUmAAAAAL4aDMPQLf0PpZ1/xLht\nUCtsfRns/A4AAAAg7tMC7fxOYgIAAAAgLp9VuUZC1g0WAQAAAHy19GpcWjMTCoXk8XjkdrtVW1ub\nsb+XX35ZFotFzzzzTNZxqZgAAAAAiOvJ41WuWCymmpoatbS0yGazqaKiQn6/X16vNy3u5z//uZYt\nW5ZzjgoVEwAAAABxvRqb1lKFw2G5XC45nU5ZrVZVVVUpGAymxW3cuFHf/e53NXPmzJzjkpgAAAAA\niDulcWktVTQalcPhiB/b7XZFo9G0mGAwqPXr10vKvVs8r3IBAAAAiDOrkKTKlWRI0m233ab7779f\nhmGov78/56tcJCYAAAAA4no0XidaX1FP6ysZY2w2myKRSPw4EonIbrcnxezdu1dVVVWSpK6uLjU3\nN8tqtcrv95v2yQaLAAAAACSdroRc1P9/aeffNuYlVTz6+vpUWlqqHTt2qKSkRAsXLlRDQ0Pa5Pcv\nrF27Vtddd52+853vZBybigkAAACAOLM5JaksFovq6upUWVmpWCym6upqeb1e1dfXS5ICgcCgx6Vi\nAgAAAEDS6YrJ9FhH2vmPz7PnnCMyVFRMAAAAAMSdOpm7YjISSEwAAAAAxPWezL0q10hgHxMAAAAA\ncb3Hx6c1M6FQSB6PR263W7W1tWnXg8GgysrK5PP5dNlll+mf//xn1nGZYwIAAABA0pn9Sd40SQ9K\njaQ5JrFYTKWlpWppaZHNZlNFRUXaqlwnTpzQhAkTJElvvPGGVq1apUOHDmUcm4oJAAAAgISTJi1F\nOByWy+WS0+mU1WpVVVWVgsFgUswXSYkkHT9+XDNmzMg6LIkJAAAAgIRPTVqKaDQqh8MRP7bb7YpG\no2lxW7duldfr1fLly/WHP/wh67AkJgAAAAASTpi0FIZh5NXVypUrdfDgQT377LP6wQ9+kDWWVbkA\nAAAAJJyU9Hqr9EZrxhCbzaZIJBI/jkQistvtGeOXLFmivr4+ffzxx5o+fbppDIkJAAAAgISTki5e\nerp94akNSSHl5eVqa2tTe3u7SkpK1NjYqIaGhqSYt99+W3PmzJFhGHr11VclKWNSIpGYAAAAABjo\neO4Qi8Wiuro6VVZWKhaLqbq6Wl6vV/X19ZKkQCCgv/3tb9q8ebOsVqsmTpyoLVu2ZO2T5YIBAAAA\nSDozd6TeJD0IJC8XPBKomAAAAABIMFke+GwgMQEAAACQYLI88NnAcsEAAAAAEvJYLliSQqGQPB6P\n3G63amtr067/9a9/VVlZmRYsWKDFixfr9ddfzzosFRMAAAAACXlUTGKxmGpqatTS0iKbzaaKigr5\n/X55vd54zJw5c/T8889r8uTJCoVCuuWWW7R79+6MfVIxAQAAAJBwyqSlCIfDcrlccjqdslqtqqqq\nUjAYTIq54oorNHnyZEnSokWL1NHRkXVYEhMAAAAACcdNWopoNCqHwxE/ttvtikajGbvctGmTVqxY\nkXVYXuUCAAAAkJDHqlyGYeTd3b/+9S899thj2rVrV9Y4EhMAAAAACSclfdQqdbVmDLHZbIpEIvHj\nSCQiu92eFvf6669r3bp1CoVCmjp1atZh2WARAAAAgKQzlZDlJulBc/IGi319fSotLdWOHTtUUlKi\nhQsXqqGhIWny+3vvvaerr75af/nLX3T55ZfnHJuKCQAAAIAEkzklqSwWi+rq6lRZWalYLKbq6mp5\nvV7V19dLkgKBgO6++24dOXJE69evlyRZrVaFw+GMfVIxAQAAACDpTMWk3CQ9eCW5YjISqJgAAAAA\nSDBZHvhsIDEBAAAAkJDHq1wjgcQEAAAAQEIeywWPBDZYBAAAAJBw0qSZCIVC8ng8crvdqq2tTbv+\nn//8R1dccYXOP/98/fa3v805LBUTAAAAAAnHcofEYjHV1NSopaVFNptNFRUV8vv9ScsFT58+XRs3\nbtTWrVvzGpaKCQAAAICEPpOWIhwOy+Vyyel0ymq1qqqqSsFgMClm5syZKi8vl9VqzWtYEhMAAAAA\ngxKNRuVwOOLHdrtd0Wh0SH3yKhcAAACAAT7LGWEYxrCPSmICAAAAYIBjkl6UtCtjhM1mUyQSiR9H\nIhHZ7fYhjUpiAgAAAGCAHkmXnmlfeCApory8XG1tbWpvb1dJSYkaGxvV0NBg2lu+O8aTmAAAAAAY\n4NOcERaLRXV1daqsrFQsFlN1dbW8Xq/q6+slSYFAQB9++KEqKip09OhRjRkzRr///e914MABTZw4\n0bRPoz/fFAYAAADAqHZ67shekyuX5V35+LKomAAAAAAYIHfFZCSQmAAAAAAYoDCJCfuYAAAAABjg\nU5OWLhQKyePxyO12q7a21jTm1ltvldvtVllZmfbt25d1VBITAAAAAAMcNWnJYrGYampqFAqFdODA\nATU0NOjgwYNJMU1NTTp06JDa2tr0yCOPaP369VlHJTEBAAAAMEDuikk4HJbL5ZLT6ZTValVVVZWC\nwWBSzLZt27RmzRpJ0qJFi9Td3a3Ozs6Mo5KYAAAAABggd2ISjUblcDjix3a7XdFoNGdMR0dHxlGZ\n/A4AAABggPRXt1KdXlY4t9QlhrN9HYkJAAAAgAF+lXYmdVNEm82mSCQSP45EIrLb7VljOjo6ZLPZ\nMo7Kq1wAAAAAJJ2ucJi1Y8eOJcWVl5erra1N7e3t6u3tVWNjo/x+f1KM3+/X5s2bJUm7d+/WlClT\nVFxcnHFsKiYAAAAABsVisaiurk6VlZWKxWKqrq6W1+tVfX29JCkQCGjFihVqamqSy+XShAkT9Pjj\nj2ft0+gf6b3lAQAAACAHXuUCAAAAUHAkJgAAAAAKjsQEAAAAQMGRmAAAAAAoOBITAAAAAAVHYgIA\nAACg4EhMAAAAABQciQkAAACAgvt/THvt5pvqBbMAAAAASUVORK5CYII=\n",
"text": [
"<matplotlib.figure.Figure at 0x6f12790>"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAyYAAACICAYAAAASyf8xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFnhJREFUeJzt3X1wVNX9x/HPxV0wPPx4NpLdbRfcNbsMEFcT0KFUdKwB\nZtxC7TihnYqY4g5tftZpZ6y1nalYR41jO1piO9HBB9oa4kwti51kcUIbVBxYRVB/hWpwjG5WzRgh\n8hAhZs3vD3A3u3v3IYawGN6vmTNy7/3u+V7DXpgv555zjP7+/n4BAAAAQAGNKvQNAAAAAACFCQAA\nAICCozABAAAAUHAUJgAAAAAKjsIEAAAAQMFRmAAAAAAoOAoTAAAAAINy8803q7i4WHPnzs0Yc+ut\nt8rtdqusrEx79uzJ2SeFCQAAAIBBWb16tUKhUMbrTU1NOnDggNra2vToo49q7dq1OfukMAEAAAAw\nKIsWLdLkyZMzXt+yZYtWrVolSVqwYIG6u7vV2dmZtU8KEwAAAACnVTQalcPhiB/b7XZ1dHRk/QyF\nCQAAAIDTrr+/P+nYMIys8RQmAAAAACRJRYYhw6RNmDBhUP3YbDZFIpH4cUdHh2w2W9bPWL7SHQMA\nAAAYcY5Lusfk/G+OHh1UP36/X3V1daqqqtLOnTs1adIkFRcXZ/0MhQkAAACAuP/JI2blypXavn27\nurq65HA4tG7dOn3++eeSpEAgoGXLlqmpqUkul0vjxo3TE088kbNPoz/15S8AAAAA5yTDMPSYyfk1\nSp8zcroxYgIAAAAgrqhAeSlMAAAAAMTl8yrXcGBVLgAAAABxRSbNTCgUksfjkdvtVm1tbdr1Q4cO\nacWKFSorK9OCBQv0n//8J2teChMAAAAAcfkUJrFYTDU1NQqFQtq3b58aGhq0f//+pJh7771Xl156\nqV5//XVt3LhRP/vZz7LmpTABAAAAEJdPYRIOh+VyueR0OmW1WlVVVaVgMJgUs3//fl111VWSpNLS\nUrW3t+vjjz/OmJfCBAAAAEDcBJOWKhqNyuFwxI/tdrui0WhSTFlZmZ599llJJwuZ9957Tx0dHRnz\nUpgAAAAAiBtr0lIZhpGznzvuuEPd3d3y+Xyqq6uTz+fTeeedlzGeVbkAAAAAxBVZpJe+kHYM3LYk\nZQsTm82mSCQSP45EIrLb7UkxEyZM0OOPPx4/njlzpmbNmpUxLxssAgAAAJB0ciSkd2L6+dGfJm+w\n2NfXp9LSUm3btk0lJSWaP3++Ghoa5PV64zGffvqpioqKNHr0aD322GPasWOHnnzyyYy5GTEBAAAA\nEGcdkzvGYrGorq5OlZWVisViqq6ultfrVX19vSQpEAho3759uummm2QYhubMmaMNGzZk7ZMREwAA\nAACSTo6Y9H/D5Pz7ySMmw4EREwAAAAAJeYyYDAcKEwAAAAAJ4wqTlsIEAAAAQML5hUnLPiYAAAAA\nEsaYNBOhUEgej0dut1u1tbVp17u6urRkyRJdcsklmjNnTtYVuSQmvwMAAAA4xTAM9S8zOd+UPPk9\nFouptLRULS0tstlsqqioSFsu+K677tKJEyd03333qaurS6Wlpers7JTFYv7SFiMmAAAAABLON2kp\nwuGwXC6XnE6nrFarqqqqFAwGk2JmzJihw4cPS5IOHz6sqVOnZixKJOaYAAAAABgoj1W5otGoHA5H\n/Nhut2vXrl1JMWvWrNHVV1+tkpISHTlyRM8880zWPhkxAQAAAJAwzqSlMAwjZzf33nuvLrnkEn3w\nwQfau3evfvrTn+rIkSMZ4xkxAQAAAJBwvtTaIbVGM4fYbDZFIpH4cSQSkd1uT4p5+eWX9etf/1qS\ndNFFF2nmzJl66623VF5ebtonhQkAAACAhDHS4otOti+teyU5pLy8XG1tbWpvb1dJSYkaGxvV0NCQ\nFOPxeNTS0qKFCxeqs7NTb731lmbNmpUxLYUJAAAAgIQ85phYLBbV1dWpsrJSsVhM1dXV8nq9qq+v\nlyQFAgHdeeedWr16tcrKyvTFF1/ogQce0JQpUzL2yXLBAAAAACSdWi74XpPzdyYvFzwcGDEBAAAA\nkJDHiMlwoDABAAAAkFCgwoTlggEAAAAk5LFcsCSFQiF5PB653W7V1tamXX/wwQfl8/nk8/k0d+5c\nWSwWdXd3Z0zLHBMAAAAAkk7NMWkwOb8yeY5JLBZTaWmpWlpaZLPZVFFRoYaGBnm9XtN+//nPf+qh\nhx5SS0tLxtyMmAAAAABIGGPSUoTDYblcLjmdTlmtVlVVVSkYDGbs8umnn9bKlSuzpqUwAQAAAJBw\nvklLEY1G5XA44sd2u13RqPmOjD09Pdq6dauuv/76rGmZ/A4AAAAgIcOckoEMw8i7u+eee07f+ta3\nNGnSpKxxFCYAAAAAEsZIrbul1tcyh9hsNkUikfhxJBKR3W43jd20aVPO17gkJr8DAAAAOMUwDPXv\nNTl/SfLk976+PpWWlmrbtm0qKSnR/PnzTSe/f/rpp5o1a5Y6OjpUVFSUNTcjJgAAAAAS8niVy2Kx\nqK6uTpWVlYrFYqqurpbX61V9fb0kKRAISJI2b96sysrKnEWJxIgJAAAAgFMMw1D/+ybnv5E8YjIc\nGDEBAAAAkGCyCteZQGECAAAAIO5EHq9yDYfTso9Jru3ogZHO6XRq3rx58vl8mj9/viTp4MGD+s53\nvqOLL75Y1157rbq7uwt8l8DwuPnmm1VcXKy5c+fGz2X7/t93331yu93yeDx6/vnnC3HLwLAxex7u\nuusu2e12+Xw++Xw+NTc3x6/xPOBsdGLM6LRmJp8aoLW1VT6fT3PmzNHixYuz5h3yHJPBbkcPjEQz\nZ87U7t27NWXKlPi522+/XdOmTdPtt9+u2tpaHTp0SPfff38B7xIYHi+++KLGjx+vG2+8UW+++aak\nzN//ffv26Qc/+IFeeeUVRaNRXXPNNXr77bc1ahT7/WJkMHse1q1bpwkTJujnP/95UizPA85GhmHo\n4/7xaeenG0eT5pjkUwN0d3dr4cKF2rp1q+x2u7q6ujRt2rSMuYf8zR/sdvTASJVa42/ZskWrVq2S\nJK1atUqbN28uxG0Bw27RokWaPHly0rlM3/9gMKiVK1fKarXK6XTK5XIpHA6f8XsGhovZ8yCZTxrm\necDZ6oTGpLVU+dQATz/9tK6//vr4/ibZihLpNBQmg9mOHhipDMPQNddco/Lycj322GOSpM7OThUX\nF0uSiouL1dnZWchbBM6oTN//Dz74IGkDLv7OwLli/fr1KisrU3V1dfzVRp4HnK0+U1FaS5VPDdDW\n1qaDBw/qqquuUnl5uf7yl79kzTvkwmQw29EDI9WOHTu0Z88eNTc365FHHtGLL76YdN0wDJ4VnLNy\nff95NjDSrV27Vu+++6727t2rGTNm6Be/+EXGWJ4HnA3yGTHJ57v6+eef67XXXlNTU5O2bt2q3/3u\nd2pra8sYP+RVuQazHT0wUs2YMUOSNH36dK1YsULhcFjFxcX66KOPdOGFF+rDDz/UBRdcUOC7BM6c\nTN//1L8zOjo6ZLPZCnWbwBkx8M//H//4x7ruuusk8Tzg7NWr0Xq19Zhebe3JGJNPDeBwODRt2jQV\nFRWpqKhI3/72t/X666/L7Xab9jnkEZPy8nK1tbWpvb1dvb29amxslN/vH2q3wNdGT0+Pjhw5Ikk6\nduyYnn/+ec2dO1d+v19PPfWUJOmpp57S8uXLC3mbwBmV6fvv9/u1adMm9fb26t1331VbW1t8JTtg\npPrwww/jv/7HP/4RX7GL5wFnqx6N1ezF03XjXd+Mt1T51ADf/e539dJLLykWi6mnp0e7du3S7Nmz\nM+Yd8ohJpu3ogXNFZ2enVqxYIUnq6+vTD3/4Q1177bUqLy/XDTfcoA0bNsjpdOqZZ54p8J0Cw2Pl\nypXavn27urq65HA4dPfdd+uOO+4w/f7Pnj1bN9xwg2bPni2LxaI//elPvLqCESX1eVi3bp1aW1u1\nd+9eGYahmTNnqr6+XhLPA85evTJfHnigTDXAl9/vQCAgj8ejJUuWaN68eRo1apTWrFmTtTAZ8nLB\nAAAAAEYGwzDU3L847fxSo9V0dbnTiZ3fAQAAAMTlM2IyHChMAAAAAMT1aGxB8uac/J7PVvMAAAAA\nRoYTGp3WzOSqE1pbWzVx4kT5fD75fD7dc889WfNmHTGJxWKqqalJ2mre7/czuR0AAAAYoXpN9i1J\nlW+dcOWVV2rLli155c06YpLPVvMAAAAARo4eFaW1VPnWCYOZMJ91xMRsq/ldu3YlxTgNQ+/lnQ4A\nAAA4V3xT/f3thb6JQctnxCSfOsEwDL388ssqKyuTzWbTgw8++NX3MclnLe33JP1W0pOSbpK0Tr/N\n+ZmBfqt1ecfm6ju1r7MtPtXAzw8m9qvcy2Dvbaj5B17P9dkzbag/y3z+356U+fMwmJ9bPtdTne58\nQ8mdy+nuL1f/Q3neBtP36YgfiuH+uX7V/E/q9DwT2WIH63T/vgzleRvuZ3uw/Q13/sEYat+n+zt/\ner6TT0q66Yz/WZTr80PJfbo/P5z9Dfefg6nZvo4yzSkZKJ864dJLL1UkEtHYsWPV3Nys5cuX6+23\n384Yn7UwyWereUlqldR96r9SuyRnzhsFAAAARpb2U+3r7TON1Xut7XqvNfN7UfnUCRMmTIj/eunS\npfrJT36igwcPasqUKaZ9Zi1MBm41X1JSosbGRjU0NKTFLdbJ34LFkrZTlAAAAOCc5FTyP9BvL8xt\nDNEJjdaFiy/WhYsvjp97ad0LSTH51AmdnZ264IILZBiGwuGw+vv7MxYlUh47vzc3N+u2226LbzX/\nq1/9Kun64sWLtX371/OHDgAAAAyXK6+8Uq2trYW+jUExDEP/2/9A2vn1xu1pE9nN6oT6+npJUiAQ\n0COPPKI///nPslgsGjt2rP7whz/o8ssvz5w7V2ECAAAA4NxgGIZu6X8o7fyjxm2DWmHrq2DndwAA\nAABxnxVo53cKEwAAAABx+azKNRyybrAIAAAA4NzSqzFpzUwoFJLH45Hb7VZtbW3G/l555RVZLBY9\n++yzWfMyYgIAAAAgriePV7lisZhqamrU0tIim82miooK+f1+eb3etLhf/vKXWrJkSc45KoyYAAAA\nAIjr1ei0liocDsvlcsnpdMpqtaqqqkrBYDAtbv369fr+97+v6dOn58xLYQIAAAAg7oTGpLVU0WhU\nDocjfmy32xWNRtNigsGg1q5dKyn3bvG8ygUAAAAgzmyEJFWuIkOSbrvtNt1///0yDEP9/f05X+Wi\nMAEAAAAQ16OxOtb6qnpaX80YY7PZFIlE4seRSER2uz0pZvfu3aqqqpIkdXV1qbm5WVarVX6/37RP\nNlgEAAAAIOnkSMhF/f+Xdv4dY07SiEdfX59KS0u1bds2lZSUaP78+WpoaEib/P6l1atX67rrrtP3\nvve9jLkZMQEAAAAQZzanJJXFYlFdXZ0qKysVi8VUXV0tr9er+vp6SVIgEBh0XkZMAAAAAEg6OWIy\nNdaRdv6T8+w554gMFSMmAAAAAOJOHM89YjIcKEwAAAAAxPUez70q13BgHxMAAAAAcb1Hx6Y1M6FQ\nSB6PR263W7W1tWnXg8GgysrK5PP5dNlll+lf//pX1rzMMQEAAAAg6dT+JG+ZlAelRtIck1gsptLS\nUrW0tMhms6mioiJtVa5jx45p3LhxkqQ333xTK1as0IEDBzLmZsQEAAAAQMJxk5YiHA7L5XLJ6XTK\narWqqqpKwWAwKebLokSSjh49qmnTpmVNS2ECAAAAIOEzk5YiGo3K4XDEj+12u6LRaFrc5s2b5fV6\ntXTpUv3xj3/MmpbCBAAAAEDCMZOWwjCMvLpavny59u/fr+eee04/+tGPssayKhcAAACAhOOS3miV\n3mzNGGKz2RSJROLHkUhEdrs9Y/yiRYvU19enTz75RFOnTjWNoTABAAAAkHBc0sWLT7YvPb0uKaS8\nvFxtbW1qb29XSUmJGhsb1dDQkBTzzjvvaNasWTIMQ6+99pokZSxKJAoTAAAAAAMdzR1isVhUV1en\nyspKxWIxVVdXy+v1qr6+XpIUCAT097//XRs3bpTVatX48eO1adOmrH2yXDAAAAAASafmjtSblAeB\n5OWChwMjJgAAAAASTJYHPhMoTAAAAAAkmCwPfCawXDAAAACAhDyWC5akUCgkj8cjt9ut2tratOt/\n+9vfVFZWpnnz5mnhwoV64403sqZlxAQAAABAQh4jJrFYTDU1NWppaZHNZlNFRYX8fr+8Xm88Ztas\nWXrhhRc0ceJEhUIh3XLLLdq5c2fGPhkxAQAAAJBwwqSlCIfDcrlccjqdslqtqqqqUjAYTIq54oor\nNHHiREnSggUL1NHRkTUthQkAAACAhKMmLUU0GpXD4Ygf2+12RaPRjF1u2LBBy5Yty5qWV7kAAAAA\nJOSxKpdhGHl39+9//1uPP/64duzYkTWOwgQAAABAwnFJH7dKXa0ZQ2w2myKRSPw4EonIbrenxb3x\nxhtas2aNQqGQJk+enDUtGywCAAAAkHRqJGSpSXnQnLzBYl9fn0pLS7Vt2zaVlJRo/vz5amhoSJr8\n/v777+vqq6/WX//6V11++eU5czNiAgAAACDBZE5JKovForq6OlVWVioWi6m6ulper1f19fWSpEAg\noLvvvluHDh3S2rVrJUlWq1XhcDhjn4yYAAAAAJB0asSk3KQ8eDV5xGQ4MGICAAAAIMFkeeAzgcIE\nAAAAQEIer3INBwoTAAAAAAl5LBc8HNhgEQAAAEDCcZNmIhQKyePxyO12q7a2Nu36f//7X11xxRU6\n//zz9fvf/z5nWkZMAAAAACQcyR0Si8VUU1OjlpYW2Ww2VVRUyO/3Jy0XPHXqVK1fv16bN2/OKy0j\nJgAAAAAS+kxainA4LJfLJafTKavVqqqqKgWDwaSY6dOnq7y8XFarNa+0FCYAAAAABiUajcrhcMSP\n7Xa7otHokPrkVS4AAAAAA3yeM8IwjNOelcIEAAAAwABHJL0kaUfGCJvNpkgkEj+ORCKy2+1Dykph\nAgAAAGCAHkmXnmpfeiApory8XG1tbWpvb1dJSYkaGxvV0NBg2lu+O8ZTmAAAAAAY4LOcERaLRXV1\ndaqsrFQsFlN1dbW8Xq/q6+slSYFAQB999JEqKip0+PBhjRo1Sg8//LD27dun8ePHm/Zp9OdbwgAA\nAAAY0U7OHdltcuWyvEc+vipGTAAAAAAMkHvEZDhQmAAAAAAYoDCFCfuYAAAAABjgM5OWLhQKyePx\nyO12q7a21jTm1ltvldvtVllZmfbs2ZM1K4UJAAAAgAEOm7RksVhMNTU1CoVC2rdvnxoaGrR///6k\nmKamJh04cEBtbW169NFHtXbt2qxZKUwAAAAADJB7xCQcDsvlcsnpdMpqtaqqqkrBYDApZsuWLVq1\napUkacGCBeru7lZnZ2fGrBQmAAAAAAbIXZhEo1E5HI74sd1uVzQazRnT0dGRMSuT3wEAAAAMkP7q\nVqqTywrnlrrEcLbPUZgAAAAAGOA3aWdSN0W02WyKRCLx40gkIrvdnjWmo6NDNpstY1Ze5QIAAAAg\n6eQIh1k7cuRIUlx5ebna2trU3t6u3t5eNTY2yu/3J8X4/X5t3LhRkrRz505NmjRJxcXFGXMzYgIA\nAABgUCwWi+rq6lRZWalYLKbq6mp5vV7V19dLkgKBgJYtW6ampia5XC6NGzdOTzzxRNY+jf7h3lse\nAAAAAHLgVS4AAAAABUdhAgAAAKDgKEwAAAAAFByFCQAAAICCozABAAAAUHAUJgAAAAAKjsIEAAAA\nQMFRmAAAAAAouP8HE7/o5pQM9Z4AAAAASUVORK5CYII=\n",
"text": [
"<matplotlib.figure.Figure at 0x3f5a710>"
]
}
],
"prompt_number": 16
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"## Compute the biased estimator of HSIC.\n",
"# @param x The data.\n",
"# @param y The labels.\n",
"# @param kernelx The kernel on the data, default to linear kernel.\n",
"# @param kernely The kernel on the labels, default to linear kernel.\n",
"#\n",
"\n",
"def BiasedHSIC(x, y, kernelx, kernely):\n",
" #nx = kernelx.input_dim\n",
" #ny = kernely.input_dim\n",
" #print x.shape\n",
" nx=float(x.shape[0])\n",
" ny=float(y.shape[0])\n",
" \n",
" assert nx == ny, \\\n",
" \"Argument 1 and 2 have different number of data points\"\n",
"\n",
" kMat=kernelx.K(x)\n",
" #if len(nx) > 1:\n",
" # kMat = kernelx.Dot(x, x)\n",
" #else:\n",
" # kMat = numpy.outerproduct(x, x)\n",
"\n",
" #print kMat\n",
" hlhMat = ComputeHLH(y, kernely)\n",
" #print hlhMat\n",
" return numpy.sum(numpy.sum(kMat * hlhMat)) / ((nx-1)*(nx-1))\n",
" \n",
"## Compute HLH give the labels.\n",
"# @param y The labels.\n",
"# @param kernely The kernel on the labels, default to linear kernel.\n",
"#\n",
"def ComputeHLH(y, kernely):\n",
" #ny = kernely.input_dim\n",
" ny=float(y.shape[0])\n",
" #if len(ny) > 1:\n",
" # lMat = kernely.Dot(y, y)\n",
" #else:\n",
" # lMat = numpy.outerproduct(y, y)\n",
" #print y.shape\n",
" lMat=kernely.K(y)\n",
" #print lMat\n",
" sL = numpy.sum(lMat, axis=1)\n",
" ssL = numpy.sum(sL)\n",
" # hlhMat\n",
" return lMat - numpy.add.outer(sL, sL)/ny + ssL/(ny*ny)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 17
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#Do GPy ML based regression\n",
"#Start with no SNPs, add SNPs by greedy forward selection - add best SNP at each step, is there a maxima?\n",
"#Y ~ N(0, s^2 XX^T + e^2 I )\n",
"\n",
"#X=X.transpose()\n",
"#Y=Y.transpose()\n",
"#print X.shape\n",
"#print Y.shape\n",
"chosenSNPS2=[]\n",
"X2=np.zeros((182,0))\n",
"#print X2.shape\n",
"#print X2.dtype\n",
"\n",
"Y=Y.transpose()\n",
"#X=X.transpose()\n",
"print X2.shape\n",
"print Y.shape\n",
"k2 = GPy.kern.linear(input_dim=0)\n",
"m2 = GPy.models.GPRegression(X2,Y,k2)\n",
"m2.optimize(messages=False)\n",
"#m['.*variance']=1\n",
"#m['linear_variance']=10000\n",
"bestll=m2.log_likelihood()\n",
"X=X.transpose()\n",
"oldHSIC=0\n",
"#print m\n",
"#print bestll\n",
"#provides initial parameters for 0 SNPS\n",
"oldparams=m['.*variance']\n",
"kernelx2 = GPy.kern.linear(input_dim=ng)\n",
"kernely2 = GPy.kern.linear(input_dim=X.shape[1])\n",
"\n",
"#print X.shape\n",
"for step in range(500):\n",
" hasbeenadded=False\n",
" #bestll=np.inf\n",
" bestSNP=np.nan\n",
" X2=np.hstack([X2, np.zeros((182, 1))])\n",
" tempbestHSIC=oldHSIC\n",
"\n",
" #for SNP in np.array(list(set(range(X.shape[1])).difference(set(chosenSNPS2)))):\n",
" for SNP in np.array(list(set(SNPsconsider).difference(set(chosenSNPS2)))):\n",
"\n",
" X2[:,-1]=X[:,SNP]\n",
" m2.X=X2\n",
" m2.kern.X=X2\n",
" m2.input_dim=step+1\n",
" m2.kern.input_dim=step+1\n",
"\n",
" m2.update_likelihood_approximation()\n",
" \n",
" m2._Xoffset = np.zeros((1, m2.input_dim))\n",
" m2._Xscale = np.ones((1, m2.input_dim))\n",
"\n",
" residuals=m2.predict(X2)[0]\n",
"\n",
" newHSIC=BiasedHSIC(residuals, X, kernelx2, kernely2)\n",
" if newHSIC > tempbestHSIC:\n",
" hasbeenadded=True\n",
" #oldHSIC=newHSIC\n",
" tempbestHSIC=newHSIC\n",
" bestSNP=SNP\n",
" if hasbeenadded==False:\n",
" print \"converged: \" + str(len(chosenSNPS2)) + \"SNPs\" \n",
" X2=X2[:,:-1]\n",
" break\n",
" X2[:,-1]=X[:,bestSNP]\n",
"\n",
" #print \"X2\" + str(X2.shape)\n",
"\n",
" k2 = GPy.kern.linear(input_dim=step+1)\n",
"\n",
" tempm2 = GPy.models.GPRegression(X2,Y,k2)\n",
" tempm2.optimize(messages=False)\n",
" residuals=tempm2.predict(X2)[0]\n",
"\n",
" newHSIC=BiasedHSIC(residuals, X, kernelx2, kernely2)\n",
" if newHSIC > oldHSIC:\n",
" #hasbeenadded=True\n",
" oldHSIC=newHSIC\n",
" m2=tempm2\n",
" print newHSIC\n",
" print m2.log_likelihood()\n",
" oldparams=m2['.*variance']\n",
" chosenSNPS2.append(bestSNP)\n",
" else:\n",
" print \"converged: \" + str(len(chosenSNPS2)) + \"SNPs\" \n",
" X2=X2[:,:-1]\n",
" break\n",
" #print X2.shape\n",
" \n",
"pb.matshow(X2[:,:].transpose())\n",
"pb.colorbar()\n",
"pb.show()\n",
"\n",
"#m.plot()\n",
"#print m\n",
"\n",
"#print m.predict(np.array([1]))\n",
"#print m.predict(np.array([2]))\n",
"\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"(182, 0)\n",
"(182, 1)\n",
"564.625529454"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-264.808464272\n",
"595.669299835"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-266.927369149\n",
"607.310018956"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-269.103409966\n",
"613.738771928"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-271.13949487\n",
"619.256789863"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-272.697020004\n",
"620.435920743"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-275.345073878\n",
"621.009692534"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-277.880692572\n",
"621.114219258"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-280.202631535\n",
"621.918054578"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-282.444740607\n",
"622.457302782"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-284.722386428\n",
"622.555324361"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-286.494271702\n",
"converged: 11SNPs"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAywAAACICAYAAAAF62/4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH+VJREFUeJzt3X9wVNXdx/HPpVnllwUUCGQ3fVZMTELFGE1AWq1YwRBm\nuoI/Q22lGDWDZNRWZ9T2man4dBSstjLEdqKjFqqGOLYlYGF1Qg0qiKkFkSFUgnXtZsVMU0TUqDHr\nPn9Qtuxmf9+7Pwjv18wZuLtnzznZ3bubb+4552sEAoGAAAAAACAHDcv2AAAAAAAgGgIWAAAAADmL\ngAUAAABAziJgAQAAAJCzCFgAAAAA5CwCFgAAAAA5i4AFAAAAgCWuv/565efna9q0aVHr3HLLLSou\nLlZ5ebl27twZt00CFgAAAACWWLx4sdxud9T7N27cqP3796urq0uPPvqolixZErdNAhYAAAAAlrjw\nwgs1bty4qPevX79eixYtkiTNmDFDhw4dUk9PT8w2CVgAAAAAZITP51NhYWHw2OFwqLu7O+ZjCFgA\nAAAAZEwgEAg5NgwjZn0CFgAAAAAxjTAMGRHKKaecklQ7drtdXq83eNzd3S273R7zMXkpjRgAAADA\nCeNzSb+IcPv/fvJJUu24XC41NjaqtrZW27dv19ixY5Wfnx/zMQQsAAAAAOL6egJ1Fi5cqC1btqi3\nt1eFhYVatmyZvvzyS0lSfX295s2bp40bN6qoqEijRo3Sk08+GbdNIxA+iQwAAAAAjmEYhh6LcPuN\nGrwmxWpcYQEAAAAQ14gs9UvAAgAAACCuRKaEpQO7hAEAAACIa0SEEonb7VZpaamKi4u1YsWKQfd/\n+OGHWrBggcrLyzVjxgzt2bMnZr8ELAAAAADiSiRg8fv9amhokNvtVmdnp5qbm7V3796QOvfdd5/O\nPfdc7dq1S2vWrNGtt94as18CFgAAAABxJRKwdHR0qKioSE6nUzabTbW1tWptbQ2ps3fvXl188cWS\npJKSEnk8Hv3rX/+K2i8BCwAAAIC4TolQwvl8PhUWFgaPHQ6HfD5fSJ3y8nL98Y9/lHQkwHnvvffU\n3d0dtV8CFgAAAABxjYxQwhmGEbedu+66S4cOHVJFRYUaGxtVUVGhr33ta1Hrs0sYAAAAgLhG5Emv\nfiVtPTbtSlgKFrvdLq/XGzz2er1yOBwhdU455RQ98cQTwePTTz9dU6ZMidoviSMBAAAAxGQYhvrH\nDL79pI9CE0cODAyopKREmzdvVkFBgaZPn67m5maVlZUF63z00UcaMWKETjrpJD322GPaunWrfve7\n30XtmyssAAAAAOKynRy/Tl5enhobG1VdXS2/36+6ujqVlZWpqalJklRfX6/Ozk796Ec/kmEYOuus\ns/T444/HbJMrLAAAAABiMgxDgW9EuP2foVdY0oErLAAAAADiS+AKSzoQsAAAAACIb1R2uiVgAQAA\nABDf8Ox0Sx4WAAAAAPGdHKFE4Ha7VVpaquLiYq1YsWLQ/b29vZo7d67OOeccnXXWWTF3CJNYdA8A\nAAAgDsMwFJgX4faNoYvu/X6/SkpK1NbWJrvdrqqqqkHbGt9zzz364osvdP/996u3t1clJSXq6elR\nXl7kyV9pvcISL7oChjqn06mzzz5bFRUVmj59uiTp4MGDmjNnjs4880xdeumlOnToUJZHCaTH9ddf\nr/z8fE2bNi14W6z3//3336/i4mKVlpbqxRdfzMaQgbSJdD7cc889cjgcqqioUEVFhTZt2hS8j/MB\nOWl4hBKmo6NDRUVFcjqdstlsqq2tVWtra0idyZMn6/Dhw5Kkw4cP67TTTosarEhpDFj8fr8aGhrk\ndrvV2dmp5uZm7d27N13dATnJMAy1t7dr586d6ujokCQtX75cc+bM0b59+3TJJZdo+fLlWR4lkB6L\nFy+W2+0OuS3a+7+zs1MtLS3q7OyU2+3WzTffrK+++iobwwbSItL5YBiGfvKTn2jnzp3auXOnampq\nJHE+IIclMCXM5/OpsLAweOxwOOTz+ULq3HjjjdqzZ48KCgpUXl6ulStXxuw2bQFLItEVcCIIn3W5\nfv16LVq0SJK0aNEirVu3LhvDAtLuwgsv1Lhx40Jui/b+b21t1cKFC2Wz2eR0OlVUVBQM8oGhINL5\nIEXOX8H5gJw1KkIJYxhG3Gbuu+8+nXPOOXr//ff15ptvaunSpfr444+j1k9bwJJIdAUMdYZhaPbs\n2aqsrNRjjz0mSerp6VF+fr4kKT8/Xz09PdkcIpBR0d7/77//vhwOR7Ae3xk4UaxatUrl5eWqq6sL\nTpHkfEDOGi6190r37PpvCWe32+X1eoPHXq835P0sSdu2bdNVV10lSTrjjDN0+umn6+23347abdoC\nlkSiK2Co27p1q3bu3KlNmzbpkUce0SuvvBJyv2EYnCs4YcV7/3NuYKhbsmSJ3n33Xb355puaPHmy\nbr/99qh1OR+QE06WZp0h3fOd/5ZwlZWV6urqksfjUX9/v1paWuRyuULqlJaWqq2tTdKRP2S9/fbb\nmjJlStRu0xawJBJdAUPd5MmTJUkTJkzQggUL1NHRofz8fH3wwQeSpAMHDmjixInZHCKQUdHe/+Hf\nGd3d3bLb7VkZI5ApEydODAbuN9xwQ3DaF+cDclYCa1jy8vLU2Nio6upqTZ06Vddcc43KysrU1NSk\npqYmSdJPf/pTvfHGGyovL9fs2bP1wAMP6NRTT43abdoClkSiK2Ao6+vrC87H/PTTT/Xiiy9q2rRp\ncrlcWr16tSRp9erVmj9/fjaHCWRUtPe/y+XS2rVr1d/fr3fffVddXV3BnfWAoerAgQPB///pT38K\n7iDG+YCcNTpCiaCmpkZvv/229u/fr7vvvluSVF9fr/r6eknS+PHjtWHDBu3atUu7d+/W97///Zjd\npi3T/bHRld/vV11dXcj+y8BQ19PTowULFkiSBgYGdO211+rSSy9VZWWlrr76aj3++ONyOp169tln\nszxSID0WLlyoLVu2qLe3V4WFhbr33nt11113RXz/T506VVdffbWmTp2qvLw8/eY3v2EKDIaU8PNh\n2bJlam9v15tvvinDMHT66acH//rM+YCcFSVRZLqROBIAAABATIZhKNAY4faGyLvdWSmtiSMBAAAA\nDBEJbGssxU8e/+CDDwYTpk6bNk15eXkxE2lzhQUAAABATIZhKNAc4faFoVdY/H6/SkpK1NbWJrvd\nrqqqKjU3N0ddGvL888/r4YcfDu4aFknKV1jiRU4AAAAAhpAEdglLNnn8M888o4ULF8bsNqWAxe/3\nq6GhQW63W52dnWpubtbevXtTaQoAAADA8WB4hBImmeTxfX19euGFF3TFFVfE7DalgCXZyAkAAADA\ncS6BNSzJ7Gi3YcMGXXDBBRo7dmzMeiltaxwpcnr99ddD6hiGU9J7qTQPAAAADGH/o0DAk+1BJO9k\nqf1vUvuO6FWSSR6/du3auNPBpBQDlsQip/ck/VxSu6RZ+rmWxay9TD+PeX/444+tH+u+SPdbLd7Y\n44k1Pqt/lmTbM/O6xJPun8Xs/WbHE6ntdkmzLG772PZTrZ9Mf8k+1sx7JBOPN/O6W/08Z/N1NNtX\nMn0fba9d/z0frP5sitVWup/3cNn8HDX7vFp9/iTz3R3OyvdEuttL9Xlq15Fzwurvt2SZ+b0k25/R\n6XwPW3sup/d307QZLs369pFy1LLHQ6scmzy+oKBALS0tam4evFr/o48+0ssvv6xnnnkmbrcpBSyJ\nR07tkjyS2uWR5EylMwAAAOC45vlPOc5F2cb4WNGSxx9NjHo02/26detUXV2tESNGxG0zpW2NBwYG\nVFJSos2bN6ugoEDTp08ftF2ZYRjHXF/J7b8UpvuvpPFk8opQrv1l0UzbyfZl9V+g4rUfua92WXHF\n0cxYIrWfzHMz1K6YxPqLb6y6icj2Xw7NSPf5eES7jl5jsfKzKN3fAen+S34yf6nP9meFlaz+/jHb\nX7KsGF+7UrsKb1Y2Pzuy/Tok8x2QrNjndvqTLVrNMAwF/hnh9m+k/2dJ6QpLtMgpEqeZ0QFDjjPb\nAwByiDPbAwByijPbAwDiibArWCakFLB4vV498MADysvLk81m06hR0a8POVMdGTAkObM9ACCHOLM9\nACCnOLM9ACCOLxKYEpYOKW1rbLPZ9Otf/1p79uzR9u3b9cgjj5CHBQAAABjCvjj5pEElkkQSzLe3\nt6uiokJnnXWWZs2aFbPflK6wTJo0SZMmTZIkjR49WmVlZXr//fejTgtLRLx5srk8X9ls/WR2QAtn\ndn5zptffxJqbbVayY0/nritW71Ji9vFW7jhjtu907wSU6XnwsdpO9j2ezvpWv07JyvTucsn0FS7d\n692S+VnT/bmW7NiSrZ9Me7m8/tOsTO9gls51V2bPh3Svlc2l30vChbafu+/XWPq/FilA6Q85Oppg\nvq2tTXa7XVVVVXK5XCFxwqFDh7R06VK98MILcjgc6u3tjdlvSldYjuXxeLRz507NmDHDbFMAAAAA\nctQXOnlQCZdIgvlnnnlGV1xxRXCX4fHjx8fs11TA8sknn+jKK6/UypUrNXr0aDNNAQAAAMhhn2nE\noBIuUoJ5n88XUqerq0sHDx7UxRdfrMrKSv3+97+P2W9KU8Ik6csvv9QVV1yhH/zgB5o/f37EOu0h\nRx6xnAwAAAAnHo+GQh6WSFdUwiWSYP7LL7/Ujh07tHnzZvX19WnmzJk6//zzVVxcHLF+SgFLIBBQ\nXV2dpk6dqttuuy1qvS0x5gLGm19pZt6r1fMrrc7NYGa9jtm2rZ4jb3bNjFV1I43FbPvpnB9tds1K\nOLPvqWR+tmznYUlWMv2ZfR7TvX7Gys8Gq/vK9Hs01jqPTK8lyqR0r6eJd3821yJlOleXmfbNfnZY\nfb6m87varFzKBZQsc5+TW6wdTIb06yS90f6p3mjvi1onkQTzhYWFGj9+vEaMGKERI0boO9/5jnbt\n2hU1YElpStjWrVv11FNP6aWXXlJFRYUqKirkdrtTaQoAAADAcaBPIzV11gRdd8//BEu4yspKdXV1\nyePxqL+/Xy0tLXK5XCF1LrvsMr366qvy+/3q6+vT66+/rqlTp0btN6UrLBdccIG++uor+f1+VVZW\nyuFwaO7cuak0BQAAAOA40K/I2xgfK1qC+aamJklSfX29SktLNXfuXJ199tkaNmyYbrzxRusDlqNW\nrlypqVOn6uOPPzbTDAAAAIAcl8gaFkmqqalRTU1NyG319fUhx3fccYfuuOOOhNpLOWDp7u7Wxo0b\n9bOf/Uy/+tWvItbJ1jxDK+frR2J2bmcy/ad7Tn02mX2dzO69Hq8/M/PezeZCSHZsZnMBJSPTOWbS\nLdYe/JnOv2F2bZ+ZvCxWr0HJ5NqHdK4LTKT9dD7e7Lls9bqpTL4H0/3Zkc5cJeEyff6YlUwelmTH\nbvX3VTrzzw3l9XCpSuQKSzqkHLD8+Mc/1i9/+UsdPnzYyvEAAAAAyEF9GpmVflNadP/8889r4sSJ\nqqioUCAQsHpMAAAAAHLMFzppUInE7XartLRUxcXFWrFixaD729vbNWbMmODmXb/4xS9i9pvSFZZt\n27Zp/fr12rhxoz7//HMdPnxY1113ndasWRM6mGP+7xRZWAAAAHAi8mgo5GHpT2ANi9/vV0NDg9ra\n2mS321VVVSWXy6WysrKQehdddJHWr1+fUL9GwOQlki1btujBBx/Uhg0bQhs2jJCZe1bPuY9V3+r8\nAVbvEZ5Mf1avRTA79zSd8zGtnqNr9RqqZF5nq9c5Wf26WDmHON1zrTO55iyTc9ojtZ/NOfXJtp3p\nzworvwOynZMmmfatzoVlZiyR2jfz3Fid1yiT60sj9ZfMd3eyfWf6OyCT78l0r9FM57oqM8/7Mum4\nm6VkGIYeDtw06PbbjEdDfpbXXntNy5YtC6Y8Wb58uSTprrvuCtZpb2/XQw89NCh+iCalKWHhEslo\nCQAAAOD41a+TB5VwPp9PhYWFwWOHwyGfzxdSxzAMbdu2TeXl5Zo3b546Oztj9pvyovtDhw7phhtu\n0J49e2QYhrZv367zzz8/1eYAAAAA5LBoa1aOlciFjHPPPVder1cjR47Upk2bNH/+fO3bty9q/ZQD\nlltvvVXz5s3Tc889p4GBAX366aepNgUAAAAgx32mkXqv3aP32t+LWsdut8vr9QaPvV6vHA5HSJ1T\nTjkl+P+amhrdfPPNOnjwoE499dSIbaYUsHz00Ud65ZVXtHr16iON5OVpzJgxqTQFAAAA4DjwhU7S\npFlnatKsM4O3vbrs5ZA6lZWV6urqksfjUUFBgVpaWtTc3BxSp6enRxMnTpRhGOro6FAgEIgarEgp\nBizvvvuuJkyYoMWLF2vXrl0677zztHLlSo0cGbo3s5WLLtO5IMzqhX/xmFk8mmxiu3QnTDO34Mxc\nkj6zSfbCWbn5QqYTPaZ7sbeVm1zEe7zZZIm5nFzV7ILlcGY+S6xOtpbNhbjJvmfM/uzxxhrvfjOf\nVVYnpc1mssV0b36QzY0iMplcNJH+zX4/munL7O8xZvtP5vGZ3dzA2o2BMiWRXcLy8vLU2Nio6upq\n+f1+1dXVqaysTE1NTZKOZLx/7rnn9Nvf/lZ5eXkaOXKk1q5dG7vNVAY7MDCgHTt2qLGxUVVVVbrt\nttu0fPly3Xvvvak0BwAAACDHJbKGRToyzaumpibktvr6+uD/ly5dqqVLlybcb0oBi8PhkMPhUFVV\nlSTpyiuvDG5ZFqr9mP87RSYWAAAAnGg8kkJ/Lz4+fZalTPcpBSyTJk1SYWGh9u3bpzPPPFNtbW36\n5je/GaHmLHOjAwAAAI5zTkmhvxdvycYwTEv0CovVUt4lbNWqVbr22mvV39+vM844Q08++eSgOmbm\nGIezcs6j2TUuZuc7x+vPTJKjZGUzUWQ4q+fsxmP2Z0vnvPN0JwEzO+c+FqvXD8RrP9n70zlXO93v\nKTP9WT1HPtPz1JNhddLYeDK9DiSXZPJnSfdaomSlM3mw2bV78Vi9jtJM21Y/3oxMris8XiWyhkWS\n3G63brvtNvn9ft1www268847I9b761//qpkzZ+rZZ5/V5ZdfHrW9lAOWjRs3qq+vT8OGDdPw4cM1\nfPjwVJsCAAAAkOP6EpgS5vf71dDQoLa2NtntdlVVVcnlcqmsrGxQvTvvvFNz585VIBCI2WZKme49\nHo8ee+wx7dixQ7t375bf74+7uh8AAADA8atfJw0q4To6OlRUVCSn0ymbzaba2lq1trYOqrdq1Spd\neeWVmjBhQtx+UwpYvv71r8tms6mvr08DAwPq6+uT3W5PpSkAAAAAx4EvdPKgEs7n86mwsDB47HA4\n5PP5BtVpbW3VkiVLJEmGYcTsN6UpYaeeeqpuv/12feMb39CIESNUXV2t2bNnJ9WG1fMxzbSd7twm\nuTS3O1w296lPltVz8OMxk/vE7Botq9fjpPN1tTrfTTir8+eYea4znStoKM1/zvS6kliynfchHivX\nMWb6syid393J1k/3ey6T/WX6dTTD6pxM8R6PzIp0RSVcvOBDUjAlimEYCgQCcaeEpRSwvPPOO3r4\n4Yfl8Xg0ZswYXXXVVXr66ad17bXXptIcAAAAgBzXp5H6tP0N9bW/EbWO3W6X1+sNHnu9XjkcjpA6\nf/vb31RbWytJ6u3t1aZNm2Sz2eRyuSK2mVLA8sYbb+hb3/qWTjvtNEnS5Zdfrm3btg0KWNpDjjwi\nDwsAAABOPB4dzcZyPOvXSbLN+pbGzPpW8LbeZY+G1KmsrFRXV5c8Ho8KCgrU0tKi5ubmkDr/+Mc/\ngv9fvHixvve970UNVqQUA5bS0lL93//9nz777DMNHz5cbW1tmj59+qB6s475/xaCFQAAAJyQnAr9\nw/3xmocl/rbGeXl5amxsVHV1tfx+v+rq6lRWVqampiZJoRnvE5VSwFJeXq7rrrtOlZWVGjZsmM49\n91zddNNNg+plal681XuzJ9u+1f2Zke75x+nMdWJ1/g6r9+Q3s94gnblBEnl8OtdNWX3+mGXluqxM\n5/aJJ51rJ7L9OmYzf0e4XF5LZDY/R7Zfx2PHk+nvynQ/N8c+Pp15ThJpL9n66cyzYvZns3odY6zH\np3sdVS59lqTqs69GJFSvpqZGNTU1IbdFC1Qi5XIMFzNguf766/XnP/9ZEydO1O7duyVJBw8e1DXX\nXKP33ntPTqdTzz77rMaOHZvQ4AEAAAAcn774PLHEkVaLua3x4sWL5Xa7Q25bvny55syZo3379umS\nSy7R8uXL0zpAAAAAANnX//lJg0omxAxYLrzwQo0bNy7ktvXr12vRokWSpEWLFmndunXpGx0AAACA\nnND/ychBJRK3263S0lIVFxdrxYoVg+5vbW1VeXm5KioqdN555+kvf/lLzH6TXsPS09Oj/Px8SVJ+\nfr56enqi1v25lqVtb7Bk5hxmeh6s1fWTke38AGbmwSbLbI4Ls/2lxiPJafl8ZbPngJXPldU5LuK1\nn+750bEemyyz52cm8+mEM5tPKvL9HlnxDWH1vPNk61ud2yvW91u41J73xPpO5PFm5NpndC7kB/HI\nmt+ZMrlGM1lWrjlJ5PHp/J3K7Pl3XPo8fujg9/vV0NCgtrY22e12VVVVyeVyqaysLFhn9uzZuuyy\nyyRJu3fv1oIFC7R///6obaaU6f4owzDiJofxmOkAGHI82R4AkEM82R4AkFM82R4AEM/nEUqYjo4O\nFRUVyel0ymazqba2Vq2trSF1Ro0aFfz/J598ovHjx8fsNukrLPn5+frggw80adIkHThwQBMnToxa\nt11HTr52Dd7MDQAAADgxeDQkQtLP4lfx+XwqLCwMHjscDr3++uuD6q1bt0533323Dhw4oBdffDFm\nm0lfYXG5XFq9erUkafXq1Zo/f37UurN0JEg5+i8AAABw4nHqyG/ER8tx6tMIJUy82VdHzZ8/X3v3\n7tWGDRv0wx/+MHblQAy1tbWByZMnB2w2W8DhcASeeOKJwL///e/AJZdcEiguLg7MmTMn8OGHH0Z8\n7EUXXRSQRKFQKBQKhUKhUI4pF110UaxfwXOSpID+HAjo/pcC+v7P/1sUGk689tprgerq6uDxfffd\nF1i+fHnMtqdMmRLo7e2Ner/xnwEAAAAAQESGYUh/iBA2XGHo2HBiYGBAJSUl2rx5swoKCjR9+nQ1\nNzeHLLp/5513NGXKFBmGoR07duiqq67SO++8E7XvlDLdAwAAADjBfBK/Sl5enhobG1VdXS2/36+6\nujqVlZWpqalJ0pGM93/4wx+0Zs0a2Ww2jR49WmvXro3ZJldYAAAAAMRkGIbUFCFsqA+9wpIOXGEB\nAAAAEF+EbYwzgYAFAAAAQHwJbGucDqYSRwIAAAA4QSSwrbEkud1ulZaWqri4WCtWrBh0/9NPP63y\n8nKdffbZ+va3v6233norZrdcYQEAAAAQXwJXWPx+vxoaGtTW1ia73a6qqiq5XK6QXcKmTJmil19+\nWWPGjJHb7dZNN92k7du3R22TKywAAAAA4vsiQgnT0dGhoqIiOZ1O2Ww21dbWqrW1NaTOzJkzNWbM\nGEnSjBkz1N3dHbNbAhYAAAAA8X0SoYTx+XwqLCwMHjscDvl8vqhNPv7445o3b17MbpkSBgAAACC+\nBHYJMwwj4eZeeuklPfHEE9q6dWvMegQsAAAAAOL7XNK/2qXe9qhV7Ha7vF5v8Njr9crhcAyq99Zb\nb+nGG2+U2+3WuHHjYnZL4kgAAAAAMRmGIdVECBs2hSaOHBgYUElJiTZv3qyCggJNnz5dzc3NIYvu\n//nPf+q73/2unnrqKZ1//vlx++YKCwAAAID4IqxZCZeXl6fGxkZVV1fL7/errq5OZWVlampqkiTV\n19fr3nvv1YcffqglS5ZIkmw2mzo6OqK2yRUWAAAAADEZhiFVRggb3gi9wpIOXGEBAAAAEF+EbYwz\ngYAFAAAAQHwJTAlLBwIWAAAAAPElsK1xOpA4EgAAAEB8n0coEbjdbpWWlqq4uFgrVqwYdP/f//53\nzZw5U8OHD9dDDz0Ut1uusAAAAACI7+P4Vfx+vxoaGtTW1ia73a6qqiq5XK6QbY1PO+00rVq1SuvW\nrUuoW66wAAAAAIhvIEIJ09HRoaKiIjmdTtlsNtXW1qq1tTWkzoQJE1RZWSmbzZZQtwQsAAAAACzh\n8/lUWFgYPHY4HPL5fKbaZEoYAAAAgAR8GbeGYRiW90rAAgAAACABH0t6VdLWqDXsdru8Xm/w2Ov1\nyuFwmOqVgAUAAABAAvoknfufctQDITUqKyvV1dUlj8ejgoICtbS0qLm5OWJrgUAgoV4JWAAAAAAk\n4LO4NfLy8tTY2Kjq6mr5/X7V1dWprKxMTU1NkqT6+np98MEHqqqq0uHDhzVs2DCtXLlSnZ2dGj16\ndMQ2jUCioQ0AAACAE9KRtSl/i3DPeQlfKUkVV1gAAAAAJCD+FZZ0IGABAAAAkIDsBCzkYQEAAACQ\ngM8ilMHcbrdKS0tVXFysFStWRKxzyy23qLi4WOXl5dq5c2fMXglYAAAAACTgcIQSyu/3q6GhQW63\nW52dnWpubtbevXtD6mzcuFH79+9XV1eXHn30US1ZsiRmrwQsAAAAABIQ/wpLR0eHioqK5HQ6ZbPZ\nVFtbq9bW1pA669ev16JFiyRJM2bM0KFDh9TT0xO1VwIWAAAAAAmIH7D4fD4VFhYGjx0Oh3w+X9w6\n3d3dUXtl0T0AAACABAyeAhbuyPbH8YVvhRzrcQQsAAAAABLwv4NuCU/2aLfb5fV6g8der1cOhyNm\nne7ubtnt9qi9MiUMAAAAQEyBQCBi+fjjj0PqVVZWqqurSx6PR/39/WppaZHL5Qqp43K5tGbNGknS\n9u3bNXbsWOXn50ftmyssAAAAACyRl5enxsZGVVdXy+/3q66uTmVlZWpqapIk1dfXa968edq4caOK\nioo0atQoPfnkkzHbNALhE8gAAAAAIEcwJQwAAABAziJgAQAAAJCzCFgAAAAA5CwCFgAAAAA5i4AF\nAAAAQM4iYAEAAACQswhYAAAAAOQsAhYAAAAAOev/Ad5Cl3cMAzVZAAAAAElFTkSuQmCC\n",
"text": [
"<matplotlib.figure.Figure at 0x8256c50>"
]
}
],
"prompt_number": 21
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print chosenSNPS2\n",
"print m2\n",
"pb.matshow(X[:,chosenSNPS2].transpose())\n",
"pb.colorbar()\n",
"pb.show()\n",
"pb.matshow(X[:,1000].reshape(1,182))\n",
"pb.colorbar()\n",
"pb.show()\n",
"#print m.plot()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[2000, 13000, 1000, 3000, 1800, 11400, 7000, 10800, 11600, 6600, 800]\n",
"\n",
"Log-likelihood: -2.888e+02\n",
"\n",
" Name | Value | Constraints | Ties | prior \n",
"-----------------------------------------------------------------\n",
" linear_variance | 2.3442 | (+ve) | | \n",
" noise_variance | 1.0246 | (+ve) | | \n",
"\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAywAAACICAYAAAAF62/4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH+VJREFUeJzt3X9wVNXdx/HPpVnllwUUCGQ3fVZMTELFGE1AWq1YwRBm\nuoI/Q22lGDWDZNRWZ9T2man4dBSstjLEdqKjFqqGOLYlYGF1Qg0qiKkFkSFUgnXtZsVMU0TUqDHr\nPn9Qtuxmf9+7Pwjv18wZuLtnzznZ3bubb+4552sEAoGAAAAAACAHDcv2AAAAAAAgGgIWAAAAADmL\ngAUAAABAziJgAQAAAJCzCFgAAAAA5CwCFgAAAAA5i4AFAAAAgCWuv/565efna9q0aVHr3HLLLSou\nLlZ5ebl27twZt00CFgAAAACWWLx4sdxud9T7N27cqP3796urq0uPPvqolixZErdNAhYAAAAAlrjw\nwgs1bty4qPevX79eixYtkiTNmDFDhw4dUk9PT8w2CVgAAAAAZITP51NhYWHw2OFwqLu7O+ZjCFgA\nAAAAZEwgEAg5NgwjZn0CFgAAAAAxjTAMGRHKKaecklQ7drtdXq83eNzd3S273R7zMXkpjRgAAADA\nCeNzSb+IcPv/fvJJUu24XC41NjaqtrZW27dv19ixY5Wfnx/zMQQsAAAAAOL6egJ1Fi5cqC1btqi3\nt1eFhYVatmyZvvzyS0lSfX295s2bp40bN6qoqEijRo3Sk08+GbdNIxA+iQwAAAAAjmEYhh6LcPuN\nGrwmxWpcYQEAAAAQ14gs9UvAAgAAACCuRKaEpQO7hAEAAACIa0SEEonb7VZpaamKi4u1YsWKQfd/\n+OGHWrBggcrLyzVjxgzt2bMnZr8ELAAAAADiSiRg8fv9amhokNvtVmdnp5qbm7V3796QOvfdd5/O\nPfdc7dq1S2vWrNGtt94as18CFgAAAABxJRKwdHR0qKioSE6nUzabTbW1tWptbQ2ps3fvXl188cWS\npJKSEnk8Hv3rX/+K2i8BCwAAAIC4TolQwvl8PhUWFgaPHQ6HfD5fSJ3y8nL98Y9/lHQkwHnvvffU\n3d0dtV8CFgAAAABxjYxQwhmGEbedu+66S4cOHVJFRYUaGxtVUVGhr33ta1Hrs0sYAAAAgLhG5Emv\nfiVtPTbtSlgKFrvdLq/XGzz2er1yOBwhdU455RQ98cQTwePTTz9dU6ZMidoviSMBAAAAxGQYhvrH\nDL79pI9CE0cODAyopKREmzdvVkFBgaZPn67m5maVlZUF63z00UcaMWKETjrpJD322GPaunWrfve7\n30XtmyssAAAAAOKynRy/Tl5enhobG1VdXS2/36+6ujqVlZWpqalJklRfX6/Ozk796Ec/kmEYOuus\ns/T444/HbJMrLAAAAABiMgxDgW9EuP2foVdY0oErLAAAAADiS+AKSzoQsAAAAACIb1R2uiVgAQAA\nABDf8Ox0Sx4WAAAAAPGdHKFE4Ha7VVpaquLiYq1YsWLQ/b29vZo7d67OOeccnXXWWTF3CJNYdA8A\nAAAgDsMwFJgX4faNoYvu/X6/SkpK1NbWJrvdrqqqqkHbGt9zzz364osvdP/996u3t1clJSXq6elR\nXl7kyV9pvcISL7oChjqn06mzzz5bFRUVmj59uiTp4MGDmjNnjs4880xdeumlOnToUJZHCaTH9ddf\nr/z8fE2bNi14W6z3//3336/i4mKVlpbqxRdfzMaQgbSJdD7cc889cjgcqqioUEVFhTZt2hS8j/MB\nOWl4hBKmo6NDRUVFcjqdstlsqq2tVWtra0idyZMn6/Dhw5Kkw4cP67TTTosarEhpDFj8fr8aGhrk\ndrvV2dmp5uZm7d27N13dATnJMAy1t7dr586d6ujokCQtX75cc+bM0b59+3TJJZdo+fLlWR4lkB6L\nFy+W2+0OuS3a+7+zs1MtLS3q7OyU2+3WzTffrK+++iobwwbSItL5YBiGfvKTn2jnzp3auXOnampq\nJHE+IIclMCXM5/OpsLAweOxwOOTz+ULq3HjjjdqzZ48KCgpUXl6ulStXxuw2bQFLItEVcCIIn3W5\nfv16LVq0SJK0aNEirVu3LhvDAtLuwgsv1Lhx40Jui/b+b21t1cKFC2Wz2eR0OlVUVBQM8oGhINL5\nIEXOX8H5gJw1KkIJYxhG3Gbuu+8+nXPOOXr//ff15ptvaunSpfr444+j1k9bwJJIdAUMdYZhaPbs\n2aqsrNRjjz0mSerp6VF+fr4kKT8/Xz09PdkcIpBR0d7/77//vhwOR7Ae3xk4UaxatUrl5eWqq6sL\nTpHkfEDOGi6190r37PpvCWe32+X1eoPHXq835P0sSdu2bdNVV10lSTrjjDN0+umn6+23347abdoC\nlkSiK2Co27p1q3bu3KlNmzbpkUce0SuvvBJyv2EYnCs4YcV7/3NuYKhbsmSJ3n33Xb355puaPHmy\nbr/99qh1OR+QE06WZp0h3fOd/5ZwlZWV6urqksfjUX9/v1paWuRyuULqlJaWqq2tTdKRP2S9/fbb\nmjJlStRu0xawJBJdAUPd5MmTJUkTJkzQggUL1NHRofz8fH3wwQeSpAMHDmjixInZHCKQUdHe/+Hf\nGd3d3bLb7VkZI5ApEydODAbuN9xwQ3DaF+cDclYCa1jy8vLU2Nio6upqTZ06Vddcc43KysrU1NSk\npqYmSdJPf/pTvfHGGyovL9fs2bP1wAMP6NRTT43abdoClkSiK2Ao6+vrC87H/PTTT/Xiiy9q2rRp\ncrlcWr16tSRp9erVmj9/fjaHCWRUtPe/y+XS2rVr1d/fr3fffVddXV3BnfWAoerAgQPB///pT38K\n7iDG+YCcNTpCiaCmpkZvv/229u/fr7vvvluSVF9fr/r6eknS+PHjtWHDBu3atUu7d+/W97///Zjd\npi3T/bHRld/vV11dXcj+y8BQ19PTowULFkiSBgYGdO211+rSSy9VZWWlrr76aj3++ONyOp169tln\nszxSID0WLlyoLVu2qLe3V4WFhbr33nt11113RXz/T506VVdffbWmTp2qvLw8/eY3v2EKDIaU8PNh\n2bJlam9v15tvvinDMHT66acH//rM+YCcFSVRZLqROBIAAABATIZhKNAY4faGyLvdWSmtiSMBAAAA\nDBEJbGssxU8e/+CDDwYTpk6bNk15eXkxE2lzhQUAAABATIZhKNAc4faFoVdY/H6/SkpK1NbWJrvd\nrqqqKjU3N0ddGvL888/r4YcfDu4aFknKV1jiRU4AAAAAhpAEdglLNnn8M888o4ULF8bsNqWAxe/3\nq6GhQW63W52dnWpubtbevXtTaQoAAADA8WB4hBImmeTxfX19euGFF3TFFVfE7DalgCXZyAkAAADA\ncS6BNSzJ7Gi3YcMGXXDBBRo7dmzMeiltaxwpcnr99ddD6hiGU9J7qTQPAAAADGH/o0DAk+1BJO9k\nqf1vUvuO6FWSSR6/du3auNPBpBQDlsQip/ck/VxSu6RZ+rmWxay9TD+PeX/444+tH+u+SPdbLd7Y\n44k1Pqt/lmTbM/O6xJPun8Xs/WbHE6ntdkmzLG772PZTrZ9Mf8k+1sx7JBOPN/O6W/08Z/N1NNtX\nMn0fba9d/z0frP5sitVWup/3cNn8HDX7vFp9/iTz3R3OyvdEuttL9Xlq15Fzwurvt2SZ+b0k25/R\n6XwPW3sup/d307QZLs369pFy1LLHQ6scmzy+oKBALS0tam4evFr/o48+0ssvv6xnnnkmbrcpBSyJ\nR07tkjyS2uWR5EylMwAAAOC45vlPOc5F2cb4WNGSxx9NjHo02/26detUXV2tESNGxG0zpW2NBwYG\nVFJSos2bN6ugoEDTp08ftF2ZYRjHXF/J7b8UpvuvpPFk8opQrv1l0UzbyfZl9V+g4rUfua92WXHF\n0cxYIrWfzHMz1K6YxPqLb6y6icj2Xw7NSPf5eES7jl5jsfKzKN3fAen+S34yf6nP9meFlaz+/jHb\nX7KsGF+7UrsKb1Y2Pzuy/Tok8x2QrNjndvqTLVrNMAwF/hnh9m+k/2dJ6QpLtMgpEqeZ0QFDjjPb\nAwByiDPbAwByijPbAwDiibArWCakFLB4vV498MADysvLk81m06hR0a8POVMdGTAkObM9ACCHOLM9\nACCnOLM9ACCOLxKYEpYOKW1rbLPZ9Otf/1p79uzR9u3b9cgjj5CHBQAAABjCvjj5pEElkkQSzLe3\nt6uiokJnnXWWZs2aFbPflK6wTJo0SZMmTZIkjR49WmVlZXr//fejTgtLRLx5srk8X9ls/WR2QAtn\ndn5zptffxJqbbVayY0/nritW71Ji9vFW7jhjtu907wSU6XnwsdpO9j2ezvpWv07JyvTucsn0FS7d\n692S+VnT/bmW7NiSrZ9Me7m8/tOsTO9gls51V2bPh3Svlc2l30vChbafu+/XWPq/FilA6Q85Oppg\nvq2tTXa7XVVVVXK5XCFxwqFDh7R06VK98MILcjgc6u3tjdlvSldYjuXxeLRz507NmDHDbFMAAAAA\nctQXOnlQCZdIgvlnnnlGV1xxRXCX4fHjx8fs11TA8sknn+jKK6/UypUrNXr0aDNNAQAAAMhhn2nE\noBIuUoJ5n88XUqerq0sHDx7UxRdfrMrKSv3+97+P2W9KU8Ik6csvv9QVV1yhH/zgB5o/f37EOu0h\nRx6xnAwAAAAnHo+GQh6WSFdUwiWSYP7LL7/Ujh07tHnzZvX19WnmzJk6//zzVVxcHLF+SgFLIBBQ\nXV2dpk6dqttuuy1qvS0x5gLGm19pZt6r1fMrrc7NYGa9jtm2rZ4jb3bNjFV1I43FbPvpnB9tds1K\nOLPvqWR+tmznYUlWMv2ZfR7TvX7Gys8Gq/vK9Hs01jqPTK8lyqR0r6eJd3821yJlOleXmfbNfnZY\nfb6m87varFzKBZQsc5+TW6wdTIb06yS90f6p3mjvi1onkQTzhYWFGj9+vEaMGKERI0boO9/5jnbt\n2hU1YElpStjWrVv11FNP6aWXXlJFRYUqKirkdrtTaQoAAADAcaBPIzV11gRdd8//BEu4yspKdXV1\nyePxqL+/Xy0tLXK5XCF1LrvsMr366qvy+/3q6+vT66+/rqlTp0btN6UrLBdccIG++uor+f1+VVZW\nyuFwaO7cuak0BQAAAOA40K/I2xgfK1qC+aamJklSfX29SktLNXfuXJ199tkaNmyYbrzxRusDlqNW\nrlypqVOn6uOPPzbTDAAAAIAcl8gaFkmqqalRTU1NyG319fUhx3fccYfuuOOOhNpLOWDp7u7Wxo0b\n9bOf/Uy/+tWvItbJ1jxDK+frR2J2bmcy/ad7Tn02mX2dzO69Hq8/M/PezeZCSHZsZnMBJSPTOWbS\nLdYe/JnOv2F2bZ+ZvCxWr0HJ5NqHdK4LTKT9dD7e7Lls9bqpTL4H0/3Zkc5cJeEyff6YlUwelmTH\nbvX3VTrzzw3l9XCpSuQKSzqkHLD8+Mc/1i9/+UsdPnzYyvEAAAAAyEF9GpmVflNadP/8889r4sSJ\nqqioUCAQsHpMAAAAAHLMFzppUInE7XartLRUxcXFWrFixaD729vbNWbMmODmXb/4xS9i9pvSFZZt\n27Zp/fr12rhxoz7//HMdPnxY1113ndasWRM6mGP+7xRZWAAAAHAi8mgo5GHpT2ANi9/vV0NDg9ra\n2mS321VVVSWXy6WysrKQehdddJHWr1+fUL9GwOQlki1btujBBx/Uhg0bQhs2jJCZe1bPuY9V3+r8\nAVbvEZ5Mf1avRTA79zSd8zGtnqNr9RqqZF5nq9c5Wf26WDmHON1zrTO55iyTc9ojtZ/NOfXJtp3p\nzworvwOynZMmmfatzoVlZiyR2jfz3Fid1yiT60sj9ZfMd3eyfWf6OyCT78l0r9FM57oqM8/7Mum4\nm6VkGIYeDtw06PbbjEdDfpbXXntNy5YtC6Y8Wb58uSTprrvuCtZpb2/XQw89NCh+iCalKWHhEslo\nCQAAAOD41a+TB5VwPp9PhYWFwWOHwyGfzxdSxzAMbdu2TeXl5Zo3b546Oztj9pvyovtDhw7phhtu\n0J49e2QYhrZv367zzz8/1eYAAAAA5LBoa1aOlciFjHPPPVder1cjR47Upk2bNH/+fO3bty9q/ZQD\nlltvvVXz5s3Tc889p4GBAX366aepNgUAAAAgx32mkXqv3aP32t+LWsdut8vr9QaPvV6vHA5HSJ1T\nTjkl+P+amhrdfPPNOnjwoE499dSIbaYUsHz00Ud65ZVXtHr16iON5OVpzJgxqTQFAAAA4DjwhU7S\npFlnatKsM4O3vbrs5ZA6lZWV6urqksfjUUFBgVpaWtTc3BxSp6enRxMnTpRhGOro6FAgEIgarEgp\nBizvvvuuJkyYoMWLF2vXrl0677zztHLlSo0cGbo3s5WLLtO5IMzqhX/xmFk8mmxiu3QnTDO34Mxc\nkj6zSfbCWbn5QqYTPaZ7sbeVm1zEe7zZZIm5nFzV7ILlcGY+S6xOtpbNhbjJvmfM/uzxxhrvfjOf\nVVYnpc1mssV0b36QzY0iMplcNJH+zX4/munL7O8xZvtP5vGZ3dzA2o2BMiWRXcLy8vLU2Nio6upq\n+f1+1dXVqaysTE1NTZKOZLx/7rnn9Nvf/lZ5eXkaOXKk1q5dG7vNVAY7MDCgHTt2qLGxUVVVVbrt\nttu0fPly3Xvvvak0BwAAACDHJbKGRToyzaumpibktvr6+uD/ly5dqqVLlybcb0oBi8PhkMPhUFVV\nlSTpyiuvDG5ZFqr9mP87RSYWAAAAnGg8kkJ/Lz4+fZalTPcpBSyTJk1SYWGh9u3bpzPPPFNtbW36\n5je/GaHmLHOjAwAAAI5zTkmhvxdvycYwTEv0CovVUt4lbNWqVbr22mvV39+vM844Q08++eSgOmbm\nGIezcs6j2TUuZuc7x+vPTJKjZGUzUWQ4q+fsxmP2Z0vnvPN0JwEzO+c+FqvXD8RrP9n70zlXO93v\nKTP9WT1HPtPz1JNhddLYeDK9DiSXZPJnSfdaomSlM3mw2bV78Vi9jtJM21Y/3oxMris8XiWyhkWS\n3G63brvtNvn9ft1www268847I9b761//qpkzZ+rZZ5/V5ZdfHrW9lAOWjRs3qq+vT8OGDdPw4cM1\nfPjwVJsCAAAAkOP6EpgS5vf71dDQoLa2NtntdlVVVcnlcqmsrGxQvTvvvFNz585VIBCI2WZKme49\nHo8ee+wx7dixQ7t375bf74+7uh8AAADA8atfJw0q4To6OlRUVCSn0ymbzaba2lq1trYOqrdq1Spd\neeWVmjBhQtx+UwpYvv71r8tms6mvr08DAwPq6+uT3W5PpSkAAAAAx4EvdPKgEs7n86mwsDB47HA4\n5PP5BtVpbW3VkiVLJEmGYcTsN6UpYaeeeqpuv/12feMb39CIESNUXV2t2bNnJ9WG1fMxzbSd7twm\nuTS3O1w296lPltVz8OMxk/vE7Botq9fjpPN1tTrfTTir8+eYea4znStoKM1/zvS6kliynfchHivX\nMWb6syid393J1k/3ey6T/WX6dTTD6pxM8R6PzIp0RSVcvOBDUjAlimEYCgQCcaeEpRSwvPPOO3r4\n4Yfl8Xg0ZswYXXXVVXr66ad17bXXptIcAAAAgBzXp5H6tP0N9bW/EbWO3W6X1+sNHnu9XjkcjpA6\nf/vb31RbWytJ6u3t1aZNm2Sz2eRyuSK2mVLA8sYbb+hb3/qWTjvtNEnS5Zdfrm3btg0KWNpDjjwi\nDwsAAABOPB4dzcZyPOvXSbLN+pbGzPpW8LbeZY+G1KmsrFRXV5c8Ho8KCgrU0tKi5ubmkDr/+Mc/\ngv9fvHixvve970UNVqQUA5bS0lL93//9nz777DMNHz5cbW1tmj59+qB6s475/xaCFQAAAJyQnAr9\nw/3xmocl/rbGeXl5amxsVHV1tfx+v+rq6lRWVqampiZJoRnvE5VSwFJeXq7rrrtOlZWVGjZsmM49\n91zddNNNg+plal681XuzJ9u+1f2Zke75x+nMdWJ1/g6r9+Q3s94gnblBEnl8OtdNWX3+mGXluqxM\n5/aJJ51rJ7L9OmYzf0e4XF5LZDY/R7Zfx2PHk+nvynQ/N8c+Pp15ThJpL9n66cyzYvZns3odY6zH\np3sdVS59lqTqs69GJFSvpqZGNTU1IbdFC1Qi5XIMFzNguf766/XnP/9ZEydO1O7duyVJBw8e1DXX\nXKP33ntPTqdTzz77rMaOHZvQ4AEAAAAcn774PLHEkVaLua3x4sWL5Xa7Q25bvny55syZo3379umS\nSy7R8uXL0zpAAAAAANnX//lJg0omxAxYLrzwQo0bNy7ktvXr12vRokWSpEWLFmndunXpGx0AAACA\nnND/ychBJRK3263S0lIVFxdrxYoVg+5vbW1VeXm5KioqdN555+kvf/lLzH6TXsPS09Oj/Px8SVJ+\nfr56enqi1v25lqVtb7Bk5hxmeh6s1fWTke38AGbmwSbLbI4Ls/2lxiPJafl8ZbPngJXPldU5LuK1\nn+750bEemyyz52cm8+mEM5tPKvL9HlnxDWH1vPNk61ud2yvW91u41J73xPpO5PFm5NpndC7kB/HI\nmt+ZMrlGM1lWrjlJ5PHp/J3K7Pl3XPo8fujg9/vV0NCgtrY22e12VVVVyeVyqaysLFhn9uzZuuyy\nyyRJu3fv1oIFC7R///6obaaU6f4owzDiJofxmOkAGHI82R4AkEM82R4AkFM82R4AEM/nEUqYjo4O\nFRUVyel0ymazqba2Vq2trSF1Ro0aFfz/J598ovHjx8fsNukrLPn5+frggw80adIkHThwQBMnToxa\nt11HTr52Dd7MDQAAADgxeDQkQtLP4lfx+XwqLCwMHjscDr3++uuD6q1bt0533323Dhw4oBdffDFm\nm0lfYXG5XFq9erUkafXq1Zo/f37UurN0JEg5+i8AAABw4nHqyG/ER8tx6tMIJUy82VdHzZ8/X3v3\n7tWGDRv0wx/+MHblQAy1tbWByZMnB2w2W8DhcASeeOKJwL///e/AJZdcEiguLg7MmTMn8OGHH0Z8\n7EUXXRSQRKFQKBQKhUKhUI4pF110UaxfwXOSpID+HAjo/pcC+v7P/1sUGk689tprgerq6uDxfffd\nF1i+fHnMtqdMmRLo7e2Ner/xnwEAAAAAQESGYUh/iBA2XGHo2HBiYGBAJSUl2rx5swoKCjR9+nQ1\nNzeHLLp/5513NGXKFBmGoR07duiqq67SO++8E7XvlDLdAwAAADjBfBK/Sl5enhobG1VdXS2/36+6\nujqVlZWpqalJ0pGM93/4wx+0Zs0a2Ww2jR49WmvXro3ZJldYAAAAAMRkGIbUFCFsqA+9wpIOXGEB\nAAAAEF+EbYwzgYAFAAAAQHwJbGucDqYSRwIAAAA4QSSwrbEkud1ulZaWqri4WCtWrBh0/9NPP63y\n8nKdffbZ+va3v6233norZrdcYQEAAAAQXwJXWPx+vxoaGtTW1ia73a6qqiq5XK6QXcKmTJmil19+\nWWPGjJHb7dZNN92k7du3R22TKywAAAAA4vsiQgnT0dGhoqIiOZ1O2Ww21dbWqrW1NaTOzJkzNWbM\nGEnSjBkz1N3dHbNbAhYAAAAA8X0SoYTx+XwqLCwMHjscDvl8vqhNPv7445o3b17MbpkSBgAAACC+\nBHYJMwwj4eZeeuklPfHEE9q6dWvMegQsAAAAAOL7XNK/2qXe9qhV7Ha7vF5v8Njr9crhcAyq99Zb\nb+nGG2+U2+3WuHHjYnZL4kgAAAAAMRmGIdVECBs2hSaOHBgYUElJiTZv3qyCggJNnz5dzc3NIYvu\n//nPf+q73/2unnrqKZ1//vlx++YKCwAAAID4IqxZCZeXl6fGxkZVV1fL7/errq5OZWVlampqkiTV\n19fr3nvv1YcffqglS5ZIkmw2mzo6OqK2yRUWAAAAADEZhiFVRggb3gi9wpIOXGEBAAAAEF+EbYwz\ngYAFAAAAQHwJTAlLBwIWAAAAAPElsK1xOpA4EgAAAEB8n0coEbjdbpWWlqq4uFgrVqwYdP/f//53\nzZw5U8OHD9dDDz0Ut1uusAAAAACI7+P4Vfx+vxoaGtTW1ia73a6qqiq5XK6QbY1PO+00rVq1SuvW\nrUuoW66wAAAAAIhvIEIJ09HRoaKiIjmdTtlsNtXW1qq1tTWkzoQJE1RZWSmbzZZQtwQsAAAAACzh\n8/lUWFgYPHY4HPL5fKbaZEoYAAAAgAR8GbeGYRiW90rAAgAAACABH0t6VdLWqDXsdru8Xm/w2Ov1\nyuFwmOqVgAUAAABAAvoknfufctQDITUqKyvV1dUlj8ejgoICtbS0qLm5OWJrgUAgoV4JWAAAAAAk\n4LO4NfLy8tTY2Kjq6mr5/X7V1dWprKxMTU1NkqT6+np98MEHqqqq0uHDhzVs2DCtXLlSnZ2dGj16\ndMQ2jUCioQ0AAACAE9KRtSl/i3DPeQlfKUkVV1gAAAAAJCD+FZZ0IGABAAAAkIDsBCzkYQEAAACQ\ngM8ilMHcbrdKS0tVXFysFStWRKxzyy23qLi4WOXl5dq5c2fMXglYAAAAACTgcIQSyu/3q6GhQW63\nW52dnWpubtbevXtD6mzcuFH79+9XV1eXHn30US1ZsiRmrwQsAAAAABIQ/wpLR0eHioqK5HQ6ZbPZ\nVFtbq9bW1pA669ev16JFiyRJM2bM0KFDh9TT0xO1VwIWAAAAAAmIH7D4fD4VFhYGjx0Oh3w+X9w6\n3d3dUXtl0T0AAACABAyeAhbuyPbH8YVvhRzrcQQsAAAAABLwv4NuCU/2aLfb5fV6g8der1cOhyNm\nne7ubtnt9qi9MiUMAAAAQEyBQCBi+fjjj0PqVVZWqqurSx6PR/39/WppaZHL5Qqp43K5tGbNGknS\n9u3bNXbsWOXn50ftmyssAAAAACyRl5enxsZGVVdXy+/3q66uTmVlZWpqapIk1dfXa968edq4caOK\nioo0atQoPfnkkzHbNALhE8gAAAAAIEcwJQwAAABAziJgAQAAAJCzCFgAAAAA5CwCFgAAAAA5i4AF\nAAAAQM4iYAEAAACQswhYAAAAAOQsAhYAAAAAOev/Ad5Cl3cMAzVZAAAAAElFTkSuQmCC\n",
"text": [
"<matplotlib.figure.Figure at 0x39b8a10>"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAyYAAACICAYAAAASyf8xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFnhJREFUeJzt3X1wVNX9x/HPxV0wPPx4NpLdbRfcNbsMEFcT0KFUdKwB\nZtxC7TihnYqY4g5tftZpZ6y1nalYR41jO1piO9HBB9oa4kwti51kcUIbVBxYRVB/hWpwjG5WzRgh\n8hAhZs3vD3A3u3v3IYawGN6vmTNy7/3u+V7DXpgv555zjP7+/n4BAAAAQAGNKvQNAAAAAACFCQAA\nAICCozABAAAAUHAUJgAAAAAKjsIEAAAAQMFRmAAAAAAoOAoTAAAAAINy8803q7i4WHPnzs0Yc+ut\nt8rtdqusrEx79uzJ2SeFCQAAAIBBWb16tUKhUMbrTU1NOnDggNra2vToo49q7dq1OfukMAEAAAAw\nKIsWLdLkyZMzXt+yZYtWrVolSVqwYIG6u7vV2dmZtU8KEwAAAACnVTQalcPhiB/b7XZ1dHRk/QyF\nCQAAAIDTrr+/P+nYMIys8RQmAAAAACRJRYYhw6RNmDBhUP3YbDZFIpH4cUdHh2w2W9bPWL7SHQMA\nAAAYcY5Lusfk/G+OHh1UP36/X3V1daqqqtLOnTs1adIkFRcXZ/0MhQkAAACAuP/JI2blypXavn27\nurq65HA4tG7dOn3++eeSpEAgoGXLlqmpqUkul0vjxo3TE088kbNPoz/15S8AAAAA5yTDMPSYyfk1\nSp8zcroxYgIAAAAgrqhAeSlMAAAAAMTl8yrXcGBVLgAAAABxRSbNTCgUksfjkdvtVm1tbdr1Q4cO\nacWKFSorK9OCBQv0n//8J2teChMAAAAAcfkUJrFYTDU1NQqFQtq3b58aGhq0f//+pJh7771Xl156\nqV5//XVt3LhRP/vZz7LmpTABAAAAEJdPYRIOh+VyueR0OmW1WlVVVaVgMJgUs3//fl111VWSpNLS\nUrW3t+vjjz/OmJfCBAAAAEDcBJOWKhqNyuFwxI/tdrui0WhSTFlZmZ599llJJwuZ9957Tx0dHRnz\nUpgAAAAAiBtr0lIZhpGznzvuuEPd3d3y+Xyqq6uTz+fTeeedlzGeVbkAAAAAxBVZpJe+kHYM3LYk\nZQsTm82mSCQSP45EIrLb7UkxEyZM0OOPPx4/njlzpmbNmpUxLxssAgAAAJB0ciSkd2L6+dGfJm+w\n2NfXp9LSUm3btk0lJSWaP3++Ghoa5PV64zGffvqpioqKNHr0aD322GPasWOHnnzyyYy5GTEBAAAA\nEGcdkzvGYrGorq5OlZWVisViqq6ultfrVX19vSQpEAho3759uummm2QYhubMmaMNGzZk7ZMREwAA\nAACSTo6Y9H/D5Pz7ySMmw4EREwAAAAAJeYyYDAcKEwAAAAAJ4wqTlsIEAAAAQML5hUnLPiYAAAAA\nEsaYNBOhUEgej0dut1u1tbVp17u6urRkyRJdcsklmjNnTtYVuSQmvwMAAAA4xTAM9S8zOd+UPPk9\nFouptLRULS0tstlsqqioSFsu+K677tKJEyd03333qaurS6Wlpers7JTFYv7SFiMmAAAAABLON2kp\nwuGwXC6XnE6nrFarqqqqFAwGk2JmzJihw4cPS5IOHz6sqVOnZixKJOaYAAAAABgoj1W5otGoHA5H\n/Nhut2vXrl1JMWvWrNHVV1+tkpISHTlyRM8880zWPhkxAQAAAJAwzqSlMAwjZzf33nuvLrnkEn3w\nwQfau3evfvrTn+rIkSMZ4xkxAQAAAJBwvtTaIbVGM4fYbDZFIpH4cSQSkd1uT4p5+eWX9etf/1qS\ndNFFF2nmzJl66623VF5ebtonhQkAAACAhDHS4otOti+teyU5pLy8XG1tbWpvb1dJSYkaGxvV0NCQ\nFOPxeNTS0qKFCxeqs7NTb731lmbNmpUxLYUJAAAAgIQ85phYLBbV1dWpsrJSsVhM1dXV8nq9qq+v\nlyQFAgHdeeedWr16tcrKyvTFF1/ogQce0JQpUzL2yXLBAAAAACSdWi74XpPzdyYvFzwcGDEBAAAA\nkJDHiMlwoDABAAAAkFCgwoTlggEAAAAk5LFcsCSFQiF5PB653W7V1tamXX/wwQfl8/nk8/k0d+5c\nWSwWdXd3Z0zLHBMAAAAAkk7NMWkwOb8yeY5JLBZTaWmpWlpaZLPZVFFRoYaGBnm9XtN+//nPf+qh\nhx5SS0tLxtyMmAAAAABIGGPSUoTDYblcLjmdTlmtVlVVVSkYDGbs8umnn9bKlSuzpqUwAQAAAJBw\nvklLEY1G5XA44sd2u13RqPmOjD09Pdq6dauuv/76rGmZ/A4AAAAgIcOckoEMw8i7u+eee07f+ta3\nNGnSpKxxFCYAAAAAEsZIrbul1tcyh9hsNkUikfhxJBKR3W43jd20aVPO17gkJr8DAAAAOMUwDPXv\nNTl/SfLk976+PpWWlmrbtm0qKSnR/PnzTSe/f/rpp5o1a5Y6OjpUVFSUNTcjJgAAAAAS8niVy2Kx\nqK6uTpWVlYrFYqqurpbX61V9fb0kKRAISJI2b96sysrKnEWJxIgJAAAAgFMMw1D/+ybnv5E8YjIc\nGDEBAAAAkGCyCteZQGECAAAAIO5EHq9yDYfTso9Jru3ogZHO6XRq3rx58vl8mj9/viTp4MGD+s53\nvqOLL75Y1157rbq7uwt8l8DwuPnmm1VcXKy5c+fGz2X7/t93331yu93yeDx6/vnnC3HLwLAxex7u\nuusu2e12+Xw++Xw+NTc3x6/xPOBsdGLM6LRmJp8aoLW1VT6fT3PmzNHixYuz5h3yHJPBbkcPjEQz\nZ87U7t27NWXKlPi522+/XdOmTdPtt9+u2tpaHTp0SPfff38B7xIYHi+++KLGjx+vG2+8UW+++aak\nzN//ffv26Qc/+IFeeeUVRaNRXXPNNXr77bc1ahT7/WJkMHse1q1bpwkTJujnP/95UizPA85GhmHo\n4/7xaeenG0eT5pjkUwN0d3dr4cKF2rp1q+x2u7q6ujRt2rSMuYf8zR/sdvTASJVa42/ZskWrVq2S\nJK1atUqbN28uxG0Bw27RokWaPHly0rlM3/9gMKiVK1fKarXK6XTK5XIpHA6f8XsGhovZ8yCZTxrm\necDZ6oTGpLVU+dQATz/9tK6//vr4/ibZihLpNBQmg9mOHhipDMPQNddco/Lycj322GOSpM7OThUX\nF0uSiouL1dnZWchbBM6oTN//Dz74IGkDLv7OwLli/fr1KisrU3V1dfzVRp4HnK0+U1FaS5VPDdDW\n1qaDBw/qqquuUnl5uf7yl79kzTvkwmQw29EDI9WOHTu0Z88eNTc365FHHtGLL76YdN0wDJ4VnLNy\nff95NjDSrV27Vu+++6727t2rGTNm6Be/+EXGWJ4HnA3yGTHJ57v6+eef67XXXlNTU5O2bt2q3/3u\nd2pra8sYP+RVuQazHT0wUs2YMUOSNH36dK1YsULhcFjFxcX66KOPdOGFF+rDDz/UBRdcUOC7BM6c\nTN//1L8zOjo6ZLPZCnWbwBkx8M//H//4x7ruuusk8Tzg7NWr0Xq19Zhebe3JGJNPDeBwODRt2jQV\nFRWpqKhI3/72t/X666/L7Xab9jnkEZPy8nK1tbWpvb1dvb29amxslN/vH2q3wNdGT0+Pjhw5Ikk6\nduyYnn/+ec2dO1d+v19PPfWUJOmpp57S8uXLC3mbwBmV6fvv9/u1adMm9fb26t1331VbW1t8JTtg\npPrwww/jv/7HP/4RX7GL5wFnqx6N1ezF03XjXd+Mt1T51ADf/e539dJLLykWi6mnp0e7du3S7Nmz\nM+Yd8ohJpu3ogXNFZ2enVqxYIUnq6+vTD3/4Q1177bUqLy/XDTfcoA0bNsjpdOqZZ54p8J0Cw2Pl\nypXavn27urq65HA4dPfdd+uOO+4w/f7Pnj1bN9xwg2bPni2LxaI//elPvLqCESX1eVi3bp1aW1u1\nd+9eGYahmTNnqr6+XhLPA85evTJfHnigTDXAl9/vQCAgj8ejJUuWaN68eRo1apTWrFmTtTAZ8nLB\nAAAAAEYGwzDU3L847fxSo9V0dbnTiZ3fAQAAAMTlM2IyHChMAAAAAMT1aGxB8uac/J7PVvMAAAAA\nRoYTGp3WzOSqE1pbWzVx4kT5fD75fD7dc889WfNmHTGJxWKqqalJ2mre7/czuR0AAAAYoXpN9i1J\nlW+dcOWVV2rLli155c06YpLPVvMAAAAARo4eFaW1VPnWCYOZMJ91xMRsq/ldu3YlxTgNQ+/lnQ4A\nAAA4V3xT/f3thb6JQctnxCSfOsEwDL388ssqKyuTzWbTgw8++NX3MclnLe33JP1W0pOSbpK0Tr/N\n+ZmBfqt1ecfm6ju1r7MtPtXAzw8m9qvcy2Dvbaj5B17P9dkzbag/y3z+356U+fMwmJ9bPtdTne58\nQ8mdy+nuL1f/Q3neBtP36YgfiuH+uX7V/E/q9DwT2WIH63T/vgzleRvuZ3uw/Q13/sEYat+n+zt/\ner6TT0q66Yz/WZTr80PJfbo/P5z9Dfefg6nZvo4yzSkZKJ864dJLL1UkEtHYsWPV3Nys5cuX6+23\n384Yn7UwyWereUlqldR96r9SuyRnzhsFAAAARpb2U+3r7TON1Xut7XqvNfN7UfnUCRMmTIj/eunS\npfrJT36igwcPasqUKaZ9Zi1MBm41X1JSosbGRjU0NKTFLdbJ34LFkrZTlAAAAOCc5FTyP9BvL8xt\nDNEJjdaFiy/WhYsvjp97ad0LSTH51AmdnZ264IILZBiGwuGw+vv7MxYlUh47vzc3N+u2226LbzX/\nq1/9Kun64sWLtX371/OHDgAAAAyXK6+8Uq2trYW+jUExDEP/2/9A2vn1xu1pE9nN6oT6+npJUiAQ\n0COPPKI///nPslgsGjt2rP7whz/o8ssvz5w7V2ECAAAA4NxgGIZu6X8o7fyjxm2DWmHrq2DndwAA\nAABxnxVo53cKEwAAAABx+azKNRyybrAIAAAA4NzSqzFpzUwoFJLH45Hb7VZtbW3G/l555RVZLBY9\n++yzWfMyYgIAAAAgriePV7lisZhqamrU0tIim82miooK+f1+eb3etLhf/vKXWrJkSc45KoyYAAAA\nAIjr1ei0liocDsvlcsnpdMpqtaqqqkrBYDAtbv369fr+97+v6dOn58xLYQIAAAAg7oTGpLVU0WhU\nDocjfmy32xWNRtNigsGg1q5dKyn3bvG8ygUAAAAgzmyEJFWuIkOSbrvtNt1///0yDEP9/f05X+Wi\nMAEAAAAQ16OxOtb6qnpaX80YY7PZFIlE4seRSER2uz0pZvfu3aqqqpIkdXV1qbm5WVarVX6/37RP\nNlgEAAAAIOnkSMhF/f+Xdv4dY07SiEdfX59KS0u1bds2lZSUaP78+WpoaEib/P6l1atX67rrrtP3\nvve9jLkZMQEAAAAQZzanJJXFYlFdXZ0qKysVi8VUXV0tr9er+vp6SVIgEBh0XkZMAAAAAEg6OWIy\nNdaRdv6T8+w554gMFSMmAAAAAOJOHM89YjIcKEwAAAAAxPUez70q13BgHxMAAAAAcb1Hx6Y1M6FQ\nSB6PR263W7W1tWnXg8GgysrK5PP5dNlll+lf//pX1rzMMQEAAAAg6dT+JG+ZlAelRtIck1gsptLS\nUrW0tMhms6mioiJtVa5jx45p3LhxkqQ333xTK1as0IEDBzLmZsQEAAAAQMJxk5YiHA7L5XLJ6XTK\narWqqqpKwWAwKebLokSSjh49qmnTpmVNS2ECAAAAIOEzk5YiGo3K4XDEj+12u6LRaFrc5s2b5fV6\ntXTpUv3xj3/MmpbCBAAAAEDCMZOWwjCMvLpavny59u/fr+eee04/+tGPssayKhcAAACAhOOS3miV\n3mzNGGKz2RSJROLHkUhEdrs9Y/yiRYvU19enTz75RFOnTjWNoTABAAAAkHBc0sWLT7YvPb0uKaS8\nvFxtbW1qb29XSUmJGhsb1dDQkBTzzjvvaNasWTIMQ6+99pokZSxKJAoTAAAAAAMdzR1isVhUV1en\nyspKxWIxVVdXy+v1qr6+XpIUCAT097//XRs3bpTVatX48eO1adOmrH2yXDAAAAAASafmjtSblAeB\n5OWChwMjJgAAAAASTJYHPhMoTAAAAAAkmCwPfCawXDAAAACAhDyWC5akUCgkj8cjt9ut2tratOt/\n+9vfVFZWpnnz5mnhwoV64403sqZlxAQAAABAQh4jJrFYTDU1NWppaZHNZlNFRYX8fr+8Xm88Ztas\nWXrhhRc0ceJEhUIh3XLLLdq5c2fGPhkxAQAAAJBwwqSlCIfDcrlccjqdslqtqqqqUjAYTIq54oor\nNHHiREnSggUL1NHRkTUthQkAAACAhKMmLUU0GpXD4Ygf2+12RaPRjF1u2LBBy5Yty5qWV7kAAAAA\nJOSxKpdhGHl39+9//1uPP/64duzYkTWOwgQAAABAwnFJH7dKXa0ZQ2w2myKRSPw4EonIbrenxb3x\nxhtas2aNQqGQJk+enDUtGywCAAAAkHRqJGSpSXnQnLzBYl9fn0pLS7Vt2zaVlJRo/vz5amhoSJr8\n/v777+vqq6/WX//6V11++eU5czNiAgAAACDBZE5JKovForq6OlVWVioWi6m6ulper1f19fWSpEAg\noLvvvluHDh3S2rVrJUlWq1XhcDhjn4yYAAAAAJB0asSk3KQ8eDV5xGQ4MGICAAAAIMFkeeAzgcIE\nAAAAQEIer3INBwoTAAAAAAl5LBc8HNhgEQAAAEDCcZNmIhQKyePxyO12q7a2Nu36f//7X11xxRU6\n//zz9fvf/z5nWkZMAAAAACQcyR0Si8VUU1OjlpYW2Ww2VVRUyO/3Jy0XPHXqVK1fv16bN2/OKy0j\nJgAAAAAS+kxainA4LJfLJafTKavVqqqqKgWDwaSY6dOnq7y8XFarNa+0FCYAAAAABiUajcrhcMSP\n7Xa7otHokPrkVS4AAAAAA3yeM8IwjNOelcIEAAAAwABHJL0kaUfGCJvNpkgkEj+ORCKy2+1Dykph\nAgAAAGCAHkmXnmpfeiApory8XG1tbWpvb1dJSYkaGxvV0NBg2lu+O8ZTmAAAAAAY4LOcERaLRXV1\ndaqsrFQsFlN1dbW8Xq/q6+slSYFAQB999JEqKip0+PBhjRo1Sg8//LD27dun8ePHm/Zp9OdbwgAA\nAAAY0U7OHdltcuWyvEc+vipGTAAAAAAMkHvEZDhQmAAAAAAYoDCFCfuYAAAAABjgM5OWLhQKyePx\nyO12q7a21jTm1ltvldvtVllZmfbs2ZM1K4UJAAAAgAEOm7RksVhMNTU1CoVC2rdvnxoaGrR///6k\nmKamJh04cEBtbW169NFHtXbt2qxZKUwAAAAADJB7xCQcDsvlcsnpdMpqtaqqqkrBYDApZsuWLVq1\napUkacGCBeru7lZnZ2fGrBQmAAAAAAbIXZhEo1E5HI74sd1uVzQazRnT0dGRMSuT3wEAAAAMkP7q\nVqqTywrnlrrEcLbPUZgAAAAAGOA3aWdSN0W02WyKRCLx40gkIrvdnjWmo6NDNpstY1Ze5QIAAAAg\n6eQIh1k7cuRIUlx5ebna2trU3t6u3t5eNTY2yu/3J8X4/X5t3LhRkrRz505NmjRJxcXFGXMzYgIA\nAABgUCwWi+rq6lRZWalYLKbq6mp5vV7V19dLkgKBgJYtW6ampia5XC6NGzdOTzzxRNY+jf7h3lse\nAAAAAHLgVS4AAAAABUdhAgAAAKDgKEwAAAAAFByFCQAAAICCozABAAAAUHAUJgAAAAAKjsIEAAAA\nQMFRmAAAAAAouP8HE7/o5pQM9Z4AAAAASUVORK5CYII=\n",
"text": [
"<matplotlib.figure.Figure at 0x6f29610>"
]
}
],
"prompt_number": 23
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment