# 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:

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

To fit this design to any data, we take the pseudoinverse:

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

Now let’s make some data to fit to:

```
>>> y_0 = np.random.normal(size=12)
>>> beta_0 = piX.dot(y_0)
>>> beta_0
array([-0.372953, 0.295955])
```

We can fit this same design to another set of data:

```
>>> y_1 = np.random.normal(size=12)
>>> beta_1 = piX.dot(y_1)
>>> beta_1
array([ 0.340509, -0.591232])
```

Now the trick. Because of the way that matrix multiplication works, we can fit
to these two sets of data with the same call to `dot`

:

```
>>> Y = np.vstack((y_0, y_1)).T
>>> betas = piX.dot(Y)
>>> betas
array([[-0.372953, 0.340509],
[ 0.295955, -0.591232]])
```

Of course this is true for any number of columns of Y.