|
19 | 19 | "from theano import tensor as tt\n",
|
20 | 20 | "\n",
|
21 | 21 | "%config InlineBackend.figure_format = 'retina'\n",
|
22 |
| - "warnings.simplefilter(action='ignore', category=(FutureWarning, UserWarning))\n", |
| 22 | + "warnings.simplefilter(action=\"ignore\", category=(FutureWarning, UserWarning))\n", |
23 | 23 | "RANDOM_SEED = 8927\n",
|
24 | 24 | "np.random.seed(RANDOM_SEED)"
|
25 | 25 | ]
|
|
111 | 111 | "\n",
|
112 | 112 | " cov = pm.Deterministic(\n",
|
113 | 113 | " \"cov\",\n",
|
114 |
| - " tt.stacklists(\n", |
115 |
| - " [[lambda1 ** -1, r * sigma1 * sigma2], [r * sigma1 * sigma2, lambda2 ** -1]]\n", |
116 |
| - " ),\n", |
| 114 | + " tt.stacklists([[lambda1 ** -1, r * sigma1 * sigma2], [r * sigma1 * sigma2, lambda2 ** -1]]),\n", |
117 | 115 | " )\n",
|
118 | 116 | "\n",
|
119 | 117 | " tau1 = pm.Deterministic(\"tau1\", tt.nlinalg.matrix_inverse(cov))\n",
|
|
166 | 164 | "my_pdf = stats.kde.gaussian_kde(r)\n",
|
167 | 165 | "x = np.linspace(-1, 1, 200)\n",
|
168 | 166 | "axes[1].plot(\n",
|
169 |
| - " x, my_pdf(x), \"r\", lw=2.5, alpha=0.6,\n", |
| 167 | + " x,\n", |
| 168 | + " my_pdf(x),\n", |
| 169 | + " \"r\",\n", |
| 170 | + " lw=2.5,\n", |
| 171 | + " alpha=0.6,\n", |
170 | 172 | ") # distribution function\n",
|
171 | 173 | "axes[1].plot(x, np.ones(x.size) * 0.5, \"k--\", lw=2.5, alpha=0.6, label=\"Prior\")\n",
|
172 | 174 | "\n",
|
|
175 | 177 | "BF01 = posterior / prior\n",
|
176 | 178 | "print(\"the Bayes Factor is %.3f\" % (1 / BF01))\n",
|
177 | 179 | "\n",
|
178 |
| - "axes[1].plot(\n", |
179 |
| - " [0, 0], [posterior, prior], \"k-\", [0, 0], [posterior, prior], \"ko\", lw=1.5, alpha=1\n", |
180 |
| - ")\n", |
| 180 | + "axes[1].plot([0, 0], [posterior, prior], \"k-\", [0, 0], [posterior, prior], \"ko\", lw=1.5, alpha=1)\n", |
181 | 181 | "# ax1.hist(r, bins=100, normed=1,alpha=.3)\n",
|
182 | 182 | "axes[1].set_xlabel(\"Correlation r\")\n",
|
183 | 183 | "axes[1].set_ylabel(\"Density\")\n",
|
|
193 | 193 | "\n",
|
194 | 194 | "# Compare to exact solution Jeffreys (numerical integration):\n",
|
195 | 195 | "BF_Jeffex = sp.integrate.quad(\n",
|
196 |
| - " lambda rho: ((1 - rho ** 2) ** ((n - 1) / 2))\n", |
197 |
| - " / ((1 - rho * freqr[0, 1]) ** ((n - 1) - 0.5)),\n", |
| 196 | + " lambda rho: ((1 - rho ** 2) ** ((n - 1) / 2)) / ((1 - rho * freqr[0, 1]) ** ((n - 1) - 0.5)),\n", |
198 | 197 | " -1,\n",
|
199 | 198 | " 1,\n",
|
200 | 199 | ")\n",
|
|
239 | 238 | "source": [
|
240 | 239 | "# Data\n",
|
241 | 240 | "# Proportion correct on erotic pictures, block 1 and block 2:\n",
|
| 241 | + "# fmt: off\n", |
242 | 242 | "prc1er=np.array([0.6000000, 0.5333333, 0.6000000, 0.6000000, 0.4666667, \n",
|
243 | 243 | " 0.6666667, 0.6666667, 0.4000000, 0.6000000, 0.6000000,\n",
|
244 | 244 | " 0.4666667, 0.6666667, 0.4666667, 0.6000000, 0.3333333,\n",
|
|
280 | 280 | " 0.5333333, 0.6000000, 0.6666667, 0.4000000, 0.4000000,\n",
|
281 | 281 | " 0.5333333, 0.8000000, 0.6000000, 0.4000000, 0.2000000,\n",
|
282 | 282 | " 0.6000000, 0.6666667, 0.4666667, 0.4666667, 0.4666667])\n",
|
| 283 | + "# fmt: on\n", |
283 | 284 | "Nt = 60\n",
|
284 | 285 | "xobs = np.vstack([prc1er, prc2er]).T * Nt\n",
|
285 | 286 | "n, n2 = np.shape(xobs) # number of participants\n",
|
|
345 | 346 | "\n",
|
346 | 347 | " cov = pm.Deterministic(\n",
|
347 | 348 | " \"cov\",\n",
|
348 |
| - " tt.stacklists(\n", |
349 |
| - " [[lambda1 ** -1, r * sigma1 * sigma2], [r * sigma1 * sigma2, lambda2 ** -1]]\n", |
350 |
| - " ),\n", |
| 349 | + " tt.stacklists([[lambda1 ** -1, r * sigma1 * sigma2], [r * sigma1 * sigma2, lambda2 ** -1]]),\n", |
351 | 350 | " )\n",
|
352 | 351 | "\n",
|
353 | 352 | " tau1 = pm.Deterministic(\"tau1\", tt.nlinalg.matrix_inverse(cov))\n",
|
|
398 | 397 | "my_pdf = stats.kde.gaussian_kde(r)\n",
|
399 | 398 | "x1 = np.linspace(0, 1, 200)\n",
|
400 | 399 | "plt.plot(\n",
|
401 |
| - " x1, my_pdf(x1), \"r\", lw=2.5, alpha=0.6,\n", |
| 400 | + " x1,\n", |
| 401 | + " my_pdf(x1),\n", |
| 402 | + " \"r\",\n", |
| 403 | + " lw=2.5,\n", |
| 404 | + " alpha=0.6,\n", |
402 | 405 | ") # distribution function\n",
|
403 | 406 | "plt.plot(x1, np.ones(x1.size), \"k--\", lw=2.5, alpha=0.6, label=\"Prior\")\n",
|
404 | 407 | "\n",
|
|
407 | 410 | "BF01 = posterior / prior\n",
|
408 | 411 | "print(\"the Bayes Factor is %.5f\" % (1 / BF01))\n",
|
409 | 412 | "\n",
|
410 |
| - "plt.plot(\n", |
411 |
| - " [0, 0], [posterior, prior], \"k-\", [0, 0], [posterior, prior], \"ko\", lw=1.5, alpha=1\n", |
412 |
| - ")\n", |
| 413 | + "plt.plot([0, 0], [posterior, prior], \"k-\", [0, 0], [posterior, prior], \"ko\", lw=1.5, alpha=1)\n", |
413 | 414 | "plt.xlim([0, 1])\n",
|
414 | 415 | "plt.xlabel(\"Correlation r\")\n",
|
415 | 416 | "plt.ylabel(\"Density\")\n",
|
416 | 417 | "\n",
|
417 | 418 | "# Compare to exact solution Jeffreys (numerical integration):\n",
|
418 | 419 | "freqr = sp.corrcoef(xobs[:, 0] / Nt, xobs[:, 1] / Nt)\n",
|
419 | 420 | "BF_Jeffex = sp.integrate.quad(\n",
|
420 |
| - " lambda rho: ((1 - rho ** 2) ** ((n - 1) / 2))\n", |
421 |
| - " / ((1 - rho * freqr[0, 1]) ** ((n - 1) - 0.5)),\n", |
| 421 | + " lambda rho: ((1 - rho ** 2) ** ((n - 1) / 2)) / ((1 - rho * freqr[0, 1]) ** ((n - 1) - 0.5)),\n", |
422 | 422 | " 0,\n",
|
423 | 423 | " 1,\n",
|
424 | 424 | ")\n",
|
|
463 | 463 | }
|
464 | 464 | ],
|
465 | 465 | "source": [
|
466 |
| - "k2 = np.array([36, 32, 36, 36, 28, 40, 40, 24, 36, 36, 28, 40, 28, \n", |
467 |
| - " 36, 20, 24, 24, 16, 20, 32, 40, 32, 36, 24, 28, 44,\n", |
468 |
| - " 40, 36, 40, 32, 32, 40, 28, 20, 24, 32, 24, 24, 20, \n", |
469 |
| - " 28, 24, 28, 28, 32, 20, 44, 16, 36, 32, 28, 24, 32,\n", |
470 |
| - " 40, 28, 32, 32, 28, 24, 28, 40, 28, 20, 20, 20, 24,\n", |
471 |
| - " 24, 36, 28, 20, 20, 40, 32, 20, 36, 28, 28, 24, 20,\n", |
472 |
| - " 28, 32, 48, 24, 32, 32, 40, 40, 40, 36, 36, 32, 20,\n", |
473 |
| - " 28, 40, 32, 20, 20, 16, 16, 28, 40])\n", |
474 |
| - " \n", |
475 |
| - "x2 = np.array([50, 80, 79, 56, 50, 80, 53, 84, 74, 67, 50, 45, 62, \n", |
476 |
| - " 65, 71, 71, 68, 63, 67, 58, 72, 73, 63, 54, 63, 70, \n", |
477 |
| - " 81, 71, 66, 74, 70, 84, 66, 73, 78, 64, 54, 74, 62, \n", |
478 |
| - " 71, 70, 79, 66, 64, 62, 63, 60, 56, 72, 72, 79, 67, \n", |
479 |
| - " 46, 67, 77, 55, 63, 44, 84, 65, 41, 62, 64, 51, 46,\n", |
480 |
| - " 53, 26, 67, 73, 39, 62, 59, 75, 65, 60, 69, 63, 69, \n", |
481 |
| - " 55, 63, 86, 70, 67, 54, 80, 71, 71, 55, 57, 41, 56, \n", |
482 |
| - " 78, 58, 76, 54, 50, 61, 60, 32, 67])\n", |
| 466 | + "k2 = np.array(\n", |
| 467 | + " [\n", |
| 468 | + " 36,\n", |
| 469 | + " 32,\n", |
| 470 | + " 36,\n", |
| 471 | + " 36,\n", |
| 472 | + " 28,\n", |
| 473 | + " 40,\n", |
| 474 | + " 40,\n", |
| 475 | + " 24,\n", |
| 476 | + " 36,\n", |
| 477 | + " 36,\n", |
| 478 | + " 28,\n", |
| 479 | + " 40,\n", |
| 480 | + " 28,\n", |
| 481 | + " 36,\n", |
| 482 | + " 20,\n", |
| 483 | + " 24,\n", |
| 484 | + " 24,\n", |
| 485 | + " 16,\n", |
| 486 | + " 20,\n", |
| 487 | + " 32,\n", |
| 488 | + " 40,\n", |
| 489 | + " 32,\n", |
| 490 | + " 36,\n", |
| 491 | + " 24,\n", |
| 492 | + " 28,\n", |
| 493 | + " 44,\n", |
| 494 | + " 40,\n", |
| 495 | + " 36,\n", |
| 496 | + " 40,\n", |
| 497 | + " 32,\n", |
| 498 | + " 32,\n", |
| 499 | + " 40,\n", |
| 500 | + " 28,\n", |
| 501 | + " 20,\n", |
| 502 | + " 24,\n", |
| 503 | + " 32,\n", |
| 504 | + " 24,\n", |
| 505 | + " 24,\n", |
| 506 | + " 20,\n", |
| 507 | + " 28,\n", |
| 508 | + " 24,\n", |
| 509 | + " 28,\n", |
| 510 | + " 28,\n", |
| 511 | + " 32,\n", |
| 512 | + " 20,\n", |
| 513 | + " 44,\n", |
| 514 | + " 16,\n", |
| 515 | + " 36,\n", |
| 516 | + " 32,\n", |
| 517 | + " 28,\n", |
| 518 | + " 24,\n", |
| 519 | + " 32,\n", |
| 520 | + " 40,\n", |
| 521 | + " 28,\n", |
| 522 | + " 32,\n", |
| 523 | + " 32,\n", |
| 524 | + " 28,\n", |
| 525 | + " 24,\n", |
| 526 | + " 28,\n", |
| 527 | + " 40,\n", |
| 528 | + " 28,\n", |
| 529 | + " 20,\n", |
| 530 | + " 20,\n", |
| 531 | + " 20,\n", |
| 532 | + " 24,\n", |
| 533 | + " 24,\n", |
| 534 | + " 36,\n", |
| 535 | + " 28,\n", |
| 536 | + " 20,\n", |
| 537 | + " 20,\n", |
| 538 | + " 40,\n", |
| 539 | + " 32,\n", |
| 540 | + " 20,\n", |
| 541 | + " 36,\n", |
| 542 | + " 28,\n", |
| 543 | + " 28,\n", |
| 544 | + " 24,\n", |
| 545 | + " 20,\n", |
| 546 | + " 28,\n", |
| 547 | + " 32,\n", |
| 548 | + " 48,\n", |
| 549 | + " 24,\n", |
| 550 | + " 32,\n", |
| 551 | + " 32,\n", |
| 552 | + " 40,\n", |
| 553 | + " 40,\n", |
| 554 | + " 40,\n", |
| 555 | + " 36,\n", |
| 556 | + " 36,\n", |
| 557 | + " 32,\n", |
| 558 | + " 20,\n", |
| 559 | + " 28,\n", |
| 560 | + " 40,\n", |
| 561 | + " 32,\n", |
| 562 | + " 20,\n", |
| 563 | + " 20,\n", |
| 564 | + " 16,\n", |
| 565 | + " 16,\n", |
| 566 | + " 28,\n", |
| 567 | + " 40,\n", |
| 568 | + " ]\n", |
| 569 | + ")\n", |
| 570 | + "\n", |
| 571 | + "x2 = np.array(\n", |
| 572 | + " [\n", |
| 573 | + " 50,\n", |
| 574 | + " 80,\n", |
| 575 | + " 79,\n", |
| 576 | + " 56,\n", |
| 577 | + " 50,\n", |
| 578 | + " 80,\n", |
| 579 | + " 53,\n", |
| 580 | + " 84,\n", |
| 581 | + " 74,\n", |
| 582 | + " 67,\n", |
| 583 | + " 50,\n", |
| 584 | + " 45,\n", |
| 585 | + " 62,\n", |
| 586 | + " 65,\n", |
| 587 | + " 71,\n", |
| 588 | + " 71,\n", |
| 589 | + " 68,\n", |
| 590 | + " 63,\n", |
| 591 | + " 67,\n", |
| 592 | + " 58,\n", |
| 593 | + " 72,\n", |
| 594 | + " 73,\n", |
| 595 | + " 63,\n", |
| 596 | + " 54,\n", |
| 597 | + " 63,\n", |
| 598 | + " 70,\n", |
| 599 | + " 81,\n", |
| 600 | + " 71,\n", |
| 601 | + " 66,\n", |
| 602 | + " 74,\n", |
| 603 | + " 70,\n", |
| 604 | + " 84,\n", |
| 605 | + " 66,\n", |
| 606 | + " 73,\n", |
| 607 | + " 78,\n", |
| 608 | + " 64,\n", |
| 609 | + " 54,\n", |
| 610 | + " 74,\n", |
| 611 | + " 62,\n", |
| 612 | + " 71,\n", |
| 613 | + " 70,\n", |
| 614 | + " 79,\n", |
| 615 | + " 66,\n", |
| 616 | + " 64,\n", |
| 617 | + " 62,\n", |
| 618 | + " 63,\n", |
| 619 | + " 60,\n", |
| 620 | + " 56,\n", |
| 621 | + " 72,\n", |
| 622 | + " 72,\n", |
| 623 | + " 79,\n", |
| 624 | + " 67,\n", |
| 625 | + " 46,\n", |
| 626 | + " 67,\n", |
| 627 | + " 77,\n", |
| 628 | + " 55,\n", |
| 629 | + " 63,\n", |
| 630 | + " 44,\n", |
| 631 | + " 84,\n", |
| 632 | + " 65,\n", |
| 633 | + " 41,\n", |
| 634 | + " 62,\n", |
| 635 | + " 64,\n", |
| 636 | + " 51,\n", |
| 637 | + " 46,\n", |
| 638 | + " 53,\n", |
| 639 | + " 26,\n", |
| 640 | + " 67,\n", |
| 641 | + " 73,\n", |
| 642 | + " 39,\n", |
| 643 | + " 62,\n", |
| 644 | + " 59,\n", |
| 645 | + " 75,\n", |
| 646 | + " 65,\n", |
| 647 | + " 60,\n", |
| 648 | + " 69,\n", |
| 649 | + " 63,\n", |
| 650 | + " 69,\n", |
| 651 | + " 55,\n", |
| 652 | + " 63,\n", |
| 653 | + " 86,\n", |
| 654 | + " 70,\n", |
| 655 | + " 67,\n", |
| 656 | + " 54,\n", |
| 657 | + " 80,\n", |
| 658 | + " 71,\n", |
| 659 | + " 71,\n", |
| 660 | + " 55,\n", |
| 661 | + " 57,\n", |
| 662 | + " 41,\n", |
| 663 | + " 56,\n", |
| 664 | + " 78,\n", |
| 665 | + " 58,\n", |
| 666 | + " 76,\n", |
| 667 | + " 54,\n", |
| 668 | + " 50,\n", |
| 669 | + " 61,\n", |
| 670 | + " 60,\n", |
| 671 | + " 32,\n", |
| 672 | + " 67,\n", |
| 673 | + " ]\n", |
| 674 | + ")\n", |
483 | 675 | "\n",
|
484 | 676 | "nsubjs = len(k2)\n",
|
485 | 677 | "ntrials = 60\n",
|
|
541 | 733 | "\n",
|
542 | 734 | " cov = pm.Deterministic(\n",
|
543 | 735 | " \"cov\",\n",
|
544 |
| - " tt.stacklists(\n", |
545 |
| - " [[lambda1 ** -1, r * sigma1 * sigma2], [r * sigma1 * sigma2, lambda2 ** -1]]\n", |
546 |
| - " ),\n", |
| 736 | + " tt.stacklists([[lambda1 ** -1, r * sigma1 * sigma2], [r * sigma1 * sigma2, lambda2 ** -1]]),\n", |
547 | 737 | " )\n",
|
548 | 738 | "\n",
|
549 | 739 | " tau1 = pm.Deterministic(\"tau1\", tt.nlinalg.matrix_inverse(cov))\n",
|
|
595 | 785 | "my_pdf = stats.kde.gaussian_kde(r)\n",
|
596 | 786 | "x1 = np.linspace(-1, 1, 200)\n",
|
597 | 787 | "plt.plot(\n",
|
598 |
| - " x1, my_pdf(x1), \"r\", lw=2.5, alpha=0.6,\n", |
| 788 | + " x1,\n", |
| 789 | + " my_pdf(x1),\n", |
| 790 | + " \"r\",\n", |
| 791 | + " lw=2.5,\n", |
| 792 | + " alpha=0.6,\n", |
599 | 793 | ") # distribution function\n",
|
600 | 794 | "plt.plot(x1, np.ones(x1.size) * 0.5, \"k--\", lw=2.5, alpha=0.6, label=\"Prior\")\n",
|
601 | 795 | "\n",
|
|
604 | 798 | "BF01 = posterior / prior\n",
|
605 | 799 | "print(\"the Bayes Factor is %.5f\" % (1 / BF01))\n",
|
606 | 800 | "\n",
|
607 |
| - "plt.plot(\n", |
608 |
| - " [0, 0], [posterior, prior], \"k-\", [0, 0], [posterior, prior], \"ko\", lw=1.5, alpha=1\n", |
609 |
| - ")\n", |
| 801 | + "plt.plot([0, 0], [posterior, prior], \"k-\", [0, 0], [posterior, prior], \"ko\", lw=1.5, alpha=1)\n", |
610 | 802 | "# ax1.hist(r, bins=100, normed=1,alpha=.3)\n",
|
611 | 803 | "plt.xlim([-1, 1])\n",
|
612 | 804 | "plt.xlabel(\"Correlation r\")\n",
|
|
0 commit comments