|
332 | 332 | }
|
333 | 333 | ],
|
334 | 334 | "source": [
|
335 |
| - "value=np.asarray(np.random.randn(N,), dtype=theano.config.floatX)\n", |
| 335 | + "value = np.asarray(np.random.randn(N,), dtype=theano.config.floatX)\n", |
336 | 336 | "\n",
|
337 | 337 | "maxwz = max([sum(w) for w in weights])\n",
|
338 | 338 | "N = len(weights)\n",
|
|
351 | 351 | "a = tt.matrix(\"a\", dtype='int32')\n",
|
352 | 352 | "a.tag.test_value = amat\n",
|
353 | 353 | "\n",
|
| 354 | + "\n", |
354 | 355 | "def get_mu(w, a):\n",
|
355 | 356 | " a1 = tt.cast(a, 'int32')\n",
|
356 | 357 | " return tt.sum(w*x[a1])/tt.sum(w)\n",
|
|
360 | 361 | "\n",
|
361 | 362 | "print(compute_elementwise(value, wmat, amat))\n",
|
362 | 363 | "\n",
|
| 364 | + "\n", |
363 | 365 | "def mu_phi(value):\n",
|
364 |
| - " N=len(weights)\n", |
365 |
| - " # Calculate mu based on average of neighbours \n", |
| 366 | + " N = len(weights)\n", |
| 367 | + " # Calculate mu based on average of neighbours\n", |
366 | 368 | " mu = np.array([np.sum(weights[i]*value[adj[i]])/Wplus[i] for i in range(N)])\n",
|
367 | 369 | " return mu\n",
|
368 | 370 | "\n",
|
|
402 | 404 | "class CAR(distribution.Continuous):\n",
|
403 | 405 | " \"\"\"\n",
|
404 | 406 | " Conditional Autoregressive (CAR) distribution\n",
|
405 |
| - " \n", |
| 407 | + "\n", |
406 | 408 | " Parameters\n",
|
407 | 409 | " ----------\n",
|
408 | 410 | " a : list of adjacency information\n",
|
|
415 | 417 | " self.w = w = tt.as_tensor_variable(w)\n",
|
416 | 418 | " self.tau = tau*tt.sum(w, axis=1)\n",
|
417 | 419 | " self.mode = 0.\n",
|
418 |
| - " \n", |
| 420 | + "\n", |
419 | 421 | " def get_mu(self, x):\n",
|
420 |
| - " \n", |
| 422 | + "\n", |
421 | 423 | " def weigth_mu(w, a):\n",
|
422 | 424 | " a1 = tt.cast(a, 'int32')\n",
|
423 | 425 | " return tt.sum(w*x[a1])/tt.sum(w)\n",
|
424 | 426 | "\n",
|
425 |
| - " mu_w, _ = scan(fn=weigth_mu, \n", |
| 427 | + " mu_w, _ = scan(fn=weigth_mu,\n", |
426 | 428 | " sequences=[self.w, self.a])\n",
|
427 |
| - " \n", |
| 429 | + "\n", |
428 | 430 | " return mu_w\n",
|
429 |
| - " \n", |
| 431 | + "\n", |
430 | 432 | " def logp(self, x):\n",
|
431 | 433 | " mu_w = self.get_mu(x)\n",
|
432 | 434 | " tau = self.tau\n",
|
|
458 | 460 | }
|
459 | 461 | ],
|
460 | 462 | "source": [
|
461 |
| - "with pm.Model() as model:\n", |
| 463 | + "with pm.Model() as model1:\n", |
462 | 464 | " # Vague prior on intercept\n",
|
463 | 465 | " beta0 = pm.Normal('beta0', mu=0.0, tau=1.0e-5)\n",
|
464 | 466 | " # Vague prior on covariate effect\n",
|
|
488 | 490 | " sd_c = pm.Deterministic('sd_c', tt.std(phi))\n",
|
489 | 491 | " # Proportion sptial variance\n",
|
490 | 492 | " alpha = pm.Deterministic('alpha', sd_c/(sd_h+sd_c))\n",
|
491 |
| - " \n", |
| 493 | + "\n", |
492 | 494 | " trace1 = pm.sample(3e3, njobs=2, tune=1000)"
|
493 | 495 | ]
|
494 | 496 | },
|
|
602 | 604 | "for i, a in enumerate(adj):\n",
|
603 | 605 | " amat2[i, a] = 1\n",
|
604 | 606 | " wmat2[i, a] = weights[i]\n",
|
605 |
| - " \n", |
606 |
| - "value=np.asarray(np.random.randn(N,), dtype=theano.config.floatX)\n", |
| 607 | + "\n", |
| 608 | + "value = np.asarray(np.random.randn(N,), dtype=theano.config.floatX)\n", |
607 | 609 | "\n",
|
608 | 610 | "print(np.sum(value*amat2, axis=1)/np.sum(wmat2, axis=1))\n",
|
609 | 611 | "\n",
|
| 612 | + "\n", |
610 | 613 | "def mu_phi(value):\n",
|
611 |
| - " N=len(weights)\n", |
612 |
| - " # Calculate mu based on average of neighbours \n", |
| 614 | + " N = len(weights)\n", |
| 615 | + " # Calculate mu based on average of neighbours\n", |
613 | 616 | " mu = np.array([np.sum(weights[i]*value[adj[i]])/Wplus[i] for i in range(N)])\n",
|
614 | 617 | " return mu\n",
|
615 | 618 | "\n",
|
|
634 | 637 | "class CAR2(distribution.Continuous):\n",
|
635 | 638 | " \"\"\"\n",
|
636 | 639 | " Conditional Autoregressive (CAR) distribution\n",
|
637 |
| - " \n", |
| 640 | + "\n", |
638 | 641 | " Parameters\n",
|
639 | 642 | " ----------\n",
|
640 | 643 | " a : adjacency matrix\n",
|
|
648 | 651 | " self.w = w = tt.as_tensor_variable(w)\n",
|
649 | 652 | " self.tau = tau*tt.sum(w, axis=1)\n",
|
650 | 653 | " self.mode = 0.\n",
|
651 |
| - " \n", |
| 654 | + "\n", |
652 | 655 | " def logp(self, x):\n",
|
653 | 656 | " tau = self.tau\n",
|
654 | 657 | " w = self.w\n",
|
655 | 658 | " a = self.a\n",
|
656 |
| - " \n", |
| 659 | + "\n", |
657 | 660 | " mu_w = tt.sum(x*a, axis=1)/tt.sum(w, axis=1)\n",
|
658 | 661 | " return tt.sum(continuous.Normal.dist(mu=mu_w, tau=tau).logp(x))"
|
659 | 662 | ]
|
|
676 | 679 | }
|
677 | 680 | ],
|
678 | 681 | "source": [
|
679 |
| - "with pm.Model() as model:\n", |
| 682 | + "with pm.Model() as model2:\n", |
680 | 683 | " # Vague prior on intercept\n",
|
681 | 684 | " beta0 = pm.Normal('beta0', mu=0.0, tau=1.0e-5)\n",
|
682 | 685 | " # Vague prior on covariate effect\n",
|
|
706 | 709 | " sd_c = pm.Deterministic('sd_c', tt.std(phi))\n",
|
707 | 710 | " # Proportion sptial variance\n",
|
708 | 711 | " alpha = pm.Deterministic('alpha', sd_c/(sd_h+sd_c))\n",
|
709 |
| - " \n", |
| 712 | + "\n", |
710 | 713 | " trace2 = pm.sample(3e3, njobs=2, tune=1000)"
|
711 | 714 | ]
|
712 | 715 | },
|
|
807 | 810 | },
|
808 | 811 | "outputs": [],
|
809 | 812 | "source": [
|
810 |
| - "X = np.hstack((np.ones((N,1)), stats.zscore(aff, ddof=1)[:,None]))\n", |
| 813 | + "X = np.hstack((np.ones((N, 1)), stats.zscore(aff, ddof=1)[:, None]))\n", |
811 | 814 | "W = wmat2\n",
|
812 | 815 | "D = np.diag(W.sum(axis=1))\n",
|
813 |
| - "log_offset = logE[:,None]" |
| 816 | + "log_offset = logE[:, None]" |
814 | 817 | ]
|
815 | 818 | },
|
816 | 819 | {
|
|
845 | 848 | }
|
846 | 849 | ],
|
847 | 850 | "source": [
|
848 |
| - "with pm.Model() as model:\n", |
| 851 | + "with pm.Model() as model3:\n", |
849 | 852 | " # Vague prior on intercept and effect\n",
|
850 | 853 | " beta = pm.Normal('beta', mu=0.0, tau=1.0, shape=(2, 1))\n",
|
851 | 854 | "\n",
|
|
859 | 862 | "\n",
|
860 | 863 | " # Likelihood\n",
|
861 | 864 | " Yi = pm.Poisson('Yi', mu=mu.ravel(), observed=O)\n",
|
862 |
| - " \n", |
| 865 | + "\n", |
863 | 866 | " trace3 = pm.sample(3e3, njobs=2, tune=1000)"
|
864 | 867 | ]
|
865 | 868 | },
|
|
1025 | 1028 | "outputs": [],
|
1026 | 1029 | "source": [
|
1027 | 1030 | "import scipy\n",
|
| 1031 | + "\n", |
| 1032 | + "\n", |
1028 | 1033 | "class Sparse_CAR(distribution.Continuous):\n",
|
1029 | 1034 | " \"\"\"\n",
|
1030 | 1035 | " Sparse Conditional Autoregressive (CAR) distribution\n",
|
1031 |
| - " \n", |
| 1036 | + "\n", |
1032 | 1037 | " Parameters\n",
|
1033 | 1038 | " ----------\n",
|
1034 | 1039 | " alpha : spatial smoothing term\n",
|
|
1044 | 1049 | " self.n = n\n",
|
1045 | 1050 | " self.median = self.mode = self.mean = 0\n",
|
1046 | 1051 | " super(Sparse_CAR, self).__init__(*args, **kwargs)\n",
|
1047 |
| - " \n", |
| 1052 | + "\n", |
1048 | 1053 | " # eigenvalues of D^−1/2 * W * D^−1/2\n",
|
1049 | 1054 | " Dinv_sqrt = np.diag(1 / np.sqrt(D))\n",
|
1050 | 1055 | " DWD = np.matmul(np.matmul(Dinv_sqrt, W), Dinv_sqrt)\n",
|
1051 | 1056 | " self.lam = scipy.linalg.eigvalsh(DWD)\n",
|
1052 |
| - " \n", |
| 1057 | + "\n", |
1053 | 1058 | " # sparse representation of W\n",
|
1054 | 1059 | " w_sparse = scipy.sparse.csr_matrix(W)\n",
|
1055 | 1060 | " self.W = theano.sparse.as_sparse_variable(w_sparse)\n",
|
1056 | 1061 | " self.D = tt.as_tensor_variable(D)\n",
|
1057 |
| - " \n", |
| 1062 | + "\n", |
1058 | 1063 | " # Presicion Matrix (inverse of Covariance matrix)\n",
|
1059 | 1064 | " # d_sparse = scipy.sparse.csr_matrix(np.diag(D))\n",
|
1060 |
| - " # self.D = theano.sparse.as_sparse_variable(d_sparse)\n", |
| 1065 | + " # self.D = theano.sparse.as_sparse_variable(d_sparse)\n", |
1061 | 1066 | " # self.Phi = self.tau * (self.D - self.alpha*self.W)\n",
|
1062 |
| - " \n", |
| 1067 | + "\n", |
1063 | 1068 | " def logp(self, x):\n",
|
1064 | 1069 | " logtau = self.n * tt.log(tau)\n",
|
1065 | 1070 | " logdet = tt.log(1 - self.alpha * self.lam).sum()\n",
|
1066 |
| - " \n", |
| 1071 | + "\n", |
1067 | 1072 | " # tau * ((phi .* D_sparse)' * phi - alpha * (phit_W * phi))\n",
|
1068 | 1073 | " Wx = theano.sparse.dot(self.W, x)\n",
|
1069 | 1074 | " tau_dot_x = self.D * x.T - self.alpha * Wx.ravel()\n",
|
1070 | 1075 | " logquad = self.tau * tt.dot(x.ravel(), tau_dot_x.ravel())\n",
|
1071 |
| - " \n", |
| 1076 | + "\n", |
1072 | 1077 | " # logquad = tt.dot(x.T, theano.sparse.dot(self.Phi, x)).sum()\n",
|
1073 | 1078 | " return 0.5*(logtau + logdet - logquad)"
|
1074 | 1079 | ]
|
|
1093 | 1098 | }
|
1094 | 1099 | ],
|
1095 | 1100 | "source": [
|
1096 |
| - "with pm.Model() as model:\n", |
| 1101 | + "with pm.Model() as model4:\n", |
1097 | 1102 | " # Vague prior on intercept and effect\n",
|
1098 | 1103 | " beta = pm.Normal('beta', mu=0.0, tau=1.0, shape=(2, 1))\n",
|
1099 | 1104 | "\n",
|
1100 | 1105 | " # Priors for spatial random effects\n",
|
1101 | 1106 | " tau = pm.Gamma('tau', alpha=2., beta=2.)\n",
|
1102 | 1107 | " alpha = pm.Uniform('alpha', lower=0, upper=1)\n",
|
1103 |
| - " phi = Sparse_CAR('phi', alpha, W, tau, shape=(N,1))\n", |
1104 |
| - " \n", |
| 1108 | + " phi = Sparse_CAR('phi', alpha, W, tau, shape=(N, 1))\n", |
| 1109 | + "\n", |
1105 | 1110 | " # Mean model\n",
|
1106 | 1111 | " mu = pm.Deterministic('mu', tt.exp(tt.dot(X, beta) + phi + log_offset))\n",
|
1107 | 1112 | "\n",
|
1108 | 1113 | " # Likelihood\n",
|
1109 | 1114 | " Yi = pm.Poisson('Yi', mu=mu.ravel(), observed=O)\n",
|
1110 |
| - " \n", |
| 1115 | + "\n", |
1111 | 1116 | " trace4 = pm.sample(3e3, njobs=2, tune=1000)"
|
1112 | 1117 | ]
|
1113 | 1118 | },
|
|
0 commit comments