diff --git a/brainiak/reprsimil/brsa.py b/brainiak/reprsimil/brsa.py index 11eec9fe0..59a466ced 100755 --- 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/examples/reprsimil/brsa_representational_similarity_estimate_example.ipynb b/examples/reprsimil/bayesian_rsa_example.ipynb similarity index 97% rename from examples/reprsimil/brsa_representational_similarity_estimate_example.ipynb rename to examples/reprsimil/bayesian_rsa_example.ipynb index f143c45e6..c6c7620ac 100644 --- a/examples/reprsimil/brsa_representational_similarity_estimate_example.ipynb +++ b/examples/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 index ad25e2d53..1d5a9a46e 100755 --- 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