mihaic · Mar 22, 2017 · Sep 2, 2016 · Sep 2, 2016 · Sep 22, 2016 · Sep 28, 2016
diff --git a/brainiak/reprsimil/brsa.py b/brainiak/reprsimil/brsa.py
@@ -418,10 +418,8 @@ def _prepare_DF(self, n_T, scan_onsets=None):
             And we denote A = I - rho1*D + rho1**2 * F"""
         if scan_onsets is None:
             # assume that all data are acquired within the same scan.
-            D = np.diag(np.ones(n_T - 1), -1) + np.diag(np.ones(n_T - 1), 1)
-            F = np.eye(n_T)
-            F[0, 0] = 0
-            F[n_T - 1, n_T - 1] = 0
+            D = self._D_gen(n_T)
+            F = self._F_gen(n_T)
             n_run = 1
             run_TRs = np.array([n_T])
         else:
@@ -439,12 +437,8 @@ def _prepare_DF(self, n_T, scan_onsets=None):
 
             D_ele = map(self._D_gen, run_TRs)
             F_ele = map(self._F_gen, run_TRs)
-            D = []
-            for d_ele in D_ele:
-                D = scipy.linalg.block_diag(D, d_ele)
-            F = []
-            for f_ele in F_ele:
-                F = scipy.linalg.block_diag(F, f_ele)
+            D = scipy.linalg.block_diag(*D_ele)
+            F = scipy.linalg.block_diag(*F_ele)
             # D and F above are templates for constructing
             # the inverse of temporal covariance matrix of noise
         return D, F, run_TRs, n_run
@@ -480,9 +474,7 @@ def _prepare_data_XYX0(self, X, Y, X0, D, F, run_TRs, no_DC=False):
             It will only take effect if X0 is not None.
         """
 
-        X_base = []
-        for r_l in run_TRs:
-            X_base = scipy.linalg.block_diag(X_base, np.ones(r_l)[:, None])
+        X_base = scipy.linalg.block_diag(*map(np.ones, run_TRs)).T
         res = np.linalg.lstsq(X_base, X)
         if np.any(np.isclose(res[1], 0)):
             raise ValueError('Your design matrix appears to have '
@@ -860,7 +852,7 @@ def _fit_RSA_UV(self, X, Y, X0,
                             X0TAX0, est_SNR_AR1_UV**2,
                             n_V, n_T, n_run, rank, n_base)
         est_sigma_AR1_UV = sigma2**0.5
-        est_beta_AR1_UV = est_sigma_AR1_UV * est_SNR_AR1_UV**2 \
+        est_beta_AR1_UV = est_SNR_AR1_UV**2 \
             * np.dot(estU_chlsk_l_AR1_UV, YTAcorrXL_LAMBDA.T)
         est_beta0_AR1_UV = np.einsum(
             'ijk,ki->ji', X0TAX0_i,
@@ -1088,8 +1080,7 @@ def _fit_diagV_noGP(
                     = self._calc_LL(rho1, LTXTAcorrXL, LTXTAcorrY, YTAcorrY,
                                     X0TAX0, current_SNR2,
                                     n_V, n_T, n_run, rank, n_base)
-                betas = current_sigma2**0.5 * current_SNR2 \
-                    * np.dot(L, YTAcorrXL_LAMBDA.T)
+                betas = current_SNR2 * np.dot(L, YTAcorrXL_LAMBDA.T)
                 residuals = Y - np.dot(X, betas)
                 u, s, v = np.linalg.svd(residuals)
                 X0 = u[:, :self.n_nureg]
@@ -1213,8 +1204,7 @@ def _fit_diagV_GP(
                     = self._calc_LL(rho1, LTXTAcorrXL, LTXTAcorrY, YTAcorrY,
                                     X0TAX0, current_SNR2,
                                     n_V, n_T, n_run, rank, n_base)
-                betas = current_sigma2**0.5 * current_SNR2 \
-                    * np.dot(L, YTAcorrXL_LAMBDA.T)
+                betas = current_SNR2 * np.dot(L, YTAcorrXL_LAMBDA.T)
                 residuals = Y - np.dot(X, betas)
                 u, s, v = np.linalg.svd(residuals)
                 X0 = u[:, :self.n_nureg]

diff --git a/...ational_similarity_estimate_example.ipynb → ...ples/reprsimil/bayesian_rsa_example.ipynb b/...ational_similarity_estimate_example.ipynb → ...ples/reprsimil/bayesian_rsa_example.ipynb
@@ -139,7 +139,7 @@
    "outputs": [],
    "source": [
     "noise_bot = 0.5\n",
-    "noise_top = 1.5\n",
+    "noise_top = 5.0\n",
     "noise_level = np.random.rand(n_V) * \\\n",
     "    (noise_top - noise_bot) + noise_bot\n",
     "# The standard deviation of the noise is in the range of [noise_bot, noise_top]\n",
@@ -413,7 +413,6 @@
    "metadata": {},
    "source": [
     "# Fit Bayesian RSA to our simulated data\n",
-    "### The data should be z-scored along time dimension. Otherwise voxels of high variance in noise can dominate the estimation of spatially shared noise component.\n",
     "\n",
     "The nuisance regressors in typical fMRI analysis (such as head motion signal) are replaced by principal components estimated from residuals after subtracting task-related response. n_nureg tells the model how many principal components to keep from the residual as nuisance regressors, in order to account for spatial correlation in noise. \n",
     "If you prefer not using this approach based on principal components of residuals, you can set auto_nuisance=False, and optionally provide your own nuisance regressors as nuisance argument to BRSA.fit()"
@@ -427,9 +426,6 @@
    },
    "outputs": [],
    "source": [
-    "Y_z = scipy.stats.zscore(Y,axis=0)\n",
-    "Y_std = np.std(Y,axis=0)\n",
-    "# Z-scoreing the data\n",
     "\n",
     "brsa = BRSA(GP_space=True, GP_inten=True,\n",
     "            n_nureg=10)\n",
@@ -438,7 +434,7 @@
     "# over both space and intensity.\n",
     "\n",
     "\n",
-    "brsa.fit(X=Y_z, design=design.design_task,\n",
+    "brsa.fit(X=Y, design=design.design_task,\n",
     "         coords=coords_flat, inten=inten, scan_onsets=scan_onsets)\n",
     "# The data to fit should be given to the argument X.\n",
     "# Design matrix goes to design. And so on.\n"
@@ -500,7 +496,7 @@
    "source": [
     "regressor = np.insert(design.design_task,\n",
     "                      0, 1, axis=1)\n",
-    "betas_point = np.linalg.lstsq(regressor, Y_z)[0]\n",
+    "betas_point = np.linalg.lstsq(regressor, Y)[0]\n",
     "point_corr = np.corrcoef(betas_point[1:, :])\n",
     "point_cov = np.cov(betas_point[1:, :])\n",
     "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n",
@@ -634,14 +630,14 @@
    },
    "outputs": [],
    "source": [
-    "plt.scatter(betas_simulated / Y_std, brsa.beta_)\n",
+    "plt.scatter(betas_simulated, brsa.beta_)\n",
     "plt.xlabel('true betas (response amplitudes)')\n",
     "plt.ylabel('recovered betas by Bayesian RSA')\n",
     "ax = plt.gca()\n",
     "plt.show()\n",
     "\n",
     "\n",
-    "plt.scatter(betas_simulated / Y_std, betas_point[1:, :])\n",
+    "plt.scatter(betas_simulated, betas_point[1:, :])\n",
     "plt.xlabel('true betas (response amplitudes)')\n",
     "plt.ylabel('recovered betas by simple regression')\n",
     "ax = plt.gca()\n",
@@ -664,7 +660,7 @@
    },
    "outputs": [],
    "source": [
-    "u, s, v = np.linalg.svd(Y_z - signal / Y_std)\n",
+    "u, s, v = np.linalg.svd(Y - signal)\n",
     "plt.plot(s)\n",
     "plt.xlabel('principal component')\n",
     "plt.ylabel('singular value of simulated noise')\n",

diff --git a/tests/reprsimil/test_brsa.py b/tests/reprsimil/test_brsa.py
@@ -42,17 +42,18 @@ def test_fit():
     design = utils.ReadDesign(fname=file_path)
 
 
-    # concatenate it by 4 times, mimicking 4 runs of itenditcal timing
-    design.design_task = np.tile(design.design_task[:,:-1],[4,1])
-    design.n_TR = design.n_TR * 4
+    # concatenate it by 2 times, mimicking 2 runs of itenditcal timing
+    n_run = 2
+    design.design_task = np.tile(design.design_task[:,:-1],[n_run,1])
+    design.n_TR = design.n_TR * n_run
 
     # start simulating some data
-    n_V = 200
+    n_V = 50
     n_C = np.size(design.design_task,axis=1)
     n_T = design.n_TR
 
     noise_bot = 0.5
-    noise_top = 1.5
+    noise_top = 5.0
     noise_level = np.random.rand(n_V)*(noise_top-noise_bot)+noise_bot
     # noise level is random.
 
@@ -110,14 +111,14 @@ def test_fit():
     signal = np.dot(design.design_task,betas_simulated)
 
     # Adding noise to signal as data
-    Y = signal + noise
+    Y = signal + noise + inten
 
 
     scan_onsets = np.linspace(0,design.n_TR,num=5)
 
 
     # Test fitting with GP prior.
-    brsa = BRSA(GP_space=True,GP_inten=True,verbose=False,n_iter = 200,auto_nuisance=False)
+    brsa = BRSA(GP_space=True,GP_inten=True,verbose=False,n_iter = 5,auto_nuisance=False)
 
     # We also test that it can detect baseline regressor included in the design matrix for task conditions
     wrong_design = np.insert(design.design_task, 0, 1, axis=1)
@@ -167,7 +168,7 @@ def test_fit():
         'the BRSA object should not have parameters of GP if GP is not requested.'
     # GP parameters are not set if not requested
     assert brsa.beta0_.shape[0] == n_nureg, 'Shape of beta0 incorrect'
-    p = scipy.stats.pearsonr(brsa.beta0_[0,:],np.mean(noise,axis=0))[1]
+    p = scipy.stats.pearsonr(brsa.beta0_[0,:],inten)[1]
     assert p < 0.05, 'recovered beta0 does not correlate with the baseline of voxels.'
 
     # Test fitting with GP over just spatial coordinates.
@@ -209,7 +210,7 @@ def test_gradient():
     design.n_TR = design.n_TR * n_run
 
     # start simulating some data
-    n_V = 200
+    n_V = 30
     n_C = np.size(design.design_task,axis=1)
     n_T = design.n_TR