Skip to content

Instantly share code, notes, and snippets.

@jamesmcm

jamesmcm/gist:35f639c4bf6c61bd4f6b

Created July 9, 2014 14:15
Show Gist options

  • Select an option

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

Select an option

Learn more about clone URLs
Save jamesmcm/35f639c4bf6c61bd4f6b to your computer and use it in GitHub Desktop.
Download ZIP
Raw
gistfile1.txt
{
"metadata": {
"name": "",
"signature": "sha256:e846f8520620b644df3be3ae6c719db97d8028c34840710ae36f7a56ebc7b724"
},
"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",
"Xbk=X.copy()\n",
"#Y = WX + N , W is sparse\n",
"SNPsconsider=np.array(range(10))*1000\n",
"\n",
"ng=1\n",
"W=np.zeros((ng, X.shape[0]))\n",
"\n",
"\n",
"W[:,[SNPsconsider[1]]]=2\n",
"W[:,[SNPsconsider[5]]]=3\n",
"W[:,[SNPsconsider[3]]]=5\n",
"#W[:,2553]=1.5\n",
"W[:,7473]=3\n",
"W[:,6849]=6\n",
"\n",
"print W.shape\n",
"print X.shape\n",
"print np.dot(W,X).shape\n",
"\n",
"Y=np.dot(W,X) + np.random.normal(size=(ng,182))\n",
"X=X[SNPsconsider,:]\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",
"print SNPsconsider\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n",
"(1, 13314)"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"(13314, 182)\n",
"(1, 182)\n",
"(1, 182)\n",
"(1, 182)"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"[ 0 1000 2000 3000 4000 5000 6000 7000 8000 9000]\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAywAAAB+CAYAAADGIJX8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFwdJREFUeJzt3X9U1XWex/HXN8ESLX+UEAJnMUHBMrxJ0bpzjzQJBhkj\n7WyJ7egElqtj5dF+zjRztN9u08yWNg25pTntIm3nGMwidzg0Xcah1Ttb0LbhUWql+KGs5dAopiTd\n/cPtXj73XuBe/HGv+nyc8zmH9/f7uV/eXr7fC2+/n+/nY7ndbrcAAAAAIAJdEO4EAAAAAKA/FCwA\nAAAAIhYFCwAAAICIRcECAAAAIGJRsAAAAACIWBQsAAAAACIWBQsAAAAASdIIy5IVoI0bNy5sOVms\nwwIAAABAkizL0hMBtj8qKVxlQ1RYvisAAACAiHRJuBPwQcECAAAAwGNEuBPwQcECAAAAwIOCBQAA\nAEDEYkgYAAAAgIjFHRYAAAAAEYuCBQAAAEDECqZgKS4uVlVVlWJjY/Xhhx9KkubPn6/du3dLkrq6\nujRmzBg1NDT4vTY5OVmXXHKJhg0bpujoaLlcrgG/F+uwAAAAAJB0Yh2WPQG2T5a5Dsv27ds1atQo\nLVy40FOw9HX//fdrzJgxevTRR/32TZw4Ue+9917Qi1FyhwUAAACAR0wQfex2u1paWgLuc7vdeuON\nN/TOO+/0+/pQ7plcEHRPAAAAAOe8EVH+LRTbt29XXFycJk2aFHC/ZVmaPXu2MjMztWHDhkGPxx0W\nAAAAAB4Xjwyw8cvgX19WVqYFCxb0u7++vl7x8fE6cOCAcnJylJaWJrvd3m9/ChYAAAAAHtEXSs6e\nEy1Ux48f19atW/X+++/32yc+Pl6SNH78eBUWFsrlcg1YsDAkDAAAAIDXRVL2JdLqy7wtWLW1tUpP\nT9eECRMC7j9y5IgOHTokSeru7lZNTY2mTZs24DEpWAAAAAB4XRig+SgqKtLMmTO1Z88eJSUlaePG\njZKk8vJyFRUVGX07Ojp08803S5L2798vu92u6dOnKysrS3PnzlVubu6A6TCtMQAAAABJJx6Id08P\nsL0xtJm9TiWeYQEAAADgdVG4EzBRsAAAAADwCjAELJwoWAAAAAB4BZrWOIwoWAAAAAB4RdiQMGYJ\nAwAAAOAVxCxhxcXFiouLM6YkXr16tRITE2Wz2WSz2eRwOAIe3uFwKC0tTampqVq7du2g6VCwAAAA\nAPAaGaD5uPPOO/0KEsuytHLlSjU0NKihoUE33XST3+t6e3u1fPlyORwONTU1qaysTLt27RowHQoW\nAAAAAF4XBWg+7Ha7xo4d67d9sKmPXS6XUlJSlJycrOjoaM2fP18VFRUDvoaCBQAAAIBXEEPC+rNu\n3TplZGSopKREXV1dfvvb29uVlJTkiRMTE9Xe3j7gMSlYAAAAAHgNsWBZunSp9u7dq8bGRsXHx2vV\nqlV+fSzLCjkdZgkDAAAA4DVKcv7PiRaK2NhYz9eLFy/WLbfc4tcnISFBra2tnri1tVWJiYkDHpeC\nBQAAAIDXhVJ2+on2rTW/H/xl+/btU3x8vCRp69atxgxi38rMzFRzc7NaWlo0YcIElZeXq6ysbMDj\nUrAAAAAA8ApiCFhRUZHq6ur0+eefKykpSWvWrJHT6VRjY6Msy9LEiRNVWloqSero6NBdd92lqqoq\nRUVFaf369ZozZ456e3tVUlKi9PT0Ab+X5R7sUX4AAAAA5wXLsuTeGGD7nYPPAHa6cIcFAAAAgFeE\nrXRPwQIAAADAK4RpjM8EChYAAAAAXtxhAQAAABCxRoY7ARMLRwIAAADwCmLhyOLiYsXFxRlTFz/w\nwANKT09XRkaGbr31Vn355ZcBD5+cnKyrr75aNptN11133aDpULAAAAAA8LooQPNx5513yuFwGNty\nc3P10Ucf6YMPPtDkyZP19NNPBzy8ZVlyOp1qaGiQy+UaNB0KFgAAAABeIwM0H3a7XWPHjjW25eTk\n6IILTpQXWVlZamtr6/dbhDJFMgULAAAAAK8ghoQN5tVXX1V+fn7AfZZlafbs2crMzNSGDRsGPRYP\n3QMAAADwOslZwp588kkNHz5cCxYsCLi/vr5e8fHxOnDggHJycpSWlia73d7v8ShYAAAAAHgcGynV\n/UH6w/bQX7tp0yZt27ZNb7/9dr994uPjJUnjx49XYWGhXC7XgAXLKRkS5nA4lJaWptTUVK1du/ZU\nHBI4qwSa7eLgwYPKycnR5MmTlZubq66urjBnCZwegWaKGej8f/rpp5Wamqq0tDTV1NSEI2XgtAl0\nPaxevVqJiYmy2Wyy2Wyqrq727ON6QCQ6duFwXZ8zXA8+5m3BcDgcevbZZ1VRUaGLLgp8m+bIkSM6\ndOiQJKm7u1s1NTXG9RLISRcsvb29Wr58uRwOh5qamlRWVqZdu3ad7GGBs0qg2S6eeeYZ5eTkaM+e\nPbrxxhv1zDPPhDlL4PQINFNMf+d/U1OTysvL1dTUJIfDoWXLlumbb74JR9rAaRHoerAsSytXrlRD\nQ4MaGhqUl5cniesBkatn2HC/5quoqEgzZ87U7t27lZSUpFdffVX33HOPDh8+rJycHNlsNi1btkyS\n1NHRoZtvvlmStH//ftntdk2fPl1ZWVmaO3eucnNzB8znpIeEuVwupaSkKDk5WZI0f/58VVRUKD09\n/WQPDZxVfGe7qKysVF1dnSRp0aJFys7OpmjBOclut6ulpcXY1t/5X1FRoaKiIkVHRys5OVkpKSly\nuVy6/vrrw5A5cOoFuh6kwDMicT0gUh0L+JT9YSMqKyvz61FcXBzweBMmTFBVVZUk6YorrlBjY2NI\n+Zz0HZb29nYlJSV54sTERLW3t5/sYYGzSqDZLjo7OxUXFydJiouLU2dnZzhTBM6o/s7/jo4OJSYm\nevrxOwPni3Xr1ikjI0MlJSWeIZJcD4hUX2mEXwunky5YLMs6FXkAZ7X6+no1NDSourpaL774orZv\nN59SsyyLawXnrcHOf64NnOuWLl2qvXv3qrGxUfHx8Vq1alW/fbkeEAmO6UK/Fk4nXbAkJCSotbXV\nE7e2thr/WwCcDwLNdhEXF6f9+/dLkvbt26fY2NhwpgicUf2d/76/M9ra2pSQkBCWHIEzJTY21lO4\nL1682POsI9cDIlWPhvu1cDrpgiUzM1PNzc1qaWlRT0+PysvLVVBQcCpyA84K/c12UVBQoNdee02S\n9Nprr2nevHnhTBM4o/o7/wsKCrRlyxb19PRo7969am5u9sysB5yr9u3b5/l669atnhmRuB4QqY4o\nxq+F00k/dB8VFaX169drzpw56u3tVUlJCQ/c47zS2dmpwsJCSdLx48d1xx13KDc3V5mZmbrtttv0\nyiuvKDk5WW+88UaYMwVOj6KiItXV1enzzz9XUlKSHnvsMT388MMBz/+pU6fqtttu09SpUxUVFaVf\n/epXDIHBOcX3elizZo2cTqcaGxtlWZYmTpyo0tJSSVwPiFzB3FEpLi5WVVWVYmNj9eGHH0o6MaX9\n7bffrk8//dTz2T9mzBi/1zocDq1YsUK9vb1avHixHnrooQG/l+UONG0FAAAAgPOOZVmqdmf7bc+z\nnMZsd9u3b9eoUaO0cOFCT8Hy4IMP6rLLLtODDz6otWvX6s9//rPfDKm9vb2aMmWKamtrlZCQoGuv\nvVZlZWUD3vA4JQtHAgAAADg3BPMMi91u19ixY41tlZWVWrRokaQTU9q/9dZbfq/ruyRKdHS0Z0mU\ngVCwAAAAAPAY6jMswSzpMJQlUQYtWBwOh9LS0pSamqq1a9cGlSwAAACAs9MxDfdroepvSvuhPKc1\n4EP3vb29Wr58uTHGrKCggIfqAQAAgHNUjy7Ubud+7XaGtuj1t1PaX3755f0u6TCUJVEGvMMylDFm\nAAAAAM5eRzRCSdkTNXv19Z4WjGCWdBjKkigD3mEJNMZs586dRh9b9ig11nUH9Y8AAAAAzh+Zcrv/\nFO4kQtYTxMr2oUxp39HRobvuuktVVVVDWhJlwIIlmDFmjXXd2u6+Rvdk79Y65xTZb3nP2F/zW7sR\n3/xFlRHnX2rG91sLjPgT922er3+4s9z85m+aYfdLZr4jd5ozNrsrfP49PsXip98db8Tlut2IH/rk\nBSNePelhI96rZCMulDkzQrL2GnFvn7d/xicfGfsen/SAEf/U+pERjz46woi7nos34rofmwtPZf+L\nWWjefcfzRvwzPWbEiX/3hRHv+bckI079hzYjTvr1HiNu253qDXx+TjU/Mc+JnE1/NGLrOZ+f2yzz\n5/bFr83jXdpoxjunmXHGYTMeMdc8/ivvmOdcyZP/asS//cmNRnzzpb/3fD22Y5+xb8VFJ34OmyT9\nUNKaP/rMGp581AhXJfzciH/+k58acdmT3zPioufNO5wXL/5fIz68wzyH3U+Z75010+e9XWrunzPB\ne85eJ/OcefzFp434ox9dYcQr9Qsj7lScEX9gzTC/d4M5HvbSaeY5df8w871p8bm+ktVixElqNeIf\nPO9z4q2o9HxZ437W2JXzkHkO2tfWGPG9Mq/9eV/+uxEPv8p8XxNbm4249bbJRjz9jf8w4jiZt9x/\noZVGfFX9J0a86W/Mz6a+//ZeDTP25f5uuxHPmuMw4mdkfo5VKd+IH68zf+5W9nNG/LK7wYjvKn3d\niOctOXE9/TH7MX3H+TMV6LfG/uHqMeIsn/Mu9XXvebH+70uMfbNVa8SbtdCIO2R+Lr5mX2rEr2w3\nr/13NdOIfd+LO2R+Njx38aNGvO7QYiN+pNt8794YZQ6NyP/4G8/X7v0+Ax4uNcMR8eZn8sOjzWdK\nff/tG3SXEd947btGPP5PnxnxO7rBiCtl/m/nj63RZkLJq40wca/POV/nPef3zzJf6/s+D9cxI577\ng98bccVvco34e9nm9fmEc5UR12q2Efue43/9U/OXRs3jPr+TnjI/D9zmKa+qL77r+dr3c2+zFhnx\n9Ts/MOIZWSeux93Z92iKc53et33HPHimGbqv8vmMvs/8u+J3NvN/r8sazN8ZC+p9ZmhK/NoIr/sr\n87y4R+s8X++U+XfEup3mOhlTs8y/9z7XZUbs0E1GfE3lLjMXn3PcOm5+jl41y/wj/8NmM5+ZqW8b\nse85v1tTjPif5b0+qz+61dj31pVzjHiedbcRv+w238e7P/qNEbsLfX7XNvf9m+5KnY2CeWalrKws\n4Pba2lq/bRMmTFBVlfdv/ry8POXl5QWdz4AFS7BjzF5d3aF9LT16dXWH9LlTuiw76AQAAACAc4Hz\niCS9GO40TtpXYV7Z3teAz7AEO8asePUExScPV/HqCRQrAAAAOC9lx0jSj/q0s9OpmCXsVBp0pfvq\n6mqtWLHCM8bskUceMfZnZ2errq7utCYJAAAAnG1mzZolp9MZ7jRCYlmW7nH/o9/2ddaDGqRsOG0G\nLVgAAAAAnB8sy9Ld7n/y2/6ytSJsBQsr3QMAAADw+Eoxfs3X7t27ZbPZPG306NF64QVzkhqn06nR\no0d7+jzxxBNDymfAh+4BAAAAnF+CeWZlypQpamg4MVPkN998o4SEBBUWFvr1mzVrliorK/22h4KC\nBQAAAIBHMOuw9FVbW6tJkyYZ6zd+61QMI2NIGAAAAACPI4rxawPZsmWLFixY4Lfdsiy9++67ysjI\nUH5+vpqamoaUDw/dAwAAAJB0osjIdlf7bXdaeQHvlvT09CghIUFNTU0aP95cwPrQoUMaNmyYYmJi\nVF1drfvuu0979uzxO8ZgGBIGAAAAwOOYLtSXzkb9xdk4aN/q6mrNmDHDr1iRpIsvvtjzdV5enpYt\nW6aDBw9q3LhxIeVDwQIAAADAo0fDNSL7Oo3Ivs6zrW3N5oB9y8rKVFRUFHBfZ2enYmNjZVmWXC6X\n3G53yMWKRMECAAAAoI/Bnln5Vnd3t2pra7VhwwbPttLSUknSkiVL9Oabb+qll15SVFSUYmJitGXL\nliHlwzMsAAAAACSdeIZlkvu//bZ/Yl0VtoUjucMCAAAAwONYiNMan24ULAAAAAA8vvpmRLhTMFCw\nAAAAAPA4dpQ7LAAAAAAiVM/R4eFOwUDBAgAAAMCj53Bws4QlJyfrkksu0bBhwxQdHS2Xy+XX5957\n71V1dbViYmK0adMm2Wy2kPOhYAEAAADgdTS4EsGyLDmdzn7XVtm2bZs+/vhjNTc3a+fOnVq6dKl2\n7NgRcjoXhPwKAAAAAOeuowFaPwaa6riyslKLFi2SJGVlZamrq0udnZ0hp0PBAgAAAMDrqwAtAMuy\nNHv2bGVmZhqLR36rvb1dSUlJnjgxMVFtbW0hp8OQMAAAAABe3cF1q6+vV3x8vA4cOKCcnBylpaXJ\nbrcbfXzvwFiWFXI6FCwAAAAAvI5K+i+n9KFzwG7x8fGSpPHjx6uwsFAul8soWBISEtTa2uqJ29ra\nlJCQEHI6DAkDAAAA4HVU0uRs6W9Xe5uPI0eO6NChQ5Kk7u5u1dTUaNq0aUafgoICbd68WZK0Y8cO\njRkzRnFxcSGnwx0WAAAAAF6HB+/S2dmpwsJCSdLx48d1xx13KDc3V6WlpZKkJUuWKD8/X9u2bVNK\nSopGjhypjRs3Dikdyz3Qo/0AAAAAzhuWZUmlAcqDJdaAM4KdTtxhAQAAAOA1wDTG4UDBAgAAAMCr\nn2mMw4WCBQAAAIBXkNManykULAAAAAC8IuwOC9MaAwAAAPA6FqD5aG1t1Q033KArr7xSV111lV54\n4QW/Pk6nU6NHj5bNZpPNZtMTTzwxpHS4wwIAAADAK4hpjaOjo/XLX/5S06dP1+HDhzVjxgzl5OQo\nPT3d6Ddr1ixVVlaeVDrcYQEAAADgdTRA83H55Zdr+vTpkqRRo0YpPT1dHR0dfv1OxVTIFCwAAAAA\nvIIoWPpqaWlRQ0ODsrKyjO2WZendd99VRkaG8vPz1dTUNKR0GBIGAAAAwCuEdVgOHz6s73//+3r+\n+ec1atQoY98111yj1tZWxcTEqLq6WvPmzdOePXtCToeV7gEAAABI+v+V7u1uqcspfen07vhsjd/w\nrq+//lpz585VXl6eVqxYMeixJ06cqPfee0/jxo0LKSfusAAAAADw+krShdlSbLZ322drjC5ut1sl\nJSWaOnVqv8VKZ2enYmNjZVmWXC6X3G53yMWKRMECAAAAoK8A0xj7qq+v1+uvv66rr75aNptNkvTU\nU0/ps88+kyQtWbJEb775pl566SVFRUUpJiZGW7ZsGVI6DAkDAAAAIOn/h4RNDFAe7LVOyYxfQ8Ed\nFgAAAABeITx0fyZQsAAAAADwomABAAAAELEOhTsBEwULAAAAAK/j4U7AxEr3AAAAAELmcDiUlpam\n1NRUrV27NmCfe++9V6mpqcrIyFBDQ8OQvg8FCwAAAIA+vg7QTL29vVq+fLkcDoeamppUVlamXbt2\nGX22bdumjz/+WM3NzXr55Ze1dOnSIWVDwQIAAACgj0MBmsnlciklJUXJycmKjo7W/PnzVVFRYfSp\nrKzUokWLJElZWVnq6upSZ2dnyNlQsAAAAADo40iAZmpvb1dSUpInTkxMVHt7+6B92traQs6Gh+4B\nAAAA9PHVoD0sywrqSL6LTQb7ur4oWAAAAAD0cUjSf0p6r98eCQkJam1t9cStra1KTEwcsE9bW5sS\nEhJCzoYhYQAAAAD6+ErSlZIW9mmmzMxMNTc3q6WlRT09PSovL1dBQYHRp6CgQJs3b5Yk7dixQ2PG\njFFcXFzI2XCHBQAAAEAfgw8Ji4qK0vr16zVnzhz19vaqpKRE6enpKi0tlSQtWbJE+fn52rZtm1JS\nUjRy5Eht3LhxSNlYbt+BZQAAAADOSyeeMakMsKfA73mUM4U7LAAAAAD6+Eu4EzBQsAAAAADoY/Ah\nYWcSBQsAAACAPihYAAAAAESsyBoSxkP3AAAAACT1v7Dj2LFjdfDgwTOczQncYQEAAAAgyX9l+kjA\nwpEAAAAAIhYFCwAAAICIRcECAAAAIGJRsAAAAACIWBQsAAAAACLW/wHrgHWalWO9ogAAAABJRU5E\nrkJggg==\n",
"text": [
"<matplotlib.figure.Figure at 0x47a6a90>"
]
}
],
"prompt_number": 444
},
{
"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.preferred_optimizer='simplex' #use this or not\n",
"\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",
" k = GPy.kern.linear(input_dim=step+1, ARD=True)\n",
" tempm = GPy.models.GPRegression(X2,Y,k)\n",
" #tempm.preferred_optimizer='simplex'\n",
"\n",
" tempm.optimize(messages=False)\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 tempm.log_likelihood() > tempbestll:\n",
" hasbeenadded=True\n",
" tempbestll=tempm.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, ARD=True)\n",
" tempm = GPy.models.GPRegression(X2,Y,k)\n",
" #tempm.preferred_optimizer='simplex'\n",
"\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",
"print chosenSNPS\n",
"print m"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"Log-likelihood: -6.827e+02\n",
"\n",
" Name | Value | Constraints | Ties | prior \n",
"-------------------------------------------------------------------\n",
" linear_variance | 1.0000 | (+ve) | | \n",
" noise_variance | 105.9333 | (+ve) | | \n",
"\n",
"-682.668583026\n",
"-601.262003795"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-549.586714841"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-532.459621967"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-522.263661533"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-520.179952453"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-518.823440096"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-518.564972765"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"-518.564253767"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"converged: 8SNPs"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAyYAAACICAYAAAASyf8xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHiBJREFUeJzt3X90k+X9//HXzRIEwS+oxUqTfAzY2oYjVlwLetTJ3KSW\nnVVR58p2lGFhHaMim+7gNp3g2QFx2x9K3U7lIJOppc45ijsQXTkLKB6oPxCd7bDsEEmDduIsIgil\nsd8/KmmTJrmT5VcLz8c594HkvnJdV9LcSd+9rvd1GT09PT0CAAAAgCwalu0OAAAAAACBCQAAAICs\nIzABAAAAkHUEJgAAAACyjsAEAAAAQNYRmAAAAADIOgITAAAAAAm54447lJubq8mTJ0cts2jRIhUU\nFKi4uFi7du0yrZPABAAAAEBC5s6dK7fbHfX8pk2btHfvXrW1tenxxx/XggULTOskMAEAAACQkKuv\nvlpnn3121PMbN27UnDlzJEnTpk1TZ2enOjo6YtZJYAIAAAAgpfx+vxwOR/C23W5Xe3t7zMcQmAAA\nAABIuZ6enpDbhmHELE9gAgAAAECSNNIwZEQ4zjrrrITqsdls8vl8wdvt7e2y2WwxH2P5n3oMAAAA\n4JRzTNKvI9x/32efJVRPRUWFamtrVVlZqR07dmjs2LHKzc2N+RgCEwAAAABB/y+OMrNnz9bWrVt1\n8OBBORwOLVu2TCdOnJAkVVdXa+bMmdq0aZPy8/M1atQorV271rROoyd88hcAAACA05JhGFod4f75\nGpgzkmqMmAAAAAAIGpmldglMAAAAAATFM5UrHViVCwAAAEDQyAhHJG63W0VFRSooKNDKlSsHnP/k\nk080a9YsFRcXa9q0aXr33XdjtktgAgAAACAonsAkEAiopqZGbrdbLS0tqq+vV2tra0iZ5cuX67LL\nLtPu3bu1bt063XXXXTHbJTABAAAAEBRPYNLc3Kz8/Hw5nU5ZrVZVVlaqsbExpExra6u+/vWvS5IK\nCwvl9Xr10UcfRW2XwAQAAABA0FkRjnB+v18OhyN42263y+/3h5QpLi7W888/L6k3kHn//ffV3t4e\ntV0CEwAAAABBZ0Y4whmGYVrPvffeq87OTk2ZMkW1tbWaMmWKvvKVr0Qtz6pcAAAAAIJGWqRXvpC2\n99+2JGwLE5vNJp/PF7zt8/lkt9tDypx11ll64okngrcnTJigiRMnRm2XDRYBAAAASOodCekaM/D+\n4YdCN1js7u5WYWGhtmzZory8PE2dOlX19fVyuVzBMocOHdLIkSM1fPhwrV69Wtu3b9cf//jHqG0z\nYgIAAAAgyHqGeRmLxaLa2lqVlZUpEAioqqpKLpdLdXV1kqTq6mq1tLToBz/4gQzD0MUXX6w1a9bE\nrJMREwAAAACSekdMev4vwv37Q0dM0oEREwAAAAB94hgxSQcCEwAAAAB9RmWnWQITAAAAAH1GZKdZ\n9jEBAAAA0OeMCEcEbrdbRUVFKigo0MqVKwecP3jwoK6//npdeumluvjii2OuyCWR/A4AAADgS4Zh\nqGdmhPs3hSa/BwIBFRYWqqmpSTabTaWlpQOWC166dKmOHz+uFStW6ODBgyosLFRHR4cslsiTthgx\nAQAAANBnRIQjTHNzs/Lz8+V0OmW1WlVZWanGxsaQMuPHj9enn34qSfr000917rnnRg1KJHJMAAAA\nAPQXx6pcfr9fDocjeNtut2vnzp0hZebPn69rr71WeXl5Onz4sJ599tmYdTJiAgAAAKDPqAhHGMMw\nTKtZvny5Lr30Uh04cEBvvfWWFi5cqMOHD0ctz4gJAAAAgD4jJE+75PFHL2Kz2eTz+YK3fT6f7HZ7\nSJlXX31Vv/zlLyVJF154oSZMmKA9e/aopKQkYp0EJgAAAAD6nCFNv7D3OGnZa6FFSkpK1NbWJq/X\nq7y8PDU0NKi+vj6kTFFRkZqamnTllVeqo6NDe/bs0cSJE6M2S2ACAAAAoE8cOSYWi0W1tbUqKytT\nIBBQVVWVXC6X6urqJEnV1dX6xS9+oblz56q4uFhffPGFHn74YZ1zzjlR60xJjonZGsbAqc7pdOqS\nSy7RlClTNHXqVEnSf//7X1133XW66KKLNGPGDHV2dma5l0B63HHHHcrNzdXkyZOD98V6/69YsUIF\nBQUqKirSSy+9lI0uA2kT6XpYunSp7Ha7pkyZoilTpmjz5s3Bc1wPGJRGRzgiKC8v1549e7R37179\n/Oc/l9QbkFRXV0uScnJy9MILL2j37t1655139L3vfS9ms0kHJoFAQDU1NXK73WppaVF9fb1aW1uT\nrRYYUgzDkMfj0a5du9Tc3CxJeuihh3Tdddfpvffe0ze+8Q099NBDWe4lkB5z586V2+0OuS/a+7+l\npUUNDQ1qaWmR2+3Wj3/8Y33xxRfZ6DaQFpGuB8Mw9NOf/lS7du3Srl27VF5eLonrAYNYnBssplrS\ngUk8axgDp4PwvUo3btyoOXPmSJLmzJmjDRs2ZKNbQNpdffXVOvvss0Pui/b+b2xs1OzZs2W1WuV0\nOpWfnx8M5oFTQaTrQRr4HSFxPWAQG6qBSaQ1jP3+GCn8wCnIMAx985vfVElJiVavXi1J6ujoUG5u\nriQpNzdXHR0d2ewikFHR3v8HDhwIWbWF7wycLlatWqXi4mJVVVUFpzZyPWDQimO5YMk8neO3v/1t\ncArj5MmTZbFYYk5tTzowiWcNY+BUt337du3atUubN2/WY489ppdffjnkvGEYXCs4bZm9/7k2cKpb\nsGCB9u3bp7feekvjx4/X3XffHbUs1wMGhTh2fo8nneOee+4JTmFcsWKFpk+frrFjx0ZtNunAJJ41\njIFT3fjx4yVJ48aN06xZs9Tc3Kzc3Fx9+OGHkqQPPvhA5513Xja7CGRUtPd/+HdGe3u7bDZbVvoI\nZMp5550XDNDnzZsXnK7F9YBBK46pXImmczzzzDOaPXt2zGaTDkz6r2Hc1dWlhoYGVVRUJFstMGQc\nPXo0uIvpkSNH9NJLL2ny5MmqqKjQk08+KUl68skndeONN2azm0BGRXv/V1RUaP369erq6tK+ffvU\n1tYWXMkOOFV98MEHwf//9a9/Da7YxfWAQSuOEZNE0jmOHj2qF198UTfffHPMZpPexyTaGsbA6aKj\no0OzZs2SJHV3d+v73/++ZsyYoZKSEt16661as2aNnE6nnn322Sz3FEiP2bNna+vWrTp48KAcDoce\nfPBB3XvvvRHf/5MmTdKtt96qSZMmyWKx6Pe//z1TV3BKCb8eli1bJo/Ho7feekuGYWjChAnBfR64\nHjBoRckp6S+R9+oLL7ygq666KuY0LkkyeiItEwEAAADgtGMYhnp2SJ43JM+bffcvWxO6utyOHTu0\ndOnS4PLYK1as0LBhw7RkyZIBdc6aNUvf/e53VVlZGbttAhMAAAAA0peByVsR7r80NDDp7u5WYWGh\ntmzZory8PE2dOlX19fUDZk4dOnRIEydOVHt7u0aOHBmz7aSncgEAAAA4hcQxlStaOsfJqYond3/f\nsGGDysrKTIMSKY4RE7fbrcWLFysQCGjevHkRh2cAAAAADH2GYahnf4T7/y/yRqEpbTtWYBIIBFRY\nWKimpibZbDaVlpZGHKIBAAAAMPQZhqGe/0S4/7z0ByYxlwtOdH1iAAAAAEPb8VEDj0yImWMSaX3i\nnTt3hpRxGobeT0/fAAAAgCHuAvX0eLPdiYQcP2N4hHu7BtwTT8qHx+PRT37yE504cUI5OTnyeDxR\n240ZmBiGoW3btik3N1fnnXdexMbel/SAJI+k6RHqWKYHQm4/oGWxmjQVXl8q606krWy0F0t4XxJ9\n3VP9XM3aj3U+1e8ZM8m87lL0/nnUe00kUn+yr7NZfYk+10TfN4mUT/Q9m+h5M6luP1ys52om2Z9r\nOn9uZqL3zaNI3xKZ/DlEksr3Wao/JxOV7DUSq/1kvyPS/X2ayfd0MnX18Sjyb03m7af7szCVbSX7\nWZTJ55revqX395h06PqKeWASCARUU1MTkvJRUVERkvLR2dmphQsX6sUXX5TdbtfBgwdjthtzKpfN\nZtPYsWOD6xP7fD7Z7fYB5TySvP3+BQAAAE5HXkm9vxWfPIae4zpjwBEunpSPZ555RjfffHMwfsjJ\nyYnZbszk95PrEz/55JP60Y9+JIvFMiD5vXfXx+hjJqmOrmPVlem/Dicqmb8Ipfov4+kubyaZ+lL9\n1+dU13/y8R7F+/ew6HWZSfeoXTpHJVJ9vWXyL2nx1J/Ovz6byfZzj1S3R4lfD5HaSvfndDIyPQqY\nyb8uD7bviGyO0CQ70rUsxu9M8Uj3d2Ai0j06m8nfIc0k97ovS3vSeCoZhqG2noEDEQVGe8jzeO65\n5/Tiiy9q9erVkqSnnnpKO3fu1KpVq4JlTk7hevfdd3X48GHddddduu2226K2HXMq18n1iW+//Xb5\n/X4tXbo0xopczlhVAacdZ7Y7AAwizmx3ABhUnNnuABBTpBGScL2DE7GdOHFCb775prZs2aKjR4/q\niiuu0OWXX66CgoKI5U03WLz44ouVm5urAwcO6Omnn9aoUaO0aNGiCCWdpp0DTifObHcAGESc2e4A\nMKg4s90BIKYuDdfrniN63XM0ahmbzSafzxe8HSnlw+FwKCcnRyNHjtTIkSP1ta99Tbt3744amMTM\nMZEkq9Wq+++/XwUFBdqxY4cee+wxtba2hpXy9Du8ZlUCAAAApySvpKGeY3JUZ2rS9HG6fekFwSNc\nSUmJ2tra5PV61dXVpYaGBlVUVISUueGGG/TKK68oEAjo6NGj2rlzpyZNmhS1XdMRk/PPPz9YwejR\no+VyuXTgwIGQKV0PaGvw/8vC5kyme7WnZGQ6LyNcMvNU0z2fN1ym50qnUrb7mswKZOmen5vu+cGZ\nnO8/2FZzSueqXGay+XPK9ns0Uan+LE2lwZa3ka660iGR6zPZusOlO781m6tyJcus79l+38SS3Gux\n1bzIINOlSKtyhTqZ8lFWVqZAIKCqqiq5XC7V1dVJkqqrq1VUVKTrr79el1xyiYYNG6b58+cnF5jM\nnj1bW7du1ccff6zx48frxIkTeuqppxJ4agAAAACGinhyTCSpvLxc5eXlIfdVV1eH3L7nnnt0zz33\nxFWf6VSutWvXyuFwqLCwUJ2dnbrmmms0evTouCoHAAAAMLR0afiAIxNMA5MRI0bopZde0vjx4/Xr\nX/9afr9fr7zySkgZj/rPovOmuo8AAADAEOHVqZBjEn5kgulUrp6eHt15552aNGmSFixYoPXr1+uc\nc84JKTO93/+3JrjSRCrneaY6N+DUX0c9+uPN+p7s+VjlM72Tc6plcr5+spLtazL5AonONR7s+1lk\nc250tnc/T+fO0em+frOZB2V2PtmfazbzZTKdz5btz/1YMp0Tmsw1lmjdZhLN9Ul1Tmis86neLyr2\n5+DQczzOERK3263FixcrEAho3rx5WrJkSch5j8ejG264QRMnTpQk3Xzzzbrvvvui1mcamGzfvl1/\n+tOfNGLECD366KPKycnR/v37YyauAAAAABiauuLIMQkEAqqpqVFTU5NsNptKS0tVUVExYM/Da665\nRhs3boyrXdPA5KqrrlJ3d7dKSkqUm5urzs5OjRgxIq7KAQAAAAwtRzXStExzc7Py8/PldDolSZWV\nlWpsbBwQmCSy671pjokkPfLII5o0aZKGDx+ub33rW3r99ddDzntEjgkAAADg1VDPMOkdMQk/wvn9\nfjkcjuBtu90uv98fUsYwDL366qsqLi7WzJkz1dLSErNd08Dk7bff1saNGzVv3jx1d3fr73//u6ZM\nmRJSZnq/g91MAQAAcLpyKvx346HnuIYPOMIZhmFaz2WXXSafz6fdu3frzjvv1I033hizvNFjMr4y\nY8YM7d+/X93d3frPf/6j+++/Xz/72c9COhUr9SjVGwElYrBtGphM+VQn9mU7kT+RhLRMJzlmc8O2\nVL/HMp0gnkj/0/1ZkOmN85J57uFSvXlrslKZOJzu92gmPx9SvXFmphO4k6k/05uAJvv4VL622dww\nOJJU/i6S6UT9RGXz52omdl+WJTSdKdsMw9Aveu7X+x6v3ve8H7z/lWXbQp7Hjh07tHTpUrndbknS\nihUrNGzYsAEJ8P1NmDBBb7zxxoCFtE6KOWLyt7/9TQUFBTp27Jh6enpkGIb+/Oc/J/TkAAAAAAwd\nxzVc50+/SNOWXhc8wpWUlKitrU1er1ddXV1qaGhQRUVFSJmOjo5gMNPc3Kyenp6oQYlkkvz+6quv\nauPGjfrwww91zjnnqLu7W0VFRQPKefr93ykmcwEAAOD05JU0dLNLesWzKpfFYlFtba3KysoUCARU\nVVUll8uluro6Sb07wD/33HP6wx/+IIvFojPPPFPr16+PXWesk8uXL9fy5cs1YcIEPfLII1q9erXW\nrVs3oNx0064DAAAApz6npPBd/oaaePcxKS8vV3l5ech91dXVwf8vXLhQCxcujLtd0+WCpd65Znff\nfbc++ugjrV69WvPnzw85n875iols3JWodPYt1Y9PZsPCeNrK9oaN6aornvrDZTI3KdObWWXyPWv2\n+GzPt8/k9Z/uDRBT3Z6ZdH7mp/r6T+ec+XTnIWY7tzDWY80kurGe2ePN+pPJ3MJ0b8yX6vOpfK1S\nnYub7vqT2eQ3mxuUDgafZ2in93Cmq3J1dnZq0qRJslgsys3N1cMPP6yXX345E30DAAAAkGHxrMqV\nDqaByV133aWbbrpJra2tevfdd3XLLbeoubk5rJSn3+FNdR8BAACAIcEraajvZBLPPiaS5Ha7VVRU\npIKCAq1cuTJqfa+99posFouef/75mO3GDEwOHTqkbdu26Tvf+Y4k6fjx49q2bZsmT54cVnJ6v8MZ\ns0EAAADgVOWUNNR3MjmqMwcc4QKBgGpqauR2u9XS0qL6+nq1trZGLLdkyRJdf/31pssmx8wx2bdv\nn8aMGaMLLrhAx44d0xlnnKHFixdrxowZIeWSmTceLp1zI7O930WicycTya9JZp5lPI9P9VrkqZzL\nmek9WNK5Rn26fw6pnlOb6Hs4kfyaRN/T4TKdnxOrP8lc6/HI9L4JsV6rVM/LTnW+XLJ5GcnkiKXz\n++1/kUj76c6DSPV3Uia/4wbbzzGTbWU6BzSVn8uZ/jkNNV1xTN1qbm5Wfn6+nE6nJKmyslKNjY1y\nuVwh5VatWqVbbrlFr732mmmdMUdMuru79c9//lPjxo1TYWGhrFarVqxYoUcffdS0YgAAAABDz3Gd\nMeAI5/f75XA4grftdrv8fv+AMo2NjVqwYIEk893iY46Y2O12ORwO7dmzR5K0bds2zZgxQ7NmzQop\n5wm55RXTuQAAAHA68koaqrklJ8UzYmIWZEjS4sWL9dBDD8kwDPX09CQ3lev888+Xw+HQe++9p4su\nukirV69WTk5OSHQkha/U7DTtJAAAAHAqckoa6vuYHNWZOuJ5XUc9r0ctY7PZ5PP5grd9Pp/sdntI\nmTfeeEOVlZWSpIMHD2rz5s2yWq0Ddog/yegxCV12796tefPmqaurS52dnVq0aJHuvvvuvgoMQ9ID\nOjlSkuo9IDIp3X1J5RxWs7qzvX5+Kn/Omc4hMRPva+FV74dTIq9dqvfSSOV83P9FOvdJSLStdL+2\nqZTp+fip/nyJ9FivIo+lZ/uzLdH+JFKXWd3JPt6sP5nM0Uz3ni5m5c1k8no1s6zf70yRHp/u557J\nfJpM5/amM+ckmWtgmWQ6UjCYGIahC3v+OeD+fxsXhzyP7u5uFRYWasuWLcrLy9PUqVNVX18/IMfk\npLlz5+rb3/62brrppqhtm26wWFxcrNdee01dXV2y2Wy6/fbbo5T0iilcQB+vuCKAk7ziegD6eMUV\ngcEsUk5JOIvFotraWpWVlSkQCKiqqkoul0t1dXWSQneAj1dcO79L0ubNm/XVr35V48aNi3DWo96L\nzMOlBgAAgNOYd4hnmEiffzEyrnLl5eUqLy8PuS9aQLJ27VrT+uIOTOrr6zV79uwoZ6erNziZLucQ\nnEcHAAAApIZziGeYSMePmY+YpINpjokkHTlyRBdccIH27duns846K+Tc9OnTtXXrUHzJAQAAgPS6\n5ppr5PF4st2NuBmGoeEfHxpwf9e5Y9KeKxNXYAIAAADg1GcYhvT+iYEnLrAOCEzcbrcWL16sQCCg\nefPmacmSJSHnGxsb9atf/UrDhg3TsGHD9Jvf/EbXXntt9LYJTAAAAABIXwYmeyKEB4VGSGASCARU\nWFiopqYm2Ww2lZaWDliV68iRIxo1apQk6Z133tGsWbO0d+/eqG3H3PkdAAAAwGnmWIQjTHNzs/Lz\n8+V0OmW1WlVZWanGxsaQMieDEkn67LPPlJOTE7NZAhMAAAAAfT6PcITx+/0hm67b7Xb5/f4B5TZs\n2CCXy6Xy8nI9+uijMZslMAEAAADQ50iEI0zvJuvmbrzxRrW2tuqFF17QbbfdFrNs3MsFAwAAADgN\nHJP0tkd6xxO1iM1mk8/nC972+Xyy2+1Ry1999dXq7u7Wxx9/rHPPPTdiGQITAAAAAH2OSbpoeu9x\n0jPLQoqUlJSora1NXq9XeXl5amhoUH19fUiZf//735o4caIMw9Cbb74pSVGDEonABAAAAEB/n5kX\nsVgsqq2tVVlZmQKBgKqqquRyuVRXVyepdwf4v/zlL1q3bp2sVqtGjx6t9evXx6yT5YIBAAAASPoy\nd6QuQnhQbaR9g0VGTAAAAAD0ibA8cCYQmAAAAADoE2F54ExguWAAAAAAfeJYLliS3G63ioqKVFBQ\noJUrVw44//TTT6u4uFiXXHKJrrzySr399tsxm2XEBAAAAECfOEZMAoGAampq1NTUJJvNptLSUlVU\nVMjlcgXLTJw4Udu2bdOYMWPkdrv1wx/+UDt27IhaJyMmAAAAAPocj3CEaW5uVn5+vpxOp6xWqyor\nK9XY2BhS5oorrtCYMWMkSdOmTVN7e3vMZglMAAAAAPT5LMIRxu/3y+FwBG/b7Xb5/f6oVa5Zs0Yz\nZ86M2SxTuQAAAAD0iWNVLsMw4q7uH//4h5544glt3749ZjkCEwAAAAB9jkn6yCMd9EQtYrPZ5PP5\ngrd9Pp/sdvuAcm+//bbmz58vt9uts88+O2azbLAIAAAAQNKXIyHlEcKDzaEbLHZ3d6uwsFBbtmxR\nXl6epk6dqvr6+pDk9/379+vaa6/VU089pcsvv9y0bUZMAAAAAPSJkFMSzmKxqLa2VmVlZQoEAqqq\nqpLL5VJdXZ0kqbq6Wg8++KA++eQTLViwQJJktVrV3NwctU5GTAAAAABI+nLEpCRCePB66IhJOjBi\nAgAAAKBPhOWBM4HABAAAAECfOKZypQOBCQAAAIA+cSwXnA5ssAgAAACgz7EIRwRut1tFRUUqKCjQ\nypUrB5z/17/+pSuuuEIjRozQ7373O9NmGTEBAAAA0OeweZFAIKCamho1NTXJZrOptLRUFRUVIcsF\nn3vuuVq1apU2bNgQV7OMmAAAAADo0x3hCNPc3Kz8/Hw5nU5ZrVZVVlaqsbExpMy4ceNUUlIiq9Ua\nV7MEJgAAAAAS4vf75XA4grftdrv8fn9SdTKVCwAAAEA/J0xLGIaR8lYJTAAAAAD0c1jSK5K2Ry1h\ns9nk8/mCt30+n+x2e1KtEpgAAAAA6OeopMu+PE56OKRESUmJ2tra5PV6lZeXp4aGBtXX10esLd4d\n4wlMAAAAAPTzuWkJi8Wi2tpalZWVKRAIqKqqSi6XS3V1dZKk6upqffjhhyotLdWnn36qYcOG6ZFH\nHlFLS4tGjx4dsU6jJ94QBgAAAMAprTd35I0IZ74a98jH/4oREwAAAAD9mI+YpAOBCQAAAIB+shOY\nsI8JAAAAgH4+j3AM5Ha7VVRUpIKCAq1cuTJimUWLFqmgoEDFxcXatWtXzFYJTAAAAAD082mEI1Qg\nEFBNTY3cbrdaWlpUX1+v1tbWkDKbNm3S3r171dbWpscff1wLFiyI2SqBCQAAAIB+zEdMmpublZ+f\nL6fTKavVqsrKSjU2NoaU2bhxo+bMmSNJmjZtmjo7O9XR0RG1VQITAAAAAP2YByZ+v18OhyN42263\ny+/3m5Zpb2+P2irJ7wAAAAD6GTh1K1zvssLmwpcYjvU4AhMAAAAA/dw34J7wTRFtNpt8Pl/wts/n\nk91uj1mmvb1dNpstaqtM5QIAAAAgqXeEI9Jx+PDhkHIlJSVqa2uT1+tVV1eXGhoaVFFREVKmoqJC\n69atkyTt2LFDY8eOVW5ubtS2GTEBAAAAkBCLxaLa2lqVlZUpEAioqqpKLpdLdXV1kqTq6mrNnDlT\nmzZtUn5+vkaNGqW1a9fGrNPoSffe8gAAAABggqlcAAAAALKOwAQAAABA1hGYAAAAAMg6AhMAAAAA\nWUdgAgAAACDrCEwAAAAAZB2BCQAAAICsIzABAAAAkHX/H9U6yQxiSz/BAAAAAElFTkSuQmCC\n",
"text": [
"<matplotlib.figure.Figure at 0xa81d510>"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[5, 3, 1, 6, 9, 7, 0, 8]\n",
"\n",
"Log-likelihood: -5.186e+02\n",
"\n",
" Name | Value | Constraints | Ties | prior \n",
"--------------------------------------------------------------------\n",
" linear_variance_0 | 27.3240 | (+ve) | | \n",
" linear_variance_1 | 24.5210 | (+ve) | | \n",
" linear_variance_2 | 7.3404 | (+ve) | | \n",
" linear_variance_3 | 5.5106 | (+ve) | | \n",
" linear_variance_4 | 1.2720 | (+ve) | | \n",
" linear_variance_5 | 1.1788 | (+ve) | | \n",
" linear_variance_6 | 0.4906 | (+ve) | | \n",
" linear_variance_7 | 0.0002 | (+ve) | | \n",
" noise_variance | 15.6334 | (+ve) | | \n",
"\n"
]
}
],
"prompt_number": 445
},
{
"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",
" #print \"---\"\n",
" #print kernelx.input_dim\n",
" #print x.shape\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": 446
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from types import MethodType\n",
"\n",
"#HSIChp=False\n",
"\n",
"def testobj(self, inpu):\n",
" global kernelx2\n",
" global kernely2\n",
" try:\n",
" self._set_params(inpu)\n",
" except:\n",
" return 10000\n",
" #print self\n",
" #print inpu\n",
" mresiduals=np.abs(self.predict(self.X)[0] - self.likelihood.Y)\n",
" #print residuals.shape\n",
" #print np.sum(mresiduals)\n",
" #print X.shape\n",
" newHSIC=1.0*BiasedHSIC(mresiduals, X, kernelx2, kernely2)\n",
" #print self\n",
" return newHSIC\n",
"\n",
"\n",
"def HSICreg(HSIChp=False):\n",
" global X\n",
" global Y\n",
"\n",
" kernelx2 = GPy.kern.linear(input_dim=1)\n",
" kernely2 = GPy.kern.linear(input_dim=10)\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",
" X=X.transpose()\n",
" #print X.shape\n",
" #print \"###\"\n",
" #print Y.shape\n",
" k2 = GPy.kern.linear(input_dim=0)\n",
" m2 = GPy.models.GPRegression(X2,Y,k2)\n",
" if HSIChp == True:\n",
" m2.objective_function = MethodType(testobj, m2)\n",
" m2.preferred_optimizer='simplex'\n",
" m2.optimize(messages=False)\n",
" #m['.*variance']=1\n",
" #m['linear_variance']=10000\n",
" bestll=m2.log_likelihood()\n",
" print m2\n",
"\n",
" oldHSIC=10000000\n",
" #print m\n",
" #print bestll\n",
" #provides initial parameters for 0 SNPS\n",
" oldparams=m['.*variance']\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",
" k2 = GPy.kern.linear(input_dim=step+1, ARD=True)\n",
"\n",
" tempm2 = GPy.models.GPRegression(X2,Y,k2)\n",
" if HSIChp == True:\n",
" tempm2.preferred_optimizer='simplex'\n",
" tempm2.objective_function = MethodType(testobj, tempm2)\n",
"\n",
" tempm2.optimize(messages=False)\n",
" residuals=np.abs(tempm2.predict(X2)[0] - Y)\n",
" #print residuals\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",
"\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, ARD=True)\n",
"\n",
"\n",
" tempm2 = GPy.models.GPRegression(X2,Y,k2)\n",
" if HSIChp == True:\n",
" tempm2.preferred_optimizer='simplex'\n",
" tempm2.objective_function = MethodType(testobj, tempm2)\n",
"\n",
" tempm2.optimize(messages=False)\n",
" #residuals=tempm2.predict(X2)[0]\n",
" residuals=np.abs(tempm2.predict(X2)[0] - Y)\n",
" newHSIC=BiasedHSIC(residuals, X, kernelx2, kernely2)\n",
" if newHSIC < oldHSIC:\n",
" #hasbeenadded=True\n",
" oldHSIC=newHSIC\n",
" m2=tempm2\n",
" print \"HSIC: \" + str(newHSIC)\n",
" print \"LogLikelihood: \" + str(tempm2.log_likelihood())\n",
" chosenSNPS2.append(bestSNP)\n",
" else:\n",
" print \"HSIC not improved\"\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",
" return (m2, chosenSNPS2)\n"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 447
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"out=HSICreg(HSIChp=True)\n",
"print out[0], out[1]"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"Log-likelihood: -9.818e+03\n",
"\n",
" Name | Value | Constraints | Ties | prior \n",
"-----------------------------------------------------------------\n",
" linear_variance | 1.0000 | (+ve) | | \n",
" noise_variance | 1.0000 | (+ve) | | \n",
"\n",
"HSIC: 1.68344184072"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"LogLikelihood: -7028.85443473\n",
"HSIC: 0.562717391998"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"LogLikelihood: -5421.45213718\n",
"HSIC: 0.280591930483"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"LogLikelihood: -4933.40407941\n",
"HSIC: 0.195837573383"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"LogLikelihood: -4251.54667043\n",
"HSIC: 0.119006412867"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"LogLikelihood: -4088.05903774\n",
"converged: 5SNPs"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAyYAAACICAYAAAASyf8xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG4BJREFUeJzt3X9wVNX9//HXtbvKzy8oYCS7+X4XzJpsBozRDUjRBm0l\nho4RqqPJdGqKgaZoaq21g9N+nIrTQdNpZ6qm7USrVqxuYm1L0IGVBrtSsJhaUfppUGPHlc2qqaki\nKApk3e8fgd3s7103mw3h+Zg5A/fu+55z98fd5J1zzzlGMBgMCgAAAADy6JR8nwAAAAAAkJgAAAAA\nyDsSEwAAAAB5R2ICAAAAIO9ITAAAAADkHYkJAAAAgLwjMQEAAACQkeuvv14FBQWaP39+wpibbrpJ\ndrtd5eXl2r17d8o6SUwAAAAAZGTlypVyu90JH9+8ebPeeOMN9fb26v7779eaNWtS1kliAgAAACAj\nF198sU4//fSEj2/atEkNDQ2SpIULF2r//v3q7+9PWieJCQAAAIAR5ff7VVRUFNq2Wq3q6+tLegyJ\nCQAAAIARFwwGI7YNw0gaT2ICAAAAQJI00TBkxClTp07NqB6LxSKfzxfa7uvrk8ViSXqM6XOdMQAA\nAIBx51NJP4mz/38++iijempra9Xa2qq6ujrt2rVL06dPV0FBQdJjSEwAAAAAhPyfNGLq6+v13HPP\naWBgQEVFRVq3bp2OHj0qSWpqatKyZcu0efNmFRcXa/LkyXr44YdT1mkEo2/+AgAAAHBSMgxDD8TZ\nv1qxY0ZGGj0mAAAAAEIm5qldEhMAAAAAIencypULzMoFAAAAIGRinBKP2+1WaWmp7Ha7WlpaYh7/\n4IMPtGLFCpWXl2vhwoX617/+lbRdEhMAAAAAIekkJoFAQM3NzXK73erp6ZHL5dLevXsjYtavX6/z\nzz9fr7zyijZs2KDvfve7SdslMQEAAAAQkk5i0t3dreLiYtlsNpnNZtXV1amzszMiZu/evbrkkksk\nSSUlJfJ6vXrvvfcStktiAgAAACBkapwSze/3q6ioKLRttVrl9/sjYsrLy/XHP/5R0lAi89Zbb6mv\nry9huyQmAAAAAEImxSnRDMNIWc9tt92m/fv3q6KiQq2traqoqNAXvvCFhPHMygUAAAAgZKJJ2vGZ\ntHP4siVRS5hYLBb5fL7Qts/nk9VqjYiZOnWqHnroodD2nDlzNHfu3ITtssAiAAAAAElDPSFHpsXu\nP/XDyAUWBwcHVVJSom3btqmwsFALFiyQy+WSw+EIxXz44YeaOHGiTj31VD3wwAPauXOnfvvb3yZs\nmx4TAAAAACHm01LHmEwmtba2qrq6WoFAQI2NjXI4HGpra5MkNTU1qaenR9/85jdlGIbmzZunBx98\nMGmd9JgAAAAAkDTUYxL8v3H274vsMckFekwAAAAAhKXRY5ILJCYAAAAAwibnp1kSEwAAAABhE/LT\nLOuYAAAAAAg7LU6Jw+12q7S0VHa7XS0tLTGPDwwM6PLLL9d5552nefPmJZ2RS2LwOwAAAIBjDMNQ\ncFmc/ZsjB78HAgGVlJSoq6tLFotFlZWVMdMF33HHHTp8+LDuuusuDQwMqKSkRP39/TKZ4t+0RY8J\nAAAAgLAJcUqU7u5uFRcXy2azyWw2q66uTp2dnRExs2fP1oEDByRJBw4c0IwZMxImJRJjTAAAAAAM\nl8asXH6/X0VFRaFtq9WqF154ISJm9erVuvTSS1VYWKiDBw/qiSeeSFonPSYAAAAAwibHKVEMw0hZ\nzfr163Xeeefp7bff1ssvv6wbb7xRBw8eTBhPjwkAAACAsAmSp0/y+BOHWCwW+Xy+0LbP55PVao2I\nef755/WjH/1IknT22Wdrzpw5eu211+R0OuPWSWICAAAAIOw0acnZQ+W4dX+PDHE6nert7ZXX61Vh\nYaE6OjrkcrkiYkpLS9XV1aXFixerv79fr732mubOnZuwWRITAAAAAGFpjDExmUxqbW1VdXW1AoGA\nGhsb5XA41NbWJklqamrSD3/4Q61cuVLl5eX67LPP9NOf/lRnnHFGwjqZLhgAAACApGPTBa+Ps/+H\nkdMF5wI9JgAAAADC0ugxyQUSEwAAAABheUpMmC4YAAAAQFga0wVLktvtVmlpqex2u1paWmIe/9nP\nfqaKigpVVFRo/vz5MplM2r9/f8JmGWMCAAAAQNKxMSauOPvrI8eYBAIBlZSUqKurSxaLRZWVlXK5\nXHI4HHHrffrpp/WLX/xCXV1dCdumxwQAAABA2GlxSpTu7m4VFxfLZrPJbDarrq5OnZ2dCat8/PHH\nVV9fn7TZEUlMUnXjAOOdzWbTueeeq4qKCi1YsECS9P777+uyyy7TOeeco6VLlybtugROZNdff70K\nCgo0f/780L5kn/+77rpLdrtdpaWl2rp1az5OGciZeNfDHXfcIavVGrqlZcuWLaHHuB4wJk2IU6L4\n/X4VFRWFtq1Wq/z++CsyHjp0SM8884yuuuqqpM1mnZgEAgE1NzfL7Xarp6dHLpdLe/fuzbZa4IRi\nGIY8Ho92796t7u5uSdLdd9+tyy67TK+//rq+/OUv6+67787zWQK5sXLlSrnd7oh9iT7/PT096ujo\nUE9Pj9xut2644QZ99tln+ThtICfiXQ+GYeiWW27R7t27tXv3btXU1EjiesAYlsYYE8Mw0q7uqaee\n0kUXXaTp06cnjcs6Mcm0GwcYr6KHa23atEkNDQ2SpIaGBm3cuDEfpwXk3MUXX6zTTz89Yl+iz39n\nZ6fq6+tlNptls9lUXFwcSuaB8SDe9SDFX/+B6wFj1mmS53+lOzaESzSLxSKfzxfa9vl8slqtcatr\nb29PeRuXNAKJSSbdOMB4ZRiGvvKVr8jpdOqBBx6QJPX396ugoECSVFBQoP7+/nyeIjCqEn3+3377\n7YgfXPzMwMnivvvuU3l5uRobG0O3NnI9YMyaIC1ZLN3xnXCJ5nQ61dvbK6/XqyNHjqijo0O1tbUx\ncR9++KG2b9+uK6+8MmWzWScmmXTjAOPVzp07tXv3bm3ZskW//OUv9de//jXiccMwuFZw0kr1+efa\nwHi3Zs0avfnmm3r55Zc1e/Zsff/7308Yy/WAMSGNW7lMJpNaW1tVXV2tsrIyXXvttXI4HGpra1Nb\nW1sobuPGjaqurtbEiRNTNpv1AouZdOMA49Xs2bMlSbNmzdKKFSvU3d2tgoICvfvuuzrrrLP0zjvv\n6Mwzz8zzWQKjJ9HnP/pnRl9fnywWS75OExgVw7//V61apSuuuEIS1wPGsDQXWKypqQmNmTquqakp\nYruhoSF0a28qWfeYpNuNA4xXhw4d0sGDByVJH3/8sbZu3ar58+ertrZWjzzyiCTpkUce0fLly/N5\nmsCoSvT5r62tVXt7u44cOaI333xTvb29oZnsgPHqnXfeCf3/T3/6U2jGLq4HjFlpzMqVC1n3mAzv\nxgkEAmpsbEy4sAowHvX392vFihWSpMHBQX3961/X0qVL5XQ6dc011+jBBx+UzWbTE088keczBXKj\nvr5ezz33nAYGBlRUVKQ777xTt912W9zPf1lZma655hqVlZXJZDLpV7/6FbeuYFyJvh7WrVsnj8ej\nl19+WYZhaM6cOaHbXLgeMFYdTrDSe66x8jsAAAAASUPjnD4cPDVm/zTTkZjZ5dxut26++WYFAgGt\nWrVKa9eujTnO4/Hoe9/7no4ePaqZM2fK4/EkbpvEBAAAAIA0lJi8F5wSs3+W8VFEYhIIBFRSUqKu\nri5ZLBZVVlbK5XJF3Dm1f/9+LV68WM8884ysVqsGBgY0c+bMhG2PyMrvAAAAAMaHwzotpkRLZy3D\nxx9/XFdddVVoYqxkSYmURmLidrtVWloqu92ulpaWTJ4TAAAAgBPMJ5oYU6Kls5Zhb2+v3n//fV1y\nySVyOp169NFHk7abdPB7IBBQc3NzRBdNbW0tg9sBAACAcSpeD0m0dCZqOHr0qF566SVt27ZNhw4d\n0qJFi3ThhRfKbrfHjU+amAzvopEU6qIhMQEAAADGpyM6VS96PtaLnkMJY9JZy7CoqEgzZ87UxIkT\nNXHiRH3pS1/SK6+88vkSk3hdNC+88EJEjM0w9FaySgAAAICT1v9TMOjN90lk5JAmqWzJJJUtCe+7\nf91ARMzwtQwLCwvV0dEhl8sVEXPllVequblZgUBAhw8f1gsvvKBbbrklYbtJExPDMLR9+3YVFBTo\nzDPPjDsF2FuSfizJI2mJpHX6cdInGu3HWpdRfLL6o+tKdS6jHR8t+vhk8aliR/t1T9V+po+Ppmxf\ny3Sfm0fxr4nhx4/06zbS71O0TNsfrbrSqT/b1zLT+j9v7EjI9Wv7edr2KPX1oDQej5brz8lI1pXv\n6zvT+nLZdqayrX8kP/Mj85n0aOiKGN3vonSOz6btkT4+l/Xl+ns4urUTzRHFThccLdFahsfX6Wlq\nalJpaakuv/xynXvuuTrllFO0evVqlZWVJa4zWYMWi0XTp09Xe3u7rrvuurhdNNLQ5eU99u/Q/2wp\nnwwAAAAw/niPlRNXOmNMJKmmpkY1NTUR+5qamiK2b731Vt16661p1Zd0HZPBwUGVlJTokUce0be/\n/W2ZTKaY+YmHBr6E+0xynfln8tfmkTbSf5lLVX82zyeffx1O5/hMjLXemnRfK4+O/z1s5Ov+vPHZ\ntpfq+Gw+JyP9GcvX+57Osfl+H0fzr5iZ9CBmWney9hIdn8vXIt/fo7nsgcl1L3u0XH8usqlvJL5r\nPAr/fMh1T3cq+bxbY6R7wkay/pG9ntfFLEw4lhmGoc7g0pj9Vxpbc/48kk4XfLyL5rrrrlNvb6+u\nvfbaJAPfbSN/dsAJzJbvEwDGEFu+TwAYQ2z5PgEghUOaFFNGQ8p1TObNm6eCggIZhqHHHntM9957\nb4JI28ieGXCCs+X7BIAxxJbvEwDGEFu+TwBI4bBOjSnxpFrv0OPxaNq0aaqoqFBFRYV+8pOfJG03\n6RgTSTKbzbr99tu1du1a/e1vf9MFF1ygyy67LKrnxBP6n1dccAAAADhZeXWijzE5ksYYk3TXO6yq\nqtKmTZvSajdlYnLWWWeFRs9PmTJFDodDb7/9dkSjP9ZzCY/P9P6/VPf0JXt8rM9glE18pueW63uR\no+VydrVs5XP2lkzry/T6yPVsT6M5e1O2bY32a5PKaL7v2X6vRsv0+yGfs+xlK5vvh2y/d1Ody2iO\n5cu0vlx/pnI9BiWXYw8ylevrK5Pvh9Ee23MijSnJ7lwS/548Vh2Ks9J7tHTXO8xkXErKW7nq6+v1\nxS9+Ua+//rpmz56tHTt2aOHChWk3AAAAAODEcUSnxZRo8dY79Pv9ETGGYej5559XeXm5li1bpp6e\nnqTtpuwxefjhh1VVVaWZM2eqt7dXy5Yt05QpU9J9XgAAAABOIInGlAw3NDNvcueff758Pp8mTZqk\nLVu2aPny5Xr99dcTxqdMTCZMmKCtW7fqmmuuUUNDg37/+99rx44duuiii0IxnmHxNjHGBAAAACcn\nr6TI345PPJ9okt7yePWW562EMRaLRT6fL7Qdb73DqVOnhv5fU1OjG264Qe+//77OOOOMuHWmTEyC\nwaC+853vqKysTGvWrFF7e3tMZUtSVQIAAACcBGySIn87PvHGmBzWqTpryTk6a8k5oX071m2PiHE6\nnert7ZXX61VhYaE6OjrkcrkiYvr7+3XmmWfKMAx1d3crGAwmTEqkNBKTnTt36tFHH9WECRN07733\naubMmdq3b1/C5eSzGbyeqVwPGM12UFW2g/eSLSaZa9kOAkslk0Ge2Q4IHemF9rIdaJzJ+zrSg1lH\neiGubN73kZ6sIlV8snPJxfEjOXA41fGZfrdFyzQ+k/ditBdoG82fQSP93MbaBA7J2s71gonZ/nxP\nJZvv1lxPrJNpfSM9YD2Z0V6oNpeTFGX6OkbL7Gf9iSedWbmOr3dYXV2tQCCgxsZGORwOtbW1SRpa\nAf7JJ5/Ur3/9a5lMJk2aNEnt7e3J60zV6EUXXaTBwUE5nU4VFBRo//79mjBhQppPCwAAAMCJJJ0x\nJtLQ7Vk1NTUR+5qamkL/v/HGG3XjjTem3W7KxESS7rnnHpWVlengwYP66le/qhdffFFLliwJPe6J\niPaKUSYAAAA4GXl1oq9iMjTGJB9SThe8Z88ebdq0SatWrdLg4KD+/Oc/q6KiIiJmybBCUgIAAICT\nlU3RvxufeNJd+X2kGcEUq54sXbpU+/bt0+DgoP7zn//o9ttv1w9+8INwBYYh6cc63lMy0vdejqZc\nn0s292aO5mKOn+fxaCP5Po/24o+ppPtaeDX05ZTN/fipjPa4jEyN5qJlo70wWD4Xn8z1/f2Znk86\nx3oV/89W+f5uy/R8MqkrVd25vKd9JNrLpK58L0icymher6msG/Y7U7zjc/3cc/k5yGXb8eqLls31\nncuxQ+uU2SKD+WYYhpYHXTH7Nxr1Mc/D7Xbr5ptvViAQ0KpVq7R27dq4df7973/XokWL9MQTT+hr\nX/tawraT9pg8/fTTstvtevXVV/Wb3/xGVVVVEUlJJG+yqoCTjjffJwCMId58nwAwpnjzfQJAUoc0\nKaZECwQCam5ultvtVk9Pj1wul/bu3Rs3bu3atbr88stTJmhJE5Pnn39emzZt0pw5c1RfX69nn31W\n1113XZxIj4YuMg+XGgAAAE5iXnmkUDkRHdGpMSVad3e3iouLZbPZZDabVVdXp87Ozpi4++67T1df\nfbVmzZqVst2kicn69evl8/n05ptvqr29XZdeeqk2bNgQJ3KJjt9RZ0vZJAAAADBe2cbBGJPTYko0\nv9+voqKi0LbVapXf74+J6ezs1Jo1aySlsVp8ME0ejyd4xRVXxOyvqqoKSqJQKBQKhUKhUChRpaqq\nKt1ft8cEScELgjtiihSZNjz55JPBVatWhbYfffTRYHNzc0TM1VdfHdy1a1cwGAwGGxoagk8++WTS\nttOaLliSqqqqVFVVFbPf4/GkWwUAAACAMe6QJuljz4s65HkxYYzFYpHP5wtt+3w+Wa3WiJh//OMf\nqqurkyQNDAxoy5YtMpvNqq2tjVtnylm5AAAAAJwcDMPQ2cH/jdn/b2NexOD1wcFBlZSUaNu2bSos\nLNSCBQvkcrnkcDji1rty5UpdccUVSWflSrvHBAAAAMD4F29MSTSTyaTW1lZVV1crEAiosbFRDodD\nbW1tkiJXgE8XPSYAAAAAJA31mMwI9MXs/+8XrDlfj4UeEwAAAAAhhz9N3WOSCyQmAAAAAEKOfBq7\nbsloSLqOCQAAAICTy5GPJsWUeNxut0pLS2W329XS0hLzeGdnp8rLy1VRUaELLrhAzz77bNJ2GWMC\nAAAAQNKxRRBfi5MelBgRY0wCgYBKSkrU1dUli8WiysrKmFm5Pv74Y02ePFmS9M9//lMrVqzQG2+8\nkbBtekwAAAAAhH0ap0Tp7u5WcXGxbDabzGaz6urq1NnZGRFzPCmRpI8++kgzZ85M2iyJCQAAAICw\nT+KUKH6/X0VFRaFtq9Uqv98fE7dx40Y5HA7V1NTo3nvvTdosiQkAAACAsI/jlCiGYaRV1fLly7V3\n71499dRT+sY3vpE0llm5AAAAAIR9KmmPR/qnJ2GIxWKRz+cLbft8Plmt1oTxF198sQYHB/Xf//5X\nM2bMiBtDYgIAAAAg7FNJ5ywZKsc9vi4ixOl0qre3V16vV4WFhero6JDL5YqI+fe//625c+fKMAy9\n9NJLkpQwKZFITAAAAAAM91HqEJPJpNbWVlVXVysQCKixsVEOh0NtbW2SpKamJv3hD3/Qhg0bZDab\nNWXKFLW3tyetk+mCAQAAAEg6NnakLU560BQ5XXAu0GMCAAAAICzO9MCjgcQEAAAAQFic6YFHA9MF\nAwAAAAhLY7pgSXK73SotLZXdbldLS0vM44899pjKy8t17rnnavHixdqzZ0/SZukxAQAAABCWRo9J\nIBBQc3Ozurq6ZLFYVFlZqdraWjkcjlDM3LlztX37dk2bNk1ut1vf+ta3tGvXroR10mMCAAAAIOxw\nnBKlu7tbxcXFstlsMpvNqqurU2dnZ0TMokWLNG3aNEnSwoUL1dfXl7RZEhMAAAAAYR/FKVH8fr+K\niopC21arVX6/P2GVDz74oJYtW5a0WW7lAgAAABCWxqxchmGkXd1f/vIXPfTQQ9q5c2fSOBITAAAA\nAGGfSnrPIw14EoZYLBb5fL7Qts/nk9VqjYnbs2ePVq9eLbfbrdNPPz1psyywCAAAAEDSsZ6Qmjjp\nwZbIBRYHBwdVUlKibdu2qbCwUAsWLJDL5YoY/L5v3z5deuml+t3vfqcLL7wwZdv0mAAAAAAIizOm\nJJrJZFJra6uqq6sVCATU2Ngoh8OhtrY2SVJTU5PuvPNOffDBB1qzZo0kyWw2q7u7O2Gd9JgAAAAA\nkHSsx8QZJz14MbLHJBfoMQEAAAAQFmd64NFAYgIAAAAgLI1buXKBxAQAAABAWBrTBecCCywCAAAA\nCPs0TonD7XartLRUdrtdLS0tMY+/+uqrWrRokSZMmKCf//znKZulxwQAAABA2MHUIYFAQM3Nzerq\n6pLFYlFlZaVqa2sjpgueMWOG7rvvPm3cuDGtZukxAQAAABA2GKdE6e7uVnFxsWw2m8xms+rq6tTZ\n2RkRM2vWLDmdTpnN5rSaJTEBAAAAkBG/36+ioqLQttVqld/vz6pObuUCAAAAMMzRlBGGYYx4qyQm\nAAAAAIY5KGmHpJ0JIywWi3w+X2jb5/PJarVm1SqJCQAAAIBhDkk6/1g57qcREU6nU729vfJ6vSos\nLFRHR4dcLlfc2tJdMZ7EBAAAAMAwn6SMMJlMam1tVXV1tQKBgBobG+VwONTW1iZJampq0rvvvqvK\nykodOHBAp5xyiu655x719PRoypQpces0gummMAAAAADGtaGxI/+I88gFafd8fF70mAAAAAAYJnWP\nSS6QmAAAAAAYJj+JCeuYAAAAABjmkzglltvtVmlpqex2u1paWuLG3HTTTbLb7SovL9fu3buTtkpi\nAgAAAGCYA3FKpEAgoObmZrndbvX09Mjlcmnv3r0RMZs3b9Ybb7yh3t5e3X///VqzZk3SVklMAAAA\nAAyTuseku7tbxcXFstlsMpvNqqurU2dnZ0TMpk2b1NDQIElauHCh9u/fr/7+/oStkpgAAAAAGCZ1\nYuL3+1VUVBTatlqt8vv9KWP6+voStsrgdwAAAADDxN66FW1oWuHUoqcYTnYciQkAAACAYf4nZk/0\noogWi0U+ny+07fP5ZLVak8b09fXJYrEkbJVbuQAAAABIGurhiFcOHjwYEed0OtXb2yuv16sjR46o\no6NDtbW1ETG1tbXasGGDJGnXrl2aPn26CgoKErZNjwkAAACAjJhMJrW2tqq6ulqBQECNjY1yOBxq\na2uTJDU1NWnZsmXavHmziouLNXnyZD388MNJ6zSCuV5bHgAAAABS4FYuAAAAAHlHYgIAAAAg70hM\nAAAAAOQdiQkAAACAvCMxAQAAAJB3JCYAAAAA8o7EBAAAAEDekZgAAAAAyLv/D+XtJx0HL6lyAAAA\nAElFTkSuQmCC\n",
"text": [
"<matplotlib.figure.Figure at 0xb360290>"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"Log-likelihood: -4.088e+03\n",
"\n",
" Name | Value | Constraints | Ties | prior \n",
"-------------------------------------------------------------------\n",
" linear_variance_0 | 1.4349 | (+ve) | | \n",
" linear_variance_1 | 0.6931 | (+ve) | | \n",
" linear_variance_2 | 2.7227 | (+ve) | | \n",
" linear_variance_3 | 0.7135 | (+ve) | | \n",
" linear_variance_4 | 0.7423 | (+ve) | | \n",
" noise_variance | 0.6932 | (+ve) | | \n",
" [1, 2, 9, 6, 8]\n"
]
}
],
"prompt_number": 449
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"out=HSICreg(HSIChp=False)\n",
"print out[0], out[1]"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"Log-likelihood: -6.827e+02\n",
"\n",
" Name | Value | Constraints | Ties | prior \n",
"-------------------------------------------------------------------\n",
" linear_variance | 1.0000 | (+ve) | | \n",
" noise_variance | 105.9333 | (+ve) | | \n",
"\n",
"HSIC: 1.68352091142"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"LogLikelihood: -621.226239734\n",
"HSIC: 0.565432510184"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"LogLikelihood: -599.253902186\n",
"HSIC: 0.291092340396"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"LogLikelihood: -592.081614024\n",
"HSIC: 0.202822719409"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"LogLikelihood: -579.821823892\n",
"HSIC: 0.126942040002"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"LogLikelihood: -577.294289371\n",
"converged: 5SNPs"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAyYAAACICAYAAAASyf8xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG4BJREFUeJzt3X9wVNX9//HXtbvKzy8oYCS7+X4XzJpsBozRDUjRBm0l\nho4RqqPJdGqKgaZoaq21g9N+nIrTQdNpZ6qm7USrVqxuYm1L0IGVBrtSsJhaUfppUGPHlc2qqaki\nKApk3e8fgd3s7103mw3h+Zg5A/fu+55z98fd5J1zzzlGMBgMCgAAAADy6JR8nwAAAAAAkJgAAAAA\nyDsSEwAAAAB5R2ICAAAAIO9ITAAAAADkHYkJAAAAgLwjMQEAAACQkeuvv14FBQWaP39+wpibbrpJ\ndrtd5eXl2r17d8o6SUwAAAAAZGTlypVyu90JH9+8ebPeeOMN9fb26v7779eaNWtS1kliAgAAACAj\nF198sU4//fSEj2/atEkNDQ2SpIULF2r//v3q7+9PWieJCQAAAIAR5ff7VVRUFNq2Wq3q6+tLegyJ\nCQAAAIARFwwGI7YNw0gaT2ICAAAAQJI00TBkxClTp07NqB6LxSKfzxfa7uvrk8ViSXqM6XOdMQAA\nAIBx51NJP4mz/38++iijempra9Xa2qq6ujrt2rVL06dPV0FBQdJjSEwAAAAAhPyfNGLq6+v13HPP\naWBgQEVFRVq3bp2OHj0qSWpqatKyZcu0efNmFRcXa/LkyXr44YdT1mkEo2/+AgAAAHBSMgxDD8TZ\nv1qxY0ZGGj0mAAAAAEIm5qldEhMAAAAAIencypULzMoFAAAAIGRinBKP2+1WaWmp7Ha7WlpaYh7/\n4IMPtGLFCpWXl2vhwoX617/+lbRdEhMAAAAAIekkJoFAQM3NzXK73erp6ZHL5dLevXsjYtavX6/z\nzz9fr7zyijZs2KDvfve7SdslMQEAAAAQkk5i0t3dreLiYtlsNpnNZtXV1amzszMiZu/evbrkkksk\nSSUlJfJ6vXrvvfcStktiAgAAACBkapwSze/3q6ioKLRttVrl9/sjYsrLy/XHP/5R0lAi89Zbb6mv\nry9huyQmAAAAAEImxSnRDMNIWc9tt92m/fv3q6KiQq2traqoqNAXvvCFhPHMygUAAAAgZKJJ2vGZ\ntHP4siVRS5hYLBb5fL7Qts/nk9VqjYiZOnWqHnroodD2nDlzNHfu3ITtssAiAAAAAElDPSFHpsXu\nP/XDyAUWBwcHVVJSom3btqmwsFALFiyQy+WSw+EIxXz44YeaOHGiTj31VD3wwAPauXOnfvvb3yZs\nmx4TAAAAACHm01LHmEwmtba2qrq6WoFAQI2NjXI4HGpra5MkNTU1qaenR9/85jdlGIbmzZunBx98\nMGmd9JgAAAAAkDTUYxL8v3H274vsMckFekwAAAAAhKXRY5ILJCYAAAAAwibnp1kSEwAAAABhE/LT\nLOuYAAAAAAg7LU6Jw+12q7S0VHa7XS0tLTGPDwwM6PLLL9d5552nefPmJZ2RS2LwOwAAAIBjDMNQ\ncFmc/ZsjB78HAgGVlJSoq6tLFotFlZWVMdMF33HHHTp8+LDuuusuDQwMqKSkRP39/TKZ4t+0RY8J\nAAAAgLAJcUqU7u5uFRcXy2azyWw2q66uTp2dnRExs2fP1oEDByRJBw4c0IwZMxImJRJjTAAAAAAM\nl8asXH6/X0VFRaFtq9WqF154ISJm9erVuvTSS1VYWKiDBw/qiSeeSFonPSYAAAAAwibHKVEMw0hZ\nzfr163Xeeefp7bff1ssvv6wbb7xRBw8eTBhPjwkAAACAsAmSp0/y+BOHWCwW+Xy+0LbP55PVao2I\nef755/WjH/1IknT22Wdrzpw5eu211+R0OuPWSWICAAAAIOw0acnZQ+W4dX+PDHE6nert7ZXX61Vh\nYaE6OjrkcrkiYkpLS9XV1aXFixerv79fr732mubOnZuwWRITAAAAAGFpjDExmUxqbW1VdXW1AoGA\nGhsb5XA41NbWJklqamrSD3/4Q61cuVLl5eX67LPP9NOf/lRnnHFGwjqZLhgAAACApGPTBa+Ps/+H\nkdMF5wI9JgAAAADC0ugxyQUSEwAAAABheUpMmC4YAAAAQFga0wVLktvtVmlpqex2u1paWmIe/9nP\nfqaKigpVVFRo/vz5MplM2r9/f8JmGWMCAAAAQNKxMSauOPvrI8eYBAIBlZSUqKurSxaLRZWVlXK5\nXHI4HHHrffrpp/WLX/xCXV1dCdumxwQAAABA2GlxSpTu7m4VFxfLZrPJbDarrq5OnZ2dCat8/PHH\nVV9fn7TZEUlMUnXjAOOdzWbTueeeq4qKCi1YsECS9P777+uyyy7TOeeco6VLlybtugROZNdff70K\nCgo0f/780L5kn/+77rpLdrtdpaWl2rp1az5OGciZeNfDHXfcIavVGrqlZcuWLaHHuB4wJk2IU6L4\n/X4VFRWFtq1Wq/z++CsyHjp0SM8884yuuuqqpM1mnZgEAgE1NzfL7Xarp6dHLpdLe/fuzbZa4IRi\nGIY8Ho92796t7u5uSdLdd9+tyy67TK+//rq+/OUv6+67787zWQK5sXLlSrnd7oh9iT7/PT096ujo\nUE9Pj9xut2644QZ99tln+ThtICfiXQ+GYeiWW27R7t27tXv3btXU1EjiesAYlsYYE8Mw0q7uqaee\n0kUXXaTp06cnjcs6Mcm0GwcYr6KHa23atEkNDQ2SpIaGBm3cuDEfpwXk3MUXX6zTTz89Yl+iz39n\nZ6fq6+tlNptls9lUXFwcSuaB8SDe9SDFX/+B6wFj1mmS53+lOzaESzSLxSKfzxfa9vl8slqtcatr\nb29PeRuXNAKJSSbdOMB4ZRiGvvKVr8jpdOqBBx6QJPX396ugoECSVFBQoP7+/nyeIjCqEn3+3377\n7YgfXPzMwMnivvvuU3l5uRobG0O3NnI9YMyaIC1ZLN3xnXCJ5nQ61dvbK6/XqyNHjqijo0O1tbUx\ncR9++KG2b9+uK6+8MmWzWScmmXTjAOPVzp07tXv3bm3ZskW//OUv9de//jXiccMwuFZw0kr1+efa\nwHi3Zs0avfnmm3r55Zc1e/Zsff/7308Yy/WAMSGNW7lMJpNaW1tVXV2tsrIyXXvttXI4HGpra1Nb\nW1sobuPGjaqurtbEiRNTNpv1AouZdOMA49Xs2bMlSbNmzdKKFSvU3d2tgoICvfvuuzrrrLP0zjvv\n6Mwzz8zzWQKjJ9HnP/pnRl9fnywWS75OExgVw7//V61apSuuuEIS1wPGsDQXWKypqQmNmTquqakp\nYruhoSF0a28qWfeYpNuNA4xXhw4d0sGDByVJH3/8sbZu3ar58+ertrZWjzzyiCTpkUce0fLly/N5\nmsCoSvT5r62tVXt7u44cOaI333xTvb29oZnsgPHqnXfeCf3/T3/6U2jGLq4HjFlpzMqVC1n3mAzv\nxgkEAmpsbEy4sAowHvX392vFihWSpMHBQX3961/X0qVL5XQ6dc011+jBBx+UzWbTE088keczBXKj\nvr5ezz33nAYGBlRUVKQ777xTt912W9zPf1lZma655hqVlZXJZDLpV7/6FbeuYFyJvh7WrVsnj8ej\nl19+WYZhaM6cOaHbXLgeMFYdTrDSe66x8jsAAAAASUPjnD4cPDVm/zTTkZjZ5dxut26++WYFAgGt\nWrVKa9eujTnO4/Hoe9/7no4ePaqZM2fK4/EkbpvEBAAAAIA0lJi8F5wSs3+W8VFEYhIIBFRSUqKu\nri5ZLBZVVlbK5XJF3Dm1f/9+LV68WM8884ysVqsGBgY0c+bMhG2PyMrvAAAAAMaHwzotpkRLZy3D\nxx9/XFdddVVoYqxkSYmURmLidrtVWloqu92ulpaWTJ4TAAAAgBPMJ5oYU6Kls5Zhb2+v3n//fV1y\nySVyOp169NFHk7abdPB7IBBQc3NzRBdNbW0tg9sBAACAcSpeD0m0dCZqOHr0qF566SVt27ZNhw4d\n0qJFi3ThhRfKbrfHjU+amAzvopEU6qIhMQEAAADGpyM6VS96PtaLnkMJY9JZy7CoqEgzZ87UxIkT\nNXHiRH3pS1/SK6+88vkSk3hdNC+88EJEjM0w9FaySgAAAICT1v9TMOjN90lk5JAmqWzJJJUtCe+7\nf91ARMzwtQwLCwvV0dEhl8sVEXPllVequblZgUBAhw8f1gsvvKBbbrklYbtJExPDMLR9+3YVFBTo\nzDPPjDsF2FuSfizJI2mJpHX6cdInGu3HWpdRfLL6o+tKdS6jHR8t+vhk8aliR/t1T9V+po+Ppmxf\ny3Sfm0fxr4nhx4/06zbS71O0TNsfrbrSqT/b1zLT+j9v7EjI9Wv7edr2KPX1oDQej5brz8lI1pXv\n6zvT+nLZdqayrX8kP/Mj85n0aOiKGN3vonSOz6btkT4+l/Xl+ns4urUTzRHFThccLdFahsfX6Wlq\nalJpaakuv/xynXvuuTrllFO0evVqlZWVJa4zWYMWi0XTp09Xe3u7rrvuurhdNNLQ5eU99u/Q/2wp\nnwwAAAAw/niPlRNXOmNMJKmmpkY1NTUR+5qamiK2b731Vt16661p1Zd0HZPBwUGVlJTokUce0be/\n/W2ZTKaY+YmHBr6E+0xynfln8tfmkTbSf5lLVX82zyeffx1O5/hMjLXemnRfK4+O/z1s5Ov+vPHZ\ntpfq+Gw+JyP9GcvX+57Osfl+H0fzr5iZ9CBmWney9hIdn8vXIt/fo7nsgcl1L3u0XH8usqlvJL5r\nPAr/fMh1T3cq+bxbY6R7wkay/pG9ntfFLEw4lhmGoc7g0pj9Vxpbc/48kk4XfLyL5rrrrlNvb6+u\nvfbaJAPfbSN/dsAJzJbvEwDGEFu+TwAYQ2z5PgEghUOaFFNGQ8p1TObNm6eCggIZhqHHHntM9957\nb4JI28ieGXCCs+X7BIAxxJbvEwDGEFu+TwBI4bBOjSnxpFrv0OPxaNq0aaqoqFBFRYV+8pOfJG03\n6RgTSTKbzbr99tu1du1a/e1vf9MFF1ygyy67LKrnxBP6n1dccAAAADhZeXWijzE5ksYYk3TXO6yq\nqtKmTZvSajdlYnLWWWeFRs9PmTJFDodDb7/9dkSjP9ZzCY/P9P6/VPf0JXt8rM9glE18pueW63uR\no+VydrVs5XP2lkzry/T6yPVsT6M5e1O2bY32a5PKaL7v2X6vRsv0+yGfs+xlK5vvh2y/d1Ody2iO\n5cu0vlx/pnI9BiWXYw8ylevrK5Pvh9Ee23MijSnJ7lwS/548Vh2Ks9J7tHTXO8xkXErKW7nq6+v1\nxS9+Ua+//rpmz56tHTt2aOHChWk3AAAAAODEcUSnxZRo8dY79Pv9ETGGYej5559XeXm5li1bpp6e\nnqTtpuwxefjhh1VVVaWZM2eqt7dXy5Yt05QpU9J9XgAAAABOIInGlAw3NDNvcueff758Pp8mTZqk\nLVu2aPny5Xr99dcTxqdMTCZMmKCtW7fqmmuuUUNDg37/+99rx44duuiii0IxnmHxNjHGBAAAACcn\nr6TI345PPJ9okt7yePWW562EMRaLRT6fL7Qdb73DqVOnhv5fU1OjG264Qe+//77OOOOMuHWmTEyC\nwaC+853vqKysTGvWrFF7e3tMZUtSVQIAAACcBGySIn87PvHGmBzWqTpryTk6a8k5oX071m2PiHE6\nnert7ZXX61VhYaE6OjrkcrkiYvr7+3XmmWfKMAx1d3crGAwmTEqkNBKTnTt36tFHH9WECRN07733\naubMmdq3b1/C5eSzGbyeqVwPGM12UFW2g/eSLSaZa9kOAkslk0Ge2Q4IHemF9rIdaJzJ+zrSg1lH\neiGubN73kZ6sIlV8snPJxfEjOXA41fGZfrdFyzQ+k/ditBdoG82fQSP93MbaBA7J2s71gonZ/nxP\nJZvv1lxPrJNpfSM9YD2Z0V6oNpeTFGX6OkbL7Gf9iSedWbmOr3dYXV2tQCCgxsZGORwOtbW1SRpa\nAf7JJ5/Ur3/9a5lMJk2aNEnt7e3J60zV6EUXXaTBwUE5nU4VFBRo//79mjBhQppPCwAAAMCJJJ0x\nJtLQ7Vk1NTUR+5qamkL/v/HGG3XjjTem3W7KxESS7rnnHpWVlengwYP66le/qhdffFFLliwJPe6J\niPaKUSYAAAA4GXl1oq9iMjTGJB9SThe8Z88ebdq0SatWrdLg4KD+/Oc/q6KiIiJmybBCUgIAAICT\nlU3RvxufeNJd+X2kGcEUq54sXbpU+/bt0+DgoP7zn//o9ttv1w9+8INwBYYh6cc63lMy0vdejqZc\nn0s292aO5mKOn+fxaCP5Po/24o+ppPtaeDX05ZTN/fipjPa4jEyN5qJlo70wWD4Xn8z1/f2Znk86\nx3oV/89W+f5uy/R8MqkrVd25vKd9JNrLpK58L0icymher6msG/Y7U7zjc/3cc/k5yGXb8eqLls31\nncuxQ+uU2SKD+WYYhpYHXTH7Nxr1Mc/D7Xbr5ptvViAQ0KpVq7R27dq4df7973/XokWL9MQTT+hr\nX/tawraT9pg8/fTTstvtevXVV/Wb3/xGVVVVEUlJJG+yqoCTjjffJwCMId58nwAwpnjzfQJAUoc0\nKaZECwQCam5ultvtVk9Pj1wul/bu3Rs3bu3atbr88stTJmhJE5Pnn39emzZt0pw5c1RfX69nn31W\n1113XZxIj4YuMg+XGgAAAE5iXnmkUDkRHdGpMSVad3e3iouLZbPZZDabVVdXp87Ozpi4++67T1df\nfbVmzZqVst2kicn69evl8/n05ptvqr29XZdeeqk2bNgQJ3KJjt9RZ0vZJAAAADBe2cbBGJPTYko0\nv9+voqKi0LbVapXf74+J6ezs1Jo1aySlsVp8ME0ejyd4xRVXxOyvqqoKSqJQKBQKhUKhUChRpaqq\nKt1ft8cEScELgjtiihSZNjz55JPBVatWhbYfffTRYHNzc0TM1VdfHdy1a1cwGAwGGxoagk8++WTS\nttOaLliSqqqqVFVVFbPf4/GkWwUAAACAMe6QJuljz4s65HkxYYzFYpHP5wtt+3w+Wa3WiJh//OMf\nqqurkyQNDAxoy5YtMpvNqq2tjVtnylm5AAAAAJwcDMPQ2cH/jdn/b2NexOD1wcFBlZSUaNu2bSos\nLNSCBQvkcrnkcDji1rty5UpdccUVSWflSrvHBAAAAMD4F29MSTSTyaTW1lZVV1crEAiosbFRDodD\nbW1tkiJXgE8XPSYAAAAAJA31mMwI9MXs/+8XrDlfj4UeEwAAAAAhhz9N3WOSCyQmAAAAAEKOfBq7\nbsloSLqOCQAAAICTy5GPJsWUeNxut0pLS2W329XS0hLzeGdnp8rLy1VRUaELLrhAzz77bNJ2GWMC\nAAAAQNKxRRBfi5MelBgRY0wCgYBKSkrU1dUli8WiysrKmFm5Pv74Y02ePFmS9M9//lMrVqzQG2+8\nkbBtekwAAAAAhH0ap0Tp7u5WcXGxbDabzGaz6urq1NnZGRFzPCmRpI8++kgzZ85M2iyJCQAAAICw\nT+KUKH6/X0VFRaFtq9Uqv98fE7dx40Y5HA7V1NTo3nvvTdosiQkAAACAsI/jlCiGYaRV1fLly7V3\n71499dRT+sY3vpE0llm5AAAAAIR9KmmPR/qnJ2GIxWKRz+cLbft8Plmt1oTxF198sQYHB/Xf//5X\nM2bMiBtDYgIAAAAg7FNJ5ywZKsc9vi4ixOl0qre3V16vV4WFhero6JDL5YqI+fe//625c+fKMAy9\n9NJLkpQwKZFITAAAAAAM91HqEJPJpNbWVlVXVysQCKixsVEOh0NtbW2SpKamJv3hD3/Qhg0bZDab\nNWXKFLW3tyetk+mCAQAAAEg6NnakLU560BQ5XXAu0GMCAAAAICzO9MCjgcQEAAAAQFic6YFHA9MF\nAwAAAAhLY7pgSXK73SotLZXdbldLS0vM44899pjKy8t17rnnavHixdqzZ0/SZukxAQAAABCWRo9J\nIBBQc3Ozurq6ZLFYVFlZqdraWjkcjlDM3LlztX37dk2bNk1ut1vf+ta3tGvXroR10mMCAAAAIOxw\nnBKlu7tbxcXFstlsMpvNqqurU2dnZ0TMokWLNG3aNEnSwoUL1dfXl7RZEhMAAAAAYR/FKVH8fr+K\niopC21arVX6/P2GVDz74oJYtW5a0WW7lAgAAABCWxqxchmGkXd1f/vIXPfTQQ9q5c2fSOBITAAAA\nAGGfSnrPIw14EoZYLBb5fL7Qts/nk9VqjYnbs2ePVq9eLbfbrdNPPz1psyywCAAAAEDSsZ6Qmjjp\nwZbIBRYHBwdVUlKibdu2qbCwUAsWLJDL5YoY/L5v3z5deuml+t3vfqcLL7wwZdv0mAAAAAAIizOm\nJJrJZFJra6uqq6sVCATU2Ngoh8OhtrY2SVJTU5PuvPNOffDBB1qzZo0kyWw2q7u7O2Gd9JgAAAAA\nkHSsx8QZJz14MbLHJBfoMQEAAAAQFmd64NFAYgIAAAAgLI1buXKBxAQAAABAWBrTBecCCywCAAAA\nCPs0TonD7XartLRUdrtdLS0tMY+/+uqrWrRokSZMmKCf//znKZulxwQAAABA2MHUIYFAQM3Nzerq\n6pLFYlFlZaVqa2sjpgueMWOG7rvvPm3cuDGtZukxAQAAABA2GKdE6e7uVnFxsWw2m8xms+rq6tTZ\n2RkRM2vWLDmdTpnN5rSaJTEBAAAAkBG/36+ioqLQttVqld/vz6pObuUCAAAAMMzRlBGGYYx4qyQm\nAAAAAIY5KGmHpJ0JIywWi3w+X2jb5/PJarVm1SqJCQAAAIBhDkk6/1g57qcREU6nU729vfJ6vSos\nLFRHR4dcLlfc2tJdMZ7EBAAAAMAwn6SMMJlMam1tVXV1tQKBgBobG+VwONTW1iZJampq0rvvvqvK\nykodOHBAp5xyiu655x719PRoypQpces0gummMAAAAADGtaGxI/+I88gFafd8fF70mAAAAAAYJnWP\nSS6QmAAAAAAYJj+JCeuYAAAAABjmkzglltvtVmlpqex2u1paWuLG3HTTTbLb7SovL9fu3buTtkpi\nAgAAAGCYA3FKpEAgoObmZrndbvX09Mjlcmnv3r0RMZs3b9Ybb7yh3t5e3X///VqzZk3SVklMAAAA\nAAyTuseku7tbxcXFstlsMpvNqqurU2dnZ0TMpk2b1NDQIElauHCh9u/fr/7+/oStkpgAAAAAGCZ1\nYuL3+1VUVBTatlqt8vv9KWP6+voStsrgdwAAAADDxN66FW1oWuHUoqcYTnYciQkAAACAYf4nZk/0\noogWi0U+ny+07fP5ZLVak8b09fXJYrEkbJVbuQAAAABIGurhiFcOHjwYEed0OtXb2yuv16sjR46o\no6NDtbW1ETG1tbXasGGDJGnXrl2aPn26CgoKErZNjwkAAACAjJhMJrW2tqq6ulqBQECNjY1yOBxq\na2uTJDU1NWnZsmXavHmziouLNXnyZD388MNJ6zSCuV5bHgAAAABS4FYuAAAAAHlHYgIAAAAg70hM\nAAAAAOQdiQkAAACAvCMxAQAAAJB3JCYAAAAA8o7EBAAAAEDekZgAAAAAyLv/D+XtJx0HL6lyAAAA\nAElFTkSuQmCC\n",
"text": [
"<matplotlib.figure.Figure at 0xb35e850>"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"Log-likelihood: -5.773e+02\n",
"\n",
" Name | Value | Constraints | Ties | prior \n",
"--------------------------------------------------------------------\n",
" linear_variance_0 | 25.7493 | (+ve) | | \n",
" linear_variance_1 | 9.3356 | (+ve) | | \n",
" linear_variance_2 | 5.1650 | (+ve) | | \n",
" linear_variance_3 | 17.2247 | (+ve) | | \n",
" linear_variance_4 | 4.1680 | (+ve) | | \n",
" noise_variance | 30.6099 | (+ve) | | \n",
" [1, 2, 9, 6, 8]\n"
]
}
],
"prompt_number": 451
},
{
"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