From e4ab872dd2f40983e1bce96904fe22624c3afe71 Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:12:26 -0400 Subject: [PATCH 01/20] a Changes to be committed: new file: brainiak/brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt new file mode 100644 index 000000000..e69de29bb From 06c869a431109610e47d20dcf8e151373aabf5bd Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:20:44 -0400 Subject: [PATCH 02/20] deleted: brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt deleted file mode 100644 index e69de29bb..000000000 From ea7b002323214d7be758de4b398f723d24cc9dc9 Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:12:26 -0400 Subject: [PATCH 03/20] a Changes to be committed: new file: brainiak/brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt new file mode 100644 index 000000000..e69de29bb From ad3af51edc4e1193894525fa88c7531e45b068a4 Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:20:44 -0400 Subject: [PATCH 04/20] deleted: brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt deleted file mode 100644 index e69de29bb..000000000 From 66277e6e91dc6bc7d94845d17cdcb2aeb0cacd1f Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:12:26 -0400 Subject: [PATCH 05/20] a Changes to be committed: new file: brainiak/brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt new file mode 100644 index 000000000..e69de29bb From 1b6127ae18d5c9206fa9f314050ce8eaec487f44 Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:20:44 -0400 Subject: [PATCH 06/20] deleted: brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt deleted file mode 100644 index e69de29bb..000000000 From 4b1a794d0ca848d49ef0736c1656539c6793077d Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:12:26 -0400 Subject: [PATCH 07/20] a Changes to be committed: new file: brainiak/brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt new file mode 100644 index 000000000..e69de29bb From 2ded8dbee5c1eabd5130bf03c48b497f7184147d Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:20:44 -0400 Subject: [PATCH 08/20] deleted: brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt deleted file mode 100644 index e69de29bb..000000000 From 66ebf5f0da4d6dabd87a6f18f9b9c4057aee5180 Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:12:26 -0400 Subject: [PATCH 09/20] a Changes to be committed: new file: brainiak/brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt new file mode 100644 index 000000000..e69de29bb From 6a243448aa23d0562dc5136e2d5454ff0a268366 Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:20:44 -0400 Subject: [PATCH 10/20] deleted: brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt deleted file mode 100644 index e69de29bb..000000000 From 746f66e274006cc3847cc7d50d41fbb6def808e2 Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:12:26 -0400 Subject: [PATCH 11/20] a Changes to be committed: new file: brainiak/brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt new file mode 100644 index 000000000..e69de29bb From 9fad12b0f263dccf38f359c2d2698e926173a2dc Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:20:44 -0400 Subject: [PATCH 12/20] deleted: brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt deleted file mode 100644 index e69de29bb..000000000 From 92a199b893f22da0c4fce8970848d04ee03365ab Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:12:26 -0400 Subject: [PATCH 13/20] a Changes to be committed: new file: brainiak/brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt new file mode 100644 index 000000000..e69de29bb From 5265476bf8b323a2fbe57fad839271cd2aa5705c Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:20:44 -0400 Subject: [PATCH 14/20] deleted: brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt deleted file mode 100644 index e69de29bb..000000000 From 3acf09fb455222ab13280bc13200b03fe1945bb1 Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:12:26 -0400 Subject: [PATCH 15/20] a Changes to be committed: new file: brainiak/brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt new file mode 100644 index 000000000..e69de29bb From bc49c16ebc4e689025803d1c81305f7ded026671 Mon Sep 17 00:00:00 2001 From: lcnature Date: Fri, 2 Sep 2016 10:20:44 -0400 Subject: [PATCH 16/20] deleted: brsa/test.txt --- brainiak/brsa/test.txt | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 brainiak/brsa/test.txt diff --git a/brainiak/brsa/test.txt b/brainiak/brsa/test.txt deleted file mode 100644 index e69de29bb..000000000 From 0661856891e7ae137c8b9e7e00ad467217f1ce93 Mon Sep 17 00:00:00 2001 From: Mingbo Cai Date: Wed, 22 Mar 2017 14:16:13 -0400 Subject: [PATCH 17/20] Changed all usage of scipy.linalg.block_diag in accordance to its change in last November. Fixed the bug in posterior estimation of betas since the review of the major pull request might take longer. --- brainiak/reprsimil/brsa.py | 26 +- ...tational_similarity_estimate_example.ipynb | 706 ------------------ 2 files changed, 8 insertions(+), 724 deletions(-) delete mode 100644 examples/reprsimil/brsa_representational_similarity_estimate_example.ipynb 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/brsa_representational_similarity_estimate_example.ipynb deleted file mode 100644 index f143c45e6..000000000 --- a/examples/reprsimil/brsa_representational_similarity_estimate_example.ipynb +++ /dev/null @@ -1,706 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": { - "collapsed": false - }, - "source": [ - "# This demo shows how to use the Bayesian Representational Similarity Analysis method in brainiak with a simulated dataset." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Load some package which we will use in this demo.\n", - "If you see error related to loading any package, you can install that package. For example, if you use Anaconda, you can use \"conda install matplotlib\" to install matplotlib." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "%matplotlib inline\n", - "import scipy.stats\n", - "import scipy.spatial.distance as spdist\n", - "import numpy as np\n", - "from brainiak.reprsimil.brsa import BRSA\n", - "import brainiak.utils.utils as utils\n", - "import matplotlib.pyplot as plt\n", - "import logging\n", - "np.random.seed(10)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### You might want to keep a log of the output." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "logging.basicConfig(\n", - " level=logging.DEBUG,\n", - " filename='brsa_example.log',\n", - " format='%(relativeCreated)6d %(threadName)s %(message)s')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# We want to simulate some data in which each voxel responds to different task conditions differently, but following a common covariance structure" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Load an example design matrix.\n", - "The user should prepare their design matrix with their favorate software, such as using 3ddeconvolve of AFNI, or using SPM or FSL.\n", - "The design matrix reflects your belief of how fMRI signal should respond to a task (if a voxel does respond).\n", - "The common assumption is that a neural event that you are interested in will elicit a slow hemodynamic response in some voxels. The response peaks around 4-6 seconds after the event onset and dies down more than 12 seconds after the event. Therefore, typically you convolve a time series A, composed of delta (stem) functions reflecting the time of each neural event belonging to the same category (e.g. all trials in which a participant sees a face), with a hemodynamic response function B, to form the hypothetic response of any voxel to such type of neural event.\n", - "For each type of event, such a convoluted time course can be generated. These time courses, put together, are called design matrix, reflecting what we believe a temporal signal would look like, if it exists in any voxel.\n", - "Our goal is to figure out how the (spatial) response pattern of a population of voxels (in an Region of Interest, ROI) are similar or disimilar to different types of tasks (e.g., watching face vs. house, watching different categories of animals, different conditions of a cognitive task). So we need the design matrix in order to estimate the similarity matrix we are interested.\n", - "\n", - "We can use the utility called ReadDesign to read a design matrix generated from AFNI. For design matrix saved as Matlab data file by SPM or or other toolbox, you can use scipy.io.loadmat('YOURFILENAME') and extract the design matrix from the dictionary returned. Basically, the Bayesian RSA in this toolkit just needs a numpy array which is in size of {time points} * {condition}\n", - "\n", - "In typical fMRI analysis, some nuisance regressors such as head motion, baseline time series and slow drift are also entered into regression. In using our method, you should not include such regressors into the design matrix, because the spatial spread of such nuisance regressors might be quite different from the spatial spread of task related signal. Including such nuisance regressors in design matrix might influence the pseudo-SNR map, which in turn influence the estimation of the shared covariance matrix. \n", - "\n", - "### We concatenate the design matrix by 2 times, mimicking 2 runs of identical timing" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "design = utils.ReadDesign(fname=\"example_design.1D\")\n", - "\n", - "n_run = 2\n", - "design.n_TR = design.n_TR * n_run\n", - "design.design_task = np.tile(design.design_task[:,:-1],\n", - " [n_run, 1])\n", - "# The last \"condition\" in design matrix\n", - "# codes for trials subjects made and error.\n", - "# We ignore it here.\n", - "\n", - "\n", - "fig = plt.figure(num=None, figsize=(12, 3),\n", - " dpi=150, facecolor='w', edgecolor='k')\n", - "plt.plot(design.design_task)\n", - "plt.ylim([-0.2, 0.4])\n", - "plt.title('hypothetic fMRI response time courses '\n", - " 'of all conditions in addition to a DC component\\n'\n", - " '(design matrix)')\n", - "plt.xlabel('time')\n", - "plt.show()\n", - "\n", - "n_C = np.size(design.design_task, axis=1)\n", - "# The total number of conditions.\n", - "ROI_edge = 20\n", - "# We simulate \"ROI\" of a square shape\n", - "n_V = ROI_edge**2\n", - "# The total number of simulated voxels\n", - "n_T = design.n_TR\n", - "# The total number of time points,\n", - "# after concatenating all fMRI runs\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## simulate data: noise + signal\n", - "### First, we start with noise, which is Gaussian Process in space and AR(1) in time" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "noise_bot = 0.5\n", - "noise_top = 1.5\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", - "# In fact, we simulate autocorrelated noise with AR(1) model. So the noise_level reflects\n", - "# the independent additive noise at each time point (the \"fresh\" noise)\n", - "\n", - "# AR(1) coefficient\n", - "rho1_top = 0.8\n", - "rho1_bot = -0.2\n", - "rho1 = np.random.rand(n_V) \\\n", - " * (rho1_top - rho1_bot) + rho1_bot\n", - "\n", - "noise_smooth_width = 10.0\n", - "coords = np.mgrid[0:ROI_edge, 0:ROI_edge, 0:1]\n", - "coords_flat = np.reshape(coords,[3, n_V]).T\n", - "dist2 = spdist.squareform(spdist.pdist(coords_flat, 'sqeuclidean'))\n", - "\n", - "# generating noise\n", - "K_noise = noise_level[:, np.newaxis] \\\n", - " * (np.exp(-dist2 / noise_smooth_width**2 / 2.0) \\\n", - " + np.eye(n_V) * 0.1) * noise_level\n", - "# We make spatially correlated noise by generating\n", - "# noise at each time point from a Gaussian Process\n", - "# defined over the coordinates.\n", - "plt.pcolor(K_noise)\n", - "plt.colorbar()\n", - "plt.xlim([0, ROI_edge * ROI_edge])\n", - "plt.ylim([0, ROI_edge * ROI_edge])\n", - "plt.title('Spatial covariance matrix of noise')\n", - "plt.show()\n", - "L_noise = np.linalg.cholesky(K_noise)\n", - "noise = np.zeros([n_T, n_V])\n", - "noise[0, :] = np.dot(L_noise, np.random.randn(n_V))\\\n", - " / np.sqrt(1 - rho1**2)\n", - "for i_t in range(1, n_T):\n", - " noise[i_t, :] = noise[i_t - 1, :] * rho1 \\\n", - " + np.dot(L_noise,np.random.randn(n_V))\n", - "# For each voxel, the noise follows AR(1) process:\n", - "# fresh noise plus a dampened version of noise at\n", - "# the previous time point.\n", - "noise = noise + np.random.randn(n_V)\n", - "fig = plt.figure(num=None, figsize=(12, 2), dpi=150,\n", - " facecolor='w', edgecolor='k')\n", - "plt.plot(noise[:, 0])\n", - "plt.title('noise in an example voxel')\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Then, we simulate signals, assuming the magnitude of response to each condition follows a common covariance matrix. \n", - "#### Our model allows to impose a Gaussian Process prior on the log(SNR) of each voxels. \n", - "What this means is that SNR turn to be smooth and local, but betas (response amplitudes of each voxel to each condition) are not necessarily correlated in space. Intuitively, this is based on the assumption that voxels coding for related aspects of a task turn to be clustered (instead of isolated)\n", - "\n", - "Our Gaussian Process are defined on both the coordinate of a voxel and its mean intensity.\n", - "This means that voxels close together AND have similar intensity should have similar SNR level. Therefore, voxels of white matter but adjacent to gray matters do not necessarily have high SNR level.\n", - "\n", - "If you have an ROI saved as a binary Nifti file, say, with name 'ROI.nii'\n", - "Then you can use nibabel package to load the ROI and the following example code to retrive the coordinates of voxels.\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "##### Note: the following code won't work if you just installed Brainiak and try this demo because ROI.nii does not exist. It just serves as an example for you to retrieve coordinates of voxels in an ROI. You can use the ROI_coords for the argument coords in BRSA.fit()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "# import nibabel\n", - "# ROI = nibabel.load('ROI.nii')\n", - "# I,J,K = ROI.shape \n", - "# all_coords = np.zeros((I, J, K, 3)) \n", - "# all_coords[...,0] = np.arange(I)[:, np.newaxis, np.newaxis] \n", - "# all_coords[...,1] = np.arange(J)[np.newaxis, :, np.newaxis] \n", - "# all_coords[...,2] = np.arange(K)[np.newaxis, np.newaxis, :] \n", - "# ROI_coords = nibabel.affines.apply_affine(\n", - "# ROI.affine, all_coords[ROI.get_data().astype(bool)])\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Let's keep in mind of the pattern of the ideal covariance / correlation below and see how well BRSA can recover their patterns." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "# ideal covariance matrix\n", - "ideal_cov = np.zeros([n_C, n_C])\n", - "ideal_cov = np.eye(n_C) * 0.6\n", - "ideal_cov[8:12, 8:12] = 0.8\n", - "for cond in range(8, 12):\n", - " ideal_cov[cond,cond] = 1\n", - "\n", - "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n", - "plt.pcolor(ideal_cov)\n", - "plt.colorbar()\n", - "plt.xlim([0, 16])\n", - "plt.ylim([0, 16])\n", - "ax = plt.gca()\n", - "ax.set_aspect(1)\n", - "plt.title('ideal covariance matrix')\n", - "plt.show()\n", - "\n", - "std_diag = np.diag(ideal_cov)**0.5\n", - "ideal_corr = ideal_cov / std_diag / std_diag[:, None]\n", - "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n", - "plt.pcolor(ideal_corr)\n", - "plt.colorbar()\n", - "plt.xlim([0, 16])\n", - "plt.ylim([0, 16])\n", - "ax = plt.gca()\n", - "ax.set_aspect(1)\n", - "plt.title('ideal correlation matrix')\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### In the following, pseudo-SNR is generated from a Gaussian Process defined on a \"square\" ROI, just for simplicity of code" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "L_full = np.linalg.cholesky(ideal_cov) \n", - "\n", - "# generating signal\n", - "snr_level = 0.6\n", - "# Notice that accurately speaking this is not SNR.\n", - "# The magnitude of signal depends not only on beta but also on x.\n", - "# (noise_level*snr_level)**2 is the factor multiplied\n", - "# with ideal_cov to form the covariance matrix from which\n", - "# the response amplitudes (beta) of a voxel are drawn from.\n", - "\n", - "tau = 0.8\n", - "# magnitude of Gaussian Process from which the log(SNR) is drawn\n", - "smooth_width = 3.0\n", - "# spatial length scale of the Gaussian Process, unit: voxel\n", - "inten_kernel = 4.0\n", - "# intensity length scale of the Gaussian Process\n", - "# Slightly counter-intuitively, if this parameter is very large,\n", - "# say, much larger than the range of intensities of the voxels,\n", - "# then the smoothness has much small dependency on the intensity.\n", - "\n", - "\n", - "inten = np.random.rand(n_V) * 20.0\n", - "# For simplicity, we just assume that the intensity\n", - "# of all voxels are uniform distributed between 0 and 20\n", - "# parameters of Gaussian process to generate pseuso SNR\n", - "# For curious user, you can also try the following commond\n", - "# to see what an example snr map might look like if the intensity\n", - "# grows linearly in one spatial direction\n", - "\n", - "# inten = coords_flat[:,0] * 2\n", - "\n", - "\n", - "inten_tile = np.tile(inten, [n_V, 1])\n", - "inten_diff2 = (inten_tile - inten_tile.T)**2\n", - "\n", - "K = np.exp(-dist2 / smooth_width**2 / 2.0 \n", - " - inten_diff2 / inten_kernel**2 / 2.0) * tau**2 \\\n", - " + np.eye(n_V) * tau**2 * 0.001\n", - "# A tiny amount is added to the diagonal of\n", - "# the GP covariance matrix to make sure it can be inverted\n", - "L = np.linalg.cholesky(K)\n", - "snr = np.exp(np.dot(L, np.random.randn(n_V))) * snr_level\n", - "sqrt_v = noise_level * snr\n", - "betas_simulated = np.dot(L_full, np.random.randn(n_C, n_V)) * sqrt_v\n", - "signal = np.dot(design.design_task, betas_simulated)\n", - "\n", - "\n", - "Y = signal + noise \n", - "# The data to be fed to the program.\n", - "\n", - "\n", - "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n", - "plt.pcolor(np.reshape(snr, [ROI_edge, ROI_edge]))\n", - "plt.colorbar()\n", - "ax = plt.gca()\n", - "ax.set_aspect(1)\n", - "plt.title('pseudo-SNR in a square \"ROI\"')\n", - "plt.show()\n", - "\n", - "idx = np.argmin(np.abs(snr - np.median(snr)))\n", - "# choose a voxel of medium level SNR.\n", - "fig = plt.figure(num=None, figsize=(12, 4), dpi=150,\n", - " facecolor='w', edgecolor='k')\n", - "noise_plot, = plt.plot(noise[:,idx],'g')\n", - "signal_plot, = plt.plot(signal[:,idx],'r')\n", - "plt.legend([noise_plot, signal_plot], ['noise', 'signal'])\n", - "plt.title('simulated data in an example voxel'\n", - " ' with pseudo-SNR of {}'.format(snr[idx]))\n", - "plt.xlabel('time')\n", - "plt.show()\n", - "\n", - "fig = plt.figure(num=None, figsize=(12, 4), dpi=150,\n", - " facecolor='w', edgecolor='k')\n", - "data_plot, = plt.plot(Y[:,idx],'b')\n", - "plt.legend([data_plot], ['observed data of the voxel'])\n", - "plt.xlabel('time')\n", - "plt.show()\n", - "\n", - "idx = np.argmin(np.abs(snr - np.max(snr)))\n", - "# display the voxel of the highest level SNR.\n", - "fig = plt.figure(num=None, figsize=(12, 4), dpi=150,\n", - " facecolor='w', edgecolor='k')\n", - "noise_plot, = plt.plot(noise[:,idx],'g')\n", - "signal_plot, = plt.plot(signal[:,idx],'r')\n", - "plt.legend([noise_plot, signal_plot], ['noise', 'signal'])\n", - "plt.title('simulated data in the voxel with the highest'\n", - " ' pseudo-SNR of {}'.format(snr[idx]))\n", - "plt.xlabel('time')\n", - "plt.show()\n", - "\n", - "fig = plt.figure(num=None, figsize=(12, 4), dpi=150,\n", - " facecolor='w', edgecolor='k')\n", - "data_plot, = plt.plot(Y[:,idx],'b')\n", - "plt.legend([data_plot], ['observed data of the voxel'])\n", - "plt.xlabel('time')\n", - "plt.show()\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### The reason that the pseudo-SNRs in the example voxels are not too small, while the signal looks much smaller is because we happen to have low amplitudes in our design matrix. The true SNR depends on both the amplitudes in design matrix and the pseudo-SNR. Therefore, be aware that pseudo-SNR does not directly reflects how much signal the data have, but rather a map indicating the relative strength of signal in differerent voxels.\n", - "#### When you have multiple runs, the noise won't be correlated between runs. Therefore, you should tell BRSA when is the onset of each scan. \n", - "#### Note that the data (variable Y above) you feed to BRSA is the concatenation of data from all runs along the time dimension, as a 2-D matrix of time x space" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "scan_onsets = np.int32(np.linspace(0, design.n_TR,num=n_run + 1)[: -1])\n", - "print('scan onsets: {}'.format(scan_onsets))" - ] - }, - { - "cell_type": "markdown", - "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()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "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", - "# Initiate an instance, telling it\n", - "# that we want to impose Gaussian Process prior\n", - "# over both space and intensity.\n", - "\n", - "\n", - "brsa.fit(X=Y_z, 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" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### We can have a look at the estimated similarity in matrix brsa.C_. \n", - "#### We can also compare the ideal covariance above with the one recovered, brsa.U_" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n", - "plt.pcolor(brsa.C_, vmin=-0.1, vmax=1)\n", - "plt.xlim([0, 16])\n", - "plt.ylim([0, 16])\n", - "plt.colorbar()\n", - "ax = plt.gca()\n", - "ax.set_aspect(1)\n", - "plt.title('Estimated correlation structure\\n shared between voxels\\n'\n", - " 'This constitutes the output of Bayesian RSA\\n')\n", - "plt.show()\n", - "\n", - "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n", - "plt.pcolor(brsa.U_)\n", - "plt.xlim([0, 16])\n", - "plt.ylim([0, 16])\n", - "plt.colorbar()\n", - "ax = plt.gca()\n", - "ax.set_aspect(1)\n", - "plt.title('Estimated covariance structure\\n shared between voxels\\n')\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### In contrast, we can have a look of the similarity matrix based on Pearson correlation between point estimates of betas of different conditions.\n", - "#### This is what vanila RSA might give" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "regressor = np.insert(design.design_task,\n", - " 0, 1, axis=1)\n", - "betas_point = np.linalg.lstsq(regressor, Y_z)[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", - "plt.pcolor(point_corr, vmin=-0.1, vmax=1)\n", - "plt.xlim([0, 16])\n", - "plt.ylim([0, 16])\n", - "plt.colorbar()\n", - "ax = plt.gca()\n", - "ax.set_aspect(1)\n", - "plt.title('Correlation structure estimated\\n'\n", - " 'based on point estimates of betas\\n')\n", - "plt.show()\n", - "\n", - "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n", - "plt.pcolor(point_cov)\n", - "plt.xlim([0, 16])\n", - "plt.ylim([0, 16])\n", - "plt.colorbar()\n", - "ax = plt.gca()\n", - "ax.set_aspect(1)\n", - "plt.title('Covariance structure of\\n'\n", - " 'point estimates of betas\\n')\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### We can make a comparison between the estimated SNR map and the true SNR map (normalized)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "fig = plt.figure(num=None, figsize=(5, 5), dpi=100)\n", - "plt.pcolor(np.reshape(brsa.nSNR_, [ROI_edge, ROI_edge]), vmin=0, vmax=5)\n", - "plt.colorbar()\n", - "ax = plt.gca()\n", - "ax.set_aspect(1)\n", - "ax.set_title('estimated pseudo-SNR')\n", - "plt.show()\n", - "\n", - "fig = plt.figure(num=None, figsize=(5, 5), dpi=100)\n", - "plt.pcolor(np.reshape(snr / np.exp(np.mean(np.log(snr))),\n", - " [ROI_edge, ROI_edge]), vmin=0, vmax=5)\n", - "plt.colorbar()\n", - "ax = plt.gca()\n", - "ax.set_aspect(1)\n", - "ax.set_title('true normalized pseudo-SNR')\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "RMS_BRSA = np.mean((brsa.C_ - ideal_corr)**2)**0.5\n", - "RMS_RSA = np.mean((point_corr - ideal_corr)**2)**0.5\n", - "print('RMS error of Bayesian RSA: {}'.format(RMS_BRSA))\n", - "print('RMS error of standard RSA: {}'.format(RMS_RSA))\n", - "print('Recovered spatial smoothness length scale: '\n", - " '{}, vs. true value: {}'.format(brsa.lGPspace_, smooth_width))\n", - "print('Recovered intensity smoothness length scale: '\n", - " '{}, vs. true value: {}'.format(brsa.lGPinten_, inten_kernel))\n", - "print('Recovered standard deviation of GP prior: '\n", - " '{}, vs. true value: {}'.format(brsa.bGP_, tau))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Empirically, the smoothness turns to be over-estimated when signal is weak." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### We can also look at how other parameters are recovered." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "\n", - "plt.scatter(rho1, brsa.rho_)\n", - "plt.xlabel('true AR(1) coefficients')\n", - "plt.ylabel('recovered AR(1) coefficients')\n", - "ax = plt.gca()\n", - "ax.set_aspect(1)\n", - "plt.show()\n", - "\n", - "plt.scatter(np.log(snr) - np.mean(np.log(snr)),\n", - " np.log(brsa.nSNR_))\n", - "plt.xlabel('true normalized log SNR')\n", - "plt.ylabel('recovered log pseudo-SNR')\n", - "ax = plt.gca()\n", - "ax.set_aspect(1)\n", - "plt.show()\n", - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Even though the variation reduced in estimated pseudo-SNR (due to overestimation of smoothness of the GP prior under low SNR situation), betas recovered by the model has higher correlation with true betas than doing simple regression, shown below. Obiously there is shrinkage of the estimated betas, as a result of variance-bias tradeoff. But we think such shrinkage does preserve the patterns of betas, and therefore the result is suitable to be further used for decoding purpose." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "plt.scatter(betas_simulated / Y_std, 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.xlabel('true betas (response amplitudes)')\n", - "plt.ylabel('recovered betas by simple regression')\n", - "ax = plt.gca()\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### The singular decomposition of noise, and the comparison between the first principal component of noise and the first principal component returned by the model.\n", - "Apparently one can imagine that the choice of the number of principal components used as nuisance regressors can influence the result. If you just choose 1 or 2, perhaps only the baseline and drift would be captured. But including too many nuisance regressors would slow the fitting speed and might have risk of overfitting. The users might consider starting in the range of 5-20. In future, we may consider using cross validation to determine the number of nuisance regressors, similar as in GLMdenoise (http://kendrickkay.net/GLMdenoise/)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "u, s, v = np.linalg.svd(Y_z - signal / Y_std)\n", - "plt.plot(s)\n", - "plt.xlabel('principal component')\n", - "plt.ylabel('singular value of simulated noise')\n", - "plt.show()\n", - "\n", - "plt.pcolor(np.reshape(v[0,:], [ROI_edge, ROI_edge]))\n", - "plt.title('Weights of the first principal component in noise')\n", - "plt.show()\n", - "\n", - "\n", - "plt.pcolor(np.reshape(brsa.beta0_[0,:], [ROI_edge, ROI_edge]))\n", - "plt.title('Weights of the first recovered principal component in noise')\n", - "plt.show()" - ] - } - ], - "metadata": { - "anaconda-cloud": {}, - "kernelspec": { - "display_name": "Python [default]", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.5.2" - } - }, - "nbformat": 4, - "nbformat_minor": 0 -} From 42e90ae17117e9c6ec669bc97c38575a8461630a Mon Sep 17 00:00:00 2001 From: Mingbo Cai Date: Wed, 22 Mar 2017 14:50:03 -0400 Subject: [PATCH 18/20] reduce testing time for brsa --- tests/reprsimil/test_brsa.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) 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 From 39a9df69bebecea383b3969a19aa3b6cd995529f Mon Sep 17 00:00:00 2001 From: Mingbo Cai Date: Wed, 22 Mar 2017 14:59:11 -0400 Subject: [PATCH 19/20] renamed BRSA example --- examples/reprsimil/Bayesian_RSA_example.ipynb | 702 ++++++++++++++++++ 1 file changed, 702 insertions(+) create mode 100644 examples/reprsimil/Bayesian_RSA_example.ipynb diff --git a/examples/reprsimil/Bayesian_RSA_example.ipynb b/examples/reprsimil/Bayesian_RSA_example.ipynb new file mode 100644 index 000000000..c6c7620ac --- /dev/null +++ b/examples/reprsimil/Bayesian_RSA_example.ipynb @@ -0,0 +1,702 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "collapsed": false + }, + "source": [ + "# This demo shows how to use the Bayesian Representational Similarity Analysis method in brainiak with a simulated dataset." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Load some package which we will use in this demo.\n", + "If you see error related to loading any package, you can install that package. For example, if you use Anaconda, you can use \"conda install matplotlib\" to install matplotlib." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import scipy.stats\n", + "import scipy.spatial.distance as spdist\n", + "import numpy as np\n", + "from brainiak.reprsimil.brsa import BRSA\n", + "import brainiak.utils.utils as utils\n", + "import matplotlib.pyplot as plt\n", + "import logging\n", + "np.random.seed(10)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### You might want to keep a log of the output." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "logging.basicConfig(\n", + " level=logging.DEBUG,\n", + " filename='brsa_example.log',\n", + " format='%(relativeCreated)6d %(threadName)s %(message)s')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# We want to simulate some data in which each voxel responds to different task conditions differently, but following a common covariance structure" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Load an example design matrix.\n", + "The user should prepare their design matrix with their favorate software, such as using 3ddeconvolve of AFNI, or using SPM or FSL.\n", + "The design matrix reflects your belief of how fMRI signal should respond to a task (if a voxel does respond).\n", + "The common assumption is that a neural event that you are interested in will elicit a slow hemodynamic response in some voxels. The response peaks around 4-6 seconds after the event onset and dies down more than 12 seconds after the event. Therefore, typically you convolve a time series A, composed of delta (stem) functions reflecting the time of each neural event belonging to the same category (e.g. all trials in which a participant sees a face), with a hemodynamic response function B, to form the hypothetic response of any voxel to such type of neural event.\n", + "For each type of event, such a convoluted time course can be generated. These time courses, put together, are called design matrix, reflecting what we believe a temporal signal would look like, if it exists in any voxel.\n", + "Our goal is to figure out how the (spatial) response pattern of a population of voxels (in an Region of Interest, ROI) are similar or disimilar to different types of tasks (e.g., watching face vs. house, watching different categories of animals, different conditions of a cognitive task). So we need the design matrix in order to estimate the similarity matrix we are interested.\n", + "\n", + "We can use the utility called ReadDesign to read a design matrix generated from AFNI. For design matrix saved as Matlab data file by SPM or or other toolbox, you can use scipy.io.loadmat('YOURFILENAME') and extract the design matrix from the dictionary returned. Basically, the Bayesian RSA in this toolkit just needs a numpy array which is in size of {time points} * {condition}\n", + "\n", + "In typical fMRI analysis, some nuisance regressors such as head motion, baseline time series and slow drift are also entered into regression. In using our method, you should not include such regressors into the design matrix, because the spatial spread of such nuisance regressors might be quite different from the spatial spread of task related signal. Including such nuisance regressors in design matrix might influence the pseudo-SNR map, which in turn influence the estimation of the shared covariance matrix. \n", + "\n", + "### We concatenate the design matrix by 2 times, mimicking 2 runs of identical timing" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "design = utils.ReadDesign(fname=\"example_design.1D\")\n", + "\n", + "n_run = 2\n", + "design.n_TR = design.n_TR * n_run\n", + "design.design_task = np.tile(design.design_task[:,:-1],\n", + " [n_run, 1])\n", + "# The last \"condition\" in design matrix\n", + "# codes for trials subjects made and error.\n", + "# We ignore it here.\n", + "\n", + "\n", + "fig = plt.figure(num=None, figsize=(12, 3),\n", + " dpi=150, facecolor='w', edgecolor='k')\n", + "plt.plot(design.design_task)\n", + "plt.ylim([-0.2, 0.4])\n", + "plt.title('hypothetic fMRI response time courses '\n", + " 'of all conditions in addition to a DC component\\n'\n", + " '(design matrix)')\n", + "plt.xlabel('time')\n", + "plt.show()\n", + "\n", + "n_C = np.size(design.design_task, axis=1)\n", + "# The total number of conditions.\n", + "ROI_edge = 20\n", + "# We simulate \"ROI\" of a square shape\n", + "n_V = ROI_edge**2\n", + "# The total number of simulated voxels\n", + "n_T = design.n_TR\n", + "# The total number of time points,\n", + "# after concatenating all fMRI runs\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## simulate data: noise + signal\n", + "### First, we start with noise, which is Gaussian Process in space and AR(1) in time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "noise_bot = 0.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", + "# In fact, we simulate autocorrelated noise with AR(1) model. So the noise_level reflects\n", + "# the independent additive noise at each time point (the \"fresh\" noise)\n", + "\n", + "# AR(1) coefficient\n", + "rho1_top = 0.8\n", + "rho1_bot = -0.2\n", + "rho1 = np.random.rand(n_V) \\\n", + " * (rho1_top - rho1_bot) + rho1_bot\n", + "\n", + "noise_smooth_width = 10.0\n", + "coords = np.mgrid[0:ROI_edge, 0:ROI_edge, 0:1]\n", + "coords_flat = np.reshape(coords,[3, n_V]).T\n", + "dist2 = spdist.squareform(spdist.pdist(coords_flat, 'sqeuclidean'))\n", + "\n", + "# generating noise\n", + "K_noise = noise_level[:, np.newaxis] \\\n", + " * (np.exp(-dist2 / noise_smooth_width**2 / 2.0) \\\n", + " + np.eye(n_V) * 0.1) * noise_level\n", + "# We make spatially correlated noise by generating\n", + "# noise at each time point from a Gaussian Process\n", + "# defined over the coordinates.\n", + "plt.pcolor(K_noise)\n", + "plt.colorbar()\n", + "plt.xlim([0, ROI_edge * ROI_edge])\n", + "plt.ylim([0, ROI_edge * ROI_edge])\n", + "plt.title('Spatial covariance matrix of noise')\n", + "plt.show()\n", + "L_noise = np.linalg.cholesky(K_noise)\n", + "noise = np.zeros([n_T, n_V])\n", + "noise[0, :] = np.dot(L_noise, np.random.randn(n_V))\\\n", + " / np.sqrt(1 - rho1**2)\n", + "for i_t in range(1, n_T):\n", + " noise[i_t, :] = noise[i_t - 1, :] * rho1 \\\n", + " + np.dot(L_noise,np.random.randn(n_V))\n", + "# For each voxel, the noise follows AR(1) process:\n", + "# fresh noise plus a dampened version of noise at\n", + "# the previous time point.\n", + "noise = noise + np.random.randn(n_V)\n", + "fig = plt.figure(num=None, figsize=(12, 2), dpi=150,\n", + " facecolor='w', edgecolor='k')\n", + "plt.plot(noise[:, 0])\n", + "plt.title('noise in an example voxel')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Then, we simulate signals, assuming the magnitude of response to each condition follows a common covariance matrix. \n", + "#### Our model allows to impose a Gaussian Process prior on the log(SNR) of each voxels. \n", + "What this means is that SNR turn to be smooth and local, but betas (response amplitudes of each voxel to each condition) are not necessarily correlated in space. Intuitively, this is based on the assumption that voxels coding for related aspects of a task turn to be clustered (instead of isolated)\n", + "\n", + "Our Gaussian Process are defined on both the coordinate of a voxel and its mean intensity.\n", + "This means that voxels close together AND have similar intensity should have similar SNR level. Therefore, voxels of white matter but adjacent to gray matters do not necessarily have high SNR level.\n", + "\n", + "If you have an ROI saved as a binary Nifti file, say, with name 'ROI.nii'\n", + "Then you can use nibabel package to load the ROI and the following example code to retrive the coordinates of voxels.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Note: the following code won't work if you just installed Brainiak and try this demo because ROI.nii does not exist. It just serves as an example for you to retrieve coordinates of voxels in an ROI. You can use the ROI_coords for the argument coords in BRSA.fit()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# import nibabel\n", + "# ROI = nibabel.load('ROI.nii')\n", + "# I,J,K = ROI.shape \n", + "# all_coords = np.zeros((I, J, K, 3)) \n", + "# all_coords[...,0] = np.arange(I)[:, np.newaxis, np.newaxis] \n", + "# all_coords[...,1] = np.arange(J)[np.newaxis, :, np.newaxis] \n", + "# all_coords[...,2] = np.arange(K)[np.newaxis, np.newaxis, :] \n", + "# ROI_coords = nibabel.affines.apply_affine(\n", + "# ROI.affine, all_coords[ROI.get_data().astype(bool)])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Let's keep in mind of the pattern of the ideal covariance / correlation below and see how well BRSA can recover their patterns." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# ideal covariance matrix\n", + "ideal_cov = np.zeros([n_C, n_C])\n", + "ideal_cov = np.eye(n_C) * 0.6\n", + "ideal_cov[8:12, 8:12] = 0.8\n", + "for cond in range(8, 12):\n", + " ideal_cov[cond,cond] = 1\n", + "\n", + "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n", + "plt.pcolor(ideal_cov)\n", + "plt.colorbar()\n", + "plt.xlim([0, 16])\n", + "plt.ylim([0, 16])\n", + "ax = plt.gca()\n", + "ax.set_aspect(1)\n", + "plt.title('ideal covariance matrix')\n", + "plt.show()\n", + "\n", + "std_diag = np.diag(ideal_cov)**0.5\n", + "ideal_corr = ideal_cov / std_diag / std_diag[:, None]\n", + "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n", + "plt.pcolor(ideal_corr)\n", + "plt.colorbar()\n", + "plt.xlim([0, 16])\n", + "plt.ylim([0, 16])\n", + "ax = plt.gca()\n", + "ax.set_aspect(1)\n", + "plt.title('ideal correlation matrix')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### In the following, pseudo-SNR is generated from a Gaussian Process defined on a \"square\" ROI, just for simplicity of code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "L_full = np.linalg.cholesky(ideal_cov) \n", + "\n", + "# generating signal\n", + "snr_level = 0.6\n", + "# Notice that accurately speaking this is not SNR.\n", + "# The magnitude of signal depends not only on beta but also on x.\n", + "# (noise_level*snr_level)**2 is the factor multiplied\n", + "# with ideal_cov to form the covariance matrix from which\n", + "# the response amplitudes (beta) of a voxel are drawn from.\n", + "\n", + "tau = 0.8\n", + "# magnitude of Gaussian Process from which the log(SNR) is drawn\n", + "smooth_width = 3.0\n", + "# spatial length scale of the Gaussian Process, unit: voxel\n", + "inten_kernel = 4.0\n", + "# intensity length scale of the Gaussian Process\n", + "# Slightly counter-intuitively, if this parameter is very large,\n", + "# say, much larger than the range of intensities of the voxels,\n", + "# then the smoothness has much small dependency on the intensity.\n", + "\n", + "\n", + "inten = np.random.rand(n_V) * 20.0\n", + "# For simplicity, we just assume that the intensity\n", + "# of all voxels are uniform distributed between 0 and 20\n", + "# parameters of Gaussian process to generate pseuso SNR\n", + "# For curious user, you can also try the following commond\n", + "# to see what an example snr map might look like if the intensity\n", + "# grows linearly in one spatial direction\n", + "\n", + "# inten = coords_flat[:,0] * 2\n", + "\n", + "\n", + "inten_tile = np.tile(inten, [n_V, 1])\n", + "inten_diff2 = (inten_tile - inten_tile.T)**2\n", + "\n", + "K = np.exp(-dist2 / smooth_width**2 / 2.0 \n", + " - inten_diff2 / inten_kernel**2 / 2.0) * tau**2 \\\n", + " + np.eye(n_V) * tau**2 * 0.001\n", + "# A tiny amount is added to the diagonal of\n", + "# the GP covariance matrix to make sure it can be inverted\n", + "L = np.linalg.cholesky(K)\n", + "snr = np.exp(np.dot(L, np.random.randn(n_V))) * snr_level\n", + "sqrt_v = noise_level * snr\n", + "betas_simulated = np.dot(L_full, np.random.randn(n_C, n_V)) * sqrt_v\n", + "signal = np.dot(design.design_task, betas_simulated)\n", + "\n", + "\n", + "Y = signal + noise \n", + "# The data to be fed to the program.\n", + "\n", + "\n", + "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n", + "plt.pcolor(np.reshape(snr, [ROI_edge, ROI_edge]))\n", + "plt.colorbar()\n", + "ax = plt.gca()\n", + "ax.set_aspect(1)\n", + "plt.title('pseudo-SNR in a square \"ROI\"')\n", + "plt.show()\n", + "\n", + "idx = np.argmin(np.abs(snr - np.median(snr)))\n", + "# choose a voxel of medium level SNR.\n", + "fig = plt.figure(num=None, figsize=(12, 4), dpi=150,\n", + " facecolor='w', edgecolor='k')\n", + "noise_plot, = plt.plot(noise[:,idx],'g')\n", + "signal_plot, = plt.plot(signal[:,idx],'r')\n", + "plt.legend([noise_plot, signal_plot], ['noise', 'signal'])\n", + "plt.title('simulated data in an example voxel'\n", + " ' with pseudo-SNR of {}'.format(snr[idx]))\n", + "plt.xlabel('time')\n", + "plt.show()\n", + "\n", + "fig = plt.figure(num=None, figsize=(12, 4), dpi=150,\n", + " facecolor='w', edgecolor='k')\n", + "data_plot, = plt.plot(Y[:,idx],'b')\n", + "plt.legend([data_plot], ['observed data of the voxel'])\n", + "plt.xlabel('time')\n", + "plt.show()\n", + "\n", + "idx = np.argmin(np.abs(snr - np.max(snr)))\n", + "# display the voxel of the highest level SNR.\n", + "fig = plt.figure(num=None, figsize=(12, 4), dpi=150,\n", + " facecolor='w', edgecolor='k')\n", + "noise_plot, = plt.plot(noise[:,idx],'g')\n", + "signal_plot, = plt.plot(signal[:,idx],'r')\n", + "plt.legend([noise_plot, signal_plot], ['noise', 'signal'])\n", + "plt.title('simulated data in the voxel with the highest'\n", + " ' pseudo-SNR of {}'.format(snr[idx]))\n", + "plt.xlabel('time')\n", + "plt.show()\n", + "\n", + "fig = plt.figure(num=None, figsize=(12, 4), dpi=150,\n", + " facecolor='w', edgecolor='k')\n", + "data_plot, = plt.plot(Y[:,idx],'b')\n", + "plt.legend([data_plot], ['observed data of the voxel'])\n", + "plt.xlabel('time')\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### The reason that the pseudo-SNRs in the example voxels are not too small, while the signal looks much smaller is because we happen to have low amplitudes in our design matrix. The true SNR depends on both the amplitudes in design matrix and the pseudo-SNR. Therefore, be aware that pseudo-SNR does not directly reflects how much signal the data have, but rather a map indicating the relative strength of signal in differerent voxels.\n", + "#### When you have multiple runs, the noise won't be correlated between runs. Therefore, you should tell BRSA when is the onset of each scan. \n", + "#### Note that the data (variable Y above) you feed to BRSA is the concatenation of data from all runs along the time dimension, as a 2-D matrix of time x space" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "scan_onsets = np.int32(np.linspace(0, design.n_TR,num=n_run + 1)[: -1])\n", + "print('scan onsets: {}'.format(scan_onsets))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Fit Bayesian RSA to our simulated data\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()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "\n", + "brsa = BRSA(GP_space=True, GP_inten=True,\n", + " n_nureg=10)\n", + "# Initiate an instance, telling it\n", + "# that we want to impose Gaussian Process prior\n", + "# over both space and intensity.\n", + "\n", + "\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" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### We can have a look at the estimated similarity in matrix brsa.C_. \n", + "#### We can also compare the ideal covariance above with the one recovered, brsa.U_" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n", + "plt.pcolor(brsa.C_, vmin=-0.1, vmax=1)\n", + "plt.xlim([0, 16])\n", + "plt.ylim([0, 16])\n", + "plt.colorbar()\n", + "ax = plt.gca()\n", + "ax.set_aspect(1)\n", + "plt.title('Estimated correlation structure\\n shared between voxels\\n'\n", + " 'This constitutes the output of Bayesian RSA\\n')\n", + "plt.show()\n", + "\n", + "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n", + "plt.pcolor(brsa.U_)\n", + "plt.xlim([0, 16])\n", + "plt.ylim([0, 16])\n", + "plt.colorbar()\n", + "ax = plt.gca()\n", + "ax.set_aspect(1)\n", + "plt.title('Estimated covariance structure\\n shared between voxels\\n')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### In contrast, we can have a look of the similarity matrix based on Pearson correlation between point estimates of betas of different conditions.\n", + "#### This is what vanila RSA might give" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "regressor = np.insert(design.design_task,\n", + " 0, 1, axis=1)\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", + "plt.pcolor(point_corr, vmin=-0.1, vmax=1)\n", + "plt.xlim([0, 16])\n", + "plt.ylim([0, 16])\n", + "plt.colorbar()\n", + "ax = plt.gca()\n", + "ax.set_aspect(1)\n", + "plt.title('Correlation structure estimated\\n'\n", + " 'based on point estimates of betas\\n')\n", + "plt.show()\n", + "\n", + "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n", + "plt.pcolor(point_cov)\n", + "plt.xlim([0, 16])\n", + "plt.ylim([0, 16])\n", + "plt.colorbar()\n", + "ax = plt.gca()\n", + "ax.set_aspect(1)\n", + "plt.title('Covariance structure of\\n'\n", + " 'point estimates of betas\\n')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### We can make a comparison between the estimated SNR map and the true SNR map (normalized)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "fig = plt.figure(num=None, figsize=(5, 5), dpi=100)\n", + "plt.pcolor(np.reshape(brsa.nSNR_, [ROI_edge, ROI_edge]), vmin=0, vmax=5)\n", + "plt.colorbar()\n", + "ax = plt.gca()\n", + "ax.set_aspect(1)\n", + "ax.set_title('estimated pseudo-SNR')\n", + "plt.show()\n", + "\n", + "fig = plt.figure(num=None, figsize=(5, 5), dpi=100)\n", + "plt.pcolor(np.reshape(snr / np.exp(np.mean(np.log(snr))),\n", + " [ROI_edge, ROI_edge]), vmin=0, vmax=5)\n", + "plt.colorbar()\n", + "ax = plt.gca()\n", + "ax.set_aspect(1)\n", + "ax.set_title('true normalized pseudo-SNR')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "RMS_BRSA = np.mean((brsa.C_ - ideal_corr)**2)**0.5\n", + "RMS_RSA = np.mean((point_corr - ideal_corr)**2)**0.5\n", + "print('RMS error of Bayesian RSA: {}'.format(RMS_BRSA))\n", + "print('RMS error of standard RSA: {}'.format(RMS_RSA))\n", + "print('Recovered spatial smoothness length scale: '\n", + " '{}, vs. true value: {}'.format(brsa.lGPspace_, smooth_width))\n", + "print('Recovered intensity smoothness length scale: '\n", + " '{}, vs. true value: {}'.format(brsa.lGPinten_, inten_kernel))\n", + "print('Recovered standard deviation of GP prior: '\n", + " '{}, vs. true value: {}'.format(brsa.bGP_, tau))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Empirically, the smoothness turns to be over-estimated when signal is weak." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### We can also look at how other parameters are recovered." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "\n", + "plt.scatter(rho1, brsa.rho_)\n", + "plt.xlabel('true AR(1) coefficients')\n", + "plt.ylabel('recovered AR(1) coefficients')\n", + "ax = plt.gca()\n", + "ax.set_aspect(1)\n", + "plt.show()\n", + "\n", + "plt.scatter(np.log(snr) - np.mean(np.log(snr)),\n", + " np.log(brsa.nSNR_))\n", + "plt.xlabel('true normalized log SNR')\n", + "plt.ylabel('recovered log pseudo-SNR')\n", + "ax = plt.gca()\n", + "ax.set_aspect(1)\n", + "plt.show()\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Even though the variation reduced in estimated pseudo-SNR (due to overestimation of smoothness of the GP prior under low SNR situation), betas recovered by the model has higher correlation with true betas than doing simple regression, shown below. Obiously there is shrinkage of the estimated betas, as a result of variance-bias tradeoff. But we think such shrinkage does preserve the patterns of betas, and therefore the result is suitable to be further used for decoding purpose." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "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, betas_point[1:, :])\n", + "plt.xlabel('true betas (response amplitudes)')\n", + "plt.ylabel('recovered betas by simple regression')\n", + "ax = plt.gca()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The singular decomposition of noise, and the comparison between the first principal component of noise and the first principal component returned by the model.\n", + "Apparently one can imagine that the choice of the number of principal components used as nuisance regressors can influence the result. If you just choose 1 or 2, perhaps only the baseline and drift would be captured. But including too many nuisance regressors would slow the fitting speed and might have risk of overfitting. The users might consider starting in the range of 5-20. In future, we may consider using cross validation to determine the number of nuisance regressors, similar as in GLMdenoise (http://kendrickkay.net/GLMdenoise/)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "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", + "plt.show()\n", + "\n", + "plt.pcolor(np.reshape(v[0,:], [ROI_edge, ROI_edge]))\n", + "plt.title('Weights of the first principal component in noise')\n", + "plt.show()\n", + "\n", + "\n", + "plt.pcolor(np.reshape(brsa.beta0_[0,:], [ROI_edge, ROI_edge]))\n", + "plt.title('Weights of the first recovered principal component in noise')\n", + "plt.show()" + ] + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python [default]", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} From 0554609ff2c98d495ec087bbe8c6c6de1b49c236 Mon Sep 17 00:00:00 2001 From: Mingbo Cai Date: Wed, 22 Mar 2017 15:16:23 -0400 Subject: [PATCH 20/20] rename brsa example --- .../{Bayesian_RSA_example.ipynb => bayesian_rsa_example.ipynb} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename examples/reprsimil/{Bayesian_RSA_example.ipynb => bayesian_rsa_example.ipynb} (100%) diff --git a/examples/reprsimil/Bayesian_RSA_example.ipynb b/examples/reprsimil/bayesian_rsa_example.ipynb similarity index 100% rename from examples/reprsimil/Bayesian_RSA_example.ipynb rename to examples/reprsimil/bayesian_rsa_example.ipynb