6464 "ax = fig.add_subplot(122, projection='3d')\n",
6565 "ax.plot_surface(X, Y, M, cmap=cm.jet, vmax=1, vmin=-.15)\n",
6666 "ax.view_init(azim=390)\n",
67- "plt.title(\"Uniform prior landscape; alternate view\")"
67+ "plt.title(\"Uniform prior landscape; alternate view\")\n "
6868 ],
6969 "language": "python",
7070 "metadata": {},
122122 "ax = fig.add_subplot(122, projection='3d')\n",
123123 "ax.plot_surface(X, Y, M, cmap=cm.jet)\n",
124124 "ax.view_init(azim=390)\n",
125- "plt.title(\"$Exp(3), Exp(10)$ prior landscape; \\nalternate view\")"
125+ "plt.title(\"$Exp(3), Exp(10)$ prior landscape; \\nalternate view\")\n "
126126 ],
127127 "language": "python",
128128 "metadata": {},
158158 "cell_type": "code",
159159 "collapsed": false,
160160 "input": [
161- "### create the observed data\n",
161+ "# create the observed data\n",
162162 "\n",
163- "#sample size of data we observe, trying varying this (keep it less than 100 ;)\n",
163+ "# sample size of data we observe, trying varying this (keep it less than 100 ;)\n",
164164 "N = 1\n",
165165 "\n",
166- "#the true parameters, but of course we do not see these values...\n",
166+ "# the true parameters, but of course we do not see these values...\n",
167167 "lambda_1_true = 1\n",
168168 "lambda_2_true = 3\n",
169169 "\n",
170170 "#...we see the data generated, dependent on the above two values.\n",
171171 "data = np.concatenate([\n",
172172 " stats.poisson.rvs(lambda_1_true, size=(N, 1)),\n",
173173 " stats.poisson.rvs(lambda_2_true, size=(N, 1))\n",
174- " ], axis=1)\n",
174+ "], axis=1)\n",
175175 "print \"observed (2-dimensional,sample size = %d):\" % N, data\n",
176176 "\n",
177- "#plotting details.\n",
177+ "# plotting details.\n",
178178 "x = y = np.linspace(.01, 5, 100)\n",
179179 "likelihood_x = np.array([stats.poisson.pmf(data[:, 0], _x)\n",
180180 " for _x in x]).prod(axis=1)\n",
181181 "likelihood_y = np.array([stats.poisson.pmf(data[:, 1], _y)\n",
182182 " for _y in y]).prod(axis=1)\n",
183- "L = np.dot(likelihood_x[:, None], likelihood_y[None, :])"
183+ "L = np.dot(likelihood_x[:, None], likelihood_y[None, :])\n "
184184 ],
185185 "language": "python",
186186 "metadata": {},
200200 "collapsed": false,
201201 "input": [
202202 "figsize(12.5, 12)\n",
203- "#matplotlib heavy lifting below, beware!\n",
203+ "# matplotlib heavy lifting below, beware!\n",
204204 "subplot(221)\n",
205205 "uni_x = stats.uniform.pdf(x, loc=0, scale=5)\n",
206206 "uni_y = stats.uniform.pdf(x, loc=0, scale=5)\n",
244244 "plt.title(\"Landscape warped by %d data observation;\\n Exponential priors on \\\n",
245245 "$p_1, p_2$.\" % N)\n",
246246 "plt.xlim(0, 5)\n",
247- "plt.ylim(0, 5)"
247+ "plt.ylim(0, 5)\n "
248248 ],
249249 "language": "python",
250250 "metadata": {},
334334 "collapsed": false,
335335 "input": [
336336 "figsize(12.5, 4)\n",
337- "data = np.loadtxt(\"data/mixture_data.csv\", delimiter=\",\")\n",
337+ "data = np.loadtxt(\"data/mixture_data.csv\", delimiter=\",\")\n",
338338 "\n",
339- "hist(data, bins=20, color=\"k\", histtype=\"stepfilled\", alpha=0.8)\n",
339+ "hist(data, bins=20, color=\"k\", histtype=\"stepfilled\", alpha=0.8)\n",
340340 "plt.title(\"Histogram of the dataset\")\n",
341341 "plt.ylim([0, None])\n",
342- "print data[:10], \"...\""
342+ "print data[:10], \"...\"\n "
343343 ],
344344 "language": "python",
345345 "metadata": {},
387387 "\n",
388388 "p = pm.Uniform(\"p\", 0, 1)\n",
389389 "\n",
390- "assignment = pm.Categorical(\"assignment\", [p, 1- p], size=data.shape[0])\n",
390+ "assignment = pm.Categorical(\"assignment\", [p, 1 - p], size=data.shape[0])\n",
391391 "print \"prior assignment, with p = %.2f:\" % p.value\n",
392- "print assignment.value[:10], \"...\""
392+ "print assignment.value[:10], \"...\"\n "
393393 ],
394394 "language": "python",
395395 "metadata": {},
426426 "cell_type": "code",
427427 "collapsed": false,
428428 "input": [
429- "taus = 1.0/ pm.Uniform(\"stds\", 0, 100, size=2) ** 2\n",
429+ "taus = 1.0 / pm.Uniform(\"stds\", 0, 100, size=2) ** 2\n",
430430 "centers = pm.Normal(\"centers\", [120, 190], [0.01, 0.01], size=2)\n",
431431 "\n",
432432 "\"\"\"\n",
447447 "\n",
448448 "print \"Random assignments: \", assignment.value[:4], \"...\"\n",
449449 "print \"Assigned center: \", center_i.value[:4], \"...\"\n",
450- "print \"Assigned precision: \", tau_i.value[:4], \"...\""
450+ "print \"Assigned precision: \", tau_i.value[:4], \"...\"\n "
451451 ],
452452 "language": "python",
453453 "metadata": {},
468468 "cell_type": "code",
469469 "collapsed": false,
470470 "input": [
471- "#and to combine it with the observations:\n",
471+ "# and to combine it with the observations:\n",
472472 "observations = pm.Normal(\"obs\", center_i, tau_i, value=data, observed=True)\n",
473473 "\n",
474- "#below we create a model class\n",
475- "model = pm.Model([p, assignment, taus, centers])"
474+ "# below we create a model class\n",
475+ "model = pm.Model([p, assignment, taus, centers])\n "
476476 ],
477477 "language": "python",
478478 "metadata": {},
495495 "collapsed": false,
496496 "input": [
497497 "mcmc = pm.MCMC(model)\n",
498- "mcmc.sample(50000)"
498+ "mcmc.sample(50000)\n "
499499 ],
500500 "language": "python",
501501 "metadata": {},
527527 "lw = 1\n",
528528 "center_trace = mcmc.trace(\"centers\")[:]\n",
529529 "\n",
530- "#for pretty colors later in the book.\n",
530+ "# for pretty colors later in the book.\n",
531531 "colors = [\"#348ABD\", \"#A60628\"] if center_trace[-1, 0] > center_trace[-1, 1] \\\n",
532- " else [\"#A60628\", \"#348ABD\"]\n",
532+ " else [\"#A60628\", \"#348ABD\"]\n",
533533 "\n",
534534 "plot(center_trace[:, 0], label=\"trace of center 0\", c=colors[0], lw=lw)\n",
535535 "plot(center_trace[:, 1], label=\"trace of center 1\", c=colors[1], lw=lw)\n",
551551 " color=\"#467821\", lw=lw)\n",
552552 "plt.xlabel(\"Steps\")\n",
553553 "plt.ylim(0, 1)\n",
554- "plt.legend()"
554+ "plt.legend()\n "
555555 ],
556556 "language": "python",
557557 "metadata": {},
583583 "cell_type": "code",
584584 "collapsed": false,
585585 "input": [
586- "mcmc.sample(100000)"
586+ "mcmc.sample(100000)\n "
587587 ],
588588 "language": "python",
589589 "metadata": {},
621621 " lw=lw, alpha=0.4, c=colors[0])\n",
622622 "\n",
623623 "x = np.arange(50000, 150000)\n",
624- "plot(x, center_trace[:, 0], label=\"new trace of center 0\", lw=lw, c=\"#348ABD\")\n",
625- "plot(x, center_trace[:, 1], label=\"new trace of center 1\", lw=lw, c=\"#A60628\")\n",
624+ "plot(x, center_trace[:, 0], label=\"new trace of center 0\", lw=lw, c=\"#348ABD\")\n",
625+ "plot(x, center_trace[:, 1], label=\"new trace of center 1\", lw=lw, c=\"#A60628\")\n",
626626 "\n",
627627 "plt.title(\"Traces of unknown center parameters\")\n",
628628 "leg = plt.legend(loc=\"upper right\")\n",
629629 "leg.get_frame().set_alpha(0.8)\n",
630- "plt.xlabel(\"Steps\")"
630+ "plt.xlabel(\"Steps\")\n "
631631 ],
632632 "language": "python",
633633 "metadata": {},
668668 " plt.title(\"Posterior of standard deviation of cluster %d\" % i)\n",
669669 " plt.hist(std_trace[:, i], color=colors[i], bins=30,\n",
670670 " histtype=\"stepfilled\")\n",
671- " #plt.autoscale(tight=True)\n",
671+ " # plt.autoscale(tight=True)\n",
672672 "\n",
673- "plt.tight_layout()"
673+ "plt.tight_layout()\n "
674674 ],
675675 "language": "python",
676676 "metadata": {},
703703 " [\"%.2f\" % s for s in np.sort(data)[::40]])\n",
704704 "plt.ylabel(\"posterior sample\")\n",
705705 "plt.xlabel(\"value of $i$th data point\")\n",
706- "plt.title(\"Posterior labels of data points\")"
706+ "plt.title(\"Posterior labels of data points\")\n "
707707 ],
708708 "language": "python",
709709 "metadata": {},
728728 "input": [
729729 "cmap = mpl.colors.LinearSegmentedColormap.from_list(\"BMH\", colors)\n",
730730 "assign_trace = mcmc.trace(\"assignment\")[:]\n",
731- "scatter(data, 1- assign_trace.mean(axis=0), cmap=cmap,\n",
731+ "scatter(data, 1 - assign_trace.mean(axis=0), cmap=cmap,\n",
732732 " c=assign_trace.mean(axis=0), s=50)\n",
733733 "plt.ylim(-0.05, 1.05)\n",
734734 "plt.xlim(35, 300)\n",
735735 "plt.title(\"Probability of data point belonging to cluster 0\")\n",
736736 "plt.ylabel(\"probability\")\n",
737- "plt.xlabel(\"value of data point\")"
737+ "plt.xlabel(\"value of data point\")\n "
738738 ],
739739 "language": "python",
740740 "metadata": {},
778778 "plt.fill_between(x, y, color=colors[0], alpha=0.3)\n",
779779 "\n",
780780 "plt.legend(loc=\"upper left\")\n",
781- "plt.title(\"Visualizing Clusters using posterior-mean parameters\")"
781+ "plt.title(\"Visualizing Clusters using posterior-mean parameters\")\n "
782782 ],
783783 "language": "python",
784784 "metadata": {},
817817 "\n",
818818 "plt.plot(ex_mcmc.trace(\"x\")[:])\n",
819819 "plt.plot(ex_mcmc.trace(\"y\")[:])\n",
820- "plt.title(\"Displaying (extreme) case of dependence between unknowns\")"
820+ "plt.title(\"Displaying (extreme) case of dependence between unknowns\")\n "
821821 ],
822822 "language": "python",
823823 "metadata": {},
891891 "p_trace = mcmc.trace(\"p\")[:]\n",
892892 "x = 175\n",
893893 "\n",
894- "v = p_trace* norm_pdf(x, loc=center_trace[:, 0], scale=std_trace[:, 0]) > \\\n",
895- " (1- p_trace)* norm_pdf(x, loc=center_trace[:, 1], scale=std_trace[:, 1])\n",
894+ "v = p_trace * norm_pdf(x, loc=center_trace[:, 0], scale=std_trace[:, 0]) > \\\n",
895+ " (1 - p_trace) * norm_pdf(x, loc=center_trace[:, 1], scale=std_trace[:, 1])\n",
896896 "\n",
897- "print \"Probability of belonging to cluster 1:\", v.mean()"
897+ "print \"Probability of belonging to cluster 1:\", v.mean()\n "
898898 ],
899899 "language": "python",
900900 "metadata": {},
988988 "plt.plot(y_t, label=\"$y_t$\", lw=3)\n",
989989 "plt.plot(x_t, label=\"$x_t$\", lw=3)\n",
990990 "plt.xlabel(\"time, $t$\")\n",
991- "plt.legend()"
991+ "plt.legend()\n "
992992 ],
993993 "language": "python",
994994 "metadata": {},
10141014 "collapsed": false,
10151015 "input": [
10161016 "def autocorr(x):\n",
1017- " #from http://tinyurl.com/afz57c4\n",
1017+ " # from http://tinyurl.com/afz57c4\n",
10181018 " result = np.correlate(x, x, mode='full')\n",
10191019 " result = result / np.max(result)\n",
10201020 " return result[result.size / 2:]\n",
10301030 "plt.legend(title=\"Autocorrelation\")\n",
10311031 "plt.ylabel(\"measured correlation \\nbetween $y_t$ and $y_{t-k}$.\")\n",
10321032 "plt.xlabel(\"k (lag)\")\n",
1033- "plt.title(\"Autocorrelation plot of $y_t$ and $x_t$ for differing $k$ lags.\")"
1033+ "plt.title(\"Autocorrelation plot of $y_t$ and $x_t$ for differing $k$ lags.\")\n "
10341034 ],
10351035 "language": "python",
10361036 "metadata": {},
10871087 "plt.ylabel(\"measured correlation \\nbetween $y_t$ and $y_{t-k}$.\")\n",
10881088 "plt.xlabel(\"k (lag)\")\n",
10891089 "plt.title(\"Autocorrelation of $y_t$ (no thinning vs. thinning) \\\n",
1090- "at differing $k$ lags.\")"
1090+ "at differing $k$ lags.\")\n "
10911091 ],
10921092 "language": "python",
10931093 "metadata": {},
11301130 "from pymc.Matplot import plot as mcplot\n",
11311131 "\n",
11321132 "mcmc.sample(25000, 0, 10)\n",
1133- "mcplot(mcmc.trace(\"centers\", 2), common_scale=False)"
1133+ "mcplot(mcmc.trace(\"centers\", 2), common_scale=False)\n "
11341134 ],
11351135 "language": "python",
11361136 "metadata": {},
12491249 "def css_styling():\n",
12501250 " styles = open(\"../styles/custom.css\", \"r\").read()\n",
12511251 " return HTML(styles)\n",
1252- "css_styling()"
1252+ "css_styling()\n "
12531253 ],
12541254 "language": "python",
12551255 "metadata": {},
13381338 "metadata": {}
13391339 }
13401340 ]
1341- }
1341+ }
0 commit comments