-
-
Save xdavio/c90347adda522ba784e4e37b31d36535 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "code", | |
"execution_count": 3, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stderr", | |
"output_type": "stream", | |
"text": [ | |
"/Users/jdavis/src/pear/examples/del/gen.py:171: RuntimeWarning: covariance is not positive-semidefinite.\n", | |
" cov=cov, size=size) > 0\n" | |
] | |
}, | |
{ | |
"data": { | |
"text/plain": [ | |
"<matplotlib.axes._subplots.AxesSubplot at 0x1130a9d68>" | |
] | |
}, | |
"execution_count": 3, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAk8AAAIJCAYAAAC8xtkfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuU5GV95/F3td09zfTUXGl6bgwjtwcTEkxkI0sQDOIF\nFhHhaI6uJyoh7hI3F3LirsRl15Pd5GTjSjbZYLKREMkasxEBRSJoDBHHWYMrJhETebg5DHNrmrnW\nzDDT3XbtH9WjnXHo7u+vu+pX1f1+nVOHrqrf9/c89avqng/P86vfU6nX60iSJGlmusrugCRJUicx\nPEmSJAUYniRJkgIMT5IkSQGGJ0mSpADDkyRJUkB3k/fvdRAkSWpvlbI70GkceZIkSQowPEmSJAUY\nniRJkgIMT5IkSQEzDk8pJYOWJEla8CpTLQycUjoduAU4HxijEbYeBW7MOT8+g/37bTtJktqb37YL\nmu5SBbcBN+WcHz72QErpAuBPgJ9sZsckSZLa0XRTcX2TgxNAzvlvm9gfSZKktjbdyNM/pJRuBx4A\n9gNV4Argm83umCRJUjuaLjz9PHA1cBGwFDgA3Afc0+R+SZIktaUpTxifA54wLklSe/OE8SAvPyBJ\nkhRgeJIkSQowPEmSJAVMd8L4rB188Muh7ZdcenGTeiJJkjR7jjxJkiQFGJ4kSZICDE+SJEkBhidJ\nkqQAw5MkSVKA4UmSJCnA8CRJkhRgeJIkSQowPEmSJAUYniRJkgIMT5IkSQGVer3ezP03deeSJGnW\nKmV3oNM48iRJkhTQ3ewG9t11b2j75ddeBUDtCw+G26q+7tJwjSRJUoQjT5IkSQGGJ0mSpADDkyRJ\nUoDhSZIkKcDwJEmSFGB4kiRJCjA8SZIkBRieJEmSAgxPkiRJAYYnSZKkAMOTJElSQKVerzdz/03d\nuSRJmrVK2R3oNI48SZIkBXQ3u4HD/+8boe0X/4sfB2D/3Z8Nt7Xsmjdy8KHN4boll/xkuEaSJC1M\njjxJkiQFGJ4kSZICDE+SJEkBhidJkqQAw5MkSVKA4UmSJCnA8CRJkhRgeJIkSQowPEmSJAUYniRJ\nkgIMT5IkSQGVer3ezP03deeSJGnWKmV3oNM48iRJkhTQ3ewGDn31a6Ht+//lTwDw/B/eHm7r5H97\nHYe//nfhusXn/xh7P3FnuG7F298SrpEkSZ3NkSdJkqQAw5MkSVKA4UmSJCnA8CRJkhRgeJIkSQow\nPEmSJAUYniRJkgIMT5IkSQGGJ0mSpADDkyRJUoDhSZIkKaBSr9ebuf+m7lySJM1apewOdBpHniRJ\nkgK6m93A2NBwaPvuwQEAjnzr2+G2+s59GbVaLVxXrVYLt3f0qe+E6xad8dJwjSRJag+OPEmSJAUY\nniRJkgIMT5IkSQEzPucppdQFrAF25pzHm9clSZKk9jXlyFNK6Y8n/vtK4HHgbuBbKaULWtA3SZKk\ntjPdtN2xr4X9BnB5zvmVwGXAf2tqryRJktrUTM95+m7O+QmAnPOOQJ0kSdK8Mt05T8tSSo8A/Sml\nnwX+DPgw8EzTeyZJktSGpgxPOedXpJQWAecBh4Fx4FHgj1vQN0mSpLYz7bftcs5Hga9NeugPm9cd\nSZKk9ua5S5IkSQGVer3ezP03deeSJGnWKmV3oNM48iRJkhQw4yuMF7X70JHQ9qv6+wCo1WrhtqrV\nauG652qHw3WnVBcXbm/kmWfDdb2nnRqukSRJc8uRJ0mSpADDkyRJUoDhSZIkKcDwJEmSFGB4kiRJ\nCjA8SZIkBRieJEmSAgxPkiRJAYYnSZKkAMOTJElSgOFJkiQpoFKv15u5/6buXJIkzVql7A50Gkee\nJEmSArqb3UCtVgttX61WC9Udqx0bfj5c1z1wcuH2RncNhet6Vg9y4P6/Ctctvfy17P/0feG6ZVdf\nGa6RJEkn5siTJElSgOFJkiQpwPAkSZIUYHiSJEkKMDxJkiQFGJ4kSZICDE+SJEkBhidJkqQAw5Mk\nSVKA4UmSJCnA8CRJkhRQqdfrzdx/U3cuSZJmrVJ2BzqNI0+SJEkB3c1u4OnhvaHtTx9YAcDQgUPh\ntgaX9lOr1cJ11Wq1cF3Rfu7/9H3humVXX1m4bt9d94brll97VbhGkqT5zpEnSZKkAMOTJElSQNOn\n7SRJklohpVQBPgKcBxwBrs85Pz3p+WuB/wCMA5/IOf9ekXYceZIkSfPF1cCinPOFwE3ALceeSCl1\nAb8JXApcCPx8SmllkUYMT5Ikab64CHgAIOf8MHD+sSdyzuPAy3LOB4GTaWSgkSKNGJ4kSdJ8sRTY\nP+n+2MSIE9AIUCmlNwN/D3wJiH9lHs95kiRJJXnioteHLqZ91lc+P90FPQ8A1Un3uyZGnL4n53wP\ncE9K6Q7gZ4A7In0AR54kSVJZKl2x2/Q2A1cApJQuAB499kRKqZpS+lJKqXfioUM0ThwPc+RJkiSV\nozLnK8PcA7w2pbR54v67U0pvA/pzzrellD4OfDmlNAJ8E/h4kUYMT5IkqRSVrrkNTznnOnDDcQ8/\nPun524DbZtuO4UmSJJVjZlNxbadSr4fO1Ypq6s4lSdKszfnc2Uw9eelVoZxw5oP3ltbXyRx5kiRJ\n5ZjjabtWaXp4qtVqoe2r1WqhumO1O/YdDNetXb6E3YeOhOtW9fcV7ufhrz0Srlv8E69gz+3xc9tW\nXvcODtz3+XDd0itfz75P3hOuW/7WN4drJEkLT2XuTxhvCUeeJElSObo685wnw5MkSSqHI0+SJEkB\nhidJkqSZqzhtJ0mSFGB4kiRJCnDaTpIkaea8VIEkSVKEF8mUJEkK6NC17QxPkiSpHI48SZIkzVyn\nnvNUqddDCxpHNXXnkiRp1kpLMFve+u5QTtj4yT9pi7TlyJMkSSqH03YnVqvVQttXq1UAtuzeF25r\n46rl4faOtbl9b7xu3Ypqy/vZ6rrdh46E61b19xVuT5K0cHiFcUmSpIgOPefJ8CRJkspheJIkSQqY\nr9N2KaU3AZcBy4B9wCbgUzlnv0knSZIK69RLFUwZnlJKtwJdwP1ADagClwOvB65veu8kSdL8NU+/\nbXduzvmS4x67N6W0uVkdkiRJC0SHLs8yXa+7UkqvmvxASuliYLR5XZIkSQtCpRK7tYnpRp7eBdyS\nUvpzGlcgHQe+Afxck/slSZLmucp8nLbLOT8FvKlFfZEkSQtJG40mRUx3wvjfAItO9FzO+cKm9EiS\nJC0M8/RSBe8HPgq8GRhrfnckSdJC0anLs1Tq9akv15RSeh/wZM75ngL791pQkiS1t9Lmzra991dD\nOWH9rf+9Leb5pr1IZs75Q63oiCRJWmDm4zlPc2F011Bo+57VgwA8smV7uK1XbFxHrVYL11WrVXbs\nOxiuW7t8Cf+0Yzhc90NrBwr3s1PqRnfsCtf1rF3NgQf+Oly39A2vCddIktpAh07bubadJEkqxbxc\nnkWSJKlpDE+SJEkBL3nJnO4upVQBPgKcBxwBrs85Pz3p+bcBv0RjpZRHc84/X6SdzpxslCRJHa/S\nVQndZuBqYNHEtShvAm459kRKqQ/4deCSnPOrgOUppSuL9NvwJEmSytHVFbtN7yLgAYCc88PA+ZOe\nOwpcmHM+OnG/m8boVJjTdpIkqRxzf87TUmD/pPtjKaWunPN4zrkODAOklH4B6M85f7FII4YnSZJU\niiZ82+4AUJ10vyvnPH7szsQ5Ub8NnAVcU7QRw5MkSSrH3F/naTNwJfCplNIFwKPHPf9HwAs556tn\n04jhSZIklWPuR57uAV6bUto8cf/dE9+w6wceAd4NbEop/Q2NJeR+N+f8mWgjhidJklSOOQ5PE+c1\n3XDcw49P+nlOco/hSZIklaLSocuzVOr10ILGUU3duSRJmrXSLvO964O/FcoJqz/4/ra4JLkjT5Ik\nqRwzu/Bl22l6eKrVaqHtq9XGNwy37N4XbmvjquXh9o61OXzwhXDdwJKT2Ln/YLhuzbIlHHns8ek3\nPE7fOWczNvx8uK574GSOPvWdcN2iM17K6Lbt4bqe9esY3TUUr1s9yMEvfSVct+TVF7H/3s+F65Zd\ndUW4RpI0h1zbTpIkaeY69Zwnw5MkSSpHxfAkSZI0c57zJEmSNHNNWJ6lJQxPkiSpHE7bSZIkBTht\nJ0mSFOC0nSRJ0sxVHHmSJEkK8JwnSZKkAKftJEmSAjp02q5Sr4cWNI5q6s4lSdKslZZgnr/1o6Gc\ncPJ7f64t0pYjT5IkqRye83RiY0PDoe27BwcAeGxnrA7gnDUD1Gq1cF21WuW52uFw3SnVxTw+tDtc\nd/bgqpb3s2h7La/7woPxutddyu4/+li4btV73sW+u+4N1y2/9qpwjSTpBDp02s6RJ0mSVAqXZ5Ek\nSYowPEmSJAV0ec6TJEnSzDnyJEmSNHOe8yRJkhThtJ0kSVKAI0+SJEkBjjxJkiTNXMWLZEqSJAU4\nbSdJkhTQoWvbVer10ILGUU3duSRJmrXShn/2/cXdoZyw/KevaYuhKkeeJElSOZy2O7FarRbavlqt\nAvD40O5wW2cPruLJoT3hujMHV/Jc7XC47pTqYp56bm+47oxTVoSPCzSOTavrxoaGw3XdgwOMbtse\nrutZv46RZ54N1/Wedip7//yucN2Kt13LkX98LFzX98PnsP/uz4brll3zxnCNJM1rHTpt58iTJEkq\nh9+2kyRJmrm5Xp4lpVQBPgKcBxwBrs85P33cNouBLwDX5ZwfL9JOZ46XSZKkztdVid2mdzWwKOd8\nIXATcMvkJ1NKrwAeAk6fVbdnUyxJklRYV1fsNr2LgAcAcs4PA+cf93wvjYAVP+F1crdnUyxJklRY\npSt2m95SYP+k+2Mppe8V5py/mnPeziwvz+A5T5IkqRRzfc4TcACoTrrflXMen+tGHHmSJEnlmPtz\nnjYDVwCklC4AHm1Gtx15kiRJ5Zj7kad7gNemlDZP3H93SultQH/O+bZJ281qBRTDkyRJKsccXyQz\n51wHbjju4R+4HEHO+dLZtGN4kiRJpah4kUxJkqSADl3brlKvz2rabzpN3bkkSZq10hLMwQe/HMoJ\nSy69uC3SliNPkiSpHDO78GXbaXp4eq52OLT9KdXFANRqtXBb1WqVJ4f2hOvOHFzJyNZt4breDesZ\neXpLvO70jRx96jvhukVnvJShA4fCdYNL+wsfz06pG3nm2XBd72mnFq8r+nkp2J4kzUsdOm3nyJMk\nSSqHJ4xLkiTNXGWOL1XQKoYnSZJUjvk4bZdSekvO+c6UUj/wQeDlwCPAf805H2xB/yRJ0nzVodN2\n042XHbtK5+8Ce4FfBLYBf9TMTkmSpAWg0hW7tYmZTtudlXO+fuLnb6eUrmlWhyRJ0sLQqVcYny7G\nnZ1SuhEYTSn9GEBK6Xygt+k9kyRJ81ulEru1ienC05XAARqL6v1oSmkZ8PvALzS7Y5IkaZ7r0PA0\n5bRdzvnvgb8H/njSwxc0tUeSJGlBqMzHK4ynlP4GWHSi53LOFzalR5IkaWGYj+EJeD/wUeDNwFjz\nuyNJkhaMNpqKi6jU61MvaJxSeh/wZM75ngL7D62WLEmSWq60BPPCo/8Yygkn/cgPt0XamjY8zZLh\nSZKk9lZaIDnyrW+HckLfuS9ri/DU9OVZarVaaPtqtQrA6Pad4bZ61q0Jt3eszZEtW8N1vRs3cOQf\nHwvX9f3wOYXbK/r6Wl235/CRcN3KxX1s2b0vXLdx1fLC/RwbGg7XdQ8O8FztcLjulOriwv3c/+n7\nwnXLrr4yXCNJLdWh03aubSdJksrRoRfJNDxJkqRyOPIkSZI0c5U2Wq8uwvAkSZLK4bSdJElSwDy9\nSKYkSVJTVDznSZIkKcCRJ0mSpABHniRJkgIMT5IkSTNX8dt2kiRJAV7nSZIkKaBDp+0q9XpoQeOo\npu5ckiTNWmkJZnTXUCgn9KwebIu01fSRp+gq8tVqtVDdsdpWr3Y/snVbuK53w/qOeX1F6/YcPhKu\nW7m4j92H4nWr+vsY3bY9XNezfh0v/MO3wnUnnXcuI1u2hut6N25gbGg4XNc9OMCBz30hXLf0itfx\n/K0fDded/N6fC9dIUhEuzyJJkhThCeOSJEkz90LfotD21WmeTylVgI8A5wFHgOtzzk9Pev6NwM3A\nKPAnOefbQh2Y0JnjZZIkST/oamBRzvlC4CbglmNPpJS6J+5fBrwaeE9KaaBII4YnSZI0X1wEPACQ\nc34YOH/Scy8Dnsg5H8g5jwJfAS4u0ojhSZIkzRdLgf2T7o+llLpe5LkasKxII4YnSZI0Xxzgn58a\n1ZVzHp/03NJJz1WBfUUaMTxJkqT5YjNwBUBK6QLg0UnPfRs4M6W0PKXUS2PK7qtFGvHbdpIkab64\nB3htSmnzxP13p5TeBvTnnG9LKf0K8AUaFwa9Lee8s0gjhidJkjQv5JzrwA3HPfz4pOf/EvjL2bbj\ntJ0kSVKA4UmSJCnAaTtJklSK0Zf0lN2FQir1emhB46im7lySJM1aaQvM7T50JJQTVvX3tcVieE0f\neXqudji0/SnVxQDUarVwW9VqleGDL4TrBpacVHi1+6L9PPrUd8J1i854KVt2xy9JsXHV8sL9HN01\nFK7rWT0Yft+h8d4X7Wfh17dte7iuZ/26lvdz75/fFa5b8bZrOfC5L4Trll7xOvZ/+r5w3bKrrwzX\nSFrYxps7gNM0TttJkqRSNHn2q2kMT5IkqRSGJ0mSpACn7SRJkgI6NDsZniRJUjmctpMkSQoY79Ar\nGhmeJElSKb47Pl52FwoxPEmSpFKMjzvyJEmSNGMdesqT4UmSJJXDE8YlSZICPGFckiQpoFNHnipN\n7nhnHhVJkhaOSlkNP7ZzOJQTzlkzUFpfJ2v6yNOew0dC269c3AdQePX5Lbv3hes2rlpeuL3RHbvC\ndT1rVzM2NByu6x4cYMe+g+G6tcuXFH59wwdfCNcNLDmpcHtDBw6F6waX9jOydVu4rnfD+sJ1o9u2\nh+t61q8r3N6u//Lb4brVN/979n7iznDdire/hX13fiZct/wtb2LPHX8erlv5zreFayTNDx36ZTun\n7SRJUjk6ddrO8CRJkkpheJIkSQoYNzxJkiTNnOFJkiQpwGk7SZKkAEeeJEmSAjo0OxmeJElSOZy2\nkyRJCnDaTpIkKcCRJ0mSpIAOzU6GJ0mSVI5OnbarNHnIrDOPiiRJC0elrIa/kreEcsJFaWNpfZ2s\n6SNPtVottH21Wi1Ud6y2aN3I01vCdb2nb2x5P4vWDR04FK4bXNpfuL3RHbvCdT1rVxdub2xoOFzX\nPThQuK7V71/R4zm6bXu8bv26wnVFXx/A6K6hWHurB8NtSWovnTry5LSdJEkqheFJkiQpoBXftksp\n9QEfB04BDgDvzDnvPsF2A8BXgB/JOY9Mtc8pw1NK6U3AZcAyYB+wCfhUzrkzo6IkSWobLbpUwQ3A\nN3POv55S+mngZuCXJ2+QUnod8FvAjM4HeNHwlFK6FegC7gdqQBW4HHg9cH2R3kuSJB0z3pqhmIuA\n/zbx8/00wtPxvgu8BnhkJjucauTp3JzzJcc9dm9KafNMdixJkjSVuR55SildB9zI97/tXwF2Afsn\n7teApcfX5Zz/eqJ+Rt/mmyo8daWUXpVz3jSpUxcDozPZsSRJ0lTmOjzlnG8Hbp/8WErpLhqzZ0z8\nd99UXZpJO1OFp3cBt6SUPkEjuY0Dfwf8wkx2LEmSNJXx1lwOcjNwBfD1if9ummLbWY88/RDwcmAE\n+EDO+f8ApJQeBC6dyc4lSZJeTItOGP8D4I6U0ibgKPB2gJTSjcATOef7JndpJjucKjx9ADgPeAlw\nZ0ppUc75Dkq8EqkkSZo/WnHCeM75BeCtJ3j8d07w2Okz2edU4Wkk57wPvnfJggdTSltxyRVJkjQH\nxlv0dbu51jXFc1tSSreklPpzzjXgGuBW4JzWdE2SJM1n9Xo9dGsXU4Wn64BvMjHSlHN+Fvgp4JMt\n6JckSZrnOjU8VZrcmfZ5pZIk6URKO5f5nq9/K5QT3nz+uW1x3nXT17Yb2bottH3vhvUALV8Nvmjd\nyDPPhut6TzuV3YeOhOtW9fcVXrW+1XUjW7aG63o3buDJoT3hujMHV7LncPx4rlzcx9Y9+6ff8Dgb\nVi5jx76D4bq1y5ewc3+8bs2yJS1//8aGhsN13YMDhdsDwrXfq/vrh2J1rzn+2r+SytJOo0kRLgws\nSZJK0aHZyfAkSZLKMd6h6cnwJEmSSuG0nSRJUoDhSZIkKcBpO0mSpADDkyRJUoDTdpIkSQEdurSd\n4UmSJJXDkSdJkqQAw5MkSVKAJ4xLkiQFdGh2otLkIbMOPSySJC0YlbIavu3Bh0M54fpLX1laXydr\n+sjTyJatoe17N25o1G3dFm6rd8N6RrdtD9f1rF9XeBX50R274u2tXV149fndh46E61b19xVur2jd\n6Pad4bqedWvYczj++lYuLv76irZX9H0oWjd88IVw3cCSk1r+vhetA8K1x+oObvq/obolr7oQgP13\nfzZUt+yaN4a2lzQ9p+0kSZICPGFckiQpYGx8vOwuFGJ4kiRJpXDkSZIkKaBDs5PhSZIklcMTxiVJ\nkgKctpMkSQowPEmSJAU4bSdJkhTQmdHJ8CRJkkriyJMkSVKA5zxJkiQFjI93ZniqNDn1deZRkSRp\n4aiU1fBv3vPFUE74tTdfVlpfJ2v6yFPRldJbvar7kW/ncF3fyxJjw8+H67oHTmZ02/ZwXc/6dYVf\n3+iuoXh7qwcLtzeyZWu4rnfjhpa/70XrtuzeF67buGo5z9UOh+tOqS4u3M+xoeFwXffgACNPbwnX\n9Z6+sfDvAxT/W7HvL+4O1S3/6WsadZ+8J1b31jc36u66N1Z37VWh7aWFxHOeJEmSAjozOhmeJElS\nSVpxwnhKqQ/4OHAKcAB4Z85593Hb3Aj8NI0897mc83+Zap9dTeqrJEnSlMbr9dCtoBuAb+acLwb+\nN3Dz5CdTSi8F3pZzviDn/C+B16eUzp1qh4YnSZJUinq9HroVdBHwwMTP9wOXHff8VuANk+73AEem\n2qHTdpIkqRRzfcJ4Suk64Ea+fzpVBdgF7J+4XwOWTq7JOX8X2DNR/yHgGznnJ6dqx/AkSZJKMden\nPOWcbwdun/xYSukuoDpxtwr8wNelU0qLJur2Az8/XTuGJ0mSVIoWXWF8M3AF8PWJ/246wTb3Al/M\nOX9oJjs0PEmSpFK06DpPfwDckVLaBBwF3g7f+4bdEzSy0KuAnpTSFTSm/G7KOT/8Yjs0PEmSpFK0\nIjzlnF8A3nqCx39n0t3FkX0aniRJUilcGFiSJCnA8CRJkhQw3pnZiUqTU1+HHhZJkhaMSlkN/8qf\nfiaUE275mTeV1tfJmj7yVHSl9KKryLd69fnRHbvCdT1rV7f89RWtG9m6LVzXu2E9I1u2xus2bmDP\n4Skv6npCKxf3dczxbHXdkX98LFzX98Pn8MLfPxquO+nlP1L49wGK/63Y/+n7QnXLrr6yUXf3Z2N1\n17wRgH133RuqW37tVQAc/NJXQnVLXn1RaHupEzltJ0mSFNCiSxXMOcOTJEkqhSNPkiRJAZ16wrjh\nSZIklWK8Pl52FwoxPEmSpFJ06Kyd4UmSJJXDc54kSZIC/LadJElSgCNPkiRJAYYnSZKkAC9VIEmS\nFODIkyRJUsA4nRmeKk1OfZ15VCRJWjgqZTX8M7f+WSgn/Ol7/3VpfZ2s6SNP2/fGVkpft6KxUvqe\nw0fCba1c3MfO/QfDdWuWLWHrnv3hug0rlzF88IVw3cCSkxjdNRSu61k9GF55Hhqrz48NDYfrugcH\nGN22PVzXs34dI09vCdf1nr6R52qHw3WnVBcXPi5F64p+Pou2V/R9aPXrK/q+A4xs2Rqr27ihUffM\ns7G6004FCB+barU6q7ro34qBJScBMDb8fKiue+Dk0PZSmcY79KQnp+0kSVIpPOdJkiQpoEMHngxP\nkiSpHPNu5Cml9Jac850ppX7gg8DLgUeA/5pzjp9YJEmSNEm9Q79X1jXFczdM/Pd3gb3ALwLbgD9q\ndqckSdL8N16vh27tYibTdmflnK+f+PnbKaVrmtkhSZK0MHTqtN1UI09np5RuBMZSSj8GkFI6H+ht\nSc8kSdK8Nl6P3drFVOHpSmA/8BjwoymlZcDvA+9rRcckSdL8Vq/XQ7d2MdW03anAfwZGgU055/3A\nBSmlB4FLW9E5SZI0f7VTIIqYauTpAzS+YfdK4D0ppXdOPN4Wl0aXJEmdbT6eMD6Sc94LkFJ6E/Bg\nSmkrrlcnSZLmQDsFooipwtOWlNItwM0559rEt+w+DyxvTdckSdJ81qlr21VebL4xpdQNvAP4ZM75\n8MRjg8BNOedfnuH+O/OoSJK0cJR2Os4bfvN/hXLCA7/2b9ri1KEXHXnKOY8BHzvusSFgpsEJKL4C\n+e5D8VXdV/X3hVcuh8bq5UVXny9aN7p9Z7iuZ92a4u3t2BVvb+3q8Iru0FjVvWg/x4aG4+0NDnD0\nqe+E6xad8dLC7RV9/4489ni4ru+cszmanwzXLUpnFn9927aH63rWr+PoE0+F6xaddQYAI1u3hep6\nN6wHCL/G7sEBoPjfplbXRX93e9auBuDAfZ8P1S298vWh7aW50KknjLu2nSRJKsV8POdJkiSpaVox\n8pRS6gM+DpwCHADemXPefdw27wXeCYwDH8453znVPqe6VIEkSVLT1OuxW0E3AN/MOV8M/G/g5slP\nppRWAf8GuAC4DPjwdDs0PEmSpFK06DpPFwEPTPx8P42A9D0To1AvzzmPA2uAaU+edtpOkiSVYq6n\n7VJK1wE38v1v+1eAXTSWmwOoAUuPr8s5j09M3X0Q+L3p2jE8SZKkUnzpg/9uTi89kHO+Hbh98mMp\npbuA6sTdKrDvRWpvTSn9L+CBlNKXc84PvVg7TttJkqT5bDNwxcTPVwCbJj+ZUjp7ImABfBc4SuPE\n8RflyJMkSZrP/gC4I6W0iUYwejtASulG4Imc830ppX9IKX2VRmi6P+e86cV3Z3iSJEnzWM75BeCt\nJ3j8dyb9/OvAr890n07bSZIkBRieJEmSAgxPkiRJAZUmXxq9MxetkSRp4ZjTywUsBE0/YXzn/oOh\n7dcsWwJdx24yAAALy0lEQVTEVyCHxirkwwenvTDoDxhYchLb98bbW7eiWrifrX59re5nq+t2HzoS\nrlvV39cx78PRJ54K1y0664yWvw9H85PhukXpTABGnt4Squs9fWOjbuu2WN2G9QDhz8yq/j4g/rep\nWm1cXmbHvtjfwrXLi/0tPNbe0G/891Dd4Ad+FYBDX/nbUF3/RReEtpfmA6ftJEmSAgxPkiRJAYYn\nSZKkAMOTJElSgOFJkiQpwPAkSZIUYHiSJEkKMDxJkiQFGJ4kSZICDE+SJEkBhidJkqQAw5MkSVJA\npV6vN3P/Td25JEmatUrZHeg03c1uYM/h2MrlKxc3Vi6PrrAOjVXWW76KfMHV7seGnw/XdQ+cXLif\nI888G67rPe1Uhg++EK4bWHJSy9+H3YdinzOAVf19hetafVyOfOvb4bq+c1/G6Lbt4bqe9eta/vsA\nhGuP1UX7Wq1WZ1U3dOBQqG5waf+s6or2c/dtfxqqW3X9zwCw/97PheqWXXUFAHtu/3iobuV17wht\nL7UTp+0kSZICDE+SJEkBhidJkqQAw5MkSVKA4UmSJCnA8CRJkhRgeJIkSQowPEmSJAUYniRJkgIM\nT5IkSQGGJ0mSpADDkyRJUkClXq83c/9N3bkkSZq1Stkd6DTdzW6g6Irgo7uGwm31rB5k+974avDr\nVlQLryI/snVbuK53w/rC7XVKXdH3r9X93H3oSLhuVX9f4fbGhobDdd2DAxx96jvhukVnvLTl70PR\n3weA0e07Q3U969YAMDb8fKiue+BkoPjfpj2HY5+ZlYv7ZtVe0boD9/9VqG7p5a8FYN9d94bqll97\nFQD7P31fqG7Z1VcCsOdjnwjVrXzX20PbS83gtJ0kSVKA4UmSJCnA8CRJkhRgeJIkSQowPEmSJAUY\nniRJkgIMT5IkSQGGJ0mSpADDkyRJUoDhSZIkKcDwJEmSFGB4kiRJCqjU6/Vm7r+pO5ckSbNWKbsD\nnaa72Q3s2HcwtP3a5UuA+Eri0FhNPLriOTRWPS/aXvT1QeM1Fm2vaF105XlorD4/umNXuK5n7WrG\nhobj7Q0OFG6v1cez1XVFP9fb98bbW7eiynO1w+G6U6qLC78+aN3fimPtRX8nugdOnlV71s1t3cgz\nz4bqek87NbS9NBWn7SRJkgIMT5IkSQGGJ0mSpIBpz3lKKb0JuAxYBuwDNgGfyjl7MrgkSVpwpgxP\nKaVbaYxO3Q/UgCpwOfB64Pqm906SJKnNTDfydG7O+ZLjHrs3pbS5WR2SJElqZ9Od89SVUnrV5AdS\nSpcAo83rkiRJUvuabuTpXcAtKaVP0LiI1knA13HKTpIkLVDTjTwtonGV8C8C7wYOAmcB5zS5X5Ik\nSW1pupGnPwRuBk4D7gTOBo7QOIH8vuZ2TZIkqf1MF566cs4PAaSULs05Pzfx81jTeyZJktSGpgtP\nOaV0G/CenPO7AFJK7wfiC5BJkiTNA5V6/cWvdZlS6gLemHP+zKTH3gHcnXOeycqhXkhTkqT2Vim7\nA51myvA0B+qjO2KDVD1rVwMUXkW+6KruY0PD4bruwQF2H4r3c1V/8X52St3orqFwXc/qQY489ni4\nru+csxnZui1c17thffHXt217uK5n/bri7W3fGW9v3ZrCv0dF64oeF4CRZ54N1fWedmqjLvje925Y\nDxB+L6rV6oKoi773Kxf3AfBcbSb/P/19p1QXz6q9Z/ceCNWdumIpAAcfil2mcMklPxnavkMZnoJc\n206SJCnA8CRJkhRgeJIkSQowPEmSJAUYniRJkgIMT5IkSQGGJ0mSpADDkyRJUoDhSZIkKcDwJEmS\nFGB4kiRJCjA8SZIkBTR9YeBm7lySJM2aCwMHdTe7gaIrew8dOBRua3Bpf+FV61tdt3XP/nDdhpXL\n2Ln/YLhuzbIlbN8b7+e6FVWeHNoTrjtzcCXDB18I1w0sOYmjTzwVrlt01hmM7tgVrutZu7pjPi9F\nj2fR36NW//4B4c/ouhWNuqJ/Y0a2bA3V9W7cMKv25nvdyDPPhup6Tzt1du1t3RZrb8N6APZ/9oFQ\n3bI3vgGAF77xD6G6k378vND26ixO20mSJAUYniRJkgIMT5IkSQGGJ0mSpADDkyRJUoDhSZIkKcDw\nJEmSFGB4kiRJCjA8SZIkBRieJEmSAgxPkiRJAYYnSZKkgEq9Xm/m/pu6c0mSNGuVsjvQabqb3UDR\nFbOf3Xsg3NapK5YWXtX96FPfCdctOuOlhdsrumr9nsNHwnUrF/ex+1C8blV/Hzv3HwzXrVm2pPBx\nia5cDo3Vy6MrukNjVfejTzwVrlt01hnhFd2hsap70eOyY1/8fVi7fEnh36Mnh/aE684cXFn49UHx\nvxXzvW50+85QXc+6NQDhz8za5UsA2L431s91K4r9zT51xdJZtff40O5Q3dmDqwA4/LVHQnWLf+IV\nAOy7695Q3fJrrwLg4Je+Eqpb8uqLQturHE7bSZIkBRieJEmSAgxPkiRJAYYnSZKkAMOTJElSgOFJ\nkiQpwPAkSZIUYHiSJEkKMDxJkiQFGJ4kSZICDE+SJEkBhidJkqSASr1eb+b+m7pzSZI0a5WyO9Bp\nHHmSJEkK6G52A7VaLbR9tVotVHestmjdyDPPhut6Tzu15f3slLqx4efDdd0DJ3fM6xvdNRSu61k9\nyOj2nfG6dWs65rgUrYPW/a04Vjfy9JZQXe/pG2fVXtG66O9S98DJs2pvvtftuf3jobqV172j0d5f\nPxRr7zWXNOq+8GCs7nWXAnDg/r8K1S29/LWh7TU7jjxJkiQFGJ4kSZICDE+SJEkBhidJkqQAw5Mk\nSVKA4UmSJCnA8CRJkhRgeJIkSQowPEmSJAUYniRJkgIMT5IkSQGVer3ezP03deeSJGnWKmV3oNM0\ne2Fg3xBJkjSvOG0nSZIUYHiSJEkKMDxJkiQFGJ4kSZICDE+SJEkBhidJkqSAZl+q4AeklCrAR4Dz\ngCPA9Tnnp1vdj3aUUnoE2D9x9zs5558tsz9lSym9EvitnPNPpZTOAD4GjAPfyjm/t9TOlei44/Jy\n4D7g8Ymn/yDnfGd5vWu9lFI3cDuwEegFfgP4J/y8vNixeRY/M13AR4FE4zPyb4Gj+JnRDJUx8nQ1\nsCjnfCFwE3BLCX1oOymlRQA550snbgs9OL2Pxh+3RRMP3QL8Ws75EqArpfSm0jpXohMcl1cAH570\nuVlQ/whOeAfwfM75YuANwO/j5+WYycfmchrH5sfxM/NGoJ5zvgi4GfhN/MwooIzwdBHwAEDO+WHg\n/BL60I7OA/pTSp9PKX1xYnRhIXsSePOk+6/IOW+a+Pl+4LLWd6kt/MBxAf5VSumhlNJtKaX+kvpV\npk/S+AcQ4CXAGPDjfl6Af35suoBRGp+ZKxfyZybn/BngPRN3TwP24mdGAWWEp6V8f2oKYGxiCHWh\nOwx8KOf8euAG4M8W8nHJOd9D4x/BYyZfrb4GLGttj9rDCY7Lw8D7Jv5v+Wngg2X0q0w558M550Mp\npSpwJ/AB/LwAJzw2/xH4GvCrC/kzA5BzHk8pfQz4PeAT+JlRQBn/OB8AqpP7kHMeL6Ef7eZx4M8A\ncs5PALuBNaX2qL1M/oxUgX1ldaTNfDrn/HcTP98DvLzMzpQlpXQq8CBwR875/+Dn5XtOcGz8zEzI\nOb8LOBu4DThp0lML+jOj6ZURnjYDVwCklC4AHi2hD+3oOuDDACmltTR+eXeW2qP28o2U0sUTP18O\nbJpq4wXk8ymlY1PfrwEeKbMzZUgpDQKfB/59zvmOiYf/zs/Lix4bPzMpvSOl9P6Ju0eA7wJfTyld\nMvHYgv3MaGZa/m07Gv+n89qU0uaJ++8uoQ/t6I+BP0kpbaLxf83XOSL3z/wq8NGUUg/wbeBTJfen\nXdwA/M+U0giwi++fx7GQ3AQsB25OKf0noA78Eo3jstA/Lyc6NjcC/2OBf2bupvH39iEa/w7+IvAY\ncJufGc1EpV6vl90HSZKkjrFgT0iWJEkqwvAkSZIUYHiSJEkKMDxJkiQFGJ4kSZICDE+SJEkBhidJ\nkqQAw5MkSVLA/weYD9T4bMXZJAAAAABJRU5ErkJggg==\n", | |
"text/plain": [ | |
"<matplotlib.figure.Figure at 0x112fee400>" | |
] | |
}, | |
"metadata": {}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"from gen import simulate_binary\n", | |
"from gen import corr_heat\n", | |
"%matplotlib inline \n", | |
"\n", | |
"data = simulate_binary(100, [5, 5, 5, 5, 5, 5, 5], .9)\n", | |
"corr_heat(data)" | |
] | |
} | |
], | |
"metadata": { | |
"kernelspec": { | |
"display_name": "Python 3", | |
"language": "python", | |
"name": "python3" | |
}, | |
"language_info": { | |
"codemirror_mode": { | |
"name": "ipython", | |
"version": 3 | |
}, | |
"file_extension": ".py", | |
"mimetype": "text/x-python", | |
"name": "python", | |
"nbconvert_exporter": "python", | |
"pygments_lexer": "ipython3", | |
"version": "3.5.2" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 0 | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
class CorrStruct(): | |
"""This class generates valid binary correlation matrices | |
""" | |
def __init__(self, K=10, n=100, mindistance=.9): | |
""" | |
Parameters | |
---------- | |
K : int | |
Dimension of simulated variable. | |
mindistance : float | |
0 < mindistance < 1 must hold. A larger value means | |
larger marginal probabilities for the K-dimensional | |
binary vector. | |
n : int | |
Number of binary vectors to simulate. | |
""" | |
self.K = K | |
self.cube = {} | |
self.marginal = {} | |
for i in range(K): | |
self.cube[i] = self.make_cube(mindistance) | |
self.marginal[i] = self.vol_cube(self.cube[i]) | |
self.sim = correlated_binary_rvs(*self.get_overlaps(), n) | |
def make_cube(self, *args): | |
cube = np.zeros((self.K, 2)) | |
for i in range(self.K): | |
cube[i, :] = self.get_pins(*args) | |
return cube | |
def get_pins(self, mindistance=.9): | |
out = np.zeros(2) | |
while out[1] - out[0] < mindistance: | |
out = np.sort(np.random.random(2)) | |
return out | |
def vol_cube(self, cube): | |
"""expects a K,2 shaped np.array | |
""" | |
return np.product(cube[:, 1] - cube[:, 0]) | |
def get_overlap(self, i, j): | |
cubea = self.cube[i] | |
cubeb = self.cube[j] | |
# the logic in the next three lines is the following: | |
# throw two red pins and two blue pins randomly on | |
# (0,1). This defines a red region and a blue region. | |
# The following lines obtain the volume of that intersection | |
# and extend it to obtaining the volume of the intersection | |
# in the hypercube space. | |
lower = np.max(np.vstack((cubea[:, 0], cubeb[:, 0])), axis=0)[:, None] | |
upper = np.min(np.vstack((cubea[:, 1], cubeb[:, 1])), axis=0)[:, None] | |
return self.vol_cube(np.hstack((lower, upper))) | |
def get_overlaps(self): | |
pAB = np.zeros((self.K, self.K)) | |
corr = np.zeros((self.K, self.K)) | |
pA = self.marginal | |
for i, j in combinations(range(self.K), 2): | |
pAB[i, j] = self.get_overlap(i, j) | |
corr[i, j] = (pAB[i, j] - pA[i]*pA[j]) / np.sqrt(pA[i]*(1-pA[i]) | |
* pA[j]*(1-pA[j])) | |
corr = corr + corr.T | |
np.fill_diagonal(corr, 1) | |
marginal = np.array([self.marginal[x] for x in self.marginal]) | |
return marginal, corr | |
def simulate_binary(n, dims, mindistances): | |
"""Horizontally stacks the structure CorrStruct and returns the output | |
as simulated random variables. | |
Examples | |
-------- | |
> simulate_binary(100, [5,5,5,5], [.9, .8, .7, .6]) | |
Parameters | |
---------- | |
n : int | |
sample size | |
dims : list of int | |
dimension of each correlation block | |
mindistances : list of floats or float | |
Larger means higher marginal probabilities. each float in (0,1). | |
If list, must match length of dims. If not list, then the same | |
marginal probability strength is used for all blocks defined by | |
dims. | |
Returns | |
------- | |
np.array | |
n rows by np.sum(dims) columns. Simulated data. | |
""" | |
if not (len(dims) == len(mindistances) or type(mindistances) is float): | |
raise Exception(('Dimension Vector should have number of ' | |
'elements equal to the number of mindistance ' | |
'parameters, or mindistance should have length 1')) | |
n_blocks = len(dims) | |
if type(mindistances) is float: | |
mindistances = np.ones(n_blocks) * mindistances | |
corrs = [] # instantiate correlation blocks | |
for i in range(n_blocks): | |
out = CorrStruct(dims[i], n, mindistances[i]) | |
corrs.append(out.sim) | |
return np.hstack(tuple(corrs)) | |
def corr_heat(df): | |
""" | |
References | |
---------- | |
https://stanford.edu/~mwaskom/software/seaborn/examples/many_pairwise_correlations.html | |
""" | |
corr = np.cov(data.T) | |
mask = np.zeros_like(corr, dtype=np.bool) | |
mask[np.triu_indices_from(mask)] = True | |
# Set up the matplotlib figure | |
f, ax = plt.subplots(figsize=(11, 9)) | |
# Generate a custom diverging colormap | |
cmap = sns.diverging_palette(220, 10, as_cmap=True) | |
# Draw the heatmap with the mask and correct aspect ratio | |
return sns.heatmap(corr, mask=mask, cmap=cmap, vmax=.3, | |
square=True, xticklabels=5, yticklabels=5, | |
linewidths=.5, cbar_kws={"shrink": .5}, ax=ax) | |
data = simulate_binary(100, [5, 5, 5, 5, 5, 5, 5], .9) | |
corr_heat(data) | |
plt.show() | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
https://gist.github.com/jseabold/b91555532584a918e51b88a818003e2e