# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Part 2 The Scientific Libraries 
# ## 2.1 [NumPy](https://lectures.quantecon.org/py/numpy.html)
# ### 2.1.1 Overview  
# [The official NumPy documentation](https://docs.scipy.org/doc/numpy/reference/)
# ### 2.1.2 Introduction to NumPy

import numpy as np
x = np.random.uniform(0, 1, size=1000000)              # 生成1000000个0~1的数
x.mean()                                               # 求均值

# ### 2.1.3 NumPy Arrays
# NumPy 中定义的最重要的对象是称为 ndarray 的N维数组类型。  
# ndarray中的每个元素在内存中使用相同大小的块。 ndarray中的每个元素是数据类型对象的对象(称为 dtype)。  
# 从ndarray对象提取的任何元素(通过切片)由一个数组标量类型的 Python 对象表示。 下图显示了ndarray，数据类型对象(dtype)和数组标量类型之间的关系。
# <img  src="https://www.yiibai.com/uploads/images/201704/2404/430090425_51678.jpg"  />
# `numpy.array(object, dtype = None, copy = True, order = None, subok = False, ndmin = 0)`  
# 上面的构造器接受以下参数：
#
#   	参数及描述
# 1.	`object` 任何暴露数组接口方法的对象都会返回一个数组或任何(嵌套)序列。
# 2.	`dtype`  数组的所需数据类型，可选。
# 3.	`copy`  可选，默认为true，对象是否被复制。
# 4.	`order` C(按行)、F(按列)或A(任意，默认)。
# 5.	`subok` 默认情况下，返回的数组被强制为基类数组。 如果为true，则返回子类。
# 6.	`ndimin` 指定返回数组的最小维数。  
#
# To create a NumPy array containing only zeros we use `np.zeros`  
# `zeros(shape, dtype=float, order='C')`

import numpy as np
a = np.zeros(3)                        # 生成3个0
a

type(a)                  

a = np.zeros(3)
type(a[0])                            # a 的元素类型为浮点型

a = np.zeros(3, dtype=int)           # a 的元素类型设置为整数型
type(a[0])

# #### Shape and Dimension

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

z = np.zeros(10)
z

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

z.shape

z.shape = (10, 1)
z

z = np.zeros(4)
z.shape = (2, 2)
z

# #### Creating Arrays
#
# Related is `np.empty`, which creates arrays in memory that can later be populated with data. The numbers you see here are garbage values `float64`  
# `empty(shape, dtype=float, order='C')`

z = np.empty(3)
z

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

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`

z = np.identity(2)
z

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

z = np.array([10, 20])              # ndarray from Python list
z

type(z)

z = np.array((10, 20), dtype=float) # Here 'float' is equivalent to 'np.float64'
z

z = np.array([[1, 2], [3, 4]])      # 2D array from a list of lists
z

na = np.linspace(10, 20, 2) 
na

np.array(na)

np.asarray(na)

na is np.asarray(na)         # Does not copy NumPy arrays

na is np.array(na)           # Does make a new copy --- perhaps unnecessarily

na == np.array(na)

print (id(na))

print (id(np.asarray(na)))

print (id(np.array(na)))

# `==` 比较操作符：用来比较两个对象是否相等，value做为判断因素；  
# `is` 同一性运算符：比较判断两个对象是否相同，id做为判断因素。

# #### Array Indexing

z = np.linspace(1,2,5)
z

z[0]

z[0:2]                # 第一个数到第二个数  [0,1）

z[-1]                # 倒数第一个数

z = np.array([[1,2],[3,4]])
z

# `d[行起始：行终止，列起始：列终止：步长]`

z[0,0]

z[0,1]

z[0:]

z[:,1]              #  第一行到第二行， 第二列

z = np.linspace(2, 4, 5)
z

indices = np.array((0, 2, 3))      # position
z[indices]

z

d = np.array([0, 1, 1, 0, 0], dtype=bool)
d

z[d]      #  position

z = np.empty(3)                 # clear z 
z

z[:] =  42                     # 给z 赋值42
z

# #### Array Methods

A = np.array((4, 3, 2, 1))
A

A.sort() # Sorts A in place 排序
A

A.sum()

A.mean()

A.max()

A.argmax()    # Returns the index of the maximal element  返回位置

A.cumsum()    # Cumulative sum of the elements of A     累计求和

A.cumprod()   # Cumulative product of the elements of A   累积相乘

A.var() # Variance   方差

A.std() # Standard deviation  标准差

A.shape = (2, 2)
A.T                    # Equivalent to A.transpose() 转置

z = np.linspace(2, 4, 5)                # [2,4] with 5 grids
z

z.searchsorted(2.2)    
#  If z is a nondecreasing array, then z.searchsorted(a) returns the index of the first element of z that is >= a
# 大于等于2.2的最小的数的位置

z.searchsorted(2.5)

z.searchsorted(2.6)

a = np.array((4, 3, 2, 1))

np.sum(a)

np.mean(a)

# ### 2.1.4 Operations on Arrays
# #### Algebraic Operations

a = np.array([1, 2, 3, 4])          # 两array相加
b = np.array([5, 6, 7, 8])
a + b

a * b                              # 对应元素相乘

a + 10                             # 元素都加上10

a * 10                             # 元素都乘上10

A = np.ones((2, 2))                  #  依据给定形状和类型(shape[, dtype, order])返回一个新的元素全部为1的数组  type of matrix（2,2）
B = np.ones((2, 2))
A + B

A + 10                              # 对应元素相加

A * B                              # 对应元素相乘

# #### Matrix Multiplication

import numpy as np                # 矩阵乘法
A = np.ones((2, 2))
B = np.ones((2, 2))
A @ B

A = np.array((1, 2))              # 矩阵乘法
B = np.array((10, 20))
A @ B

A = np.array(((1, 2), (3, 4)))    # 赋值
A

A @ (0, 1)                       # 矩阵乘法

# #### Mutability and Copying Arrays

a = np.array([42, 44])
a

a[-1] = 0         # Change last element to 0
a

a = np.random.randn(3)
a

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

b = a
b[0] = 0.0
a

# #### Making Copies

# It is of course possible to make b an independent copy of a when required
#
# With recent versions of NumPy, this is best done using `np.copyto`  
# 产生一个跟 `a` 大小一样，类型一样的对应数组。  
# `np.empty_like(a)`  
# `np.zeros_like(a)`  
# `np.ones_like(a)`

a = np.random.randn(3)
a

b = np.empty_like(a)  # empty array with same shape as a
np.copyto(b, a)       # copy to b from a
b

b[:] = 1
b

a

# ### 2.1.5 Additional Functionality
# #### Vectorized Functions

z = np.array([1, 2, 3])
np.sin(z)

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

n

z

y

# caculating $\frac{1}{\sqrt{2\pi}}{e^{-\frac{1}{2}z^2}}$

(1 / np.sqrt(2 * np.pi)) * np.exp(- 0.5 * z**2)


def f(x):
    return 1 if x > 0 else 0


# The NumPy function `np.where` provides a vectorized alternative:

# +
import numpy as np

x = np.random.randn(4)
x
# -

np.where(x > 0, 1, 0)  # Insert 1 if x > 0 true, otherwise 0


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

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

# #### Comparisons

z = np.array([2, 3])
y = np.array([2, 3])
z == y

y[0] = 5
z == y

z != y

z = np.linspace(0, 10, 5)
z

z > 3

b = z > 3   # 将 z > 3 赋值给 b
b

z[b]

z[z > 3]

# #### Subpackages

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

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()

# Another commonly used subpackage is `np.linalg`
# - `np.linalg.inv()`：矩阵求逆
# - `np.linalg.det()`：矩阵求行列式（标量）

# +
A = np.array([[1, 2], [3, 4]])

np.linalg.det(A)           # Compute the determinant
# -

np.linalg.inv(A)           # Compute the inverse

# ### 2.1.6 Exercises
# #### Exercise 1
# 考虑多项式
#
# $p\left( x \right)={{a}_{0}}+{{a}_{1}}x+{{a}_{2}}{{x}^{2}}+...{{a}_{N}}{{x}^{N}}=\sum\limits_{n=0}^{N}{{{a}_{n}}{{x}^{n}}}$
#
# 使用NumPy arrays和array operations，编写函数p(x，coeff)计算上式

import numpy as np
import matplotlib.pyplot as plt


def p(x, coef):
    X = np.empty(len(coef))
    X[0] = 1
    X[1:] = x
    y = np.cumprod(X)   # y = [1, x, x**2,...]    cumprod为累积相乘
    return coef @ y


# `np.poly1d(c_or_r, r=0, variable=None)`  
# `p=np.poly1d(1,2,3,4)`  
# $1*x**3+2*x**2+3*x+4$  
# `p(2)`

coef = np.ones(3)
print(coef)
print(p(1, coef))
# For comparison
q = np.poly1d(coef)
print(q(1))

coef = [6,8,9,10]
print(coef)
print(p(2, coef))
# For comparison
q = np.poly1d(coef)
print(q(2))

# #### Exercise 2
# 令q是长度为n的NumPy数组，其中q.sum（）== 1
#
# 假设q代表概率密度函数
#
# 生成一个离散的随机变量x，x取值范围range(len(q))且x = i的概率为q[i]
#
# 标准（逆变换）算法如下：
# - 将单位区间[0,1]划分为n个子区间${{I}_{0}},{{I}_{1}},...,{{I}_{n-1}}$，使得${{I}_{i}}$的长度为${{q}_{i}}$
# - 在[0,1]绘制一个均匀的随机变量U,并返回i使$U\in {{I}_{i}}$
#
# 绘制i的概率是${{I}_{i}}$的长度，等于${{q}_{i}}$
#
# 我们可以如下实现算法

# +
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]


# +
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))    #  返回大于（）的最小的数的位置   


# -

# 假设q在`discreteRV`创建实例之后修改了该实例，例如通过

q = (0.1, 0.9)
d = discreteRV(q)
d.q = (0.5, 0.5)

# 问题是Q不会相应地改变，并且Q是`draw`方法中使用的数据
#
# 为了解决这个问题，一种选择是在Q每次调用绘制方法时进行计算
#
# 但相对于Q一次计算而言，这是低效的
#
# 更好的选择是使用描述符
#
# Quantecon库中 使用描述符的解决方案可以在[这里](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/discrete_rv.py)找到

# #### Exercise 3
# 经验累积分布函数
#
# 1.使用NumPy使__call__方法更高效
#
# 2.在[a,b]绘制ECDF，其中a和b是参数

"""
Modifies ecdf.py from QuantEcon to add in a plot method

"""
import matplotlib.pyplot as plt
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)        #  对于小于等于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__)               # 向量化 self.__call__
        plt.plot(x_vals, f(x_vals) )                  # 画图
        plt.show()


X = np.random.randn(1000)             # 正太分布1000个值
F = ECDF(X)
F.plot()

# ## 2.2 [Matplotlib](https://lectures.quantecon.org/py/matplotlib.html)
# ### 2.2.1 Overview
# ### 2.2.2 The MATLAB-Style API

from pylab import *  # Deprecated
x = linspace(0, 10, 200)
y = sin(x)
plot(x, y, 'b-', linewidth=2)              # 'b-'颜色 linewidth 线的宽度
show()

import matplotlib.pyplot as plt
import numpy as np
x = np.linspace(0, 10, 200)
y = np.sin(x)
plt.plot(x, y, 'b-', linewidth=2)
plt.show()

# ### 2.2.3 The Object-Oriented API

# Incidentally, regarding the above lines of code,  
# - the form of the import statement import `matplotlib.pyplot` as plt is standard  
# - Here the call `fig, ax = plt.subplots()` returns a pair, where  
#     `fig` is a `Figure` instance—like a blank canvas  
#    ` ax` is an `AxesSubplot` instance—think of a frame for plotting in
# - The `plot()` function is actually a method of `ax`

import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots()
x = np.linspace(0, 10, 200)
y = np.sin(x)
ax.plot(x, y, 'b-', linewidth=2)
plt.show()

# #### Tweaks

import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots()
x = np.linspace(0, 10, 200)
y = np.sin(x)
ax.plot(x, y, 'r-', linewidth=2, label='sine function', alpha=0.6)    # alpha 透明度
ax.legend()
plt.show()

import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots()
x = np.linspace(0, 10, 200)
y = np.sin(x)
ax.plot(x, y, 'r-', linewidth=2, label='sine function', alpha=0.6)
ax.legend(loc='upper center')
plt.show()

# The `r` in front of the label string tells Python that this is a raw string

import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots()
x = np.linspace(0, 10, 200)
y = np.sin(x)
ax.plot(x, y, 'r-', linewidth=2, label=r'$y=\sin(x)$', alpha=0.6)
ax.legend(loc='upper center')
plt.show()

import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots()
x = np.linspace(0, 10, 200)
y = np.sin(x)
ax.plot(x, y, 'r-', linewidth=2, label=r'$y=\sin(x)$', alpha=0.6)
ax.legend(loc='upper center')
ax.set_yticks([-1, 0, 1])                   # 设置 y 的值
ax.set_title('Test plot')                   # 设置标题
plt.show()

# ### 2.2.4 More Features
# #### Multiple Plots on One Axis

# +
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm
from random import uniform

fig, ax = plt.subplots()
x = np.linspace(-4, 4, 150)
for i in range(3):
    m, s = uniform(-1, 1), uniform(1, 2)        # 设置随机数，作为均值和方差
    y = norm.pdf(x, loc=m, scale=s)             # 正态概率密度函数 loc: mean  scale: deviation
    current_label = r'$\mu = {0:.2f}$'.format(m)     # {0:.2f} 格式参数
    ax.plot(x, y, linewidth=2, alpha=0.6, label=current_label)
ax.legend()
plt.show()
# -

# #### Multiple Subplots

# - 正态分布，`scipy.stats.norm.pdf`
# - 正态直方图，`scipy.stats.norm.rvs `
#
# `matplotlib.pyplot.hist(  
#     x, bins=10, range=None, normed=False,   
#     weights=None, cumulative=False, bottom=None,   
#     histtype=u'bar', align=u'mid', orientation=u'vertical',   
#     rwidth=None, log=False, color=None, label=None, stacked=False,   
#     hold=None, **kwargs) ` 
#     
# `x : (n,) array or sequence of (n,) arrays`
#
# 这个参数是指定每个bin(箱子)分布的数据,对应x轴
#
# `bins : integer or array_like, optional`
#
# 这个参数指定bin(箱子)的个数,也就是总共有几条条状图
#
# `color : color or array_like of colors or None, optional`
#
# 这个指定条状图的颜色

import matplotlib.pyplot as plt
from scipy.stats import norm            # 正态分布
from random import uniform
num_rows, num_cols = 3, 2               # 设置[3,2]图阵
fig, axes = plt.subplots(num_rows, num_cols, figsize=(10, 12))
for i in range(num_rows):               # 在每个图阵内
    for j in range(num_cols):
        m, s = uniform(-1, 1), uniform(1, 2)
        x = norm.rvs(loc=m, scale=s, size=100)
        axes[i, j].hist(x, alpha=0.6, bins=20)    # 画柱状图
        t = r'$\mu = {0:.1f}, \quad \sigma = {1:.1f}$'.format(m, s)   # 保留一位小数
        axes[i, j].set_title(t)
        axes[i, j].set_xticks([-4, 0, 4])
        axes[i, j].set_yticks([])
plt.show()

# In fact the preceding figure was generated by the code above preceded by the following three lines

from matplotlib import rc
rc('font',**{'family':'serif','serif':['Palatino']})
rc('text', usetex=True)

# #### 3D Plots
# [more examples](http://blog.csdn.net/guduruyu/article/details/78050268)  
# [details](http://python.jobbole.com/87471/)

# +
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
import numpy as np
from matplotlib import cm


def f(x, y):
    return np.cos(x**2 + y**2) / (1 + x**2 + y**2)

xgrid = np.linspace(-3, 3, 50)
ygrid = xgrid
x, y = np.meshgrid(xgrid, ygrid)  # 绘制三维图形,首先,需要XY平面的网格数据,这就是meshgrid函数所实现的内容

fig = plt.figure(figsize=(12, 10))   # 图片大小
ax = fig.add_subplot(111, projection='3d')  #加入图像  参数349的意思是：将画布分割成3行4列，图像画在从左到右从上到下的第9块
ax.plot_surface(x,
                y,
                f(x, y),
                rstride=2, cstride=2,  #rstride和cstride是两个方向上的采样，越小越精细，lw是线宽
                cmap=cm.jet,
                alpha=0.7,
                linewidth=0.25)
ax.set_zlim(-0.5, 1.0)                #  z 轴的大小
plt.show()
# -

# #### A Customizing Function

# The custom `subplots` function
# 1. calls the standard `plt.subplots` function internally to generate the `fig, ax` pair,
# 2. makes the desired customizations to `ax`, and
# 3. passes the `fig, ax` pair back to the calling code

# +
import matplotlib.pyplot as plt
import numpy as np


def subplots():
    "Custom subplots with axes throught the origin"
    fig, ax = plt.subplots()

    # Set the axes through the origin
    # 将 坐标轴重新定位
    for spine in ['left', 'bottom']:
        ax.spines[spine].set_position('zero')
    # 将右方和上方的边框设置为不显示
    for spine in ['right', 'top']:
        ax.spines[spine].set_color('none')

    ax.grid()
    return fig, ax


fig, ax = subplots()  # Call the local version, not plt.subplots()
x = np.linspace(-2, 10, 200)
y = np.sin(x)
ax.plot(x, y, 'r-', linewidth=2, label='sine function', alpha=0.6)
ax.legend(loc='lower right')
plt.show()
# -

# ### 2.2.5 Further Reading
# - The [Matplotlib gallery](http://matplotlib.org/gallery.html) provides many examples
# - A nice [Matplotlib tutorial](http://www.scipy-lectures.org/intro/matplotlib/matplotlib.html) by <NAME>, <NAME> and <NAME>
# - [mpltools](http://tonysyu.github.io/mpltools/index.html) allows easy switching between plot styles
# - [Seaborn](https://github.com/mwaskom/seaborn) facilitates common statistics plots in Matplotlib
#
# ## 2.3 SciPy
# SciPy builds on top of NumPy to provide common tools for scientific programming, such as
#
# - [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)
# - etc., etc
#
# ### 2.3.1 SciPy versus NumPry
# Tips for importing SciPy
#

# +
# Import numpy symbols to scipy name space
import numpy as _num
linalg = None
from numpy import *
from numpy.random import rand, randn
from numpy.fft import fft, ifft
from numpy.lib.scimath import *

__all__  = []
__all__ += _num.__all__
__all__ += ['randn', 'rand', 'fft', 'ifft']
del _num
# Remove the linalg imported from numpy so that the scipy.linalg package can be imported.
del linalg
__all__.remove('linalg')

# +
import numpy as np

a = np.identity(3)           # 单位矩阵
# -

# What is useful in SciPy is the functionality in its subpackages
#
# - `scipy.optimize`, `scipy.integrate`, `scipy.stats`, etc.

from scipy.integrate import quad
from scipy.optimize import brentq
# etc

# ### 2.3.2 Statistics
# The `scipy.stats` subpackage supplies
# - numerous random variable objects (densities, cumulative distributions, random sampling, etc.)
# - some estimation procedures
# - some statistical tests
#
# #### Random Variables and Distributions
# $Beta(\alpha,\beta): prob(x|\alpha,\beta) =\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)}$  
# 其Beta函数$B(\alpha,\beta)=\int_0^1{x^{\alpha-1}(1-x)^{\beta-1}}dx$

import numpy as np
np.random.beta(5, 5, size=3)

# This generates a draw from the distribution below when `a, b = 5, 5`  
# $f(x; a, b) = \frac{x^{(a - 1)} (1 - x)^{(b - 1)}}{\int_0^1 u^{(a - 1)} u^{(b - 1)} du}\qquad (0 \leq x \leq 1)$
#
# In this code we created a so-called `rv_frozen object`, via the call `q = beta(5, 5)`
# 画出 beta（5,5）的概率密度图

# +
import numpy as np
from scipy.stats import beta
import matplotlib.pyplot as plt

q = beta(5, 5)      # Beta(a, b), with a = b = 5
obs = q.rvs(2000)   # 2000 observations
grid = np.linspace(0.01, 0.99, 100)

fig, ax = plt.subplots()
ax.hist(obs, bins=40, normed=True)
ax.plot(grid, q.pdf(grid), 'k-', linewidth=2)
plt.show()
# -

q.cdf(0.4)             # Cumulative distribution function

q.pdf(0.4)      # Density function

q.ppf(0.8)      # Quantile (inverse cdf) function

q.mean()

# The general syntax for creating these objects is
#
# `identifier = scipy.stats.distribution_name(shape_parameters)`
#
# There are also two keyword arguments, `loc` and `scale`, which following our example above, are called as
#
# `identifier = scipy.stats.distribution_name(shape_parameters, loc=c, scale=d)`

# +
import numpy as np
from scipy.stats import beta
import matplotlib.pyplot as plt

obs = beta.rvs(5, 5, size=2000)
grid = np.linspace(0.01, 0.99, 100)

fig, ax = plt.subplots()
ax.hist(obs, bins=40, normed=True)
ax.plot(grid, beta.pdf(grid, 5, 5), 'k-', linewidth=2)
plt.show()
# -

# #### Other Goodies in scipy.stats
# There are a variety statistical functions in `scipy.stats`
#
# For example, `scipy.stats.linregress` implements simple linear regression

# +
from scipy.stats import linregress

x = np.random.randn(200)
y = 2 * x + 0.1 * np.random.randn(200)
gradient, intercept, r_value, p_value, std_err = linregress(x, y)
gradient, intercept


# -

# ### 2.3.3 Roots and Fixed Points
# <p>A <em>root</em> of a real function <em>f</em> on <span class="math">[a,b]</span> is an $x \in [a, b]$ such that $f(x)=0$
# <p>For example, if we plot the function</p>
#
# $f(x) = \sin(4 (x - 1/4)) + x + x^{20} - 1$
# with $x \in [0,1]$ we get</p>
# <img  src="https://lectures.quantecon.org/_images/root.png" style="width: 568.4px; height: 428.4px;" />
# <p>The unique root is approximately 0.408</p>
# <p>Let&#8217;s consider some numerical techniques for finding roots</p>
# #### Bisection 二分法  
# It works for all sufficiently well behaved increasing continuous functions with $f(a)<0<f(b)$

def bisect(f, a, b, tol=10e-5):
    """
    Implements the bisection root finding algorithm, assuming that f is a
    real-valued function on [a, b] satisfying f(a) < 0 < f(b).
    """
    lower, upper = a, b

    while upper - lower > tol:
        middle = 0.5 * (upper + lower)
        # === if root is between lower and middle === #
        if f(middle) > 0:  
            lower, upper = lower, middle
        # === if root is between middle and upper  === #
        else:              
            lower, upper = middle, upper

    return 0.5 * (upper + lower)


# +
from scipy.optimize import bisect

f = lambda x: np.sin(4 * (x - 0.25)) + x + x**20 - 1
bisect(f, 0, 1)
# -

# #### The Newton-Raphson Method （寻根算法）

# +
from scipy.optimize import newton

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

# %timeit bisect(f, 0, 1)              # 运算速度

# %timeit newton(f, 0.2)

# #### Hybrid Methods
# In scipy.optimize, the function `brentq` is such a hybrid method, and a good default

brentq(f, 0, 1)

# %timeit brentq(f, 0, 1)

# #### Multivariate Root Finding
# Use `scipy.optimize.fsolve`, a wrapper for a hybrid method in MINPACK

# #### Fixed Points
# $f(x)= x $

# +
from scipy.optimize import fixed_point

fixed_point(lambda x: x**2, 10)  # 10.0 is an initial guess

# +
from scipy.optimize import fixed_point

fixed_point(lambda x: x**2, 0.2) # 10.0 is an initial guess

# +
from scipy.optimize import fixed_point

fixed_point(lambda x: x**3, 7) 
# -

 from scipy.optimize import fixed_point
# 修改迭代次数，默认值为500
fixed_point(lambda x: x**3, 7, maxiter=10000)

# ### 2.3.4 Optimization

# +
# Most numerical packages provide only functions for minimization
from scipy.optimize import fminbound

fminbound(lambda x: x**2, -1, 2)  # Search in [-1, 2]
# -

# #### Multivariate Optimization
# Multivariate local optimizers include `minimize`, `fmin`, `fmin_powell`, `fmin_cg`, `fmin_bfgs`, and `fmin_ncg`
#
# Constrained multivariate local optimizers include `fmin_l_bfgs_b`, `fmin_tnc`, `fmin_cobyla`
#
# ### 2.3.5 Integration 积分
#
# In SciPy, the relevant module for numerical integration is `scipy.integrate`
#
# A good default for univariate integration is `quad`

# +
from scipy.integrate import quad

integral, error = quad(lambda x: x**2, 0, 1)
integral


# -

# There are other options for univariate integration—a useful one is `fixed_quad`, which is fast and hence works well inside `for` loops
# ### 2.3.6 Linear Algebra
# We saw that NumPy provides a module for linear algebra called `linalg`
# ### 2.3.7 Exercises

# +
def bisect(f, a, b, tol=10e-5):
    """
    Implements the bisection root finding algorithm, assuming that f is a
    real-valued function on [a, b] satisfying f(a) < 0 < f(b).
    """
    lower, upper = a, b

    while upper - lower > tol:
        middle = 0.5 * (upper + lower)
        # === if root is between lower and middle === #
        if f(middle) > 0:  
            lower, upper = lower, middle
        # === if root is between middle and upper  === #
        else:              
            lower, upper = middle, upper

    return 0.5 * (upper + lower)


# -

def bisect(f, a, b, tol=10e-5):
    """
    Implements the bisection root finding algorithm, assuming that f is a
    real-valued function on [a, b] satisfying f(a) < 0 < f(b).
    """
    lower, upper = a, b
    if upper - lower < tol:
        return 0.5 * (upper + lower)
    else:
        middle = 0.5 * (upper + lower)
        print('Current mid point = {}'.format(middle))
        if f(middle) > 0:   # Implies root is between lower and middle
            bisect(f, lower, middle)
        else:               # Implies root is between middle and upper
            bisect(f, middle, upper)


# +
import numpy as np
f = lambda x: np.sin(4 * (x - 0.25)) + x + x**20 - 1

bisect(f, 0, 1)

StarcoderdataJupyterScriptsDedupFiltered