# Modeling a single voxelΒΆ

A long time ago (*Voxel time courses*), we were looking at a single voxel
time course.

Let’s get that same voxel time course back again:

```
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> import nibabel as nib
```

```
>>> img = nib.load('ds114_sub009_t2r1.nii')
>>> data = img.get_data()
>>> data = data[..., 4:]
```

The voxel coordinate (3D coordinate) that we were looking at before was at (42, 32, 19):

```
>>> voxel_time_course = data[42, 32, 19]
>>> plt.plot(voxel_time_course)
[...]
```

Now we are going to use our new convolved regressor to do a simple regression on this voxel time course.

If you don’t have it already, you will need to download
`ds114_sub009_t2r1_conv.txt`

.

```
>>> convolved = np.loadtxt('ds114_sub009_t2r1_conv.txt')
>>> # Knock off first 4 elements to match data
>>> convolved = convolved[4:]
>>> plt.plot(convolved)
[...]
```

First we make our *design matrix*. It has a column for the convolved
regressor, and a column of ones:

```
>>> N = len(convolved)
>>> X = np.ones((N, 2))
>>> X[:, 0] = convolved
>>> plt.imshow(X, interpolation='nearest', cmap='gray', aspect=0.1)
<...>
```

\(\newcommand{\yvec}{\vec{y}}\) \(\newcommand{\xvec}{\vec{x}}\) \(\newcommand{\evec}{\vec{\varepsilon}}\) \(\newcommand{Xmat}{\boldsymbol X} \newcommand{\bvec}{\vec{\beta}}\) \(\newcommand{\bhat}{\hat{\bvec}} \newcommand{\yhat}{\hat{\yvec}}\)

As you will remember from *Introduction to the general linear model*, our model is:

We can get our least squares parameter *estimates* for \(\bvec\) with:

where \(\Xmat^+\) is the *pseudoinverse* of \(\Xmat\). When \(\Xmat\) is
invertible, the pseudoinverse is given by:

Let’s calculate the pseudoinverse for our design:

```
>>> import numpy.linalg as npl
>>> Xp = npl.pinv(X)
>>> Xp.shape
(2, 169)
```

We calculate \(\bhat\):

```
>>> beta_hat = Xp.dot(voxel_time_course)
>>> beta_hat
array([ 31.185514, 2029.367685])
```

We can then calculate \(\yhat\) (also called the *fitted data*):

```
>>> y_hat = X.dot(beta_hat)
>>> e_vec = voxel_time_course - y_hat
>>> print(np.sum(e_vec ** 2))
41405.5727762
>>> plt.plot(voxel_time_course)
[...]
>>> plt.plot(y_hat)
[...]
```