2015_01_18a_pymc_sklearn_gp_mashup.ipynb

{"nbformat": 3, "worksheets": [{"cells": [{"cell_type": "code", "language": "python", "outputs": [{"output_type": "stream", "stream": "stdout", "text": "Sun Jan 18 14:37:28 PST 2015\r\n"}], "collapsed": false, "prompt_number": 8, "input": "!date\nimport matplotlib.pyplot as plt, seaborn as sns\n%matplotlib inline", "metadata": {"trusted": true}}, {"source": "Experiment to use the PyMC2 Matern covariance function in the sklearn GP model:", "cell_type": "markdown", "metadata": {}}, {"cell_type": "code", "language": "python", "outputs": [], "collapsed": true, "prompt_number": 9, "input": "import pymc.gp, sklearn", "metadata": {"trusted": true}}, {"cell_type": "code", "language": "python", "outputs": [{"output_type": "pyout", "prompt_number": 20, "metadata": {}, "text": "[<matplotlib.lines.Line2D at 0x7fbd55606210>]"}, {"output_type": "display_data", "metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEDCAYAAAA2k7/eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYHFW5x/HvZCYLu2yGQMIictjRcJFNRC4gBIksYkEB\nCgJ6LbiK3qDIjspFRS+boBSoLAJSUCBhX4IIQkAkEDZZDgECCUsI+xaSSabvH6fGDiGz9XT1qe76\nfZ4nT83WVW/O1Lxddeqc94CIiIiIiIiIiIiIiIiIiIiIiIiIiNRB22B3YIw5H9gVeNVau3H2tRWA\ny4E1gOnA3tbatwZ7LBERqd2QOuzjAmDcIl87CphkrTXAX7PPRUTEo0EnfGvtXcCbi3x5N+Ci7OOL\ngD0GexwRERmcelzhL85Ia+2s7ONZwMicjiMiIv2UV8L/N2ttBajkfRwREeldR077nWWMWcVa+4ox\nZhTwam8/XKlUPgSG5xSLiEhLamtrG9DAm7wS/rXAgcAp2XZiHz8/fKCBt6pKpVJRWzhqiyq1RZXa\nonb1GJZ5GfBFYCVcf/0JwDXAFcDq9GNYpn6BVWqLKrVFldqiSm3R5CqVivr4M2qLKrVFldqiSm1R\nu9wf2oqISDEo4YuIlIQSvohISSjhi4iUhBK+iEhJKOGLiJSEEr6ISEko4YuIlIQSvohISSjhi4iU\nRF7F00REchUk0RLA5sAGwKeBZYERwNu4Cr2PAfenYTzDW5AFU4gCRCqGVKW2qFJbVJW9LcZPmOjq\n5wyZT/sKr9C+4ssMWeZN2oZ09eflj+GKOZ6fhvGLOYZZeIU4gcp+Mi9MbVGltqgqe1uM//HllY5V\nnqPjkzNoa18AQNcHy7Dg7RUZOmr6AYAF3gDmAcsBo4CxwOeBnYBhQCdwMXBSGsbTG/+/8K8QJ1DZ\nT+aFqS2q1BZVZW2LIImWBo6tdA05qm1IF5V5w5n/6hgWvLYqlXlLAnD9aXv02i5BEi0H7A38EDDA\nHOBk4P/SMJ6b7/+gWApxApX1ZF4ctUWV2qKqjG0RJNGewG+A0V1zRzD/5U+xYPZqUGn/yM/1lfAX\n2l87sB/wa9w62w8C+6RhPK2+kRdXIU6gMp7MPVFbVKktqsrUFkESLQucBRyA66L51ZwpXzqOrvbF\n/nx/E/5C+18OOB04CHgX+HoaxtcOKugmoWGZIlIYQRKNBabikv0UYJM0jI/vKdnXIg3jt9MwPjg7\nRjtwdZBEUd0OUGBK+CJSCEES7Q1MBtYCfgFsnYbxU3kdLw3ji4HtgNeBc4IkOjavYxVFIW4Ry3S7\n2he1RZXaoqqV2yJIojbgp8DxwHvAfmkYX9fTz9e7LYIkWhu4DVgT+GEaxqfWa99FU4gTqJVP5oFS\nW1SpLapatS2yB6nnAN8GngV2S8P4X729Jo+2CJJoLeAuYDXgO2kYn1fP/RdFIU6gVj2Za6G2qFJb\nVLViWwRJNAw3Ln5v3IiZcWkYz+7rdXm1RZBEBrgbWCGL5bZ6H8M39eGLSMMFSTQCmIhL9ncB2/cn\n2ecpDWML7AksAK4Mkmg9n/HkQQlfRBoqSKKhuFIHuwA34a6m3/YblZOG8WTgYNxs3auCJFrKc0h1\npYQvIg2T9dlfDHwFmAR8NQ3jD/xG9VFpGF+Km/C1AXC253DqSglfRBoiG43ze2AfXF/5nmkYf+g3\nqh4dCTwAfDNIom/4DqZelPBFpFFOws1uvR/YNQ3j9z3H06Osxs7euGGiZwVJtJrnkOpCCV9Echck\n0SHAscAzuGT/jueQ+pSG8bO4gmvLAXF2h9LUlPBFJFdBEu0EnIub0bqL79E4A3QecDswHtjfcyyD\nVoh3rFYcY1wrtUWV2qKqWdsiSKKNgHtw9eh3yEbBDEqj2yJIojVxi6jMBUwaxq836tj1pit8EclF\nkEQrANcAywAH1iPZ+5AtlnICbkLW//qNZnCU8EWk7rLhl38GPgWcnIbx5Z5DGqyzgCeA72QVPZuS\nEr6I5OEkYGfgRuBEz7EMWhrGncDhuG7ws5r1Aa4SvojUVZBEXwOOBqYB+6dhvMBzSHWR1db5C26d\n3L09h1MTJXwRqZus/syFwPu4iVVv+Y2o7o4E5gMnZSUimooSvojURZBESwCXA0sBB6dh/JjnkOou\nDeNncEM118FNImsqSvgiUi+nAZsAcRrGV/gOJkf/C8wBTsze5JpGrgnfGDPdGPOIMWaqMeafeR5L\nRPzJlieMgEeACZ7DyVUaxi8DZwCrAt/1HM6A5H2FXwG2s9aOtdZunvOxRMSDbInA3+P67fdJw3iO\n55Aa4VfA28CPgiRa0ncw/dWILp2mHL4kIn3LVq1KgGWBQ9MwftJzSA2RPYw+G1gZ+JbncPqtEVf4\ntxljphhjvp3zsUSk8X4CbAZclIbxxZ5jabQzgA9wV/nDfAfTH7lefRtjRllrXzbGrIxb7OB71tq7\nFv25SqVSyTMOEam/J2dP48S/ncbKS67Ar3c+jiWGjvAdUsNdNPVKbrB/5Tub7c8Oa2/T8OMPtKZQ\nw7pbjDEnAu9Za09d9HvNWhgqD2qLKrVFlc+2GD9h4scvyIbMZ/hGkxkyYk4XsG0j6+QU6bwIkmhV\n4DlgBrBu0SeZ5dalY4xZ0hizTPbxUsBOwKN5HU9EGmfo6k8yZMQcgFOatShaPaRh/BJwEbA2roRy\noeXZhz8SuMsY8xBwH3C9tfbWHI8nIg0w5BOz6PjkTLreXxZcH37ZnZltf+A1in4oxG1RkW7RfFNb\nVKktqgrTpdMxlxEb3w3tC5j72NZc9/OvNzymIp4XQRLdCnwJGJuG8UO+4+mJZtqKSD9VGLbWY7QN\n7aRzhqHy4dK+AyqSM7Lt971G0QclfBHpl/aVXqR9+dkseHtFFsxaw3c4RXMzYIH9giQa6TuYnijh\ni0jfhn7I0NWfpLKgnc7nNqIgvcGFkYZxF/Ab3FKOkedwelSI31oR++R8UVtUqS2qfLZFttjHtbhR\nKFEaxuf6iKNbUc+LIImWBl4E3gXWTMN4vueQPkZX+CLSl/1xyf52XGlgWYw0jN8DLgVWA8Z5Dmex\nlPBFpEdBEq2C66p4H/hWGsaaFd+77jfE//IaRQ+U8EVksbKunN8BywNHpWH8nOeQCi8bkjkF2DVI\notG+41mUEr6I9GRvYE/g77jEL/1zHi63Huw7kEV1+A6g1QVJtBKuP28r3GpAa+GumEbgVs15C3gW\neAK4e/b7r3uKVKQqSKKVceV/5wCHZKNQpH8S3OpfhwRJdHKR6usU4kl3UZ+61ypb3Hh33PCs7am2\ncxfwAvA6MBeX9FcCxvDR38VU3ElzSVaro5Ra7bwYjEa3RZBEl+Ou8CekYXx6o47bH81wXgRJdC6u\nH3+XNIxv9h1PN13h11GQREOAfYGTcFfyAPfihrTdBvxrcasBBUk0HNgI+MLYURudPvXlxzYGxgIn\nB0mUAKcWebq2tJYgib6KS/b34h7YysCdj0v4B+AmZRVCId4lm+Eduy9BEm0M/BH4HNCJW/Ltt2kY\nPz6Q/VQqlcrelx+6Iu4P7nvABtm3rgCOTcN4Wv2iLrZWOC/qpVFtESTRCsDjwCeAzxZxBatmOC+y\nB95P4e7eR6Zh/I7nkAAl/EHLrup/DPwUGApcBhyThvH0Wva3cFtkJ83OwM9wbyTzgbOA49Mwfn/w\n0RdbM58X9dbAhH8hcCBuVM4peR+vFs1yXgRJdDzub/fgNIwv8B0PaJTOoARJ9AlgIvBz4FVgfBrG\n+9Wa7BeVhnEl6//bAnfF/zzwP8BjQRLtVI9jiHTLzqkDgQeBjy1UJAN2Sbb9htcoFlKId8lmecde\nWJBEY4BbgfWAvwJhGsavDXa/vbVFkERLAMcDRwLtQIx7qPax5wKtoBnPi7zk3RZBEi0FPIbrgtis\nyM+Mmum8CJLoLuDzwBppGM/wHY+u8GsQJNF6wGRcsj8DGFePZN+XNIznpGF8DG7R6Edxo4DuC5Jo\n/byPLS3vJGBN4NdFTvZN6GLchfX+vgMBJfwBC5Lo08AduCuho3BX2A0tkpT9QW4BnANsDEwJkmi/\nRsYgrSNIos1xddyfxvU5S/2kwDzgG9kzOa+8BwDNc4uWTZW+G1gDODwN47PqfYyBtkWQRHvhhoAt\nC/wCOK5VJsk0y3nRCHm1RZBEw4AHcMOCt0vD+M56H6Pemu28CJLoKuCrwCZpGHtd11tX+P2UDVeb\nhEv2x+WR7GuRhvFVwJbAM8DRwNVBEi3jNyppIj/GJfvzmiHZN6nLs+3eXqNACb9fgiTqwM18XQ83\nZfrnfiP6qDSMnwA2xz083g24PZsaL9Kj7NnPccBLuIEAko8bcCUq9vbdraOE3z+/xC1QfAPwoyKW\niE3D+A1gF1z3zmbA5CCJ1vQalBRWNn/kD7gVmg5Lw/htzyG1rGzOzPWAwdXT8kYJvw9BEu0LHIGb\nNbd/kfvH0zDuBL6Fe4NaB7gnSCKvJ5gU1qHA1kCahvE1voMpgSuyrdduHSX8XgRJtBZwLvAesHsz\nXAVlk7WOxk3QGgXcESTRpp7DkgIJkmh13EXBm7jyHZK/G4EPgMBnt44Sfg+yfvtLgGWA76Zh/JTn\nkAYkDeMzgG/iaqL8NUiizfxGJEWQJZtzgKWBI9IwnuU5pFJIw/gD4DrcnfdnfMWhhN+zY3C3vJcD\nf/IcS03SML4IV61vWeC2bLy1lNu+wJdx1Vsv9BtK6XR36+zjKwAl/MXI+r2PB2YChxbxIW1/pWF8\nCa6WxzLApCCJtvQckniSLcZzJq5r4TvNfF43qZtwbb+nrwCU8BcRJFE7bomyDtwfxZueQxq0NIz/\njJvavRRwc5BEYz2HJH6cgVtw5/g0jJ/1HUzZZDWvbgHWzcqzNJwS/scdiitbkKRhfKPvYOolDeOE\navfOLUESres5JGmgIIl2wb3p34+7yhc/ukdE7e7j4Er4CwmSaDVceYK3gB94Dqfusiv9Q4GVcX36\na3gOSRogm3kd49ZT+FaR1lgtoRtwS50q4fswfsLESve/+a+tOhNYet5zG35izj/HveI7tjykYXwu\nbjr9aFzSH+k5JMnfycDqwC/TMH7EdzBlllXVvRvY0sffXukTfre2pd6iY6WX6Hp/GRbMHu07nFyl\nYfwr3J3Mp4FbgyRa3nNIkpMgibYGvoubOHiy53DEuQZXuPIrjT6wEj4AFYau7pbu7HxhPQpSRDRv\nxwK/w031vjZIohGe45E6yxbM6V5a75A0jD/0GY/8m7d+fCV8oH2FV2hf5i0WvDGSrndX9B1OQ2RD\n8r6Hm2ewDXBJNkJJWsfPcPVbzkzDeLLvYMRJw/gZ4F/Al4IkWrqRx+5o5MGKqYuO0U9T6Wqjc0a5\nBq6kYdwVJNGBwEhgr/mvrD5//ISrWdwdzvWn7VGK255WESTRVsAEXNnsYz2HIx93DW5y507AXxp1\n0NJf4bev9BJDRnzAgtmjqcxd0nc4DZeG8VzcRJDHOlZ5gY5VpnuOSAYr6547H/fOfVA2rV+KZWK2\nbWi3Tm4J3xgzzhjzpDHmaWPMj/M6zmAESTSsY7VpVLqG0PnS2r7D8SYN47eAXbrmjmDo6k/RvuJL\nvkOSwfkJbu2Gs9MwvstzLLJ4DwCvALtkpaobIpcuHWNMO3A2sCPwInC/MeZaa+0TeRxvEA4ZMvxD\ngDOvPyVsuXH3A5GG8cyvHHsxw9e/j6FrPUqlczhd75TjeUYryeol/Qh4FrcCmhRQ1p16M67A4abA\nlEYcN693ls2Badba6dbaTtxqUV4mGvQku+09Flfb4peewymEypxlmPe0q7owbJ0HaVviHc8RyUBk\n5/QFuL/rQ7KFN6S4bsq2uzTqgHkl/NWAGQt9PjP7WpEchIvp7DSMW3KSVS263l2Rzmc3oa19AcPX\nfYC2YXN8hyT9dwKwAfDbNIzv8ByL9G0SbtZt0yf8Qlfhy4Yf/hCYi1ujVhay4I1RdL6wLm3D5jLM\nPADtnb5Dkj4ESfQ53Lq004Gj/EYj/ZEVZrwX2CJIoob0n+Y1LPNFYMxCn4/BXeX3qFKpNOxN4p4X\nHuCMe//Ajmt/gf/abL9X2OecRh26XxrZFr3EwAVTr+Dmp+9g8y/P5OpfdFY62hs/ircIbVEUPbXF\n3PnzOPLWk3n53Vc5Ybvvr7nRyPXeLdo5XW+tcl785fGbSB69lsO3PPi1KyoD/521tbUNaLh0LmOr\njTEduKncOwAvAf8E9u3poW2lUqkMNPBaZSv+3I97ULJuGsZPN+K4/dXItuhLdid0JbAHcDFwYCNr\nqBepLXzrrS2CJPotcBhwehrGExobWeO10nmRlSp/ELg4DeMD8j5eLl061tr5uPodtwCPA5cXaITO\n9sB/AFcVLdkXTVZVcX/gH7hFVH7mNyJZVFb2+DDc39kxnsORgXsINzxzXCOGZ+Z2AGvtTdbada21\nn7bW/iKv49TgyGz7K69RNIls0s5uwDTguCCJvu05JMlk/b7nA53A11Urp/lkd8w340qWb5r38Uo1\n0zZIovVxU5nvTMP4ft/xNIs0jGfjRhK8BpyTXVWKR1nX5LnAKsAJaRhP9RyS1K5hwzNLlfCB/862\nv/EaRRNKw3garpxrJ5AGSZT71Yj06uvAXrja6r/2HIsMTvfwzHF5H6g0CT9IomWBA3Gjha71HE5T\nSsP4H8B+wJLADUESrek3onLK2v1s4D3gAK1g1dyy4Zn344ZnLpvnsUqT8HEPHZcG4jSM5/sOplml\nYXw18H1cV8KNWjylsYIkGgpchlub+PA0jJ/zHJLUx21AO7BtngcpRcLP+ju/C8wDfu85nKaXhvFZ\nuAlr6wMTgyQa7jmkMvkZsCUu6V/oNxSpo9uy7Y55HqQUCR83FHM9IE3D+FXfwbSIHwEp7orkwkZW\n/CurIIl2ws2ifQaIGjknQnJ3LzAHJfy6ODTb/tZrFC0kDeMu4ABgMhDi1siVnLw1521wk986gTAN\nY1W2ayHZuhR/BzYMkmhUXsdp+YQfJNHKuHHkj+ImEEmdZOO+d8fNqj4ySKLDPIfUkoIkGnLWfRcA\nfBI4Mg3jhpTSlYbr7tbZPq8DtHzCxw1fGwqcr1vg+kvD+HXc+OFXgbOCJNrNc0gtY/yEiZXxEyZW\nOmess+DRWU+x4M2VmfPPnU8fP2GizuPWlHs/fksn/Oxh7SG42+BLPIfTsrKRIrsCHwJJtgiH1MGQ\nZV+jY/TTdM0dwbznNian8ldSDI/gJjfumOWuumvphI9biGVD4Jo0jF/zHUwry7oZ9gGGA9cHSfQp\nzyE1vbZhHzDs0w9DpY1508bC/GG+Q5IcZc/FbgdGAyaPY7R6wj842/7RaxQlkYbx9bjZzCsDNzWq\nxncrCpJoiWHrTKWto5PO5zeg8v5yvkOSxsi1W6dlE36QREsB++Jm1k7yHE5ppGEc45aMNMC1QRIt\n4TmkppPdzsdDlnqX+a+OZsHsMX2+RlqGEn6N9gSWAS7S1POGOxb4M7A1cGWQROqLGJjDgAO63luO\nzufX9x2LNFD2PGw6sG0ec1taOeHvl20v9hpFCWV9kQfhqgB+Gbg0SKLGL5fVhIIk2hk4E5g9b9pn\nodLuOyRpvDuBFYCN6r3jlvwjzMbe7wQ8kIbxU77jKaM0jOcFSbQXcCPwNeCDIIkOyt4MZDGCJNoQ\nuAKYD+xx3S/3vaf7e620ypP06U5cocftcCN36qZVr/D3xhUiutR3IGWWhvEc3KS3+3Czcs/Oa7hZ\nswuSaCRwA64o2kFpGN/Tx0ukdd2Rbb9Y7x23asLfD6gAl/sOpOzSMH4XNzHrYVyJi1OU9D8qe7A9\nEVgDODEN48s8hyR+TQdmkEM/fssl/CCJ1sI9LPxbGsYv+Y5H/l3veyfgSVzRtV8p6TvZs43LcBUw\nLwVO8huR+JZVBLgDWAnYoJ77brmEjxuKCerOKZSsSun2uKT/Q+CMsif97P9/Hq4e0e3AISr/IZk7\ns21du3VaKuFnf0D74+re/8VzOLKINIxfxj2Iegw4HPhdycsqn4IbzTQF2COrmCgCSvj9sjHuFuiG\nNIzf8h2MfFwaxrOA/wQeAiLg90ESlW7sYZBER+O6t54CdsmedYh0ewZ4EfhiPe+EWy3hB9k28RqF\n9Cqra7QD8ACu/EUSJNEIv1E1TpBExwA/xz2Y20l1nmRRWdfenbiS2OvVa7+tlvD3wlVsvNF3INK7\nNIzfwE0fvxM3Tv/mIIk+4Teq/AVJdBxwMvACsF0axi94DkmK645sW7dunZZJ+EESrY9bY/XmNIzf\n8x2P9C3rdhsHXIk7qe8Kkmg1v1HlI0iitiCJTsSNwnke+GIaxs96DkuKrbsff7t67bBlEj7u6h7g\nKq9RyIBkq2aFwNm4qeT/CJJoU79R1Vf2jOJM4Ce4MdZfTMN4useQpDk8DcwCtqlXP36rJfxO4Hrf\ngcjAZMXtDgd+DKwG3B0kUeg3qvrInk0kwPeAfwFfSMP4eb9RSTPI+vEn4/4mVq/HPlsi4WeLbXwW\nuE2jc5pTGsaVNIx/hRuTPh+47M+PTKSZR/BkNZ1uwT2j+DuwTRrGM/1GJU1mcrb9fD121hIJn2p3\njsbeN7k0jK/DzTp9ZuITtwDcGiTRqn6jGrisW2oKsC2QAjvrYkRqcHe2VcJfyFeBLuAa34HI4KVh\n/Diw+WarbgJudu7DQRLt4jeq/guS6Ou4K7MxwPFAmD2rEBmoqcAcYJt67KwQU9sHU/o1SKLRuPHM\nf0vDePv6RtZ4KoNb1dXVVdnnisMOB/4PGAb8Bji2qKOwgiRaFjgLVxn0HWC/NIxvqMe+dV5Ula0t\ngiS6A3enuHwaxm8PZl+tcIW/W7ZVd06LaWtrIw3js4AtcDNSDwceDZJoJ7+RfVyQRF/AVQQ9ANeV\ns1m9kr2U3mTcxfmWg91RKyT8r2Tba71GIblJw/gh3EP5X+C6SW4JkujPQRLVZeTCYARJtFKQRH/E\nPZRdHTepaus0jJ/2G5m0kLo9uC3EbVGtt2hBEi0NvA48mYbxZ+ofWeOV7Xa1N4triyCJPgP8AdgM\nmAucDvxysLe6AxUk0XDg27ix9SviVib6ThrG/8jjeDovqsrWFkESLQ+8AdyehvEOg9lXs1/h74Dr\n273OdyDSGGkYP4zr4jkAmA0cBUwPkuh/gyT6ZN7HD5JoeJBEB+O6mM4ChgNHAP+RV7KXcsvWk/gX\nsEWQREMHs69mX9O2uztHk61KJFsX9+Igia7E9esfARwLHBEk0RXAH4G76llbPiv58G1chc+RfPTu\n4tV6HUekB5OBDYHP4J4R1SSX2yJjzE+Ab+GuwACOttbe3NPP13KLltVRfxG3du2obLZm0yvb7Wpv\n+tsW2RKBBwETgLWzLz+LG6Z7A3D3QGvNj59wdaVt+AcMWe412ld4hSHLvElbG1Tmd9DWMf9U4Mw0\njGcM7H9UO50XVWVsiyCJDgAuAn6QhvGZte4nryv8CnCatfa0nPYP8B/AKsCFrZLspTbZYum/C5Io\nxg1fOwTYA/if7N+8IIkeBu4HpuEqVb6Cu0qfh+uWWR535b4usH5l3h20DXPvEZUKdL27PAteH8WC\n11fl+v/72g8b+h8U+egErMIlfMj/gfD4bKvuHAH+3dVzB3BH9lB1W9x58nncKJ/PDWR/819fha53\nV2DBm5+EztKU65dieg53kfL5IInaau2uzDPhf88Y0z0m+Qhrbb2nlY/HFUu7tc77lRaQdeFMyv51\nj6rZCFgDN3xyJO6B/zDcVf4bwGu4CoVPffjQdjMLMohNhDSMK0ESTcaVkVkDV3V1wGpO+MaYSbgu\nlUUdC5wD/Cz7/CTgVNxtdl1kD9A2BSZpaTjpj+wN4IHsX5/GTxiXb0AiA/cPXMLfghoTfu6XMMaY\nNYHrrLUb9/QzlUplQLcntz1zN+dNuZRvjg34smn6agpSQF85oueyTNedunsDIxFxnpj9NCfefhq7\nmh04cOzXABjow+tcunSMMaOstS9nn+4JPNrXawYSeDYcb68Lp6br7rruDrbGMAupjCMQeuKzLcZP\nmNjjRYiPmHReVJW1LYIkWhJ45wb71/u+uWlQ06zbvIZl/gn3kKyCe9jwHWvtrJ5+fiC/wCCJOnB9\nrW8Cn6rnWOsiKOvJvDhqiyq1RVWZ2yJIogeBDYBl0zCeN9DX53KFb609II/9ZjYHlgOSVkv2IiJ9\nuA8YC2xCDROwmrG0QnelxFu8RiEi0nj3ZdstanlxMyb8nYEFwO2+AxERabDuek2tn/CzqnGbA/c2\nujqiiEgBWOBtaqyN31QJH9gRF7O6c0SkdLLZ5P8E1gmSaIWBvr7ZEv7O2VYJX0TKqrsff/OBvrBp\nEn6QRG24hP868KDncEREfOlO+APu1mmahA+sD4zGlVNQdUwRKauaR+o0U8LfMduqWJqIlFYaxrNx\n6z20dMLvLpqj4ZgiUnb34dZwGJCmSPhBErUD2wHPpGH8vOdwRER8u6/vH/m4pkj4uKnEy6GrexER\naPGE392d8zevUYiIFMNDwPyBvkgJX0SkyaRh/CHwyEBfV/iEHyTRMGAb4PE0jF/xHY+ISEFMGugL\nCp/wcQtPL4X670VEFnbcQF/QDAlf3TkiIotIw7gl+/C3x62cdafvQEREmlmhE36QREsAWwMPpWH8\nuu94RESaWaETPi7ZD0P99yIig1b0hP+f2VYJX0RkkJoh4XcBd/sORESk2RU24Wf9958DHkzD+B3f\n8YiINLvCJnzcai5Dgbt8ByIi0gqKnPC3zbZK+CIidVDkhP+FbKv+exGROihkwg+SqAPYCngiW91F\nREQGqZAJH/gssDTqzhERqZuiJvzu7hwlfBGROilqwtcDWxGROmvzHQBApVKptLW1tQEESdQGvAp8\nkIbxGn4ja7yF26Ls1BZVaosqtUXtiniFvx6wErq6FxGpqyImfPXfi4jkQAlfRKQkipjwtwVeA57w\nHYiISCspVMIPkmg0sDowOQ3jiu94RERaSaESPm7BE4B7vEYhItKCOmp9oTEmAH6CG1XzOWvtgwt9\n72jgYGABcLi19tZ+7narbHtvrXGJiMjiDeYK/1FgT+DvC3/RGLMBsA+wATAO+J0xpr/H2QqYD0wZ\nRFwiIrLKOwGEAAAGFklEQVQYNSd8a+2T1lq7mG/tDlxmre201k4HpuFq2/cqSKIRwKbA1DSM59Qa\nl4iILF4effirAjMX+nwmsFo/XrcpbsETdeeIiOSg1z58Y8wkYJXFfOsYa+11AzhOnyNuvvGZvSZf\n/PBVfH+rgw+/onLO4QPYd8upVCoaoZRRW1SpLarUFs5AS0z0mvCttV+qIYYXgTELfT46+1qvLn74\nqr8AXz3z3vPX3GaNzZ+v4bgtQXVCqtQWVWqLKrVF7WoepbOIhRv/WuDPxpjTcF056wD/7O3F2Zv1\nVsDLwAt1iklERBZScx++MWZPY8wMYEvgBmPMTQDW2seBK4DHgZuAw6y1vd5+vfbBGwCjgHs14UpE\nJB81X+Fba68Gru7hez8Hft7vfb3+bPeHemArIpKTQsy0ta891/2hEr6ISE6KkfDdFX4n8IDnUERE\nWlYhEv70N2cAPJiG8Ye+YxERaVWFSPgLKl2g7hwRkVwVIuFnbvMdgIiI5GzGWy9VssXLS08zCKvU\nFlVqiyq1RZPTL7BKbVGltqhSW1SpLWpXpC4dERHJkRK+iEhJKOGLiJSEEr6ISEko4YuIlIQSvohI\nSSjhi4iUhBK+iEhJKOGLiJSEEr6ISEko4YuIlIQSvohISSjhi4iUhBK+iEhJKOGLiJSEEr6ISEko\n4YuIlIQSvohISSjhi4iUhBK+iEhJKOGLiJSEEr6ISEko4YuIlIQSvohISSjhi4iUhBK+iEhJKOGL\niJSEEr6ISEl01PpCY0wA/ARYD/ictfbB7OtrAk8AT2Y/eq+19rDBhSkiIoNVc8IHHgX2BM5dzPem\nWWvHDmLfIiJSZzUnfGvtkwDGmPpFIyIiuRnMFX5v1jLGTAXeBo6z1t6d03FERKSfek34xphJwCqL\n+dYx1trrenjZS8AYa+2bxphNgYnGmA2tte8OMlYREfHJGPO3LLHX9H0REWmMeg3LbOv+wBizkjGm\nPfv4U8A6wLN1Oo6IiNSore8fWTxjzJ7Ab4CVcH31U621uxhj9gJ+CnQCXcAJ1tob6hGsiIiIiIiI\niIiIiIiIiEirqfmhbb0YY8YBZwDtwB+stad4DskLY8wY4E/AJ4EKcJ619jd+o/IrG+01BZhprf2K\n73h8McZ8AvgDsCHu3DjYWvsPv1H5YYw5Gvg6bkDIo8BB1tq5fqNqDGPM+cCuwKvW2o2zr60AXA6s\nAUwH9rbWvtXTPrxWy8z+oM8GxgEbAPsaY9b3GZNHncD/WGs3BLYE/rvEbdHt+8DjuCRXZmcCN1pr\n1wc2wRUnLJ2sMOO3gU2zhNcOhF6DaqwLcLlyYUcBk6y1Bvhr9nmPfJdH3hxXaG26tbYTSIDdPcfk\nhbX2FWvtQ9nH7+H+qFf1G5U/xpjRwJdxV7be70R9McYsB3zBWns+gLV2vrX2bc9h+fIO7sJoSWNM\nB7Ak8KLfkBrHWnsX8OYiX94NuCj7+CJgj9724TvhrwbMWOjzmdnXSi27khkL3Oc5FJ9OB36Eu3Uv\ns7WA2caYC4wxDxpjfm+MWdJ3UD5Ya98ATgVewJVwectae5vfqLwbaa2dlX08CxjZ2w/7Tvhlv1X/\nGGPM0sCVwPezK/3SMcaMx/VTTqXEV/eZDmBT4HfW2k2B9+njtr1VGWPWBn4ArIm7+13aGLO/16AK\nxFpboY+c6jvhvwiMWejzMbir/FIyxgwFrgIusdZO9B2PR1sDuxljngMuA7Y3xvzJc0y+zMQ9tL4/\n+/xK3BtAGW0G3GOtfd1aOx/4C+5cKbNZxphVAIwxo4BXe/th3wl/CrCOMWZNY8wwYB/gWs8xeWGM\naQP+CDxurT3Ddzw+WWuPsdaOsdauhXsod7u19gDfcflgrX0FmGGqC0/sCPzLY0g+PQlsaYxZIvt7\n2RH3UL/MrgUOzD4+EOj1QjGvevj9Yq2db4z5LnAL7on7H621pRyBAHweN9zskWwtAYCjrbU3e4yp\nKMre9fc94NLsougZ4CDP8XhhrX04u9Obgnu28yBwnt+oGscYcxnwRWAlY8wM4ATgl8AVxphDyIZl\n+otQRERERERERERERERERERERERERERERERERKQJ/T8VsNXWd03wTQAAAABJRU5ErkJggg==\n", "text": "<matplotlib.figure.Figure at 0x7fbd557a9490>"}], "collapsed": false, "prompt_number": 20, "input": "import numpy as np\nfrom sklearn import gaussian_process\n\ndef f(x):\n    return x * np.sin(x)\nX = np.atleast_2d([1., 3., 5., 6., 7., 8.]).T\ny = f(X).ravel()\n\ngp = gaussian_process.GaussianProcess(theta0=1e-2, thetaL=1e-4, thetaU=1e-1, nugget=.00001)\ngp.fit(X, y)  \n\nx = np.atleast_2d(np.linspace(0, 10, 1000)).T\ny_pred, sigma2_pred = gp.predict(x, eval_MSE=True)\nsigma_pred = np.sqrt(sigma2_pred)\n\nplt.plot(X.ravel(), y, 's')\nplt.plot(x, y_pred)", "metadata": {"trusted": true}}, {"source": "sklearn has a nice gaussian process implementation, but doesn't have the fanciest Matern covariance functions.", "cell_type": "markdown", "metadata": {}}, {"cell_type": "code", "language": "python", "outputs": [{"output_type": "pyout", "prompt_number": 58, "metadata": {}, "text": "[<matplotlib.lines.Line2D at 0x7fbd4f36e3d0>]"}, {"output_type": "display_data", "metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEDCAYAAAAhsS8XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYHFXVx/HvZIEIISBbEkgEFK7GsEhAQBRxCZhg2Cko\nQGRRoBBkVwgI8iqyCQEEoV5kURZTUCwBY1BWFWVRIOzLZRUCJIQtrAkJ6fePW0OHvJnMTHd136ru\n3+d5eHpm0l11OOk5qb5177kgIiIiIiIiIiIiIiIiIiIiIiIiIlJoHT15kjHmYuC7wKvW2nWyny0P\nXAmsBjwP7GytfatBcYqISA/16eHzLgHGLPSzo4GbrbUGuDX7XkREPOtRYbfW3gG8udCPtwH+kH39\nB2C7HOMSEZEa9fSKfVEGW2tnZF/PAAbnEI+IiNSpnsL+MWttBajkcSwREalPvzpeO8MYM8RaO90Y\nMxR4tbsXVCqV2cCSdZxTRKTtdHR09GiiS6d6CvsNwJ7AqdnjpB68ZsneBtiqKpVKRblwlIsq5aJK\nuahdT6c7TgQ2B1bEjacfD1wPXAV8hh5Od9RfVJVyUaVcVCkXVcpFSVQqFY3DZ5SLKuWiSrmoUi5q\nl8vNUxERKQ4VdhGRFlPPzVMREQGCJFoWGA7MBp5Lw/gjn/E09caEboZUKRdVykWVclFV9FwESdSB\nW3F/JPAVqvX0bWAicEYaxk/5iE2F3RPlokq5qFIuqoqciyCJhgKXAd8G5gP/Ah4FBuJmEA4HPgRO\nAE5Nw3h+M+NTYfdEuahSLqqUi6qi5iJIonWAm4AhwBTgiDSMn1jgz/sCOwFnZc+5Bvh+GsbvNytG\nFXZPlIsq5aJKuagqYi6CJPoScBvwaeAnuOGWRU7LDJJoBVxR3xz3D8E2aRjPaUacKuyeKBdVykWV\nclFVtFwESbQacBfuKnzvNIz/0M1LCJJoCVxxHwdcDezSjGEZTXcUEelGkERLAX8GhuKGXrot6gBp\nGH8IBMA/cMMzTdm3QoVdRKR7pwEjgfPSMD6zNy9Mw3g2rri/CJwYJNF3GhDfJ6iwi4gsRpBEWwIH\nAo8BR9RyjDSMXwV2AOYBFwdJ9On8Ivz/VNhFRLoQJNHyuK1B5wHfy66+a5KG8b246Y+rAL/JJcAu\nqLCLiHTtt7hCfEIaxlNzON5pwL+B7wVJtPA+0rlRYRcRWYQgiXYFQtxMmFPzOGYaxvOAfXGLms7O\nZs3kToVdRGQhQRINA84D3sMtLpqX17HTMH4oO7YBDs3ruAtSYRcRWUCQRH2Ai4HlcFMbn27AaY4H\nXgOOz9oT5EqFXUTkkw4EtsC1C7igESdIw/hN4DhgaeCYvI+vwi4ikgmS6Au4G5yvAz/sql1ATi4C\nngX2z1a15kaFXUQECJKoP65j4wBg/zSMX2nk+dIwnoub/tgfNzSTm7oLuzFmvDHmUWPMw8aYPxpj\nlswjMBGRJjsW2BC4LA3ja5p0zj/iFj7tFSTRmnkdtK7CboxZHTd1Z5S1dh2gL256kIhIaQRJtDHw\nM9yy/x8367zZTku/wNXiI/M6br1X7G8Dc4GljDH9gKWAl+qOSkSkSYIkWga4AlcP90rDeFaTQ7ga\neAZ31T4kjwPW3RLTGLMfcAbwAfBXa+0eXT23aG04fVIuqpSLqnbPxbjDJ3V5s3LyhO0akpcgiS4B\n9gJOS8P4qEacowcxRMD5wClpGI+v93j1DsV8DjfBfnXcstuBxpjd6w1KRKQZgiTaGVfU78NNP/Tl\nD8CrwAFBEg2q92D96nz9hsCd1trXAYwx1wKb4j7WLFKlUmnk9KFSUS6qlIuqds7F1kdc3+Wf5Z2X\nGe/O5OibTmbe/I84dcvxG6wyaMgcdjk/z1P0ynWP/YWJD1/P99bbYdZVlU/G0dtPcfUW9ieA44wx\nnwJmA6NxDW661M4fMxfU7h+5F6RcVLV7LhY3FJNnXrKNM+4C1gX2XnXZob/P69i1ylr5Trv8wWtf\nvfzBa9fMbqzWpK6hGGvtg8ClwL3AQ9mPG7JSS0QkD0ESdQC/wxX1OA3j3/uNyMlWo16GG9oeV8+x\n6r1ix1p7Gm6llohIGRwN7AbcTYOacNXht8D+wEFA1+NS3dDKUxFpG9nsk5Nw89V3SsN4jueQPiEN\n44eBvwGjgyQaUetx6r5iFxHJy4JTGvO+3xAkUYhrlzsTGJ2GcVHX3JwDfAN31X5gLQfQFbuItLwg\nibbCjV+/A3wnDWPrOaTFuQH3iWLPIImWreUAKuwi0tKCJNoMuAa3b+m4nLa4a5hsU4/zcS1996zl\nGCrsItKygiRaH5iMG3beMQ3jOzyH1FMXAh8CUTaLp1dU2EWkJQVJZIC/AssAe6RhPMVzSD2WhvFM\nXA+ZEcBmvX29CruItJwgiVYBbgJWAn6UhnHiOaRaxNlj1NsXqrCLSEvJbjjeCKwGHJeGcdzNS4rq\nn7he7Tv19oUq7CLSMoIk6gukuFWl5wG/8htR7bJt+f4Xt8NSr6iwi0gr+R/cRtSTgYMbvGdpM1yG\na4neKyrsItISgiT6Dm57u2dxN0trbqJVFFn/mN16+7qmdpFr9851C1IuqpSLKuWiqje5yHqYPwoM\nBjYu+lz1RtMVu4i0gpOBYcDJ7V7UQYVdREouSKINgB8Bj+MafLU9FXYRKbtTsseDitat0RcVdhEp\nrSCJRuN2brspDePbfMdTFCrsIlJKWQ+VznnqR/uMpWhU2EWkrDYDNgKu0w3TT1JhF5GyOiJ7PMNr\nFAVU9w5KxpjlcC0mRwIVYB9r7d31HldEpCtZ58atgXuAOz2HUzh5XLGfDUyx1o7A9Wd4PIdjiogs\nzoG4BZYTWqBtQO7qWuFmjFkWmGqt/WxPnq9VdVXKRZVyUaVcVHWViyCJlgRexu2INCwN47lND67g\n6h2KWQOYaYy5BFgPuA84xFr7ft2RiYgs2rbA8sDpKuqLVu8V+4bAXcCm1tr/GGPOAt621h6/qOdX\nKhV9ZBKRuvzq7+fw4PTHmDD2eIYNGuo7nKbo7ae4eq/YpwHTrLX/yb6/mm7mk+pjpqOP3FXKRZVy\nUbWoXARJNAx4Abh7+LKrbOonsuKr6+aptXY68KIxxmQ/Go3rsCYi0ggBbqThUt+BFFnd0x2BHwNX\nGGOWAJ4B9s7hmCIii7IzMB+4xncgRaZ+7J4oF1XKRZVyUbVwLoIkWg14Hrg1DePR3gIrAa08FZGy\n6NzU+SqvUZSACruIlMXOwEfAdb4DKToVdhEpvGwYZiPgtjSMZ/qOp+hU2EWkDLbKHnW13gMq7CJS\nBmOzx794jaIkNCvGE+WiSrmoUi6qOnOR9YZ5HZiWhvEXfMdVBrpiF5Gi2wxYGrjRdyBlocIuIkXX\nOQyjwt5DKuwiUnRjgPeBf/gOpCxU2EWksIIk+gzwReD2NIxn+46nLFTYRaTINAxTAxV2ESkyTXOs\nQR7dHUWkBYw7fFKXG+FMnrBd06dgzvtoHsC3gafSMH6m2ecvM12xi0ghPfHa0wAD0TBMr6mwi0gh\nPTD9sc4vVdh7SYVdRApp6iuPAswG/u45lNJRYReRwgmSaPiLs14GN83xA9/xlI0Ku4gU0ZjsUbNh\naqDCLiJFpPnrdchlCpMxpi9wLzDNWrt1V89T57oq5aJKuahSLiBIov7A64MHrrTMueN+2da5qFVe\nV+yHAI8BXc6DFRHpoU2BZdYfMtJ3HKVVd2E3xgzD7W5yIU3u7y4iLWkswJeGqrDXKo8r9jOBnwDz\ncziWiMhYYM7IlY3vOEqrrpYCxphxwKvW2qnGmG/05DWVSkXDNRnlokq5qGrnXLzx/ltEfxrPekNG\nsGS/Jdo6Fwvq7X2XenvFbApsY4zZChgADDLGXGqt/X5eAbYq3SSrUi6q2j0XQRLtA1z04PTHDwPO\nbOdc1KOuoRhr7THW2uHW2jWAELhtcUVdRKQbmuaYg7znsetjk4jUJEiifsAWwHOA9RxOqeXWttda\n+3dK1NMhSKIlgFWAVbP/lgeWBQZljwMWePpHwEzgJeB+4IE0jOc0NWCR1vcV3O/eFWkYV9jlfN/x\nlFZb9GMPkmgV3E7nGwBrAyOBz9RxyLeDJJoInJqG8XM5hCgi2lQjNy1Z2IMk6gNsAgTAOGDNhZ7y\nCu7TxYvANNyV+OvALODt7HE21aGl/sBKwOrAxsA2wP7AXkES/QI4JQ1jTfcUqc8Y4EPgNt+BlF1T\n7zg3+o5/kESDgH2AHwOfzX78Dm53838AdwOPpGH8Rp3n6QvsDPwaN4yTArunYTy3p8do99kPC1Iu\nqto1F9mn6peAW9Iw3gLaNxd5aIkr9iCJPgUcChyNGyOfDVwKXIl7o3yY5/nSMP4ImBgk0S3A1bhP\nBnODJNpDV+4iNekchvmz1yhaROkLe5BEY4ALgOG44ZRjgQvSMH6t0edOw3hmkERbATcDu+H65fyq\n0eeVYgiSaCSwNfA53HDdU8DkNIwf9BpYOX03e1Rhz0Fph2Kyq/QJQATMzb4+OQ3jWXkcv5exrAA8\nCAwBNk3D+N/dvUYfM6vKlosgidYBTge27OIpNwOHpmH8WBd/3qWy5SIPQRItCbwGTE/DeK3On7dj\nLvJSyn7sQRINBf6GK+qPABulYXy0j6IOkIbx68D3cfm8OJuPKy0oSKIf4VpUb4m7ybc7MAJYC9gF\nuBU3F/u+IIl+4CvOkvkabtNqXa3npHQFKEiizwO3AMOAy4D90jCe7TcqSMP4tiCJLgZ+AOwLaBJu\nCwmSqAM4Bfgp7upynzSM/7TQ054GrgqSaFvgIuDCIIlWSsP4lOZGWzoahslZqa7Ys6L+N1xRPwbY\nswhFfQHH4mbh/DJIomV8ByO5OhFX1C2wwSKK+sfSML4e10fpReDk7CpfuvZd4D3czDXJQWkKe5BE\nqwO348axD0nD+OQ0jAvVwiAN4xm4KZAr4Oa5SwsIkmg/3IXE08BmaRi/0N1r0jC2wDeAV4FzgiTq\najy+rQVJtCZgcLPXtJo7J6Uo7NnV75+AocARaRj/xnNIi3Mu7qr9iCCJBnT3ZCm2IIm+ApyHG34Z\nm4bxqz19bRrGz+IWs30EXBYk0eDGRFlqW2WPU7xG0WIKX9izxUB/xLUCODcN4wmeQ1qsNIzfxBWC\nIcAensOROgRJtBwwETd7bKc0jJ/u7THSML4Ht75iZeCSbKxeqrbNHlXYc1T4wo4b1xwH3AQc5jmW\nnjoHd5UW+Q5EapMV4BhYDTgxDeN6GtydhXv/jkX/2H8smya8OfDvNIyn+Y6nlRS6sAdJ9GXgF8DL\nwG5pGM/zHFKPpGH8EjAZGBUk0Ya+45GabI+bvvgv4Jf1HChbjbwf8D5wepBEy9cfXkvYGugLXOs7\nkFZT2MIeJNHSwBW4KZnfz+aKl0mcPeqqvWSynkPn4BpS7ZPHBUUaxv/FXaSshFYnd9ohe7zOaxQt\nqLCFHfg5btHHGWkY3+o7mBrcBDwP7Kqpj6VzIq5X/6+y2S15mQA8AewXJNEXcjxu6QRJNBC3yOuR\nnHMsFLSwZ0u2DweeBY7zHE5Nso/fvweWArbzG430VDZ0dhDwJHBqnsfOun+Ox/3etftV+1hgSTQM\n0xCFK+xZL/UYN/Z2UBrGH3gOqR4Ts8ddvUYhPZLdMJ2AmwVzQIPmVV8P3AXsECTRxg04fll0DsOo\nsDdA4Qo7rkvipsA1aRiXekPb7CPmvcCWQRKt5Dse6dZ2uJ22rk/D+PZGnCBbVHd09u3JjThH0QVJ\ntBRupttzwEOew2lJdRd2Y8xwY8ztxphHjTGPGGMOrvVYWZe3E3E3rY6oN7aCmIj79LGT70Cka9ke\nuKcB84CjGnmuNIz/AfwV+GaQRJs28lwFtTWu6dfEoq0ebxV5XLHPBQ6z1o7EbUd3oDFmRI3HOgA3\nb/jcbBZBK7gSt8Ve6DsQWaz9cVsoxmkYP9mE83WOsR/bhHMVze7Z4xVeo2hhdRd2a+10a+0D2dfv\nAo/jZhT0SpBEywI/w+05elK9cRVFNqf9TuBrGo4ppmyF6c9x771fNOOcaRjfAdwBbBUk0frNOGcR\nBEm0Iu7G6QO19KuXnsl1jN0YszqwPnBPDS8/Atc865QSzlnvziRcrrf2HYgs0njce++kNIxnNvG8\nnVftxzTxnL4FuLUpl/sOpJXl1rfCGDMQ11L3RGvtpEU9p1KpLHI87f0PP+CAycewRJ/+nDvuRJbs\nt0ReYRXC9Hde5eApP2eDVdbhqM3UwbVIXn3vdQ6bcgLLDhjEWWN/zhJNfO9VKhXG33wKz735IhPG\nHs+qg4Y07dy+HHfr6djXnuX8rU9i+aWW8x1OafR2J6lcrtiNMf2Ba4DLuyrqnToWYa/rDj/2g7mz\nmTXnnaMH9F9yUU8ptaGDBncAj9z38sNzdr7ygGU6Otxfku+4isJnLg6a/LM/zp0/j9fef+N7Szb5\nvdenT5+OZ998YccKFQ678X8u7uho7ffFzlceMOLJ156hQuWWFZb+dLfPb+Vc9FZva3Ies2I6cLvF\nPGatPau3r89aBxwGvEVr7zo0Cbcg4zu+AxEnW4y0G3A/1TUHzXY9bvOOPYIkWtVTDM3SuUfBBV6j\naAN5XLF/Ffge8E1jzNTsvzG9eP2+wIrAOWkYv51DPEXV2Q9je69RCPDxYqTTs2+PzFYKN10axh/h\nNmfpDxzqI4ZmyDaf3xOYgfvHTBqo7j1PrbX/pMZ/IIIk6g8cidsW6+x6Yym4qcA0YGzWY1782hrX\nMnZyoxYj9cJluNk40bsfvuc5lIYJgE8DJ6dh/KHvYFqd75WnOwKrAhe14EyYT8gWYkwBlgc28hxO\nW8suKE7D9cz/qedwyFoXnAkMvOnplt32c3/ceo7f+Q6kHfgu7Afj/rLP8RxHs3TuErPVYp8ljfZD\n4PPAhWkYP+47mMz/ArNutLd3Dlu0jCCJNsG1CbkxDePnfMfTDrwV9mwTja8AU2rZcqykbsWt1FVh\n9yTrtX4CbvjvBK/BLCC7v3T+rDnvgBuLbiU/yR5/7TWKNuLzir2zp0yrj61/LA3jd4F/AKPe/GCW\n73Da1Qm4/UdPScN4uudYFnZ2/z79AI5slfswQRIZ3ISBe4F6theUXvBS2IMkGoLbduxx4BYfMXg0\nBeCBVx71HUfbyfr8Hww8Q3VGTGGkYTx989U3Afgc7v5TK/gpbiHkaWr41Ty+rtj3x03vOqcN/7Kn\nAExVYW+qbHrjubhOmwenYTzbc0iLtM0XtgCYDxydxVxa2dX6Xrhdo9R3vYmaXtiz9qgHALNw07za\nzZPAcw/OeKxzdoY0x+7A13G91qd092RfhiyzMrhV3OsDo/1GU7f/wf1D+rNsvr40iY8r9p2Bwbgp\nju96OL9XndMeP5g7G9xMAWmwbOjvbNwN0zIsAurckq+hfeEbKUiiUbhW1fehq/Wmq3uBUm9kPcAO\nwU1xPLeZ5x53+KQuh3wmT9iu2R95pwAH4mbH6IZSA2XDGTFu/cBBaRg/7zei7qVhfF+QRLcC3w6S\naIM0jO/zHVNvZNtbdv5+H92Gw63eNfWK/anXnwPYELihzeez/q1/3/6gaY/NsBuwLa7zaJl6EZX5\nqn0v3FTmNA3jdpscUQhNLew3PvXxyu3fNPO8RZOG8ftrr2wA1g6S6DO+42lVQRKtBZyHG4LZx1c/\nmBrdgmtOtmP2/1EK2bDXabicH+45nLbV1MJ+94v3AzwC+O7N4d36Q9fu/HKszzhaVdY19BpgEBCV\n7RNiNnxxKu539EjP4fRINux1MW7TkvFpGE/zHFLbamph/6gyH+A3GnOD9YeO7PxSwzE5W2BcfR3g\nvDSMy7pbzzXAU8A+QRJ93ncwPXAg7kLlJuC3nmNpa00t7Guv/HnQBrYADB64Erj5vaODJBrgOZxW\ncwCulfQ9lHg4IJsieBRukkPhFlQtKEiib+Aamb0G7F2yYa+W09RZMcd/81A6Ojreb+Y5C24KrvB8\nHXeV0xQFmyGUq6zh1Fm4AhNknRPLbBLuxu+4IIm2TMO4ae+TnsruAVydfbtjGsYv+4xHmlzYfSpo\nweos7FvRxMLeqoIkWglIcYtiwjSMX/QcUt3SMK4ESXQY7kbqWUESjSrSqtkgidYAbsONq++XhnHL\n9h0uE99te9vdHcC7wHd9B1J2WdOsicAw3ErHWz2HlJs0jB/AjVmPwG3IUQhBEn0R92liGHBUGsbq\ntV4QKuweZTvJ3AysWaYpbQV1AvBt4Aaqc8BbydHA07jOj5v5DiZIoi2BO4HP4GbAnOY5JFmACrt/\n2nyjTlmRORZ4FtizFW/cpWH8Hq5PewVIgyQa5iOOIIkGBEk0AfgrsBSwRxrGp/iIRbpWd2E3xowx\nxjxhjHnKGFPGVXK+dRZ2DcfUIEiiVYDLgXnALmkYv+U5pIZJw/hO3D2ZwcDkIIlWaNa5gyTqCJJo\ne+Bh4DDAAl8t8VTSllZXYTfG9MX1hBgDfBHY1RgzIo/A2kU2g+ABYPMgiQb6jqdMsp4klwErAUek\nYXyv55Ca4Te41bTrAbc1+so9SKL+QRLtips6ei2wOq6h2qg0jP/TyHNL7eqdFbMR8LS19nkAY0yC\n68tRlH0ky+LPwJeAb+HGiBuqoDOEarEfLmd/oslN5XzJZsn8GDckcyAwNUiiCLg2r4V/2QXG13Cb\nfeyAa6BWwRX2Y9IwfjKP80jj1FvYVwUWnFI2Ddi4zmO2oym4MeLv0oTC3gqCJBqO60kyC9cyoG1W\nM6dhPD8r7o/i5uxfDdwZJNH5wF/SMH6tJ8fJZhINxu3YNBL3qXsjYAOqtWF6do7fttHexKVXb2Hv\n9S9TJevdK9VczJ8/nx9e/1OW7LvEfvPnz9+vo6NVLqh7rrfviwn/+h13T7uf6Mt78K3PbvoSu5Sp\ncePi9SYXL78zg8sfuJZ7X35oU7L+/gdNPo6Vl16BQUsOpF+ffvTt05e5H83lg3lzmDNvNu9/OJs3\nZr/FW7PfZuFT9e3ow2eXX40RK63F+kNHMmLFNYf06dPnUOBQHzlWvXA6elkU6i3sLwHDF/h+OO6q\nvUu9DbBVVSqVyoK5CJLoind5b7ddrvrRumkYP+wztmZbOBfdCZLoa7g1AHfH/7ls029/7qst88vf\n21x0CpJoTdw+wt+a8e7MtWe8O3PlxTz9A+AV4GXc7/ALuKv/Rz+qzH/8pC2Oeq+G0HNXay6k/sJ+\nL7CWMWZ13JtkF2DXeoNqU1NwvcPH4WYeyCJkN0zPzL49rJ2GYBYnGyb5VfYfQRItBXwat7dwf2A2\nbjHcu2kYz/UVpzRHXYXdWjvPGHMQbk5rX+Aia61unNbmRtyUve2Bkz3HUmTb4zZruTIN47t9B1NU\naRi/D6gvU5tq6sccfbSqWlQugiS6CdgCWC0N4xf8RNZ8PX1fZO1478PNIPpCGsa24cE1mX5HqpSL\n2mnlabFclz1u5zWK4hoLrA9c1YpFXSQvKuzFMgk302gH34EU1LHZ40leoxApOBX2AknD+BXgLmCz\nrAWtZIIk2gA3nW9KGsYP+Y5HpMhU2IvnOtzfyza+AymYH2WPbbHCVKQeKuzF0znOvqPXKAokSKLl\ncVNBn8HNwBKRxVBhL5g0jJ8BpgJbBEm0ou94CmIvYABwfiu25BXJmwp7Mf0Rt8Yg8B1IQewJzAV+\n7zkOkVJQYS+mBDc7ZjffgfgWJNG6wLrAn9Mwft13PCJloMJeQGkYTwP+DnwtSKLVfMfj2R7Z42Ve\noxApERX24roie2zbq/asrexuwFu4nvUi0gMq7MV1DfAhsEe2lL4dfRNYBbfSdI7vYETKQoW9oNIw\nfhO3EnUE8BXP4fjS+WlF+2qK9IIKe7FdkD3u5zUKD4Ik6ofbZvEV4F+ewxEpFRX2Yrsdtyhn5yCJ\nlvMdTJN9HbfX5nWauy7SOyrsBZYVtN8BnwJ29xxOs3U2QrvWaxQiJaR+7J70ogf5YNx2g08A67bi\njkGL2CawD26T9AHAkFbf8Wfc4ZO6/DudPGG7tv19Ub2ona7YCy4N4xlACqyN24SjHWyEmw1zQ6sX\ndZFGUGEvhzOyxyO8RtE822ePGoYRqYEKewmkYXwf8Ddgy2yJfasbi9t8+RbfgYiUUV2bWRtjfg2M\nwy2keQbY21o7K4/A5P85HfgGcCTwfb+hNE6QRKsC6wB/ScP4A9/xiJRRvVfsNwEjrbXrARYYX39I\n0oUbgceA3YIkWtN3MA00Jnv8i9coREqsrsJurb3ZWts5x/geYFj9IcmiZFMfTwD6Aj/3G01DqbCL\n1KmuoZiF7ANMzPF48v9dAzwI7B4k0clpGD/mO6A8ZatNtwCex30CbAsLTmnUFD/JQ7dvIGPMzcCQ\nRfzRMdbaP2XPORYYZa1d7HZulUql5eZgN9u9Lz3Iaf+M2WTYKA7/6r6+w8nVEzOf4fjbTmfLz32d\nH264q+9wRAqjt//Y131lYIzZC9gX+La1dvbinqurkapac5F1erwL2Bj4ehrGd+QeXJN15iJIohOB\nY4Ft0zC+wXdcPuh3pEq5qF1dY+zGmDHAT4Btuyvqko9s5ekh2bdnZz3LW8UY3BZ4t/sORKTM6p0V\ncw4wELjZGDPVGHNeDjFJN9Iwvge4FFgfd2+j9IIkWhnYAPhnGsbv+I5HpMzqunlqrV0rr0Ck18YD\nOwInBUl0dda/vcy2zB5v9BqFSAvQytOSSsP4ZeCXwIrAqZ7DycPY7FHTHEXqpMJebhOAh4B9gyTa\n3HcwtZo/fz7Ad4CXgUf8RiNSfirsJZZ1PtwXqAAXBEk0wHNINXn2zRcAVsC1EdCUWJE6qbCXXBrG\n/wbOBgzwM8/h1OSB6R+vs9L4ukgOVNhbw3HAf4Gjytj98YFXHgX4CHVzFMmFCnsLSMP4XeAA3Cyn\nC8o0tz1IouWfeuM5gLvTMH7LdzwirUCFvUWkYXwjrlfPxrgiXxajs04Tmg0jkhMV9tZyGPAmcHKQ\nRMN9B9NDndMcNb4ukhMV9haS7Y/6E9xq4HM8h9OtrO/NmEFLDgSY6jkckZahwt56Lgb+DmwbJNH2\n3T3Zs3WBIesN+WJnv3kRyYEKe4vJ5oHvj9uu8NwgiZb1HNLijAH40pCRvuMQaSkq7C0oDeMngROB\nVYCTPYfibGLJAAAIj0lEQVSzOGOBynpDRviOQ6SlqLC3rlOBx4EoSKJRvoNZWPZJ4mvAvwcNWMZ3\nOCItRYW9RaVh/CFwMG4zlbOzG5VFsgVu/1bNhhHJmQp7C0vD+BZgEu7KeBfP4Sxsq+xxitcoRFqQ\nCnvrOwJ3I/XXQRIt7TsYgCCJ+uDG12cC93kOR6TlqLC3uDSMnwXOAIYBP/UcTqcv4TZIv1HTHEXy\np8LeHk7C9Tr/aZBEq/kOhupqUw3DiDSACnsbyJqEHQUMoBi7LW0FzAdu8h2ISCuqu7AbY44wxsw3\nxiyfR0DSMH8E/g3sEiTRpr6CCJJoBWAT4K4W2KdVpJDqKuzGmOG4aWv/zSccaZRsLPvQ7NuzshuY\nPozFve80DCPSIPX+ck+gODfkpBtpGN8FJMCXgd09hbFD9nidp/OLtLyaC7sxZltgmrX2oRzjkcY7\nGpgNnNLs6Y/Z+cYAT6Rh/Hgzzy3STvot7g+NMTfjpqUt7FhgPLDlAj/r0crGSrargvjLRfLwDVz7\n2I2r7DRyq3evqpzftPPeM20qZ/zrArYfMeYLV1XO/8T/u94XVcpFlXLhdHR09GrleE3LzI0xawO3\nAu9nPxoGvARsZK19tavXVSqVSm8DbFU+cxEk0UDAAssBn0/D+MUmnfdy3BDQBmkY39/5c70vqpSL\nKuWidjUNxVhrH7HWDrbWrmGtXQOYBoxaXFGX4simP44HPkWTuj8GSbQEsDXuRrs21RBpoLxmRujj\nUvlchlvOv3uQRJs04XyjgUHAdVnPeBFpkFwKu7X2s9baN/I4ljRHNv3xsOzbM5vQ/XGP7HFig88j\n0va08rSNpWF8B5DiFgyFjTpP1nt9O+BJ4D+NOo+IOCrschQwBzg1SKKlGnSOnXDtDC7VMIxI46mw\nt7k0jJ/DLTQbjmvx2wh7Zo+XN+j4IrKApk4l0vSlqiLlIkiiZYCncDc310nD+Jkcj/1F4FHg9jSM\nv7Wo5/jMxbjDJ3X5CWLyhO2aHlOR3he+KRe10xW7kIbxO7gbqZ8CLsq5j8xB2eNvczymiCyGCrt0\nSoDrgc2BKI8DZjdNvw+8mB1bRJpAhV0AyG5qHgC8CZwWJNFnczjsj4ClgfPTMJ6Xw/FEpAdU2OVj\naRi/AhyMK8ZpkEQDaj1WkESDgCOBt4Dz8olQRHpChV0+IQ3jy4FLgFHAmXUc6mBgeeCMNIxn5RGb\niPSMCrssykHAQ0AUJNF+vX1xkESrA8cAM4Hf5BuaiHRH0x09KXougiRaE7gLd9W9YxrGk3r4ug5g\nMm5f0z2yTwCLVfRcNJNyUaVc1E5X7LJIaRg/DXwX+AC4MkiinXr40kNxRf0W4IoGhScii6Erdk/K\nkosgib6Jm6o4ENd+4IysgdiinrsLbtPsmcD62c3YbpUlF82gXFQpF7VTYfekTLkIkmh93ObTQ4C/\nA0elYXzPAn++NK6/+zHAO8DoNIx73OyrTLloNOWiSrmonQq7J2XLRZBEKwMXANtmP3oSuB9YCreo\naTncJho7pmF8X2+OXbZcNJJyUaVc1E6F3ZOy5iJIos1x7Qe2wBV1cDtoXUKNUxvLmotGUC6qlIva\nqbB7UvZcBEnUHzc0MweYWU873rLnIk/KRZVyURLacbxKuahSLqqUiyrlonb96nmxMebHuH4gHwF/\nttYelUtUIiJSs5oLuzHmm8A2wLrW2rnGmJXyC0tERJrOGHOVMWaRGyd0RR+tqpSLKuWiSrmoUi5q\nV8/K07WArxtj7jbG/M0Ys2FeQYmISO0WOxRjjLkZN/NhYcdmr/20tXYTY8yXgauAPHp4i4iID8aY\nG40xmy/w/dPGmBV8xiQiIvUNxUwCvgVgjDHAEtba13OJSkREalbPdMeLgYuNMQ8DH+L2thQRERER\nERERERERERERKaWmdU4zxowBzgL6Ahdaa09t1rmLxBgzHLgUWBmoABdYa9t2w2djTF/gXmCatXZr\n3/H4ZIxZDrgQGIl7b+xjrb3bb1R+GGPGA98D5gMPA3tba+f4jao5jDEX47alfNVau072s+WBK4HV\ngOeBna21b3V1jKbseZr98p4LjAG+COxqjBnRjHMX0FzgMGvtSGAT4MA2zgXAIcBjuELW7s4Gplhr\nRwDrAo97jscLY8zqwL7AqKyw9QVCr0E11yW4Wrmgo4GbrbUGuDX7vkvN2sx6I+Bpa+3z1tq5QEJ1\nJ562Yq2dbq19IPv6Xdwv7yp+o/LDGDMMt/H1hTR5b4CiMcYsC2xmrb0YwFo7z1rb601LWsTbuAug\npYwx/XAburzkN6TmsdbeAby50I+3Af6Qff0HYLvFHaNZhX1V4MUFvp+W/aytZVcm6wP3dPPUVnUm\n8BPcx+12twYw0xhziTHmfmPM74wxS3X7qhZkrX0DOAN4AXgZeMtae4vfqLwbbK2dkX09Axi8uCc3\nq7DrY/ZCjDEDgauBQ7Ir97ZijBmHG0OcSptfrWf6AaOA86y1o4D36ObjdqsyxnwOOBRYHfdpdqAx\nZnevQRWItbZCNzW1WYX9JWD4At8Px121tyVjTH/gGuBya+0k3/F4simwjTHmOWAi8C1jzKWeY/Jp\nGu4G8n+y76/GFfp2tCFwp7X2dWvtPOBa3Pulnc0wxgwBMMYMBV5d3JObVdjvBdYyxqxujFkC2AW4\noUnnLhRjTAdwEfCYtfYs3/H4Yq09xlo73Fq7Bu7G2G3W2rZtS2GtnQ68mPVdAhgNPOoxJJ+eADYx\nxnwq+30ZjbvB3s5uAPbMvt4T16urS3VtjddT1tp5xpiDgL/i7nBfZK1tyzv+wFdx07geMsZMzX42\n3lr7F48xFYGG6+DHwBXZxc8zwN6e4/HCWvtg9untXtz9l/uBC/xG1TzGmInA5sCKxpgXgeOBU4Cr\njDE/IJvu6C9CEREREREREREREREREREREREREREREREREWlp/wc+Xq0gpSeWZgAAAABJRU5ErkJg\ngg==\n", "text": "<matplotlib.figure.Figure at 0x7fbd4f3af0d0>"}], "collapsed": false, "prompt_number": 58, "input": "gp = gaussian_process.GaussianProcess(corr=sklearn.gaussian_process.correlation_models.generalized_exponential,\n    theta0=[1e-2,10], thetaL=[1e-4,10], thetaU=[1e-1,10], nugget=1)\ngp.fit(X, y)\n\nx = np.atleast_2d(np.linspace(0, 10, 1000)).T\ny_pred, sigma2_pred = gp.predict(x, eval_MSE=True)\nsigma_pred = np.sqrt(sigma2_pred)\n\nplt.plot(X.ravel(), y, 's')\nplt.plot(x, y_pred)", "metadata": {"trusted": true}}, {"source": "PyMC has a very flexible GP implementation, but it has very spartan documentation, and is a bit too much for many common situations.\n\nI will use the PyMC Matern covariance in the sklearn GP code:", "cell_type": "markdown", "metadata": {}}, {"cell_type": "code", "language": "python", "outputs": [{"output_type": "pyout", "prompt_number": 59, "metadata": {}, "text": "matrix([[ 0.15555892],\n        [ 0.13196298],\n        [ 0.10330798],\n        [ 0.08970638],\n        [ 0.0771805 ],\n        [ 0.06590167]])"}], "collapsed": false, "prompt_number": 59, "input": "pymc.gp.cov_funs.matern.euclidean(X, np.array([0]), diff_degree=1.4, amp=0.4, scale=10)", "metadata": {"trusted": true}}, {"cell_type": "code", "language": "python", "outputs": [{"output_type": "stream", "stream": "stdout", "text": "[1.4, 0.4, 10]\n"}, {"output_type": "pyout", "prompt_number": 100, "metadata": {}, "text": "(6,)"}], "collapsed": false, "prompt_number": 100, "input": "def iso_matern(theta, d):\n    \"\"\"\n    Matern correlation model.\n\n        theta, d --> r(diff_degree=theta[0], amp=theta[1], scale=theta[2]; d)\n\n    Parameters\n    ----------\n    theta : array_like\n        An array with shape 3 (isotropic) giving the\n        autocorrelation parameters (diff_degree, amp, scale)\n    d : array_like\n        An array with shape (n_eval, n_features) giving the componentwise\n        distances between locations x and x' at which the correlation model\n        should be evaluated.\n    Returns\n    -------\n    r : array_like\n        An array with shape (n_eval, ) with the values of the autocorrelation\n        model.\n    \"\"\"\n    print theta\n    #return sklearn.gaussian_process.correlation_models.generalized_exponential(theta[:2], d)\n    theta = np.asarray(theta, dtype=np.float).reshape(3)\n    #theta[0] = 2\n    #theta[1] = .4\n    d = np.asarray(d, dtype=np.float)\n\n    assert theta.shape == (3,), 'expected theta.shape == (3,),  with diff_degree=theta[0], amp=theta[1], scale=theta[2]'\n\n    r = pymc.gp.cov_funs.matern.euclidean(\n        d, np.array([0]),\n        diff_degree=theta[0], amp=theta[1], scale=theta[2])\n    return np.array(r).reshape(d.shape[0],)\niso_matern([1.4,.4,10], X).shape", "metadata": {"trusted": true}}, {"cell_type": "code", "language": "python", "outputs": [{"output_type": "stream", "stream": "stdout", "text": "[[  1.5   1.   10. ]]\n[[  1.5   1.   10. ]]\n"}, {"output_type": "pyout", "prompt_number": 153, "metadata": {}, "text": "[<matplotlib.lines.Line2D at 0x7fbd4cf99250>]"}, {"output_type": "display_data", "metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEDCAYAAAAhsS8XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHtpJREFUeJzt3Xe8XkWdx/HPTTeA0gNIkAgMSwlNmtllKRukLkU5yUGy\nsDQ5IiEQFxdQqgVBibDAcoRABCmHHEpoCwQpipFqEFkQZkMoIUAaVQxJLnn2jznx3IR7b54+T/m+\nX6/74rlPnuc8P36Z/O7cmTkzICIiIiIiIiIiIiIiIiIiIiIiIiLS0DqKeZEx5lrgAGCutXZ49ty5\nwHHAvOxlZ1hr769FkCIiUrx+Rb5uEnAZcH2X5wrABGvthKpHJSIiZetTzIustY8B73XzR0X1+EVE\npH6KKuy9GGuMec4Yc40xZvWqRCQiIhWppLBfCQwDtgPeBi6uSkQiIlKRYsfYP8NaO3fZY2PMRODu\nlb2nUCh8Agws9zNFRNpRR0dHScPeZRd2Y8z61tq3s28PBZ4v4m0DSw2wVRUKhYJy4SgXOeUip1yU\nr9jljjcDuwNrA3OAc4A9cMMwBeBV4ARr7ZzerqO/qJxykVMucspFTrloEoVCoeA7hkahXOSUi5xy\nkVMuylfpqhgREWkwKuwiIg0qSKI1giS6otT3qbCLiDSgIIm+BPweOLHU96qwi4g0mCCJdgCeALYE\nLin1/SrsIiINJEii/YDfAUOAU9MwPrXUa5S9jl1ERKorSKKjgGuAJUCQhvFt5VxHhV1ExLMgiTqA\n04ALcRsu/msaxtPKvZ4Ku4iIR0ES9QF+BowHZgP7pGH8QiXXVGEXEfEkSKL+wLXAGOAlXFF/o9Lr\nqrCLiHgQJNEqQArsBzwJHJiG8fxqXFuFXUSkzoIkWgu4F9gFuB84LA3jj6t1fS13FBGpoyCJhgKP\n4Yr6DcBB1SzqoMIuIlI3QRJtirubdAtgAnBUGsZLqv05KuwiInUQJNGWuBuPNgJ+APxHGsZLa/FZ\nGmMXEamxbIuAqcBauLtJS94moBQq7CIiNRQk0QjgPmA14Pg0jCfW+jNV2EVEaiRIor2Au4BBwBFp\nGN9cj89VYRcRqYEgiQ4AbsMdQfqNNIzvrNdna/JURKTKgiQKgCnAUty+L3Ur6qDCLiJSVUEShUAC\nLMRtETC13jFoKEZEpEqCJDocd9PRX4GvpWH8pI841GMXEamCIIm+SQMUdVBhFxGpWJBERwC/Bj4C\n9vZZ1EGFXUSkIkESjQGuJy/qT3kOSYVdRKRcWVG/jryoP+05JECFXUSkLF166h8CIxulqINWxYiI\nlCxIotG4nvoHuKL+R88hLUc9dhGREgRJdDDLr35pqKIOKuwiIkULkmgfYDKwGNi/kYZfulJhFxEp\nQpBEe7D8NgHT/EbUM42xi4isRJBEXwXuAfoCB6dh/LDnkHqlwi4i0osgib6CO3B6EBCkYXyf55BW\nSoVdRKQHQRJtjTv5aDXcfup3eA6pKCrsIiLdCJJoE+A3wJrAMfU6JKMaNHkqIrKCIInWBx4EhgBj\n0zCe5Dmkkqiwi4h0ESTRGsADwDDg3DSML/ccUslU2EVEMkESDcatfhkOXA6c7zei8qiwi4gAQRL1\nB24FRgA3A+PSMC74jao8Kuwi0vaCJOoD/ArYD7gP+Pc0jJd6DaoCKuwi0taCJOoALgW+CTyOW6u+\n2G9UlVFhF5F2dxZwEvC/wIFpGH/sOZ6KqbCLSNsKkug44DzgNWCfNIzf9RtRdaiwi0hbCpJofyAG\nFuCK+lueQ6oaFXYRaTtBEu0IpMAS3E6N1nNIVaUtBUSkrQRJ9GXgXuBzwNfTMH7cc0hVp8IuIm0j\nSKK1cMsZ1wVOSsN4iueQakJDMSLSFoIk+hxwF2CAi9IwvsJzSDVTVI/dGHMtcAAw11o7PHtuTeAW\n4Eu4GeVR1tr3axSniEjZgiTqC9xIflfpGX4jqq1ie+yTgH1XeO504EFrrQEeyr4XEWlEE4BDgUeB\no5v5rtJiFFXYrbWPAe+t8PRBwHXZ4+uAQ6oYl4hIVQRJNBY4GXgBODQN40WeQ6q5SsbYh1hr52SP\n5+D2LRYRaRjZWvVLcDVq/zSM22K4uCqrYqy1BWNMUbugFQqFptwtrRaUi5xykVMucpXk4o33Z3PW\nQz+ns/Ap5+05fsima238OqOvrGZ4ddPR0dFRyusrKexzjDHrWWvfMcasD8wt5k2lBtiqCoVCQblw\nlIuccpGrJBdBEg0BngI2AkZttvawtKrBNbhKhmLuAo7KHh8FtOR6UBFpLtmyxjtxRf0HaRi3VVGH\n4pc73gzsDqxtjJkFnA38FJhsjDmWbLljrYIUESlGtq/6JGAX4NfAT/xG5EdRhd1ae3gPfzSyirGI\niFTqXGA08Hvg+GY9AalSdR3L0/hhTrnIKRc55SJXai6CJBqD66XPBHZJw3h+zYJrcCrsnigXOeUi\np1zkSslFkES7Ar8FFgK7pmH8Uk2Da3DaK0ZEmlqQRBsAt+OGlke1e1EH7e4oIk0sSKJBwB3A+sD4\nNIyneg6pIajHLiJNKTuEOgZ2Bq7H3WEqqLCLSPMah7uH5mnghHZdAdMdTZ56olzklIuccpHrLRdB\nEo0EHsDd8b5jGsaz6xpcg1OPXUSaSpBEmwCTgU7c0XYq6ivQ5KmINI0giVbDbRewBnBMK55XWg3q\nsYtIU8i2C7ge2Aq4LA3jSZ5Dalgq7CLSLE7HHejzCPBdz7E0NE2eeqJc5JSLnHKR65qLIIn2xk2W\nzgZ2SMN4ntfgGpx67CLS0IIk2gh3AHUncJiK+spp8lREGlaQRAOBW4G1gG+nYfyk55CagnrsItLI\nLgV2Aq4Dfuk5lqahwi4iDemRmX8AOAF4DjhRd5YWT5OnnigXOeUi1+65OHD8lAJAx+APGbjlE7C0\nD4teGEFh0WDumXBI2+alVOqxi0hj6buYAZs+S0efpSx+ZRsKiwb7jqjpqLCLSAMpMGCT5+kzaCFL\nZm/C0g/W9R1QU1JhF5GG0W/9mfRdfR6fvr82nbM39R1O09JyRxFpCEES7VYo/B+FxQNZPHMb6jwF\n2FLUYxcR74IkWge4GTpYPGM76BzgO6SmpsIuIl512dzri51vbsbSv67hO6Smp6EYEfHte8C+wP39\nh9oDppw2filo6WcltI7dE+Uip1zk2i0XQRL9E/Ao8A6wfdd9YNotF9WkoRgR8SJIorWBJPs21OZe\n1aOhGBGpu67j6sCZaRj/3nNILUU9dhHx4TRgP9we6xd6jqXlaIzdE+Uip1zk2iEXQRLtAkwD5gDb\n9TQE0w65qBX12EWkbrLDqG/C1Z4xGlevDRV2Eamny4EvAxemYfyI72BalQq7iNRFkETfBI4EngbO\n9hxOS9MYuyfKRU65yLVqLoIkGgb8CeiLG1efsbL3tGou6kHLHUWkpoIk6gfcCHwe+PdiirpURkMx\nIlJrZwFfxd2MdL3nWNqChmI8US5yykWu1XIRJNFuuC0DZuGGYN4v9r2tlot6Uo9dRGoiSKI1cEMw\nAEeUUtSlMirsIlIrVwBDgfPTMJ7mO5h2osIuIlUXJNFo4HDgCeDHnsNpOxpj90S5yCkXuVbIRZBE\nGwDPA4Nw4+r/V851WiEXvmi5o4hUTZBEHcBEYE3gO+UWdamMhmJEpJq+Rb5r45WeY2lbGorxRLnI\nKRe5Zs5FkESbAs8Bi4Gt0zCeXcn1mjkXvmkoRkQqFiRRX+A6YDBwXKVFXSqjoRgRqYbTgBHAZPLj\n7sQTDcV4olzklItcM+YiSKJtcTs2zgeGp2G8oBrXbcZcNAr12EWkbEESDQB+DfQHjq1WUZfKVDzG\nbox5DfgQ+BRYYq3dudJrikjT+D4wHLg6DeP7fAcjTjUmTwvAHtbad6twLRFpEkESbQecCbyJG2OX\nBlGtoRiNg4m0kSCJ+gPX4jqHx6dh/IHnkKSLahT2AvAbY8wzxpjjq3A9EWl83wO2B36VhvH9voOR\n5VXc0zbGrG+tfdsYsw7wIDDWWvtYd68tFAqFSj9PRPya9cFb/OfUC1htwCpcvN9ZrDpgFd8htbxS\nVwdVPMZurX07++88Y8wdwM5At4UdSg+wVWkpV065yDV6LrJj7v4A7PTeJx8ctNrAVe+u1Wc1ei4a\nWUWF3RgzGOhrrf3IGLMK8DXgvKpEJiJ1c+D4KT3+Nn3PhEO6FtfxwE7AjWkY16yoS2Uq7bEPAe4w\nxiy71o3W2qkVRyUiDSdIos2B84E5wDjP4UgvKirs1tpXge2qFIuINKhsL5hrgYHAiboRqbHpzlMR\nKcaJZHvBpGF8u+9gpHcq7CLSqyCJhgI/Ad4DxnoOR4qgbXtFpBcFgMuBVXF7wcz1G48UQz12EelR\nnzXmABwEPApM8hqMFE3b9nqiXOSUi1wj5SJIoi8Af8GdX7pNGsa2np/fSLloNuqxi0hPLgDWB35c\n76IulVFhF5HPCJJoBBABLwIXeg5HSqTCLiLLyQ7PuAo3VPutNIwXew5JSqTCLiIrOg3YCojTMJ7m\nOxgpnQq7iPxdkESbAWcB7wBneA5HyqTCLiIABEnUgVuzPhA4OQ3j9z2HJGVSYReRZb6O26F1KnCr\n51ikAlrH7olykVMucr5yESTRKsBLwLrA8EZY3qh2UT712EUE3Lj6hsDPGqGoS2VU2EXaXJBEWwDf\nBV7HbfYlTU6FXaSNdZkw7QeMS8P4b55DkipQYRdpb6OBvYB7gbs8xyJVoslTT5SLnHKRq2cugiRa\nDXgZt8nXVmkYv1KPzy2W2kX51GMXaV/n4jb5uqDRirpURoVdpA0FSbQ17kDqmcBFnsORKlNhF2kz\n2YTpFUBf3B2mCz2HJFWmwi7Sfo4A/hm4Mw3je30HI9Wnwi7SRoIkWh34ObAQNxQjLUiHWYu0l/OB\nIcD30zB+3XcwUhvqsYu0iSCJtgO+A1jgYs/hSA2psIu0gSCJ+uAmTPsAY9MwXuQ5JKkhFXaR9nAk\nMAK4NQ3jqb6DkdpSYRdpcUESrYFbq/43YLzncKQONHkq0vrOB9YBzkjDeJbvYKT21GMXaWHZhOmJ\nuD1hJngOR+pEhV2kRXUzYbrYc0hSJyrsIq2r64Tpg76DkfpRYRdpQdkdppowbVOaPBVpTT9EE6Zt\nSz12kRbTZcLUognTtlTX00mueSYp3D/j0R8A72Rfb2f/nZuGcWc9Y/FNp8PklItcpbnIJkwfw42t\n79PMNyOpXZSvrkMx9894FOBH3fxRIUiieXy24Hf3+KM0jAv1iFekCS2bML2tmYu6VKauPw1ffXdW\n4XtTf7wvsB7uSK71unn8+ZVcZiF5kZ8LzOvyNXeFx/MbdYmXeiM55SJXSS6yCVMLrAJskYbxG1UN\nrs7ULsrXcIdZB0m0Cm5b0d6K/7KvvkV87Ad8tuh390NgLjAvDeMlJf+PlUGNNqdc5Cos7JcBJwFn\npmF8QXUjqz+1i/I1XGEvVjaWuAZu5n8dYN0uj5f7vrB44HD6L6ajo6gRnAW4Ij9nha+uz80G3knD\n+NNy41ejzSkXuXJzkU2Y/hGYAWzTCrs3ql2Ur2kLeykOHD+lAAXou4SO/ovp6LcE+i+io99iBgx7\n8SzcD4F1cb8pLPtacyWX7QRmAW9kX69n/30Fd/v2W73NBajR5pSLXDm5yDo5vwe+SpNPmHaldlG+\nNlrH3gGfDqDw6QC6Vtv0jJO7m8wlSKL+uB7/skLftfBvCGwEfAl3dmR3je+vQRK9hCvyLwEvANOB\nNzT5K1V2HK6oa0teAdqqsJcmG2t/K/vqUZBEA1i+0G8GbA78AzAc2HGFt7wbJNH0G567gyCJRgHT\n0jCeXe34pT0ESTQEuBD4CJ1hKhkV9gplq25mZl/LCZKoL67YLyvyO2RfI+96aSrALdnrXgV+h1t/\n/DtgRq169W5Yqnv3TDhEv/Y2n4uB1YGT0zDutRMi7UOFvYayydVlRf9/lj0fJNEXzt3z1PfPfeQX\n/wnsBvwTcFT2BTAzSKL7gPuAR9Iw/lt9I5dmECTRSOAI4Bngvz2HIw2kjSZPu+erl9o1F9nk11a4\nIr8XsDf5ev5FwEPAZODONIzfr+RzGz0X7a7YXARJNAj4M7AJsFMaxtNrHlydqV2Ury0KeyPqLRfZ\nxO0IYD9gf9wwDsAS4AHcEM4daRh/XOrnqrA3thIK+7nAOcClaRifUvPAPFC7KF/FSTPG7AtcgrtZ\naKK19sKeXqu/qFwpuQiSaDMgAEYB22ZPfwjcBEwEphc7Jq/C3tiKvIlvc1xvfR7uDtOP6hJcnald\nlK+ipBlj+uKW843E3bTzNHC4tfYv3b1ef1G5Cm5E2RwYAxwNfDF7+k+4MdYb0jBe2Nv7Vdgb28py\nESRRB/Ab3JDdN9Iwvr1uwdWZ2kX5Kt22d2dghrX2NWvtEiABDq48LOlJGsYvp2F8Fm61zQHAHcDW\nwFXAG0ESnZctgZPW9G+4on4P7u9e5DMq7bEfBuxjrT0++34MsIu1dmx3r9dP4FyVt1fYALdHSITb\nZmERMAn4STMcsqB2kVvJ3Mt6wIvAAGDrNIxfq2ds9aZ2Ub5KlzuWvNa6UCjorstMtXPxSeciHn31\nce59+aGBcz6eH/Xr0y+a+MzNHLrFvqw5ePVqflTVqV3kesrFxdOu4sk3n+Xo7Uexn9nzVUZfWe/Q\n6k7twin1B1ylhX02MLTL90OBN3t7g34CO7XsjQRJ1A8Y07m08+wHZvx22AMzfrsIuBz4UaXLJWtB\nPbNcT7kIkugbwK3AtEnPTv7n/Tffa2n9o6svtYvyVToU0w83efovuFvvn0KTp0WpRy6yZZNHAsvG\n5OcDZwNXN9KJVWoXue5yESTRmrghmNWBbdMwftlLcHWmdlG+iiZPrbWduLHdB3AN75aeirrUXxrG\nS9Iwvga3pcEZwCDc6pk/BUm0l9fgpBQTcJvPndsuRV0qoxuUPPGRi2zy7YfAsbi/++uA76ZhvKCe\ncaxI7SK3Yi6CJNoXt7XEdGCXRvpNq9bULsqnwu6Jz1wESfQV4Gpge9xNLqcAN/vaTljtIrfCVhOf\nB54HNgB2TMP4Oa/B1ZnaRfkqXccuTSgN4z/i7kE4DVgVuBG4O+vRS+O4FLcd9AXtVtSlMuqxe9Io\nuQiSaBju5qaRuMnV49MwnlLPGBolF41gWS6CJDoEdwPSdGDXep3F20jULsqnwu5JI+Ui213yO8BF\nuAnWScC4eu1B0ki58K1QKBRG3fLt9YD/BVYDdkjD+EXPYXmhdlE+DcUIaRgvTcP4MuArwLO4fWie\nC5JoJ7+RtZ/sfpyJwNrA6e1a1KUyKuzyd1kR2RW4ANgYmBYk0cnZxlNSBw/PnAZwIPAw8F9+o5Fm\npaEYTxo9F0ES7Y2bVF0HN9Z7TK3uWm30XNRLkESbDeo30H7SuegDYHgz7PNTS2oX5VOPXbqVhvGD\nwHbAb4FDgekamqmdIIkGArd80rkI4MR2L+pSGRV26VF2OPJI3E1NG+OGZk7U0ExN/AzYfq9hI0jD\n+CbfwUhz01CMJ82Wi2xo5ibcpN51wLdXdqhHsZotF9UWJNGhwO3Ai7/+xqVbDuo/sG1z0VW7t4tK\nqMcuRcmGZnbAnZJ1FK73vrHXoFpAlsNrgYXAqIH9BvgNSFqCeuyeNGsugiQahFutcTzwLnB4GsZT\nK7lms+aiUtm4+u+AnRfP3JpP52/4mdf4Oq6wEbRru6gGFXZPmj0XQRIdB1wB9Ae+D/y03L1mmj0X\n5cjmKa4Gju2cvwFLZg6nu3+OKuzt1S6qRUMxUpY0jCcCu+EOW/kJcFu2aZUU5wTcLpvTl7y6FXXu\nY0mLU2GXsqVh/BTubtVHcEsinwqSaAu/UTW+IIn+ETecNR/4OoW+niOSVqPCLhVJw3gu8DXg58Dm\nuOJ+mN+oGleQRENxR9z1AUalYfy655CkBamwS8XSMO5Mw/g0YDRuTCENkuji7Gg+yWRDVfcC6+EO\nOHnEc0jSolTYpWrSMJ4M7AS8BIwHHgqSaH2/UTWG7IdcCgzHTTprHxipGa2K8aSVcxEk0WrANUAA\nvIMbcnisp9e3ci7g7ytgrgKOA+4BDu3piLtWz0UplIvyqccuVZft4z4a12tfB3gkSKLxbbwVwVm4\noj4dt+6/bc4tFT/UY/ekXXIRJNFuwGTcuHIKHLviAR6tnIsgiU4BfgG8BoxIw/jt3l7fyrkolXJR\nPvXYpaayIZgdgMdwQzNPB0m0ld+o6iO7iesXwNvAv6ysqItUi3rsnrRbLrLJw5/ihmc+AU4FfpmG\ncaEVcxEk0RjgemABsHuxJyG1Yi7KpVyUT4Xdk3bNRZBEB+M2vVoTmAIcN3n0lfNbKRdBEh0P/BL4\nENgrDePpxb63XdtFd5SL8mkoRuoqDeM7gW2BR4FDgOdemGu9xlRNQRKNw62AWQDsWUpRF6kW9dg9\nafdcBEnUFzgdOK+Djr4FChcB56Rh/Inn0MqSrfg5J/t6GxhZzkHU7d4uulIuyqfC7oly4QRJ9NV1\nV1nrD3M/XgDuxqZj0jB+3HNYJcm2350IjMGtftk7DeMZ5VxL7SKnXJRPhd0T5SK3cPHCwpG3n3oZ\nMBYoABOAs6p1QlMtBUm0Fu6w792AJ4CDs/1zyqJ2kVMuyqcxdvFuUP9BpGF8MrA7MBP4LvB8kET7\n+42sd0ES7Qo8iyvqk3ETpWUXdZFqUY/dE+Ui1zUXQRINxh2ePQ7oC9wNnJqG8SseQ1xOkER9gFOA\nC3ExngP8OA3jpZVeW+0ip1yUT4XdE+Ui110ugiTaGrgM2ANYBFwKXJSG8YL6R7hcXF/GjafvCcwB\nvpmG8cPVur7aRU65KJ8KuyfKRa6nXGQrTUbh9nrfELcu/OfAJStuS1BrQRINAE7C/TYxGPebxLfS\nMH6nmp+jdpFTLsqnwu6JcpFbWS6yA7S/DZwJrA28B8TA5WkYv1XL2LIfLofihl02xa1PHwsk5Z7x\n2hu1i5xyUT4Vdk+Ui1yxuci2Ax6Xfa0NLMFtLPYr4OE0jD+tVkxBEvUDDgP+A3f8Xyfuh8n5aRjP\nq9bnrEjtIqdclE+F3RPlIldqLoIk+hzwb7j9Zv4he3o2kOBOKJqWhvHiUuPIeudbA0dkXxvill/e\nBvwgDeOXS71mqdQucspF+VTYPVEucuXmIivEI4Ajcfu/fyH7o4+Ax4E/4vZAfw2YBcxbtnIl++Gw\nJrAJsBWwMzASV8zBjeffAEyo54octYucclE+FXZPlItcNXKRjcPvAewH7AuYHl76KbAU6O481gXA\nw7jhnbt9bG+gdpFTLsqnwu6JcpGrRS6CJFoDtw/8tsBGwFBgLVxB7wt8gJuEfQ14Efgz8OdqrEWv\nhNpFTrloEoVCoeqrCJqVcpFTLnLKRU65KJ+2FBARaTEq7CIiLUaFXUSkxaiwi4i0GBV2EZEWo8Iu\nItJiVNhFRFpMv3LfaIw5FzgOWLYh0hnW2vurEZRIPRw4fkqP66TvmXCIboyRplV2YSc7m9JaO6Fa\nwYiISOUqHYpRr0ZEpMFUWtjHGmOeM8ZcY4xZvSoRiYhIRXodijHGPAis180ffR+4Ejg/+/6HwMXA\nsVWNTkREStZrYbfW7l3MRYwxE3FnQPZKO7XllIuccpFTLnLKhQfGmPW7PD7VGHOTz3hERMSpZFXM\nhcaY7XCrY14FTqhOSCIiIiIiIiIiIiIiIiIi9VeX5UTGmH2BS3CHCE+01l5Yj89tRMaYocD1wLq4\nieerrLX/5Tcqf4wxfYFngDettf/qOx5fshv8JgJb4drFMdbaJ/xG5Ycx5gxgDLAUeB442lq7yG9U\n9WGMuRY4AJhrrR2ePbcmcAvwJdzh66Oste/3dp2a7+6Y/cO9HNgX2BI43BizRa0/t4EtAU611m4F\n7Ap8p83zMQ54EVfM2tmlwP9Ya7cAtgH+4jkeL4wxGwPHAztkha0vEHoNqr4m4WplV6cDD1prDfBQ\n9n2v6rFt787ADGvta9baJUACHFyHz21I1tp3rLV/yh7/FfcPeAO/UflhjNkQ2B/XU23bm1GMMV8A\ndrPWXgtgre201n7gOSxfPsR1fgYbY/oBg4HZfkOqH2vtY8B7Kzx9EHBd9vg64JCVXacehf2LwKwu\n37+ZPdf2st7J9sCTnkPx5RfAabhfudvZMGCeMWaSMWa6MeZqY8xg30H5YK19F7c9yRvAW8D71trf\n+I3KuyHW2jnZ4znAkJW9oR6Fvd1/xe6WMWZV4FZgXNZzbyvGmANx44jP0sa99Uw/YAfgv621OwAf\nU8Sv263IGLMJcAqwMe432VWNMUd4DaqBWGsLFFFT61HYZwNDu3w/FNdrb1vGmP7AbcAN1topvuPx\nZARwkDHmVeBmYC9jzPWeY/LlTdzk8dPZ97fiCn072hH4g7V2gbW2E7gd11ba2RxjzHrw961c5q7s\nDfUo7M8AmxljNjbGDABGA3fV4XMbkjGmA7gGeNFae4nveHyx1p5prR1qrR2Gmxx72Fp7pO+4fLDW\nvgPMMsaY7KmRwAseQ/LpJWBXY8znsn8rI3GT6+3sLuCo7PFRwEo7g5XsFVMUa22nMeYk4AHcDPc1\n1tq2nPHP/CNuKdefjTHPZs/pWEEN2Y0Fbsw6P68AR3uOxwtr7XPZb27P4OZepgNX+Y2qfowxNwO7\nA2sbY2YBZwM/BSYbY44lW+7oL0IRERERERERERERERERERERERERERERERERaXn/Dzt/g2VWKGjB\nAAAAAElFTkSuQmCC\n", "text": "<matplotlib.figure.Figure at 0x7fbd4cf99210>"}], "collapsed": false, "prompt_number": 153, "input": "gp = gaussian_process.GaussianProcess(corr=iso_matern,\n    theta0=[1.5,1,10], nugget=.001)\ngp.fit(X, y)\n\nx = np.atleast_2d(np.linspace(0, 10, 1000)).T\ny_pred, sigma2_pred = gp.predict(x, eval_MSE=True)\nsigma_pred = np.sqrt(sigma2_pred)\n\nplt.plot(X.ravel(), y, 's')\n\nplt.plot(x, y_pred)", "metadata": {"trusted": true}}, {"cell_type": "code", "language": "python", "outputs": [{"output_type": "stream", "stream": "stdout", "text": "[ 2.5  0.5  1. ]\n[ 25.    0.5   1. ]\n[ 2.5  5.   1. ]\n[  2.5   0.5  10. ]\n[ nan  nan  nan]\n"}, {"ename": "error", "evalue": "(diff_degree>0) failed for 2nd argument diff_degree: matern:diff_degree=nan", "traceback": ["\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31merror\u001b[0m                                     Traceback (most recent call last)", "\u001b[1;32m<ipython-input-155-9965b7b21862>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m gp = gaussian_process.GaussianProcess(corr=iso_matern,\n\u001b[0;32m      2\u001b[0m     theta0=[2.5,.5,1], thetaL=[1., .1, 1], thetaU=[10., 1., 1], nugget=y/10.)\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0mgp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      4\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0matleast_2d\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m10\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1000\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;32m/homes/abie/anaconda/lib/python2.7/site-packages/sklearn/gaussian_process/gaussian_process.pyc\u001b[0m in \u001b[0;36mfit\u001b[1;34m(self, X, y)\u001b[0m\n\u001b[0;32m    337\u001b[0m                       \"autocorrelation parameters...\")\n\u001b[0;32m    338\u001b[0m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtheta_\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreduced_likelihood_function_value_\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpar\u001b[0m \u001b[1;33m=\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 339\u001b[1;33m                 \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_arg_max_reduced_likelihood_function\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    340\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0misinf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreduced_likelihood_function_value_\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    341\u001b[0m                 raise Exception(\"Bad parameter region. \"\n", "\u001b[1;32m/homes/abie/anaconda/lib/python2.7/site-packages/sklearn/gaussian_process/gaussian_process.pyc\u001b[0m in \u001b[0;36m_arg_max_reduced_likelihood_function\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    726\u001b[0m                         optimize.fmin_cobyla(minus_reduced_likelihood_function,\n\u001b[0;32m    727\u001b[0m                                              \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlog10\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtheta0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mconstraints\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 728\u001b[1;33m                                              iprint=0)\n\u001b[0m\u001b[0;32m    729\u001b[0m                 \u001b[1;32mexcept\u001b[0m \u001b[0mValueError\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mve\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    730\u001b[0m                     \u001b[1;32mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Optimization failed. Try increasing the ``nugget``\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;32m/homes/abie/anaconda/lib/python2.7/site-packages/scipy/optimize/cobyla.pyc\u001b[0m in \u001b[0;36mfmin_cobyla\u001b[1;34m(func, x0, cons, args, consargs, rhobeg, rhoend, iprint, maxfun, disp, catol)\u001b[0m\n\u001b[0;32m    169\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    170\u001b[0m     sol = _minimize_cobyla(func, x0, args, constraints=con,\n\u001b[1;32m--> 171\u001b[1;33m                            **opts)\n\u001b[0m\u001b[0;32m    172\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0miprint\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m0\u001b[0m \u001b[1;32mand\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0msol\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'success'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    173\u001b[0m         \u001b[1;32mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"COBYLA failed to find a solution: %s\"\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0msol\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmessage\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;32m/homes/abie/anaconda/lib/python2.7/site-packages/scipy/optimize/cobyla.pyc\u001b[0m in \u001b[0;36m_minimize_cobyla\u001b[1;34m(fun, x0, args, constraints, rhobeg, tol, iprint, maxiter, disp, catol, **unknown_options)\u001b[0m\n\u001b[0;32m    244\u001b[0m     xopt, info = _cobyla.minimize(calcfc, m=m, x=np.copy(x0), rhobeg=rhobeg,\n\u001b[0;32m    245\u001b[0m                                   \u001b[0mrhoend\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mrhoend\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0miprint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0miprint\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmaxfun\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mmaxfun\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 246\u001b[1;33m                                   dinfo=info)\n\u001b[0m\u001b[0;32m    247\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    248\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0minfo\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m>\u001b[0m \u001b[0mcatol\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;32m/homes/abie/anaconda/lib/python2.7/site-packages/scipy/optimize/cobyla.pyc\u001b[0m in \u001b[0;36mcalcfc\u001b[1;34m(x, con)\u001b[0m\n\u001b[0;32m    236\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    237\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mcalcfc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcon\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 238\u001b[1;33m         \u001b[0mf\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    239\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mk\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mc\u001b[0m \u001b[1;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mconstraints\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    240\u001b[0m             \u001b[0mcon\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mc\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'fun'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0mc\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'args'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;32m/homes/abie/anaconda/lib/python2.7/site-packages/sklearn/gaussian_process/gaussian_process.pyc\u001b[0m in \u001b[0;36mminus_reduced_likelihood_function\u001b[1;34m(log10t)\u001b[0m\n\u001b[0;32m    698\u001b[0m             \u001b[1;32mdef\u001b[0m \u001b[0mminus_reduced_likelihood_function\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlog10t\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    699\u001b[0m                 return - self.reduced_likelihood_function(\n\u001b[1;32m--> 700\u001b[1;33m                     theta=10. ** log10t)[0]\n\u001b[0m\u001b[0;32m    701\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    702\u001b[0m             \u001b[0mconstraints\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;32m/homes/abie/anaconda/lib/python2.7/site-packages/sklearn/gaussian_process/gaussian_process.pyc\u001b[0m in \u001b[0;36mreduced_likelihood_function\u001b[1;34m(self, theta)\u001b[0m\n\u001b[0;32m    593\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    594\u001b[0m         \u001b[1;31m# Set up R\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 595\u001b[1;33m         \u001b[0mr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcorr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtheta\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mD\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    596\u001b[0m         \u001b[0mR\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0meye\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn_samples\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m*\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnugget\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    597\u001b[0m         \u001b[0mR\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mij\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mij\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mr\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;32m<ipython-input-100-f61aaf39ca0a>\u001b[0m in \u001b[0;36miso_matern\u001b[1;34m(theta, d)\u001b[0m\n\u001b[0;32m     31\u001b[0m     r = pymc.gp.cov_funs.matern.euclidean(\n\u001b[0;32m     32\u001b[0m         \u001b[0md\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 33\u001b[1;33m         diff_degree=theta[0], amp=theta[1], scale=theta[2])\n\u001b[0m\u001b[0;32m     34\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0md\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     35\u001b[0m \u001b[0miso_matern\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1.4\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m.4\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m10\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mX\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;32m/homes/abie/anaconda/lib/python2.7/site-packages/pymc/gp/cov_funs/cov_utils.pyc\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, x, y, amp, scale, symm, *args, **kwargs)\u001b[0m\n\u001b[0;32m    170\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    171\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mn_threads\u001b[0m \u001b[1;33m<=\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 172\u001b[1;33m             \u001b[0mtarg\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mC\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0msymm\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    173\u001b[0m         \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    174\u001b[0m             \u001b[0mthread_args\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mC\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mbounds\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mbounds\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0msymm\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mxrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn_threads\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;32m/homes/abie/anaconda/lib/python2.7/site-packages/pymc/gp/cov_funs/cov_utils.pyc\u001b[0m in \u001b[0;36mtarg\u001b[1;34m(C, x, y, cmin, cmax, symm, d_kwargs, c_args, c_kwargs)\u001b[0m\n\u001b[0;32m    166\u001b[0m                 \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcov_fun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mC\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcmin\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcmin\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcmax\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcmax\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0msymm\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msymm\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mc_args\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m**\u001b[0m\u001b[0mc_kwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    167\u001b[0m             \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 168\u001b[1;33m                 \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcov_fun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mC\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcmin\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcmin\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcmax\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcmax\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0msymm\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msymm\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0mc_args\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mc_kwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    169\u001b[0m             \u001b[0mimul\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mC\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mamp\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mamp\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcmin\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcmin\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcmax\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcmax\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msymm\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msymm\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    170\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;31merror\u001b[0m: (diff_degree>0) failed for 2nd argument diff_degree: matern:diff_degree=nan"], "output_type": "pyerr"}], "collapsed": false, "prompt_number": 155, "input": "gp = gaussian_process.GaussianProcess(corr=iso_matern,\n    theta0=[2.5,.5,1], thetaL=[1., .1, 1], thetaU=[10., 1., 1], nugget=y/10.)\ngp.fit(X, y)\n\nx = np.atleast_2d(np.linspace(0, 10, 1000)).T\ny_pred, sigma2_pred = gp.predict(x, eval_MSE=True)\nsigma_pred = np.sqrt(sigma2_pred)\n\nplt.plot(X.ravel(), y, 's')\n\nplt.plot(x, y_pred)", "metadata": {"trusted": true}}, {"cell_type": "code", "language": "python", "outputs": [], "collapsed": true, "input": "", "metadata": {"trusted": true}}], "metadata": {}}], "metadata": {"name": "", "signature": "sha256:6a6cca6d7efb2fd262107abdb260934c49367912b66bfdafd43003cf95946f26"}, "nbformat_minor": 0}