# Numbers

The `numpy` array is the foundation of essentially all numerical computing in Python, so it is important to understand the array and how to use it well.

## Learning objectives

1. Attributes of an array
2. How to create vectors, matrices, tensors
3. How to index and slice arrays
4. Generating random arrays and sampling
5. Universal functions, vectorization and matrix multiplication
6. Array axes and marginal calculations
7. Broadcasting
8. Masking
9. Combining and splitting arrays
10. Vectorizing loops


```python
import numpy as np
```

## The `ndarray`: Vectors, matrices and tenosrs

dtype, shape, strides

### Vector


```python
x = np.array([1,2,3])
x
```




    array([1, 2, 3])




```python
type(x)
```




    numpy.ndarray




```python
x.dtype
```




    dtype('int64')




```python
x.shape
```




    (3,)




```python
x.strides
```




    (8,)



### Matrix


```python
x = np.array([[1,2,3], [4,5,6]], dtype=np.int32)
x
```




    array([[1, 2, 3],
           [4, 5, 6]], dtype=int32)




```python
x.dtype
```




    dtype('int32')




```python
x.shape
```




    (2, 3)




```python
x.strides
```




    (12, 4)



### Tensor


```python
x = np.arange(24).reshape((2,3,4))
```


```python
x
```




    array([[[ 0,  1,  2,  3],
            [ 4,  5,  6,  7],
            [ 8,  9, 10, 11]],
    
           [[12, 13, 14, 15],
            [16, 17, 18, 19],
            [20, 21, 22, 23]]])



## Creating `ndarray`s

### From a file


```python
%%file numbers.txt
a,b,c # can also skip headers
1,2,3
4,5,6
```

    Writing numbers.txt



```python
np.loadtxt('numbers.txt', dtype='int', delimiter=',',
           skiprows=1, comments='#')
```




    array([[1, 2, 3],
           [4, 5, 6]])



### From Python lists or tuples


```python
np.array([
    [1,2,3],
    [4,5,6]
])
```




    array([[1, 2, 3],
           [4, 5, 6]])



### From ranges

arange, linspace, logspace


```python
np.arange(1, 7).reshape((2,3))
```




    array([[1, 2, 3],
           [4, 5, 6]])




```python
np.linspace(1, 10, 4)
```




    array([ 1.,  4.,  7., 10.])




```python
np.logspace(0, 4, 5, dtype='int')
```




    array([    1,    10,   100,  1000, 10000])



### From a function

`fromfunciton`


```python
np.fromfunction(lambda i, j: i*3 + j + 1, (2,3))
```




    array([[1., 2., 3.],
           [4., 5., 6.]])




```python
np.fromfunction(lambda i, j: (i-2)**2 + (j-2)**2, (5,5), dtype='int')
```




    array([[8, 5, 4, 5, 8],
           [5, 2, 1, 2, 5],
           [4, 1, 0, 1, 4],
           [5, 2, 1, 2, 5],
           [8, 5, 4, 5, 8]])



#### How to visualize `fromfunction` 


```python
j = np.repeat([np.arange(5)], 5, axis=0)
i = j.T
```


```python
i
```




    array([[0, 0, 0, 0, 0],
           [1, 1, 1, 1, 1],
           [2, 2, 2, 2, 2],
           [3, 3, 3, 3, 3],
           [4, 4, 4, 4, 4]])




```python
j
```




    array([[0, 1, 2, 3, 4],
           [0, 1, 2, 3, 4],
           [0, 1, 2, 3, 4],
           [0, 1, 2, 3, 4],
           [0, 1, 2, 3, 4]])




```python
(i-2)**2 + (j-2)**2
```




    array([[8, 5, 4, 5, 8],
           [5, 2, 1, 2, 5],
           [4, 1, 0, 1, 4],
           [5, 2, 1, 2, 5],
           [8, 5, 4, 5, 8]])



#### Using element-wise functions in `fromfunction`


```python
np.fromfunction(lambda i, j: np.where(i==j,0, -1), (5,5))
```




    array([[ 0, -1, -1, -1, -1],
           [-1,  0, -1, -1, -1],
           [-1, -1,  0, -1, -1],
           [-1, -1, -1,  0, -1],
           [-1, -1, -1, -1,  0]])




```python
np.fromfunction(lambda i, j: np.where(i<j, 1, np.where(i==j,0, -1)), (5,5))
```




    array([[ 0,  1,  1,  1,  1],
           [-1,  0,  1,  1,  1],
           [-1, -1,  0,  1,  1],
           [-1, -1, -1,  0,  1],
           [-1, -1, -1, -1,  0]])




```python
np.fromfunction(lambda i, j: np.minimum(i,j), (5,5), dtype='int')
```




    array([[0, 0, 0, 0, 0],
           [0, 1, 1, 1, 1],
           [0, 1, 2, 2, 2],
           [0, 1, 2, 3, 3],
           [0, 1, 2, 3, 4]])




```python
np.fromfunction(lambda i, j: np.maximum(i,j), (5,5), dtype='int')
```




    array([[0, 1, 2, 3, 4],
           [1, 1, 2, 3, 4],
           [2, 2, 2, 3, 4],
           [3, 3, 3, 3, 4],
           [4, 4, 4, 4, 4]])



### From special constructors

zeros, ones, eye, diag


```python
np.zeros((2,3))
```




    array([[0., 0., 0.],
           [0., 0., 0.]])




```python
np.ones((2,3))
```




    array([[1., 1., 1.],
           [1., 1., 1.]])




```python
np.eye(3)
```




    array([[1., 0., 0.],
           [0., 1., 0.],
           [0., 0., 1.]])




```python
np.eye(3, 4)
```




    array([[1., 0., 0., 0.],
           [0., 1., 0., 0.],
           [0., 0., 1., 0.]])




```python
np.eye(4, k=-1)
```




    array([[0., 0., 0., 0.],
           [1., 0., 0., 0.],
           [0., 1., 0., 0.],
           [0., 0., 1., 0.]])




```python
np.diag([1,2,3,4])
```




    array([[1, 0, 0, 0],
           [0, 2, 0, 0],
           [0, 0, 3, 0],
           [0, 0, 0, 4]])




```python
np.diag([1,2,3,4], k=1)
```




    array([[0, 1, 0, 0, 0],
           [0, 0, 2, 0, 0],
           [0, 0, 0, 3, 0],
           [0, 0, 0, 0, 4],
           [0, 0, 0, 0, 0]])



### From random variables

#### Convenience functions

rand, randn


```python
np.random.rand(2,3)
```




    array([[0.96634527, 0.19003724, 0.64836211],
           [0.16154332, 0.37613518, 0.56387698]])




```python
np.random.randn(2,3)
```




    array([[-0.68305569,  1.02565577,  1.58315507],
           [ 0.40577645,  0.35897126, -1.26921157]])



#### Distributions

uniform, normal, randint, poisson, multinomial, multivariate_ normal


```python
np.random.uniform(0, 1, (2,3))
```




    array([[0.57022181, 0.27519768, 0.61746728],
           [0.46176807, 0.07558838, 0.49391296]])




```python
np.random.normal(0, 1, (2,3))
```




    array([[-0.86289848,  0.22753939, -0.84960146],
           [-1.28280192,  0.538776  , -0.6454793 ]])




```python
np.random.randint(0, 10, (4,5))
```




    array([[5, 1, 6, 1, 1],
           [3, 2, 4, 2, 3],
           [1, 4, 1, 8, 7],
           [6, 2, 6, 2, 2]])




```python
np.random.poisson(10, (4,5))
```




    array([[ 5, 12, 11, 10, 15],
           [15,  7,  9, 11, 13],
           [13,  6, 11,  7, 15],
           [12,  9, 17, 10,  9]])




```python
np.random.multinomial(n=5, pvals=np.ones(5)/5, size=8)
```




    array([[0, 1, 2, 1, 1],
           [1, 1, 0, 3, 0],
           [1, 1, 1, 1, 1],
           [1, 2, 0, 0, 2],
           [2, 0, 0, 2, 1],
           [3, 0, 1, 1, 0],
           [1, 0, 1, 2, 1],
           [2, 1, 2, 0, 0]])




```python
np.random.multivariate_normal(mean=[10,20,30], cov=np.eye(3), size=4)
```




    array([[10.17827907, 21.16944403, 29.31747505],
           [ 9.93282935, 18.76152496, 31.14436234],
           [ 9.45692559, 19.1736249 , 32.1318146 ],
           [ 9.63006312, 20.2535023 , 29.52264229]])



### Sampling using `choice`

Works much like the R `sample` function.


```python
x = np.random.permutation(list('ABCDEF'))
```


```python
x
```




    array(['E', 'F', 'A', 'B', 'D', 'C'], dtype='<U1')




```python
np.random.choice(x, 3)
```




    array(['E', 'F', 'E'], dtype='<U1')




```python
np.random.choice(x, 10)
```




    array(['A', 'B', 'C', 'F', 'D', 'E', 'E', 'D', 'F', 'D'], dtype='<U1')




```python
try:
    np.random.choice(x, 10, replace=False)
except ValueError as e:
    print(e)
```

    Cannot take a larger sample than population when 'replace=False'


## Indexing 


```python
x = np.arange(20).reshape((4,5))
x
```




    array([[ 0,  1,  2,  3,  4],
           [ 5,  6,  7,  8,  9],
           [10, 11, 12, 13, 14],
           [15, 16, 17, 18, 19]])



### Extracing a scalar


```python
x[1,1]
```




    6



### Extracting a vector


```python
x[1]
```




    array([5, 6, 7, 8, 9])



### Using slices


```python
x[1,:]
```




    array([5, 6, 7, 8, 9])




```python
x[:,1]
```




    array([ 1,  6, 11, 16])




```python
x[1:3,1:3]
```




    array([[ 6,  7],
           [11, 12]])



### Using slices with strides


```python
x[::2,::2]
```




    array([[ 0,  2,  4],
           [10, 12, 14]])



### Extrcting blocks with arbitrary row and column lists (fancy indexing)

`np.ix_`


```python
x[:, [0,3]]
```




    array([[ 0,  3],
           [ 5,  8],
           [10, 13],
           [15, 18]])



Warning: Fancy indexing can only be used for 1 dimension at a time.

In the example below, `numpy` treats the arguments as *paired* coordinates, and returns the values at (0,0) and (2,3).


```python
x[[0,2],[0,3]]
```




    array([ 0, 13])



Use the helper `np.ix_` to extract arbitrary blocks.


```python
x[np.ix_([0,2], [0,3])]
```




    array([[ 0,  3],
           [10, 13]])



### A slice is a view, not a copy

**Warning**

```python
b = a[:]
```

makes a copy if `a` is a list but not if `a` is a numpy array


```python
a1 = list(range(3))
a2 = np.arange(3)
```


```python
b = a1[:]
b[1] = 9
a1
```




    [0, 1, 2]




```python
b = a2[:]
b[1] = 9
a2
```




    array([0, 9, 2])




```python
x
```




    array([[ 0,  1,  2,  3,  4],
           [ 5,  6,  7,  8,  9],
           [10, 11, 12, 13, 14],
           [15, 16, 17, 18, 19]])




```python
y = x[1:-1, 1:-1]
y
```




    array([[ 6,  7,  8],
           [11, 12, 13]])




```python
y *= 10
```


```python
y
```




    array([[ 60,  70,  80],
           [110, 120, 130]])




```python
x
```




    array([[  0,   1,   2,   3,   4],
           [  5,  60,  70,  80,   9],
           [ 10, 110, 120, 130,  14],
           [ 15,  16,  17,  18,  19]])



Use the copy method to convert a view to a copy


```python
z = x[1:-1, 1:-1].copy()
```


```python
z
```




    array([[ 60,  70,  80],
           [110, 120, 130]])




```python
z[:] = 0
```


```python
z
```




    array([[0, 0, 0],
           [0, 0, 0]])




```python
x
```




    array([[  0,   1,   2,   3,   4],
           [  5,  60,  70,  80,   9],
           [ 10, 110, 120, 130,  14],
           [ 15,  16,  17,  18,  19]])



### Boolean indexing


```python
x[x % 2 == 0]
```




    array([  0,   2,   4,  60,  70,  80,  10, 110, 120, 130,  14,  16,  18])




```python
x [x > 3]
```




    array([  4,   5,  60,  70,  80,   9,  10, 110, 120, 130,  14,  15,  16,
            17,  18,  19])



### Functions that return indexes


```python
idx = np.nonzero(x)
idx
```




    (array([0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]),
     array([1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4]))




```python
x[idx]
```




    array([  1,   2,   3,   4,   5,  60,  70,  80,   9,  10, 110, 120, 130,
            14,  15,  16,  17,  18,  19])




```python
idx = np.where(x > 3)
idx
```




    (array([0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]),
     array([4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4]))




```python
x[idx]
```




    array([  4,   5,  60,  70,  80,   9,  10, 110, 120, 130,  14,  15,  16,
            17,  18,  19])



## Universal functions


```python
x
```




    array([[  0,   1,   2,   3,   4],
           [  5,  60,  70,  80,   9],
           [ 10, 110, 120, 130,  14],
           [ 15,  16,  17,  18,  19]])



Operations


```python
x + x
```




    array([[  0,   2,   4,   6,   8],
           [ 10, 120, 140, 160,  18],
           [ 20, 220, 240, 260,  28],
           [ 30,  32,  34,  36,  38]])



Element-wise functions


```python
np.log1p(x)
```




    array([[0.        , 0.69314718, 1.09861229, 1.38629436, 1.60943791],
           [1.79175947, 4.11087386, 4.26267988, 4.39444915, 2.30258509],
           [2.39789527, 4.7095302 , 4.79579055, 4.87519732, 2.7080502 ],
           [2.77258872, 2.83321334, 2.89037176, 2.94443898, 2.99573227]])




```python
x.clip(10, 100)
```




    array([[ 10,  10,  10,  10,  10],
           [ 10,  60,  70,  80,  10],
           [ 10, 100, 100, 100,  14],
           [ 15,  16,  17,  18,  19]])



Scans


```python
np.cumsum(x, axis=1)
```




    array([[  0,   1,   3,   6,  10],
           [  5,  65, 135, 215, 224],
           [ 10, 120, 240, 370, 384],
           [ 15,  31,  48,  66,  85]])



Reductions


```python
np.sum(x)
```




    703




```python
x.prod()
```




    0



## Margins and the `axis` argument


```python
x
```




    array([[  0,   1,   2,   3,   4],
           [  5,  60,  70,  80,   9],
           [ 10, 110, 120, 130,  14],
           [ 15,  16,  17,  18,  19]])



The 0th axis has 4 items, the 1st axis has 5 items.


```python
x.shape
```




    (4, 5)




```python
x.mean()
```




    35.15



### Marginalizing out the 0th axis = column summaries


```python
x.mean(axis=0)
```




    array([ 7.5 , 46.75, 52.25, 57.75, 11.5 ])



### Marginalizing out the 1st axis = row summaries


```python
x.mean(axis=1)
```




    array([ 2. , 44.8, 76.8, 17. ])



Note marginalizing out the last axis is a common default.


```python
x.mean(axis=-1)
```




    array([ 2. , 44.8, 76.8, 17. ])



### Marginalization works for higher dimensions in the same way


```python
x = np.random.random((2,3,4))
x
```




    array([[[0.53438657, 0.9188977 , 0.63964362, 0.98093193],
            [0.85858239, 0.15679935, 0.20372358, 0.60655892],
            [0.65779681, 0.84197462, 0.85369405, 0.96133565]],
    
           [[0.95375215, 0.95867659, 0.74715273, 0.4698344 ],
            [0.98018627, 0.07420097, 0.5479176 , 0.14725169],
            [0.05951145, 0.86547273, 0.38129592, 0.96394799]]])




```python
x.shape
```




    (2, 3, 4)




```python
x.mean(axis=0).shape
```




    (3, 4)




```python
x.mean(axis=1).shape
```




    (2, 4)




```python
x.mean(axis=2).shape
```




    (2, 3)




```python
x.mean(axis=(0,1)).shape
```




    (4,)




```python
x.mean(axis=(0,2)).shape
```




    (3,)




```python
x.mean(axis=(1,2)).shape
```




    (2,)



## Broadcasting

Broadcasting is what happens when `numpy` tries to perform binary operations on two arrays with different shapes. In general, shapes are *promoted* to make the arrays compatible using the following rule

- For each axis from highest to lowest
    - If both dimensions are the same, do nothing
    - If one of the dimensions is 1 or None and the other is $k$, promote to $k$
    - Otherwise print error message


```python
x = np.zeros((3,2))
x.shape
```




    (3, 2)




```python
x
```




    array([[0., 0.],
           [0., 0.],
           [0., 0.]])



Shapes are compatible


```python
y = np.ones(2)
y.shape
```




    (2,)




```python
x + y
```




    array([[1., 1.],
           [1., 1.],
           [1., 1.]])



Shapes are compatible


```python
y = np.ones((1,2))
y.shape
```




    (1, 2)




```python
x + y
```




    array([[1., 1.],
           [1., 1.],
           [1., 1.]])



Shapes are incompatible but can be made compaible by adding empty dimension


```python
y = np.ones(3)
y.shape
```




    (3,)




```python
try:
    x + y
except ValueError as e:
    print(e)
```

    operands could not be broadcast together with shapes (3,2) (3,) 



```python
y[:, None].shape
```




    (3, 1)




```python
x + y[:, None]
```




    array([[1., 1.],
           [1., 1.],
           [1., 1.]])



Shapes are incompatible


```python
y = np.ones((2,2))
y.shape
```




    (2, 2)




```python
try:
    x + y
except ValueError as e:
    print(e)
```

    operands could not be broadcast together with shapes (3,2) (2,2) 


### More examples of broadcasting


```python
x1 = np.arange(12)
```


```python
x1
```




    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])




```python
x1 * 10
```




    array([  0,  10,  20,  30,  40,  50,  60,  70,  80,  90, 100, 110])




```python
x2 = np.random.randint(0,10,(3,4))
```


```python
x2
```




    array([[1, 1, 0, 6],
           [6, 9, 3, 7],
           [5, 4, 7, 5]])




```python
x2 * 10
```




    array([[10, 10,  0, 60],
           [60, 90, 30, 70],
           [50, 40, 70, 50]])




```python
x2.shape
```




    (3, 4)



### Column-wise broadcasting


```python
mu = np.mean(x2, axis=0)
mu.shape
```




    (4,)




```python
x2 - mu
```


```python
(x2 - mu).mean(axis=0)
```

### Row wise broadcasting


```python
mu = np.mean(x2, axis=1)
mu.shape
```


```python
try:
    x2 - mu
except ValueError as e:
    print(e)
```

### We can add a "dummy" axis using None or `np.newaxis`


```python
mu[:, None].shape
```


```python
x2 - mu[:, None]
```


```python
x2 - mu[:, np.newaxis]
```


```python
np.mean(x2 - mu[:, None], axis=1)
```

#### Reshaping works too


```python
x2 - mu.reshape((-1,1))
```

#### Exercise in broadcasting

Creating a 12 by 12 multiplication table


```python
x = np.arange(1, 13)
x[:,None] * x[None,:]
```

Scaling to have zero mean and unit standard devation for each feature.


```python
x = np.random.normal(10, 5,(3,4))
x
```

Scaling column-wise


```python
(x - x.mean(axis=0))/x.std(axis=0)
```

Scaling row-wise


```python
(x - x.mean(axis=1)[:, None])/x.std(axis=1)[:, None]
```

## Masking

- [Ref](https://docs.scipy.org/doc/numpy/reference/maskedarray.generic.html)


```python
import numpy.ma as ma
```


```python
x = np.arange(20).reshape(4,5)
```


```python
x
```


```python
mask = x % 2 == 0
mask
```

- Note that values that are True in the mask are not used in the array
- Values that are False are *not* masked and so remain
- So the above mask keeps only the *odd* numbers in the array `x`


```python
m = ma.masked_array(x, mask)
```


```python
m
```


```python
m.data
```


```python
m.mask
```


```python
m.sum(axis=0).data
```


```python
m.sum(axis=1).data
```

### Often used with missing value sentinels


```python
import warnings

with warnings.catch_warnings():
    warnings.simplefilter('ignore', RuntimeWarning)
    x1 = x / mask
```


```python
x1
```


```python
x1.sum()
```


```python
x2 = ma.masked_invalid(x1)
x2
```


```python
x2.data
```


```python
x2.mask
```


```python
x2.sum()
```


```python
x2.filled(0)
```

## Combining `ndarray`s


```python
x1 = np.zeros((3,4))
x2 = np.ones((3,5))
x3 = np.eye(4)
```


```python
x1
```


```python
x2
```


```python
x3
```

### Binding rows when number of columns is the same


```python
np.r_[x1, x3]
```

### Binding columns when number of rows is the same


```python
np.c_[x1, x2]
```

### You can combine more than 2 at a time


```python
np.c_[x1, x2, x1]
```

### Stacking


```python
np.vstack([x1, x3])
```


```python
np.hstack([x1, x2])
```


```python
np.dstack([x2, 2*x2, 3*x2])
```

### Generic stack with axis argument


```python
np.stack([x2, 2*x2, 3*x2], axis=0)
```


```python
np.stack([x2, 2*x2, 3*x2], axis=1)
```


```python
np.stack([x2, 2*x2, 3*x2], axis=2)
```

### Repetition and tiling

#### For a vector


```python
x = np.array([1,2,3])
```


```python
np.repeat(x, 3)
```


```python
np.tile(x, 3)
```

#### For a matrix


```python
x = np.arange(6).reshape((2,3))
x
```


```python
np.repeat(x, 3)
```


```python
np.repeat(x, 3, axis=0)
```


```python
np.repeat(x, 3, axis=1)
```


```python
np.tile(x, (3,2))
```

## Splitting `ndarray`s


```python
x = np.arange(32).reshape((4,8))
```


```python
x
```


```python
np.split(x, 4)
```


```python
np.split(x, 4, axis=1)
```

## Saving and loading arrays


```python
x = np.arange(16).reshape(4,4)
y = np.arange(20).reshape(-1,4)
```


```python
np.save('x.npy', x)
```


```python
np.load('x.npy')
```


```python
np.savez('xy.npz', x=x, y=y)
```


```python
arr = np.load('xy.npz')
```


```python
arr['x']
```


```python
arr['y']
```


```python
import os
```


```python
os.remove('x.npy')
```


```python
os.remove('xy.npz')
```

## Vectorization

### Example 1

The operators and functions (ufuncs) in Python are vectorized, and will work element-wise over all entries in an `ndarray`.


```python
xs = np.zeros(10, dtype='int')
for i in range(10):
    xs[i] = i**2
xs
```


```python
xs = np.arange(10)**2
xs
```

Using ufuncs


```python
np.sqrt(xs)
```


```python
np.log1p(xs)
```

### Example 2

Scalar product.


```python
n = 10

xs = np.random.rand(n)
ys = np.random.rand(n)

s = 0
for i in range(n):
    s += xs[i] * ys[i]
s
```


```python
np.dot(xs, ys)
```


```python
xs @ ys
```

### Example 3

\begin{align}
y_0 &= \alpha + \beta_1 x_1 + \beta_2 x_2 \\
y_1 &= \alpha + \beta_1 x_1 + \beta_2 x_2 \\
y_2 &= \alpha + \beta_1 x_1 + \beta_2 x_2 \\
\end{align}





```python
m = 3
n = 2

alpha = np.random.rand(1)
betas = np.random.rand(n,1)
xs = np.random.rand(m,n)
```


```python
alpha
```


```python
betas
```


```python
xs
```

### Using loops


```python
ys = np.zeros((m,1))
for i in range(m):
    ys[i] = alpha
    for j in range(n):
        ys[i] += betas[j] * xs[i,j]
ys
```

### Removing inner loop


```python
ys = np.zeros((m,1))
for i in range(m):
    ys[i] = alpha + xs[i,:].T @ betas
ys
```

### Removing all loops


```python
ys = alpha + xs @ betas
ys
```

### Alternative approach

The calculaiton with explicit intercepts and coefficients is common in deep learning, where $\alpha$ is called the bias ($b$) and $\beta$ are called the weights ($w$), and each equation is $y[i] = b + w[i]*x[i]$.

It is common in statisiics to use an augmented matrix in which the first column is all ones, so that all that is needed is a single matrix multiplicaiotn.


```python
X = np.c_[np.ones(m), xs]
X
```


```python
alpha
```


```python
betas
```


```python
betas_ = np.concatenate([[alpha], betas])
betas_
```


```python
ys = X @ betas_
ys
```

### Simulating diffusion


```python
%matplotlib inline
import matplotlib.pyplot as plt
```


```python
w = 100
h = 100
x = np.zeros((w+2,h+2), dtype='float')
x[(w//2-1):(w//2+2), (h//2-1):(h//2+2)] = 1

wts = np.ones(5)/5

for i in range(41):
    if i % 10 == 0:    
        plt.figure()
        plt.imshow(x[1:-1, 1:-1], interpolation='nearest')
        
    center = x[1:-1, 1:-1]
    left = x[:-2, 1:-1]
    right = x[2:, 1:-1]
    bottom = x[1:-1, :-2]
    top = x[1:-1, 2:]
    nbrs = np.dstack([center, left, right, bottom, top])
    x = np.sum(wts * nbrs, axis=-1)
```


```python

```

a93c9bc08d4f9c2abf5a9f012c7cb73d7fd0967f
ipynb
Jupyter Notebook
notebooks/S03_Numpy_Annotated.ipynb
Beckkmint/sta-663-2020
acbb3c87f97e0ac5a159d1608af98f8dc81b57cb
MIT
notebooks/S03_Numpy_Annotated.ipynb
Beckkmint/sta-663-2020
acbb3c87f97e0ac5a159d1608af98f8dc81b57cb
MIT
notebooks/S03_Numpy_Annotated.ipynb
Beckkmint/sta-663-2020
acbb3c87f97e0ac5a159d1608af98f8dc81b57cb
MIT
Qwen/Qwen-72B
1. YES
2. YES