Skip to content

Instantly share code, notes, and snippets.

@ricardoV94

ricardoV94/untitled0.ipynb Secret

Created November 13, 2019 10:56
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ricardoV94/916c0138cc9f8902d5a479be44af4198 to your computer and use it in GitHub Desktop.
Save ricardoV94/916c0138cc9f8902d5a479be44af4198 to your computer and use it in GitHub Desktop.
Download ZIP
Untitled0.ipynb
Display the source blob
Display the rendered blob
Raw
untitled0.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "Untitled0.ipynb",
"provenance": [],
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/ricardoV94/916c0138cc9f8902d5a479be44af4198/untitled0.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "code",
"metadata": {
"id": "q24VF7dtYPuZ",
"colab_type": "code",
"colab": {}
},
"source": [
"import pymc3 as pm\n",
"import numpy as np\n",
"import theano.tensor as tt"
],
"execution_count": 0,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "KCFfR5zCYspI",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 35
},
"outputId": "503ab7fd-3a55-4ff7-a085-3224f2bfc053"
},
"source": [
"# simulate career choices among 500 individuals\n",
"N = 500 # number of individuals\n",
"income = np.arange(3) + 1 # expected income of each career\n",
"score = 0.5 * income # scores for each career, based on income\n",
"\n",
"# next line converts scores to probabilities\n",
"def softmax(w):\n",
" e = np.exp(w)\n",
" return e/np.sum(e, axis=0)\n",
"p = softmax(score)\n",
"\n",
"# now simulate choice\n",
"# outcome career holds event type values, not counts\n",
"career = np.random.multinomial(1, p, size=N)\n",
"career = np.where(career==1)[1]\n",
"career[:11]"
],
"execution_count": 2,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([0, 1, 0, 1, 1, 1, 2, 2, 1, 2, 0])"
]
},
"metadata": {
"tags": []
},
"execution_count": 2
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "5OJiM6OQYtXV",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 285
},
"outputId": "063c8b64-c4d5-4764-9f82-ba2e5e7425d7"
},
"source": [
"with pm.Model() as m_10_16:\n",
" b = pm.Normal('b', 0., 5.)\n",
" s2 = b*2\n",
" s3 = b*3\n",
" \n",
" p_ = pm.math.stack([0, s2, s3])\n",
" obs = pm.Categorical('career', p=tt.nnet.softmax(p_), observed=career)\n",
" \n",
" trace_10_16 = pm.sample(1000, tune=2000)\n",
" \n",
"pm.summary(trace_10_16, credible_interval=.89, round_to=2)"
],
"execution_count": 8,
"outputs": [
{
"output_type": "stream",
"text": [
"/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:7: UserWarning: DEPRECATION: If x is a vector, Softmax will not automatically pad x anymore in next releases. If you need it, please do it manually. The vector case is gonna be supported soon and the output will be a vector.\n",
" import sys\n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Sequential sampling (2 chains in 1 job)\n",
"NUTS: [b]\n",
"100%|██████████| 3000/3000 [00:17<00:00, 167.55it/s]\n",
"100%|██████████| 3000/3000 [00:18<00:00, 162.85it/s]\n",
"/usr/local/lib/python3.6/dist-packages/pymc3/stats.py:991: FutureWarning: The join_axes-keyword is deprecated. Use .reindex or .reindex_like on the result to achieve the same functionality.\n",
" axis=1, join_axes=[dforg.index])\n"
],
"name": "stderr"
},
{
"output_type": "execute_result",
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>mean</th>\n",
" <th>sd</th>\n",
" <th>mc_error</th>\n",
" <th>hpd_2.5</th>\n",
" <th>hpd_97.5</th>\n",
" <th>n_eff</th>\n",
" <th>Rhat</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>b</th>\n",
" <td>0.320809</td>\n",
" <td>0.001863</td>\n",
" <td>0.000057</td>\n",
" <td>0.317203</td>\n",
" <td>0.324328</td>\n",
" <td>970.991239</td>\n",
" <td>0.999697</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat\n",
"b 0.320809 0.001863 0.000057 0.317203 0.324328 970.991239 0.999697"
]
},
"metadata": {
"tags": []
},
"execution_count": 8
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "-84KJ_S-Yw34",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 304
},
"outputId": "6cee3a8e-a8aa-4ffd-f37c-64deb3e3056e"
},
"source": [
"with pm.Model() as m_10_16_2:\n",
" b = pm.Normal('b', 0., 5.)\n",
" s1 = b*1\n",
" s2 = b*2\n",
" s3 = b*3\n",
" \n",
" p_ = pm.math.stack([s1, s2, s3])\n",
" obs = pm.Categorical('career', p=tt.nnet.softmax(p_), observed=career)\n",
" \n",
" trace_10_16_2 = pm.sample(1000, tune=2000)\n",
" \n",
"pm.summary(trace_10_16_2, credible_interval=.89, round_to=2)"
],
"execution_count": 9,
"outputs": [
{
"output_type": "stream",
"text": [
"INFO (theano.gof.compilelock): Refreshing lock /root/.theano/compiledir_Linux-4.14.137+-x86_64-with-Ubuntu-18.04-bionic-x86_64-3.6.8-64/lock_dir/lock\n",
"/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:8: UserWarning: DEPRECATION: If x is a vector, Softmax will not automatically pad x anymore in next releases. If you need it, please do it manually. The vector case is gonna be supported soon and the output will be a vector.\n",
" \n",
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Sequential sampling (2 chains in 1 job)\n",
"NUTS: [b]\n",
"100%|██████████| 3000/3000 [00:17<00:00, 174.93it/s]\n",
"100%|██████████| 3000/3000 [00:16<00:00, 184.86it/s]\n",
"/usr/local/lib/python3.6/dist-packages/pymc3/stats.py:991: FutureWarning: The join_axes-keyword is deprecated. Use .reindex or .reindex_like on the result to achieve the same functionality.\n",
" axis=1, join_axes=[dforg.index])\n"
],
"name": "stderr"
},
{
"output_type": "execute_result",
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>mean</th>\n",
" <th>sd</th>\n",
" <th>mc_error</th>\n",
" <th>hpd_2.5</th>\n",
" <th>hpd_97.5</th>\n",
" <th>n_eff</th>\n",
" <th>Rhat</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>b</th>\n",
" <td>0.492978</td>\n",
" <td>0.002616</td>\n",
" <td>0.000085</td>\n",
" <td>0.488229</td>\n",
" <td>0.498476</td>\n",
" <td>933.254766</td>\n",
" <td>1.002209</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat\n",
"b 0.492978 0.002616 0.000085 0.488229 0.498476 933.254766 1.002209"
]
},
"metadata": {
"tags": []
},
"execution_count": 9
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "DSRGXoHjZrEI",
"colab_type": "code",
"colab": {}
},
"source": [
""
],
"execution_count": 0,
"outputs": []
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment