# Estimation for many voxels at the same time

# Estimation for many voxels at the same time#

We often want to fit the same design to many different voxels.

Let’s make a design with a linear trend and a constant term:

```
import numpy as np
import matplotlib.pyplot as plt
# Print arrays to 2 decimal places
np.set_printoptions(precision=2)
```

```
X = np.ones((12, 2))
X[:, 0] = np.linspace(-1, 1, 12)
plt.imshow(X, cmap='gray')
```

```
<matplotlib.image.AxesImage at 0x7fb26497c0d0>
```

To fit this design to any data, we take the pseudo-inverse:

```
import numpy.linalg as npl
piX = npl.pinv(X)
piX
```

```
array([[-0.21, -0.17, -0.13, -0.1 , -0.06, -0.02, 0.02, 0.06, 0.1 ,
0.13, 0.17, 0.21],
[ 0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08,
0.08, 0.08, 0.08]])
```

Notice the shape of the pseudo-inverse:

```
piX.shape
```

```
(2, 12)
```

Now let’s make some data to fit to. We will draw some samples from the standard normal distribution.

```
# Random number generator.
rng = np.random.default_rng()
# 12 random numbers from normal distribution, mean 0, std 1.
y_0 = rng.normal(size=12)
y_0
```

```
array([-0.27, -1.37, -0.84, 0.53, 0.42, 1.5 , -1.23, -0.14, -0.64,
1.54, -1.94, -1.33])
```

```
beta_0 = piX @ y_0
beta_0
```

```
array([-0.2 , -0.31])
```

We can fit this same design to another set of data, using our already-calculated pseudo-inverse.

```
y_1 = rng.normal(size=12)
y_1
```

```
array([ 0.81, -0.44, 0.34, -0.25, -1. , -1.52, 1.38, -1.84, -0.32,
-0.81, -0.82, 1.28])
```

```
beta_1 = piX @ y_1
beta_1
```

```
array([-0.12, -0.27])
```

And another!:

```
y_2 = rng.normal(size=12)
beta_2 = piX @ y_2
beta_2
```

```
array([-0.42, -0.59])
```

Now the trick. Because of the way that matrix multiplication works, we can fit to these three sets of data with the one single matrix multiply:

```
# Stack the data vectors into columns in a 2D array.
Y = np.vstack((y_0, y_1, y_2)).T
Y
```

```
array([[-0.27, 0.81, 0.87],
[-1.37, -0.44, -0.81],
[-0.84, 0.34, -2.19],
[ 0.53, -0.25, -0.33],
[ 0.42, -1. , -1.35],
[ 1.5 , -1.52, 1.35],
[-1.23, 1.38, -0.12],
[-0.14, -1.84, 1.02],
[-0.64, -0.32, -1.82],
[ 1.54, -0.81, -1.57],
[-1.94, -0.82, -0.67],
[-1.33, 1.28, -1.47]])
```

```
betas = piX @ Y
betas
```

```
array([[-0.2 , -0.12, -0.42],
[-0.31, -0.27, -0.59]])
```

Notice some features of the `betas`

array:

There is one

*row*per*column in the design X*.There is one

*column*per*column in the data Y*.The first

*row*of`betas`

contains the parameters corresponding to the first*column*of the design. The first row therefore contains the*slope*parameters.The second row contains the parameters corresponding to the second

*column*of the design. The second row therefore contains the*intercept*parameters.The first

*column*of`betas`

contains the slope, intercept for the first data vector`y_0`

.

Of course this trick will work for any number of columns of Y.

This trick is important because it allows us to estimate the same design on huge number of data vectors very efficiently. In imaging, this is useful because we can arrange our image data in a one row per volume, one column per voxel array, and use technique here for very fast estimation.