<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Scientific Python Quickstart
John Stachurski
ANU
This is a fast-paced, hands-on introduction to scientific computing with Python, contained in a Jupyter notebook.  The main focus will be on introducing Python's four most important scientific libraries
Step1:  Basic NumPy
Perhaps the single most important scientific library for Python is NumPy.  NumPy provides foundational data structures and routines on which many other libraries rely.
Step2:  NumPy defines a basic data type called an array (actually a numpy.ndarray)
Step3:  Note that array data must be homogeneous
The most important data types are
Step4:  When we create an array such as
Step5:  z is a "flat" array with no dimension--- neither row nor column vector
Step6:  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 
For example, let's make it a column vector
Step7:  Creating arrays
Creating empty arrays --- initializing memory
Step8:  These are just garbage numbers --- whatever was in those memory slots
Here's how to make a regular gird sequence
Step9:  Creating an array of ones
Step10:  Arrays can be made from Python lists or tuples
Step11:  Array indexing
Step12:  Array methods
Step13:  Operations on arrays
Standard arithmetic operations on arrays act elementwise
Step14:  For Python $\geq 3.5$ and NumPy $\geq 1.1$ the @ operator also works.
Step15:  I'll continue to use np.dot below for the benefit of those who are using older versions.  But in my opinion the @ operator is much nicer.
Comparisons
Step16:  Matplotlib
Matplotlib is an outstanding plotting and visualization library for Python that interacts nicely with NumPy.  Here are a few quick examples.  We'll see more below when we discuss the SciPy library.
Step17:  Display figures in this browser window rather than having them open up separately
Step18:  Create something to plot
Step19:  Here's a slightly more complex plot
Step20:  SciPy
Let's just cover some simple examples --- references for further reading are below
Statistics and distributions
Let's use scipy.stats to generate some data from the Beta distribution
Step21:  Now let's histogram it and compare it to the original density
Step22:  Other methods
Step23:  Basic linear regression
Step24:  Let's plot this with data and line of best fit
Step25:  Roots and fixed points
Let's choose an arbitrary function to work with
Step26:  Here we see that the algorithm gets it wrong --- newton is fast but not robust
Let's try a hybrid method
Step27:  Note that the hybrid method is robust but still quite fast...
Numerical optimization and integration
Step28:  Linear Algebra
Let's look at some of the most common routines from linear and matrix algebra
Step29:  We'll experiment with matrices
$$ 
  A = \begin{bmatrix} 2 & -1 \ 3 & 0 \end{bmatrix} 
  \quad \text{and} \quad
  b = \begin{bmatrix} 1 \ 1 \end{bmatrix}
$$
Step30:  Let's check that $Ax = b$
Step31:  We can also invert directly
Step32:  Let's compute the eigenvalues and eigenvectors
Step33:  More information
linear algebra
Step34:  Let's start by writing a test data set to the present working directory, so we can read it back in as a dataframe using pandas.  We use an IPython magic to write the data from a cell to a file
Step35:  Let's try that again but this time using the country as the index column
Step36:  Let's drop the year since it's not very informative
Step37:  Let's add a column for GDP per capita
Step38:  Let's sort the whole data frame by GDP per capita
Step39:  Now we'll plot per capital GDP using the dataframe's plot method
Step40:  Exercises
Here are two exercises.  Feel free to consult documentation such as can be found here.  The solutions are below.  The cell with "solution below" is mean to push them below your line of sight and save you from temptation.
Exercise 1
Generate 10000 data points from the exponential distribution with density
$$
f(x; \alpha) = \alpha \exp(-\alpha x)
\qquad
(x > 0, \alpha > 0)
$$
using scipy.stats and taking $\alpha = 0.5$.  Then, after looking up the maximum likelihood estimator of $\alpha$, compute the estimate given your data and check that it is in fact close to $\alpha$.
Step41:  Exercise 2
Using the same data set, implement maximum likelihood again, but this time pretending that you don't know the analytical expression for the maximum likelihood estimator.  Set up the log likelihood function and maximize it numerically using a routine from scipy.optimize.
Step42:  Solutions
Step43:  Solution to Exercise 1
After checking the docs for the exponential distribution we proceed as follows
Step44:  Let's check we've got the right distribution here
Step45:  It's well-known that the MLE of $\alpha$ is $1/\bar x$ where $\bar x$ is the mean of the sample.  Let's check that it is indeed close to $\alpha$.
Step46:  Minimize over a reasonable parameter space
<ASSISTANT_TASK:>
Python Code:
import sys
print(sys.version)
Explanation: Scientific Python Quickstart
John Stachurski
ANU
This is a fast-paced, hands-on introduction to scientific computing with Python, contained in a Jupyter notebook.  The main focus will be on introducing Python's four most important scientific libraries: NumPy, Scipy, Pandas and Matplotlib.
If you don't know how to use this notebook you need to first work through this page.
A slower, more detailed and more systematic treatment of Python for scientific applications can be found at quant-econ.net.  But this notebook is a good place to start for those who like to learn by doing.
Here's some information on the version of Python that I'm using:
End of explanation
import numpy as np     # Import library and give it alias np
print(np.__version__)  #  The version I'm using
Explanation: Basic NumPy
Perhaps the single most important scientific library for Python is NumPy.  NumPy provides foundational data structures and routines on which many other libraries rely.
End of explanation
a = np.zeros(3)     # Create an array of zeros
a                   # Print a
type(a)
Explanation: NumPy defines a basic data type called an array (actually a numpy.ndarray)
End of explanation
a = np.zeros(3)
type(a[1])
Explanation: Note that array data must be homogeneous
The most important data types are:
float64: 64 bit floating point number
float32: 32 bit floating point number
int64: 64 bit integer
int32: 32 bit integer
bool:  8 bit True or False
There are also dtypes to represent complex numbers, unsigned integers, etc
On most machines, the default dtype for arrays is float64
End of explanation
z = np.zeros(10)
Explanation: When we create an array such as
End of explanation
 z.shape
Explanation: z is a "flat" array with no dimension--- neither row nor column vector:
End of explanation
z.shape = (10, 1)
z
z = np.zeros(4)
z.shape = (2, 2)
z
Explanation: 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 
For example, let's make it a column vector
End of explanation
z = np.empty(3)
z
Explanation: Creating arrays
Creating empty arrays --- initializing memory:
End of explanation
z = np.linspace(2, 4, 5)  # From 2 to 4, with 5 elements
z
Explanation: These are just garbage numbers --- whatever was in those memory slots
Here's how to make a regular gird sequence
End of explanation
z = np.ones(3)
z
z = np.identity(2)
z
Explanation: Creating an array of ones
End of explanation
z = np.array([10, 20]) 
z
z = np.array((10, 20), dtype=float) 
z
z = np.array([[1, 2], [3, 4]])         # 2D array from a list of lists
z
Explanation: Arrays can be made from Python lists or tuples
End of explanation
z = np.linspace(1, 2, 5)
z
z[0]  # First element --- Python sequences are zero based, like C, Java, etc.
z[-1]  # Special syntax for last element
z[0:2]  # Meaning: Two elements, starting from element 0
z = np.array([[1, 2], [3, 4]])
z
z[0, 0]
z[0,:]  # First row
z[:,0]  # First column
z = np.linspace(2, 4, 5)
z
d = np.array([0, 1, 1, 0, 0], dtype=bool)
d
z[d]
Explanation: Array indexing
End of explanation
A = np.array((4, 3, 2, 1))
A
A.sort()
A
A.mean()
A.sum()
A.max()
A.cumsum()
A.var()
A.shape = (2, 2)
A
A.T  # Transpose, equivalent to A.transpose()
Explanation: Array methods
End of explanation
a = np.array([1, 2, 3, 4])
b = np.array([5, 6, 7, 8])
a + b
a - b
a + 10
a.shape = 2, 2
b.shape = 2, 2
a
b
a * b # Pointwise multiplication!!
np.dot(a, b) # Matrix multiplication
Explanation: Operations on arrays
Standard arithmetic operations on arrays act elementwise
End of explanation
a @ b
Explanation: For Python $\geq 3.5$ and NumPy $\geq 1.1$ the @ operator also works.
End of explanation
z = np.array([2, 3])
y = np.array([2, 3])
z == y
y[0] = 3
z == y
z = np.linspace(0, 10, 5)
z
z > 3
z[z > 3]  # Conditional extraction
Explanation: I'll continue to use np.dot below for the benefit of those who are using older versions.  But in my opinion the @ operator is much nicer.
Comparisons
End of explanation
import matplotlib.pyplot as plt  # Import main functionality
Explanation: Matplotlib
Matplotlib is an outstanding plotting and visualization library for Python that interacts nicely with NumPy.  Here are a few quick examples.  We'll see more below when we discuss the SciPy library.
End of explanation
%matplotlib inline 
Explanation: Display figures in this browser window rather than having them open up separately:
End of explanation
x = np.linspace(-2, 2, 100)
y = x**2
fig, ax = plt.subplots()  # Create axes and figure window
ax.plot(x, y, 'b-')
Explanation: Create something to plot
End of explanation
y3 = x**3
fig, ax = plt.subplots()  # Create axes and figure window
ax.plot(x, y, 'b-', lw=2, alpha=0.8, label='$x^2$')
ax.plot(x, y3, 'g-', lw=2, alpha=0.8, label='$x^3$')
ax.legend(loc='lower right')
Explanation: Here's a slightly more complex plot
End of explanation
from scipy.stats import beta
q = beta(5, 5)      # Beta(a, b), with a = b = 5
obs = q.rvs(2000)   # 2000 observations
Explanation: SciPy
Let's just cover some simple examples --- references for further reading are below
Statistics and distributions
Let's use scipy.stats to generate some data from the Beta distribution
End of explanation
fig, ax = plt.subplots()
ax.hist(obs, bins=40, normed=True)
grid = np.linspace(0.01, 0.99, 100)
ax.plot(grid, q.pdf(grid), 'k-', linewidth=2)
Explanation: Now let's histogram it and compare it to the original density
End of explanation
type(q)
dir(q)  # Let's see all its methods
q.cdf(0.5)
q.pdf(0.5)
q.mean()
Explanation: Other methods
End of explanation
from scipy.stats import linregress
n = 100
alpha, beta, sigma = 1, 2, 1.5
x = np.random.randn(n)  # n standard normals
y = alpha + beta * x + sigma * np.random.randn(n)
beta_hat, alpha_hat, r_value, p_value, std_err = linregress(x, y)
print("gradient = {}".format(beta_hat))
print("intercept = {}".format(alpha_hat))
Explanation: Basic linear regression:
End of explanation
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(x, y, 'bo', alpha=0.6, label='observations')
xgrid = np.linspace(-3, 3, 2)
ax.plot(xgrid, alpha_hat + beta_hat * xgrid, 'k-', lw=2, alpha=0.8, label='best fit')
ax.grid()
ax.legend(loc='upper left')
Explanation: Let's plot this with data and line of best fit
End of explanation
fig, ax = plt.subplots()
def f(x):
     return np.sin(4 * (x - 0.25)) + x + x**20 - 1
x = np.linspace(0, 1, 100)
ax.plot(x, f(x))
ax.plot(x, 0 * x)
from scipy.optimize import bisect  # Bisection algorithm --- slow but robust
bisect(f, 0, 1)
from scipy.optimize import newton  #  Newton's method --- fast but less robust
newton(f, 0.2)   # Start the search at initial condition x = 0.2
newton(f, 0.7)   # Start the search at x = 0.7 instead 
Explanation: Roots and fixed points
Let's choose an arbitrary function to work with
End of explanation
from scipy.optimize import brentq
brentq(f, 0, 1) # Hybrid method
timeit bisect(f, 0, 1)
timeit newton(f, 0.2)
timeit brentq(f, 0, 1)
Explanation: Here we see that the algorithm gets it wrong --- newton is fast but not robust
Let's try a hybrid method
End of explanation
from scipy.optimize import fminbound
fminbound(lambda x: x**2, -1, 2)  # Search in [-1, 2]
from scipy.integrate import quad
integral, error = quad(lambda x: x**2, 0, 1)
integral
Explanation: Note that the hybrid method is robust but still quite fast...
Numerical optimization and integration
End of explanation
import scipy.linalg as la
Explanation: Linear Algebra
Let's look at some of the most common routines from linear and matrix algebra
End of explanation
A = [[2, -1],
     [3, 0]]
A = np.array(A) # Convert from list to NumPy array
b = np.ones((2, 1))  # Shape is 2 x 1
A
b
x = la.solve(A, b)  # Solve for x in Ax = b
print(x)
Explanation: We'll experiment with matrices
$$ 
  A = \begin{bmatrix} 2 & -1 \ 3 & 0 \end{bmatrix} 
  \quad \text{and} \quad
  b = \begin{bmatrix} 1 \ 1 \end{bmatrix}
$$
End of explanation
np.dot(A, x)
Explanation: Let's check that $Ax = b$
End of explanation
la.inv(A)
np.dot(A, la.inv(A))  # Should be the identity
Explanation: We can also invert directly
End of explanation
eigvals, eigvecs = la.eig(A)
print("eigenvalues = {}".format(eigvals))
print("first eigenvector = {}".format(eigvecs[:, 0]))
Explanation: Let's compute the eigenvalues and eigenvectors
End of explanation
import pandas as pd
Explanation: More information
linear algebra: http://docs.scipy.org/doc/scipy/reference/linalg.html
numerical integration: http://docs.scipy.org/doc/scipy/reference/integrate.html
interpolation: http://docs.scipy.org/doc/scipy/reference/interpolate.html
optimization: http://docs.scipy.org/doc/scipy/reference/optimize.html
distributions and random number generation: http://docs.scipy.org/doc/scipy/reference/stats.html
signal processing: http://docs.scipy.org/doc/scipy/reference/signal.html
Pandas
Pandas is a very popular library for working with data sets.  In pandas, data is held in a dataframe, which is kind of like a spread sheet
End of explanation
%%file test_data.csv
"country","country isocode","year","POP","XRAT","tcgdp","cc","cg"
"Argentina","ARG","2000","37335.653","0.9995","295072.21869","75.716805379","5.5788042896"
"Australia","AUS","2000","19053.186","1.72483","541804.6521","67.759025993","6.7200975332"
"India","IND","2000","1006300.297","44.9416","1728144.3748","64.575551328","14.072205773"
"Israel","ISR","2000","6114.57","4.07733","129253.89423","64.436450847","10.266688415"
"Malawi","MWI","2000","11801.505","59.543808333","5026.2217836","74.707624181","11.658954494"
"South Africa","ZAF","2000","45064.098","6.93983","227242.36949","72.718710427","5.7265463933"
"United States","USA","2000","282171.957","1","9898700","72.347054303","6.0324539789"
"Uruguay","URY","2000","3219.793","12.099591667","25255.961693","78.978740282","5.108067988"
%ls ./*.csv # Check it's there
df = pd.read_csv('./test_data.csv')
df
Explanation: Let's start by writing a test data set to the present working directory, so we can read it back in as a dataframe using pandas.  We use an IPython magic to write the data from a cell to a file:
End of explanation
df = pd.read_csv('./test_data.csv', index_col='country')
df
Explanation: Let's try that again but this time using the country as the index column
End of explanation
df.drop(['year'], axis=1, inplace=True)
df
Explanation: Let's drop the year since it's not very informative
End of explanation
df['GDP percap'] = df['tcgdp'] / df['POP']
df
Explanation: Let's add a column for GDP per capita
End of explanation
df.sort_values(by='GDP percap', inplace=True)
df
Explanation: Let's sort the whole data frame by GDP per capita
End of explanation
df['GDP percap'].plot(kind='bar')
Explanation: Now we'll plot per capital GDP using the dataframe's plot method
End of explanation
# Put your solution here
Explanation: Exercises
Here are two exercises.  Feel free to consult documentation such as can be found here.  The solutions are below.  The cell with "solution below" is mean to push them below your line of sight and save you from temptation.
Exercise 1
Generate 10000 data points from the exponential distribution with density
$$
f(x; \alpha) = \alpha \exp(-\alpha x)
\qquad
(x > 0, \alpha > 0)
$$
using scipy.stats and taking $\alpha = 0.5$.  Then, after looking up the maximum likelihood estimator of $\alpha$, compute the estimate given your data and check that it is in fact close to $\alpha$.
End of explanation
# Put your solution here
Explanation: Exercise 2
Using the same data set, implement maximum likelihood again, but this time pretending that you don't know the analytical expression for the maximum likelihood estimator.  Set up the log likelihood function and maximize it numerically using a routine from scipy.optimize.
End of explanation
# Print some nonsense to partially hide solutions
filler_text = "solution below\n" * 25
print(filler_text)
Explanation: Solutions
End of explanation
from scipy.stats import expon
alpha = 0.5
n = 10000
ep = expon(scale=1.0/alpha)  # scale controls the exponential parameter
x = ep.rvs(n)
Explanation: Solution to Exercise 1
After checking the docs for the exponential distribution we proceed as follows
End of explanation
fig, ax = plt.subplots(figsize=(8, 5))
xmin, xmax = 0.001, 10.0
ax.set_xlim(xmin, xmax)
ax.hist(x, normed=True, bins=40, alpha=0.3)
grid = np.linspace(xmin, xmax, 200)
ax.plot(grid, ep.pdf(grid), 'g-', lw=2, label='true density')
ax.legend()
Explanation: Let's check we've got the right distribution here
End of explanation
alpha_mle = 1.0 / x.mean()
print("max likelihood estimate of alpha is {}".format(alpha_mle))
s = x.sum()
def neg_loglike(a):
    "Minus the log likelihood function for exponential"
    return - n * np.log(a) + a * s
Explanation: It's well-known that the MLE of $\alpha$ is $1/\bar x$ where $\bar x$ is the mean of the sample.  Let's check that it is indeed close to $\alpha$.
End of explanation
from scipy.optimize import fminbound
fminbound(neg_loglike, 0.01, 10.0)
Explanation: Minimize over a reasonable parameter space
End of explanation
<END_TASK>