Download PDF

We are working to support a site-wide PDF but it is not yet available. You can download PDFs for individual lectures through the download badge on each lecture page.

Update: New build system
QuantEcon is migrating to a new build system - please report any errors to contact@quantecon.org

How to read this lecture...

Code should execute sequentially if run in a Jupyter notebook

NumPy

“Let’s be clear: the work of science has nothing whatever to do with consensus. Consensus is the business of politics. Science, on the contrary, requires only one investigator who happens to be right, which means that he or she has results that are verifiable by reference to the real world. In science consensus is irrelevant. What is relevant is reproducible results.” – Michael Crichton

Overview

NumPy is a first-rate library for numerical programming

  • Widely used in academia, finance and industry
  • Mature, fast, stable and under continuous development

In this lecture we introduce NumPy arrays and the fundamental array processing operations provided by NumPy

Introduction to NumPy

The essential problem that NumPy solves is fast array processing

For example, suppose we want to create an array of 1 million random draws from a uniform distribution and compute the mean

If we did this in pure Python it would be orders of magnitude slower than C or Fortran

This is because

  • Loops in Python over Python data types like lists carry significant overhead
  • C and Fortran code contains a lot of type information that can be used for optimization
  • Various optimizations can be carried out during compilation, when the compiler sees the instructions as a whole

However, for a task like the one described above there’s no need to switch back to C or Fortran

Instead we can use NumPy, where the instructions look like this:

In [1]:
import numpy as np

x = np.random.uniform(0, 1, size=1000000)
x.mean()
Out[1]:
0.5002105347549897

The operations of creating the array and computing its mean are both passed out to carefully optimized machine code compiled from C

More generally, NumPy sends operations in batches to optimized C and Fortran code

This is similar in spirit to Matlab, which provides an interface to fast Fortran routines

A Comment on Vectorization

NumPy is great for operations that are naturally vectorized

Vectorized operations are precompiled routines that can be sent in batches, like

  • matrix multiplication and other linear algebra routines
  • generating a vector of random numbers
  • applying a fixed transformation (e.g., sine or cosine) to an entire array

In a later lecture we’ll discuss code that isn’t easy to vectorize and how such routines can also be optimized

NumPy Arrays

The most important thing that NumPy defines is an array data type formally called a numpy.ndarray

NumPy arrays power a large proportion of the scientific Python ecosystem

To create a NumPy array containing only zeros we use np.zeros

In [2]:
a = np.zeros(3)
a
Out[2]:
array([0., 0., 0.])
In [3]:
type(a)
Out[3]:
numpy.ndarray

NumPy arrays are somewhat like native Python lists, except that

  • Data must be homogeneous (all elements of the same type)
  • These types must be one of the data types (dtypes) provided by NumPy

The most important of these dtypes are:

  • float64: 64 bit floating point number
  • int64: 64 bit integer
  • bool: 8 bit True or False

There are also dtypes to represent complex numbers, unsigned integers, etc

On modern machines, the default dtype for arrays is float64

In [4]:
a = np.zeros(3)
type(a[0])
Out[4]:
numpy.float64

If we want to use integers we can specify as follows:

In [5]:
a = np.zeros(3, dtype=int)
type(a[0])
Out[5]:
numpy.int64

Shape and Dimension

Consider the following assignment

In [6]:
z = np.zeros(10)

Here z is a flat array with no dimension — neither row nor column vector

The dimension is recorded in the shape attribute, which is a tuple

In [7]:
z.shape
Out[7]:
(10,)

Here the shape tuple has only one element, which is the length of the array (tuples with one element end with a comma)

To give it dimension, we can change the shape attribute

In [8]:
z.shape = (10, 1)
z
Out[8]:
array([[0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.]])
In [9]:
z = np.zeros(4)
z.shape = (2, 2)
z
Out[9]:
array([[0., 0.],
       [0., 0.]])

In the last case, to make the 2 by 2 array, we could also pass a tuple to the zeros() function, as in z = np.zeros((2, 2))

Creating Arrays

As we’ve seen, the np.zeros function creates an array of zeros

You can probably guess what np.ones creates

Related is np.empty, which creates arrays in memory that can later be populated with data

In [10]:
z = np.empty(3)
z
Out[10]:
array([0., 0., 0.])

The numbers you see here are garbage values

(Python allocates 3 contiguous 64 bit pieces of memory, and the existing contents of those memory slots are interpreted as float64 values)

To set up a grid of evenly spaced numbers use np.linspace

In [11]:
z = np.linspace(2, 4, 5)  # From 2 to 4, with 5 elements

To create an identity matrix use either np.identity or np.eye

In [12]:
z = np.identity(2)
z
Out[12]:
array([[1., 0.],
       [0., 1.]])

In addition, NumPy arrays can be created from Python lists, tuples, etc. using np.array

In [13]:
z = np.array([10, 20])                 # ndarray from Python list
z
Out[13]:
array([10, 20])
In [14]:
type(z)
Out[14]:
numpy.ndarray
In [15]:
z = np.array((10, 20), dtype=float)    # Here 'float' is equivalent to 'np.float64'
z
Out[15]:
array([10., 20.])
In [16]:
z = np.array([[1, 2], [3, 4]])         # 2D array from a list of lists
z
Out[16]:
array([[1, 2],
       [3, 4]])

See also np.asarray, which performs a similar function, but does not make a distinct copy of data already in a NumPy array

In [17]:
na = np.linspace(10, 20, 2)
na is np.asarray(na)   # Does not copy NumPy arrays
Out[17]:
True
In [18]:
na is np.array(na)     # Does make a new copy --- perhaps unnecessarily
Out[18]:
False

To read in the array data from a text file containing numeric data use np.loadtxt or np.genfromtxt—see the documentation for details

Array Indexing

For a flat array, indexing is the same as Python sequences:

In [19]:
z = np.linspace(1, 2, 5)
z
Out[19]:
array([1.  , 1.25, 1.5 , 1.75, 2.  ])
In [20]:
z[0]
Out[20]:
1.0
In [21]:
z[0:2]  # Two elements, starting at element 0
Out[21]:
array([1.  , 1.25])
In [22]:
z[-1]
Out[22]:
2.0

For 2D arrays the index syntax is as follows:

In [23]:
z = np.array([[1, 2], [3, 4]])
z
Out[23]:
array([[1, 2],
       [3, 4]])
In [24]:
z[0, 0]
Out[24]:
1
In [25]:
z[0, 1]
Out[25]:
2

And so on

Note that indices are still zero-based, to maintain compatibility with Python sequences

Columns and rows can be extracted as follows

In [26]:
z[0, :]
Out[26]:
array([1, 2])
In [27]:
z[:, 1]
Out[27]:
array([2, 4])

NumPy arrays of integers can also be used to extract elements

In [28]:
z = np.linspace(2, 4, 5)
z
Out[28]:
array([2. , 2.5, 3. , 3.5, 4. ])
In [29]:
indices = np.array((0, 2, 3))
z[indices]
Out[29]:
array([2. , 3. , 3.5])

Finally, an array of dtype bool can be used to extract elements

In [30]:
z
Out[30]:
array([2. , 2.5, 3. , 3.5, 4. ])
In [31]:
d = np.array([0, 1, 1, 0, 0], dtype=bool)
d
Out[31]:
array([False,  True,  True, False, False])
In [32]:
z[d]
Out[32]:
array([2.5, 3. ])

We’ll see why this is useful below

An aside: all elements of an array can be set equal to one number using slice notation

In [33]:
z = np.empty(3)
z
Out[33]:
array([2. , 3. , 3.5])
In [34]:
z[:] = 42
z
Out[34]:
array([42., 42., 42.])

Array Methods

Arrays have useful methods, all of which are carefully optimized

In [35]:
a = np.array((4, 3, 2, 1))
a
Out[35]:
array([4, 3, 2, 1])
In [36]:
a.sort()              # Sorts a in place
a
Out[36]:
array([1, 2, 3, 4])
In [37]:
a.sum()               # Sum
Out[37]:
10
In [38]:
a.mean()              # Mean
Out[38]:
2.5
In [39]:
a.max()               # Max
Out[39]:
4
In [40]:
a.argmax()            # Returns the index of the maximal element
Out[40]:
3
In [41]:
a.cumsum()            # Cumulative sum of the elements of a
Out[41]:
array([ 1,  3,  6, 10])
In [42]:
a.cumprod()           # Cumulative product of the elements of a
Out[42]:
array([ 1,  2,  6, 24])
In [43]:
a.var()               # Variance
Out[43]:
1.25
In [44]:
a.std()               # Standard deviation
Out[44]:
1.118033988749895
In [45]:
a.shape = (2, 2)
a.T                   # Equivalent to a.transpose()
Out[45]:
array([[1, 3],
       [2, 4]])

Another method worth knowing is searchsorted()

If z is a nondecreasing array, then z.searchsorted(a) returns the index of the first element of z that is >= a

In [46]:
z = np.linspace(2, 4, 5)
z
Out[46]:
array([2. , 2.5, 3. , 3.5, 4. ])
In [47]:
z.searchsorted(2.2)
Out[47]:
1

Many of the methods discussed above have equivalent functions in the NumPy namespace

In [48]:
a = np.array((4, 3, 2, 1))
In [49]:
np.sum(a)
Out[49]:
10
In [50]:
np.mean(a)
Out[50]:
2.5

Operations on Arrays

Arithmetic Operations

The operators +, -, *, / and ** all act elementwise on arrays

In [51]:
a = np.array([1, 2, 3, 4])
b = np.array([5, 6, 7, 8])
a + b
Out[51]:
array([ 6,  8, 10, 12])
In [52]:
a * b
Out[52]:
array([ 5, 12, 21, 32])

We can add a scalar to each element as follows

In [53]:
a + 10
Out[53]:
array([11, 12, 13, 14])

Scalar multiplication is similar

In [54]:
a * 10
Out[54]:
array([10, 20, 30, 40])

The two dimensional arrays follow the same general rules

In [55]:
A = np.ones((2, 2))
B = np.ones((2, 2))
A + B
Out[55]:
array([[2., 2.],
       [2., 2.]])
In [56]:
A + 10
Out[56]:
array([[11., 11.],
       [11., 11.]])
In [57]:
A * B
Out[57]:
array([[1., 1.],
       [1., 1.]])

In particular, A * B is not the matrix product, it is an element-wise product

Matrix Multiplication

With Anaconda’s scientific Python package based around Python 3.5 and above, one can use the @ symbol for matrix multiplication, as follows:

In [58]:
A = np.ones((2, 2))
B = np.ones((2, 2))
A @ B
Out[58]:
array([[2., 2.],
       [2., 2.]])

(For older versions of Python and NumPy you need to use the np.dot function)

We can also use @ to take the inner product of two flat arrays

In [59]:
A = np.array((1, 2))
B = np.array((10, 20))
A @ B
Out[59]:
50

In fact, we can use @ when one element is a Python list or tuple

In [60]:
A = np.array(((1, 2), (3, 4)))
A
Out[60]:
array([[1, 2],
       [3, 4]])
In [61]:
A @ (0, 1)
Out[61]:
array([2, 4])

Since we are postmultiplying, the tuple is treated as a column vector

Mutability and Copying Arrays

NumPy arrays are mutable data types, like Python lists

In other words, their contents can be altered (mutated) in memory after initialization

We already saw examples above

Here’s another example:

In [62]:
a = np.array([42, 44])
a
Out[62]:
array([42, 44])
In [63]:
a[-1] = 0  # Change last element to 0
a
Out[63]:
array([42,  0])

Mutability leads to the following behavior (which can be shocking to MATLAB programmers…)

In [64]:
a = np.random.randn(3)
a
Out[64]:
array([-0.3980519 , -1.78898247,  1.03612545])
In [65]:
b = a
b[0] = 0.0
a
Out[65]:
array([ 0.        , -1.78898247,  1.03612545])

What’s happened is that we have changed a by changing b

The name b is bound to a and becomes just another reference to the array (the Python assignment model is described in more detail later in the course)

Hence, it has equal rights to make changes to that array

This is in fact the most sensible default behavior!

It means that we pass around only pointers to data, rather than making copies

Making copies is expensive in terms of both speed and memory

Making Copies

It is of course possible to make b an independent copy of a when required

This can be done using np.copy

In [66]:
a = np.random.randn(3)
a
Out[66]:
array([-0.15839771,  0.06738944,  1.03912633])
In [67]:
b = np.copy(a)
b
Out[67]:
array([-0.15839771,  0.06738944,  1.03912633])

Now b is an independent copy (called a deep copy)

In [68]:
b[:] = 1
b
Out[68]:
array([1., 1., 1.])
In [69]:
a
Out[69]:
array([-0.15839771,  0.06738944,  1.03912633])

Note that the change to b has not affected a

Additional Functionality

Let’s look at some other useful things we can do with NumPy

Vectorized Functions

NumPy provides versions of the standard functions log, exp, sin, etc. that act element-wise on arrays

In [70]:
z = np.array([1, 2, 3])
np.sin(z)
Out[70]:
array([0.84147098, 0.90929743, 0.14112001])

This eliminates the need for explicit element-by-element loops such as

In [71]:
n = len(z)
y = np.empty(n)
for i in range(n):
    y[i] = np.sin(z[i])

Because they act element-wise on arrays, these functions are called vectorized functions

In NumPy-speak, they are also called ufuncs, which stands for “universal functions”

As we saw above, the usual arithmetic operations (+, *, etc.) also work element-wise, and combining these with the ufuncs gives a very large set of fast element-wise functions

In [72]:
z
Out[72]:
array([1, 2, 3])
In [73]:
(1 / np.sqrt(2 * np.pi)) * np.exp(- 0.5 * z**2)
Out[73]:
array([0.24197072, 0.05399097, 0.00443185])

Not all user defined functions will act element-wise

For example, passing the function f defined below a NumPy array causes a ValueError

In [74]:
def f(x):
    return 1 if x > 0 else 0

The NumPy function np.where provides a vectorized alternative:

In [75]:
x = np.random.randn(4)
x
Out[75]:
array([-1.11617263, -0.65557289, -0.2993439 , -0.08697371])
In [76]:
np.where(x > 0, 1, 0)  # Insert 1 if x > 0 true, otherwise 0
Out[76]:
array([0, 0, 0, 0])

You can also use np.vectorize to vectorize a given function

In [77]:
def f(x): return 1 if x > 0 else 0

f = np.vectorize(f)
f(x)                # Passing the same vector x as in the previous example
Out[77]:
array([0, 0, 0, 0])

However, this approach doesn’t always obtain the same speed as a more carefully crafted vectorized function

Comparisons

As a rule, comparisons on arrays are done element-wise

In [78]:
z = np.array([2, 3])
y = np.array([2, 3])
z == y
Out[78]:
array([ True,  True])
In [79]:
y[0] = 5
z == y
Out[79]:
array([False,  True])
In [80]:
z != y
Out[80]:
array([ True, False])

The situation is similar for >, <, >= and <=

We can also do comparisons against scalars

In [81]:
z = np.linspace(0, 10, 5)
z
Out[81]:
array([ 0. ,  2.5,  5. ,  7.5, 10. ])
In [82]:
z > 3
Out[82]:
array([False, False,  True,  True,  True])

This is particularly useful for conditional extraction

In [83]:
b = z > 3
b
Out[83]:
array([False, False,  True,  True,  True])
In [84]:
z[b]
Out[84]:
array([ 5. ,  7.5, 10. ])

Of course we can—and frequently do—perform this in one step

In [85]:
z[z > 3]
Out[85]:
array([ 5. ,  7.5, 10. ])

Subpackages

NumPy provides some additional functionality related to scientific programming through its subpackages

We’ve already seen how we can generate random variables using np.random

In [86]:
z = np.random.randn(10000)  # Generate standard normals
y = np.random.binomial(10, 0.5, size=1000)    # 1,000 draws from Bin(10, 0.5)
y.mean()
Out[86]:
4.928

Another commonly used subpackage is np.linalg

In [87]:
A = np.array([[1, 2], [3, 4]])

np.linalg.det(A)           # Compute the determinant
Out[87]:
-2.0000000000000004
In [88]:
np.linalg.inv(A)           # Compute the inverse
Out[88]:
array([[-2. ,  1. ],
       [ 1.5, -0.5]])

Much of this functionality is also available in SciPy, a collection of modules that are built on top of NumPy

We’ll cover the SciPy versions in more detail soon

For a comprehensive list of what’s available in NumPy see this documentation

Exercises

Exercise 1

Consider the polynomial expression

$$ p(x) = a_0 + a_1 x + a_2 x^2 + \cdots a_N x^N = \sum_{n=0}^N a_n x^n \tag{1} $$

Earlier, you wrote a simple function p(x, coeff) to evaluate (1) without considering efficiency

Now write a new function that does the same job, but uses NumPy arrays and array operations for its computations, rather than any form of Python loop

(Such functionality is already implemented as np.poly1d, but for the sake of the exercise don’t use this class)

  • Hint: Use np.cumprod()

Exercise 2

Let q be a NumPy array of length n with q.sum() == 1

Suppose that q represents a probability mass function

We wish to generate a discrete random variable $ x $ such that $ \mathbb P\{x = i\} = q_i $

In other words, x takes values in range(len(q)) and x = i with probability q[i]

The standard (inverse transform) algorithm is as follows:

  • Divide the unit interval $ [0, 1] $ into $ n $ subintervals $ I_0, I_1, \ldots, I_{n-1} $ such that the length of $ I_i $ is $ q_i $
  • Draw a uniform random variable $ U $ on $ [0, 1] $ and return the $ i $ such that $ U \in I_i $

The probability of drawing $ i $ is the length of $ I_i $, which is equal to $ q_i $

We can implement the algorithm as follows

In [89]:
from random import uniform

def sample(q):
    a = 0.0
    U = uniform(0, 1)
    for i in range(len(q)):
        if a < U <= a + q[i]:
            return i
        a = a + q[i]

If you can’t see how this works, try thinking through the flow for a simple example, such as q = [0.25, 0.75] It helps to sketch the intervals on paper

Your exercise is to speed it up using NumPy, avoiding explicit loops

  • Hint: Use np.searchsorted and np.cumsum

If you can, implement the functionality as a class called discreteRV, where

  • the data for an instance of the class is the vector of probabilities q
  • the class has a draw() method, which returns one draw according to the algorithm described above

If you can, write the method so that draw(k) returns k draws from q

Exercise 3

Recall our earlier discussion of the empirical cumulative distribution function

Your task is to

  1. Make the __call__ method more efficient using NumPy
  2. Add a method that plots the ECDF over $ [a, b] $, where $ a $ and $ b $ are method parameters

Solutions

In [90]:
import matplotlib.pyplot as plt
%matplotlib inline

Exercise 1

This code does the job

In [91]:
def p(x, coef):
    X = np.empty(len(coef))
    X[0] = 1
    X[1:] = x
    y = np.cumprod(X)   # y = [1, x, x**2,...]
    return coef @ y

Let’s test it

In [92]:
coef = np.ones(3)
print(coef)
print(p(1, coef))
# For comparison
q = np.poly1d(coef)
print(q(1))
[1. 1. 1.]
3.0
3.0

Exercise 2

Here’s our first pass at a solution:

In [93]:
from numpy import cumsum
from numpy.random import uniform

class DiscreteRV:
    """
    Generates an array of draws from a discrete random variable with vector of
    probabilities given by q.
    """

    def __init__(self, q):
        """
        The argument q is a NumPy array, or array like, nonnegative and sums
        to 1
        """
        self.q = q
        self.Q = cumsum(q)

    def draw(self, k=1):
        """
        Returns k draws from q. For each such draw, the value i is returned
        with probability q[i].
        """
        return self.Q.searchsorted(uniform(0, 1, size=k))

The logic is not obvious, but if you take your time and read it slowly, you will understand

There is a problem here, however

Suppose that q is altered after an instance of discreteRV is created, for example by

In [94]:
q = (0.1, 0.9)
d = DiscreteRV(q)
d.q = (0.5, 0.5)

The problem is that Q does not change accordingly, and Q is the data used in the draw method

To deal with this, one option is to compute Q every time the draw method is called

But this is inefficient relative to computing Q once off

A better option is to use descriptors

A solution from the quantecon library using descriptors that behaves as we desire can be found here

Exercise 3

An example solution is given below

In essence we’ve just taken this code from QuantEcon and added in a plot method

In [95]:
"""
Modifies ecdf.py from QuantEcon to add in a plot method

"""

class ECDF:
    """
    One-dimensional empirical distribution function given a vector of
    observations.

    Parameters
    ----------
    observations : array_like
        An array of observations

    Attributes
    ----------
    observations : array_like
        An array of observations

    """

    def __init__(self, observations):
        self.observations = np.asarray(observations)

    def __call__(self, x):
        """
        Evaluates the ecdf at x

        Parameters
        ----------
        x : scalar(float)
            The x at which the ecdf is evaluated

        Returns
        -------
        scalar(float)
            Fraction of the sample less than x

        """
        return np.mean(self.observations <= x)

    def plot(self, a=None, b=None):
        """
        Plot the ecdf on the interval [a, b].

        Parameters
        ----------
        a : scalar(float), optional(default=None)
            Lower end point of the plot interval
        b : scalar(float), optional(default=None)
            Upper end point of the plot interval

        """

        # === choose reasonable interval if [a, b] not specified === #
        if a is None:
            a = self.observations.min() - self.observations.std()
        if b is None:
            b = self.observations.max() + self.observations.std()

        # === generate plot === #
        x_vals = np.linspace(a, b, num=100)
        f = np.vectorize(self.__call__)
        plt.plot(x_vals, f(x_vals))
        plt.show()

Here’s an example of usage

In [96]:
X = np.random.randn(1000)
F = ECDF(X)
F.plot()