# - compatibility with Python 3
from __future__ import print_function # print('me') instead of print 'me'
from __future__ import division # 1/2 == 0.5, not 0
# - show figures inside the notebook
%matplotlib inline
# - import common modules
import numpy as np # the Python array package
import matplotlib.pyplot as plt # the Python plotting package
range is a built-in function in Python.
"built-in" means you can use range without having to import it.
In Python 2, range returns a list:
range(10)
type(range(10))
In Python 3, range returns a special object that behaves rather like a list, but doesn't create the values until you need them. This is more efficient in memory.
Here is range in Python 3:
In [1]: range(10)
Out[1]: range(0, 10)
In [2]: type(range(10))
Out[2]: range
In [3]: list(range(10))
Out[3]: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]There are three ways of calling range - with 1, 2 or 3 arguments.
With just one argument, the argument value is the length of the returned list, or, equivalently, the largest value in the list, plus 1 (as when you do slicing), or, equivalently, the value stop such that no list element is >= stop:
range(5)
range(2)
If you pass 2 arguments, the first value is the first value in the returned list, and the second value is the stop value:
range(4, 10)
If you pass 3 arguments, the first value is the first returned value, the second is the stop value, and the third is the step size.
range(4, 10, 2)
range(4, 11, 2)
arange in numpy is very like range with two important differences:
arange returns an array rather than a list (in Python 2 or Python 3);arange arguments can be floating point values.np.arange(4, 11, 2)
np.arange(4, 11, 0.5)
Because arange returns arrays, you can also multiply by the step size and maybe add the start value to create equally spaced vectors:
np.arange(10) * 0.5 + 4
We have so far seen floating point numbers in this kind of format in text files : 314.15926.
It is common for numpy and MATLAB to write floating point numbers in exponential format, like this 3.1415926e02. You can read this as meaning 3.1415926 * 10**2 == 3.1415926 * 100 == 314.15926. Both np.loadtxt and the float function will happily convert these to floating point values in memory:
float('2.2e02')
We have already seen that we can plot two things onto the same plot:
x = np.arange(0, 2*np.pi, 0.1)
y1 = np.sin(x)
y2 = np.cos(x)
plt.plot(x, y1)
plt.plot(x, y2)
It can be hard to remember which plot is which (is green the sin or the cos)? label and legend are useful for that:
plt.plot(x, y1, label='sin')
plt.plot(x, y2, label='cos')
plt.legend()
from scipy.stats import gamma
Scipy is a large library of scientific routines that builds on top of numpy.
You can think of numpy as being a subset of MATLAB, and numpy + scipy as being as being roughly equivalent to MATLAB plus the MATLAB toolboxes.
Scipy has many sub-packages, for doing things like reading MATLAB .mat files (scipy.io) or working with sparse matrices (scipy.sparse). We are going to be using the functions and objects for working with statistical distributions in scipy.stats:
import scipy.stats
scipy.stats contains objects for working with many different distributions. We are going to be working with scipy.stats.gamma, which implements the gamma distribution.
In particular we are interested in the probability density function (PDF) of the gamma distribution.
Because this is a function, we need to pass it an array of values at which it will evaluate.
We can also pass various parameters which change the shape, location and width of the gamma PDF. The most important is the first parameter (after the input array) known as the shape parameter ($k$ in the wikipedia page on gamma distributions).
First we chose some x values at which to sample from the gamma PDF:
x = np.arange(0, 25, 0.1)
Next we plot the gamma PDF for shape values of 2, 4, 8, 12.
plt.plot(x, gamma.pdf(x, 2), label='k=2')
plt.plot(x, gamma.pdf(x, 4), label='k=4')
plt.plot(x, gamma.pdf(x, 8), label='k=8')
plt.plot(x, gamma.pdf(x, 12), label='k=12')
plt.legend()
We will see more from the statistics distributions next week.
By "correlation coefficient", we mean the Pearson product moment correlation coefficient.
Let's say I have two arrays, x and y:
x = np.array([ 20.57, 18.33, 20.8 , 18.77, 18.46, 21.09, 19.96, 20.61,
18.73, 19.6 , 18.15, 19.7 , 20.36, 20.39, 19.67, 20.73])
y = np.array([ 28.98, 30.13, 29.64, 29.87, 28.78, 33.48, 31.36, 31.04,
30.43, 31.69, 27.08, 29.25, 29.5 , 30.04, 29.74, 30.06])
plt.plot(x, y, 'o')
The numpy routine np.corrcoef gives me the correlation matrix between x and y:
corr_matrix = np.corrcoef(x, y)
corr_matrix
corr_matrix[0, 0] is the correlation of x with itself, corr_matrix[1, 0] (and corr_matrix[0, 1]) is the correlation of x with y.
See the correlation exercises for more on correlation.
Be careful - np.correlate is something else entirely.
If I have two vectors $\mathbf{a}$ with elements $a_0, a_1, ... a_{n-1}$, and $\mathbf{b}$ with elements $b_0, b_1, ... b_{n-1}$ then the dot product is defined as:
$$ \mathbf{a}\cdot \mathbf{b} = \sum_{i=0}^{n-1} a_ib_i = a_0b_0 + a_1b_1 + \cdots + a_{n-1}b_{n-1} $$In code:
a = np.arange(5)
b = np.arange(10, 15)
np.dot(a, b)
# The same thing as
np.sum(a * b) # Elementwise multiplication
Matrix multiplication operates by taking dot products of the rows of the first array (matrix) with the columns of the second.
Let's say I have a matrix $\mathbf{X}$, and $X_{i,:}$ is row $i$ in $\mathbf{X}$. I have a matrix $\mathbf{Y}$, and $Y_{:,j}$ is column $j$ in $\mathbf{Y}$. The output matrix $\mathbf{Z} = \mathbf{X} \mathbf{Y}$ has entry $Z_{i,j} = X_{i,:} \cdot Y_{:, j}$. We will see this often over the next few weeks.
X = np.array([[0, 1, 2], [3, 4, 5]])
X
Y = np.array([[7, 8], [9, 10], [11, 12]])
Y
X.dot(Y)
X[0, :].dot(Y[:, 0])
X[1, :].dot(Y[:, 0])
Let's say I have a neural vector with a couple of spikes:
times = np.arange(0, 60, 0.5) # samples every 0.5 seconds
neural_vector = np.zeros(times.shape)
neural_vector[10] = 1 # At 5 seconds
neural_vector[20] = 1 # At 10 seconds
plt.plot(times, neural_vector)
Then I have my HRF function, sampled every half second, to match:
def hrf(t):
"A hemodynamic response function"
values = t ** 8.6 * np.exp(-t / 0.547)
# Scale max to 1
return values / np.max(values)
hrf_times = np.arange(0, 20, 0.5)
hrf_samples = hrf(hrf_times)
plt.plot(hrf_times, hrf_samples)
The input neural vector is length 120, and the HRF vector is length 40:
print(len(neural_vector))
print(len(hrf_samples))
For reasons that may be familiar to you now, when we convolve the neural vector with the hrf signal, we get an output that is length 120 + 40 - 1 = 159:
hemodynamic_prediction = np.convolve(neural_vector, hrf_samples)
len(hemodynamic_prediction)
This is because of the HRF falling off the end of the input vector. The value at index 120 in the new vector refers to time 60, and value 121 refers to time 60.5 seconds. To retain only the values in the new hemodynamic vector that refer to times up to (not including) 60s, we can just drop the last len(hrf_signal) - 1 values:
hemodynamic_for_60s = hemodynamic_prediction[:len(neural_vector)]
plt.plot(times, neural_vector, label='neural vector')
plt.plot(times, hemodynamic_for_60s, label='hemodynamic prediction')
plt.legend()
Let's say I have a 4D image:
import nibabel as nib
img = nib.load('ds114_sub009_t2r1.nii', mmap=False)
data = img.get_data()
data.shape
Each element slicing over the last dimension, is a 3D volume:
vol0 = data[..., 0]
vol1 = data[..., 1]
vol1.shape
The term voxel refers to a 3 dimensional pixel - a position in three dimensional space.
These 3D volumes contain voxel values, for all of the 64, 64, 30 voxels, for a single time point.
For example, let's say I was interested at the voxel at (i, j, k) coordinates (10, 23, 9):
vox_value = vol0[10, 23, 9]
vox_value
This value is the value for this voxel (voxel at 10, 23, 9) for this time point (first time point). This would be the value for the same voxel at the second time point (volume index 1):
vol1[10, 23, 9]
Of course I could get the value for this voxel, at every time point, by doing a loop:
vox_values = []
for i in range(data.shape[-1]):
vox_values.append(data[10, 23, 9, i])
len(vox_values)
vox_values[:15]
I can get the same effect by slicing over the last dimension:
sliced_vox_values = data[10, 23, 9, :]
sliced_vox_values[:15]
np.all(vox_values == sliced_vox_values)