# Numbers


```python
%matplotlib inline
```


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


```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
```

    Overwriting 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))
```




    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.]])




```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]])



Adapted from [Circle Squared](http://igpphome.ucsd.edu/~shearer/COMP233/SciAm_Mandel.pdf)


```python
xs = np.fromfunction(lambda i, j:  0.27**2*(i**2 + j**2) % 1.5, 
                     (400, 400),)
plt.figure(figsize=(8,8))
plt.imshow(xs, interpolation='nearest', cmap=plt.cm.YlOrBr)
plt.xticks([])
plt.yticks([])
pass
```

### From special construcotrs

zeoros, ones, eye


```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.]])



### From random variables

#### Convenience functions

rand, randn


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




    array([[ 0.04737547,  0.82228651,  0.22224204],
           [ 0.96886745,  0.54149529,  0.94853004]])




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




    array([[-0.80516331,  0.9086527 , -1.21134755],
           [-0.08444392,  1.00515737,  0.94762797]])



#### Distributions

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


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




    array([[ 0.12880127,  0.32643222,  0.52082391],
           [ 0.35404701,  0.54964578,  0.456206  ]])




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




    array([[ 0.8195548 ,  0.43940477, -0.24901744],
           [ 2.6248388 ,  0.16150201,  0.12905971]])




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




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




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




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




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




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




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




    array([[  9.34907429,  20.28742704,  31.20954084],
           [ 10.61473257,  19.60538925,  29.09199082],
           [  9.10517086,  19.82157961,  28.57954932],
           [  8.35947889,  22.60010629,  29.05471801]])



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


```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


```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])



## 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.149999999999999



### 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.25091919,  0.97511897,  0.6520031 ,  0.56832053],
            [ 0.97444408,  0.72376855,  0.15067567,  0.10752009],
            [ 0.20627991,  0.0234196 ,  0.70143528,  0.87408988]],
    
           [[ 0.95505841,  0.79364654,  0.76668484,  0.47971008],
            [ 0.34780788,  0.44169709,  0.63645331,  0.96215426],
            [ 0.92889169,  0.52710155,  0.22280321,  0.04262088]]])




```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


```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([[3, 0, 1, 2],
           [2, 4, 0, 7],
           [9, 9, 3, 0]])




```python
x2 * 10
```




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




```python
x2.shape
```




    (3, 4)



### Column-wise broadcasting


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




    (4,)




```python
x2 - mu
```




    array([[-1.66666667, -4.33333333, -0.33333333, -1.        ],
           [-2.66666667, -0.33333333, -1.33333333,  4.        ],
           [ 4.33333333,  4.66666667,  1.66666667, -3.        ]])




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




    array([ -2.96059473e-16,   2.96059473e-16,   7.40148683e-17,
             0.00000000e+00])



### Row wise broadcasting


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




    (3,)




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

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


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


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




    (3, 1)




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




    array([[ 1.5 , -1.5 , -0.5 ,  0.5 ],
           [-1.25,  0.75, -3.25,  3.75],
           [ 3.75,  3.75, -2.25, -5.25]])




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




    array([[ 1.5 , -1.5 , -0.5 ,  0.5 ],
           [-1.25,  0.75, -3.25,  3.75],
           [ 3.75,  3.75, -2.25, -5.25]])




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




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



#### Reshaping works too


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




    array([[ 1.5 , -1.5 , -0.5 ,  0.5 ],
           [-1.25,  0.75, -3.25,  3.75],
           [ 3.75,  3.75, -2.25, -5.25]])



#### Broadcasting examples

Creating a 12 by 12 multiplication table


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




    array([[  1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12],
           [  2,   4,   6,   8,  10,  12,  14,  16,  18,  20,  22,  24],
           [  3,   6,   9,  12,  15,  18,  21,  24,  27,  30,  33,  36],
           [  4,   8,  12,  16,  20,  24,  28,  32,  36,  40,  44,  48],
           [  5,  10,  15,  20,  25,  30,  35,  40,  45,  50,  55,  60],
           [  6,  12,  18,  24,  30,  36,  42,  48,  54,  60,  66,  72],
           [  7,  14,  21,  28,  35,  42,  49,  56,  63,  70,  77,  84],
           [  8,  16,  24,  32,  40,  48,  56,  64,  72,  80,  88,  96],
           [  9,  18,  27,  36,  45,  54,  63,  72,  81,  90,  99, 108],
           [ 10,  20,  30,  40,  50,  60,  70,  80,  90, 100, 110, 120],
           [ 11,  22,  33,  44,  55,  66,  77,  88,  99, 110, 121, 132],
           [ 12,  24,  36,  48,  60,  72,  84,  96, 108, 120, 132, 144]])



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


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




    array([[  7.65396872,   8.74106757,   7.26659687,   9.68096491],
           [  9.27132075,   7.89349589,   9.47813473,   7.04088855],
           [ 13.78187972,   7.39243968,  10.90034325,   4.27477767]])



Standardize column-wise


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




    array([[-0.99566783,  1.31524702, -1.30321654,  1.21511731],
           [-0.37192708, -0.20751913,  0.17598238,  0.01903327],
           [ 1.36759491, -1.10772789,  1.12723416, -1.23415058]])



Standardize row-wise


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




    array([[-0.72037796,  0.42843253, -1.12973985,  1.42168528],
           [ 0.84786879, -0.52591836,  1.0540767 , -1.37602714],
           [ 1.31012381, -0.47301018,  0.50595856, -1.34307219]])



## Combining `ndarray`s


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


```python
x1
```




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




```python
x2
```




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




```python
x3
```




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



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


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




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



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


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




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



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


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




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



### Stacking


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




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




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




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




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




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



### Generic stack with axis argument


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




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




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




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




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




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



### Repetition and tiling

#### For a vector


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


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




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




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




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



#### For a matrix


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




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




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




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




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




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




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




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




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




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



## Splitting `ndarray`s


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


```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],
           [24, 25, 26, 27, 28, 29, 30, 31]])




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




    [array([[0, 1, 2, 3, 4, 5, 6, 7]]),
     array([[ 8,  9, 10, 11, 12, 13, 14, 15]]),
     array([[16, 17, 18, 19, 20, 21, 22, 23]]),
     array([[24, 25, 26, 27, 28, 29, 30, 31]])]




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




    [array([[ 0,  1],
            [ 8,  9],
            [16, 17],
            [24, 25]]), array([[ 2,  3],
            [10, 11],
            [18, 19],
            [26, 27]]), array([[ 4,  5],
            [12, 13],
            [20, 21],
            [28, 29]]), array([[ 6,  7],
            [14, 15],
            [22, 23],
            [30, 31]])]



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




    array([ 0,  1,  4,  9, 16, 25, 36, 49, 64, 81])




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




    array([ 0,  1,  4,  9, 16, 25, 36, 49, 64, 81])



Using ufuncs


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




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




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




    array([ 0.        ,  0.69314718,  1.60943791,  2.30258509,  2.83321334,
            3.25809654,  3.61091791,  3.91202301,  4.17438727,  4.40671925])



### Example 2


```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
```




    2.6868481407430282




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




    2.6868481407430282




```python
xs @ ys
```




    2.6868481407430282



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




    array([ 0.05007867])




```python
betas
```




    array([[ 0.06399532],
           [ 0.37993848]])




```python
xs
```




    array([[ 0.79258489,  0.20951971],
           [ 0.89978267,  0.044122  ],
           [ 0.34222118,  0.06629664]])



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




    array([[ 0.180405  ],
           [ 0.1244242 ],
           [ 0.09716787]])



### Removing inner loop


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




    array([[ 0.180405  ],
           [ 0.1244242 ],
           [ 0.09716787]])



### Removing all loops


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




    array([[ 0.180405  ],
           [ 0.1244242 ],
           [ 0.09716787]])



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




    array([[ 1.        ,  0.79258489,  0.20951971],
           [ 1.        ,  0.89978267,  0.044122  ],
           [ 1.        ,  0.34222118,  0.06629664]])




```python
alpha
```




    array([ 0.05007867])




```python
betas
```




    array([[ 0.06399532],
           [ 0.37993848]])




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




    array([[ 0.05007867],
           [ 0.06399532],
           [ 0.37993848]])




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




    array([[ 0.180405  ],
           [ 0.1244242 ],
           [ 0.09716787]])



### Simulating diffusion


```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

```

f5e622fd2023850b995e77cf30fd79c7be3d8a1d
ipynb
Jupyter Notebook
notebooks/S03_Numpy_Annotated.ipynb
taotangtt/sta-663-2018
67dac909477f81d83ebe61e0753de2328af1be9c
BSD-3-Clause
2018-01-20T20:50:22.000Z
2022-02-27T23:24:21.000Z
notebooks/S03_Numpy_Annotated.ipynb
taotangtt/sta-663-2018
67dac909477f81d83ebe61e0753de2328af1be9c
BSD-3-Clause
2020-02-03T13:43:46.000Z
2020-02-03T13:43:46.000Z
notebooks/S03_Numpy_Annotated.ipynb
taotangtt/sta-663-2018
67dac909477f81d83ebe61e0753de2328af1be9c
BSD-3-Clause
2018-01-12T17:13:14.000Z
2022-03-14T20:22:46.000Z
Qwen/Qwen-72B
1. YES
2. YES