2.数学基础

import numpy as np
import math

a = [5, 6, 6, 5, 9, 0, 7, 2, 8, 7, 7, 9, 4, 5, 6, 1, 8, 6, 8, 5]
b = [5, 6, 6, 5, 9, 0, 7, 2, 8, 7, 7, 9, 4, 5, 6, 1, 8, 6, 8, 5]

# 协方差
def covariance(x, y):
    if len(x) != len(y):
        raise Exception('两个列表的长度必须相同')
    n = len(x)
    avg_a = sum(x) / n #均值
    avg_b = sum(y) / n #均值
    mean_a = [round(i - avg_a, 5) for i in x] # 均差值,精确到小数点3位
    mean_b = [round(i - avg_b, 5) for i in y] # 均差值,精确到小数点3位
    return np.dot(mean_a, mean_b)/(n - 1) # dot 点乘，表示，对应索引的两个值相乘，最后在相加

print("协方差", covariance(a, b)) #6.115789473684211

import numpy as np
import math

a = [5, 6, 6, 5, 9, 0, 7, 2, 8, 7, 7, 9, 4, 5, 6, 1, 8, 6, 8, 5]
b = [5, 6, 6, 5, 9, 0, 7, 2, 8, 7, 7, 9, 4, 5, 6, 1, 8, 6, 8, 5]

# 协方差
def covariance(x, y):
    if len(x) != len(y):
        raise Exception('两个列表的长度必须相同')
    n = len(x)
    avg_a = sum(x) / n #均值
    avg_b = sum(y) / n #均值
    mean_a = [round(i - avg_a, 5) for i in x] # 均差值,精确到小数点3位
    mean_b = [round(i - avg_b, 5) for i in y] # 均差值,精确到小数点3位
    return np.dot(mean_a, mean_b)/(n - 1) # dot 点乘，表示，对应索引的两个值相乘，最后在相加

print("协方差", covariance(a, b))
# 标准差
def standard_deviation(x):
    avg = sum(x) / len(x)  # 样本期望值
    sumValue = sum([(avg - i) ** 2 for i in x])
    return math.sqrt(sumValue / (len(x) - 1))

# 相关系数
def correlation(x, y):
    stdev_x = standard_deviation(x) # 标准差
    stdev_y = standard_deviation(y) # 标准差
    if stdev_x > 0 and stdev_y > 0:
        return round(covariance(x, y) / stdev_x / stdev_y, 3)
    else:
        return 0  # 如果没有变动，相关系数为零

print(correlation(a, b))

# 使用numpy计算,返回结果均为(2, 2)矩阵result[i][j]，表示第i个变量与第j个变量的协方差和相关系数。我们一般取result[1][1]即可。
x = np.vstack((a, b))
print("协方差", np.cov(x))
print("协方差", np.cov(a, b))
print("协方差", np.cov(a, b, bias=True))
print("相关系数", np.corrcoef(a, b))
print("相关系数", np.corrcoef(x))

import math
from matplotlib import pyplot

# 均值 μ（mu）和标准差 σ （sigma）
def normal_pdf(x, mu = 0, sigma = 1):
    return 1 / (math.sqrt(2 * math.pi) * sigma) * math.exp(
        -( (x - mu) ** 2 ) / (2 * sigma ** 2)
    )

xs = [x / 10 for x in range(-50, 50)]

pyplot.plot(xs, [normal_pdf(x, mu=0,sigma=1) for x in xs],'-',label='mu=0,sigma=1')
pyplot.plot(xs, [normal_pdf(x, mu=0,sigma=2) for x in xs],'--',label='mu=0,sigma=2')
pyplot.plot(xs, [normal_pdf(x, mu=0,sigma=0.5) for x in xs],':',label='mu=0,sigma=0.5')
pyplot.plot(xs, [normal_pdf(x, mu=-1,sigma=1) for x in xs],'-.',label='mu=-1,sigma=1')

pyplot.legend()
pyplot.title("多个正态分布的概率密度函数")
pyplot.show()

如果 μ=0 并且 σ =1，这个分布称为标准正态分布。
如果 Z 是服从标准正态分布的随机变量，则有如下转换式：
X= σZ+μ
其中 X 也是正态分布，但均值是 μ，标准差是 σ 。

相反，如果 X 是均值为 μ 标准差为 σ 的正态分布，那么：
Z=(X-μ)/ σ
是标准正态分布的随机变量。

累积分
def normal_cdf(x, mu=0,sigma=1):
    return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2
# 我们再绘出一系列概率累积分布函数（如图 6-3）：
xs = [x / 10.0 for x in range(-50, 50)]
plt.plot(xs,[normal_cdf(x,sigma=1) for x in xs],'-',label='mu=0,sigma=1')
plt.plot(xs,[normal_cdf(x,sigma=2) for x in xs],'--',label='mu=0,sigma=2')
plt.plot(xs,[normal_cdf(x,sigma=0.5) for x in xs],':',label='mu=0,sigma=0.5')
plt.plot(xs,[normal_cdf(x,mu=-1) for x in xs],'-.',label='mu=-1,sigma=1')
plt.legend(loc=4) # 底部右边
plt.title("多个正态分布的累积分布函数")
plt.show()

我们有时会需要对 normal_cdf 取逆，从而可以求出特定的概率的相应值。不存在计算逆函数的简便方法，
但由于 normal_cdf 连续且严格递增，因而我们可以使用二分查找（https://en.wikipedia.org/wiki/Binary_search_algorithm）的方法：

def inverse_normal_cdf(p, mu=0, sigma=1, tolerance=0.00001):
"""find approximate inverse using binary search"""
# 如果非标准型，先调整单位使之服从标准型
if mu != 0 or sigma != 1:
    return mu + sigma * inverse_normal_cdf(p, tolerance=tolerance)
low_z, low_p = -10.0, 0 # normal_cdf(-10)是（非常接近）0
hi_z, hi_p = 10.0, 1 # normal_cdf(10)是（非常接近）1
while hi_z - low_z > tolerance:
    mid_z = (low_z + hi_z) / 2 # 考虑中点
    mid_p = normal_cdf(mid_z) # 和cdf在那里的值
    if mid_p < p:
        # midpoint仍然太低，搜索比它大的值
        low_z, low_p = mid_z, mid_p
    elif mid_p > p:
        # midpoint仍然太高，搜索比它小的值
        hi_z, hi_p = mid_z, mid_p
    else:
        break
return mid_z
这个函数反复分割区间，直到分割到一个足够接近于期望概率的精细的 Z 值。

from functools import partial
import inspect

from matplotlib import pyplot as plt

# 求差商
def difference_quotient(f, x, h):
    if inspect.isfunction(f) :
        return (f(x + h) - f(x)) / h
    else:
        return None

# 函数 f(x) = x * x
def square(x):
    return x * x

# 我们已知f(x)的导数
def derivative(x):
    return 2 * x

derivative_estimate = partial(difference_quotient, square, h=0.000001)

x = range(-10,10)
plt.title("精确的导数值与估计值")
plt.plot(x, [y for y in map(derivative_estimate, x)], 'b+', label='差商值') # 用 + 表示
plt.plot(x, [y for y in map(derivative, x)], 'rx', label='导数值') # 用 x 表示
plt.legend(loc=9)
plt.show()   

通过生产的plot可知，差商与导数值，拟合度高，当 h 值越小， 两者的拟合度越高。

当 f 是一个多变量函数时，它有多个偏导数，每个偏导数表示仅有一个输入变量发生微小变化时函数 f 的变化。
我们把导数看成是其第 i 个变量的函数，其他变量保持不变，以此来计算它第 i 个偏导数：

def partial_difference_quotient(f, v, i, h):
    """compute the ith partial difference quotient of f at v"""
    w = [v_j + (h if j == i else 0) # 只对v的第i个元素增加h
        for j, v_j in enumerate(v)]
    return (f(w) - f(v)) / h

再以同样的方法估算它的梯度函数：

def estimate_gradient(f, v, h=0.00001):
    return [partial_difference_quotient(f, v, i, h)
        for i, _ in enumerate(v)]    

# 1. 这个是常见的计算相关系数的方法
import numpy
import math
import random
import matplotlib.pyplot as plt

# 协方差
def covariance(x, y):
    if len(x) != len(y):
        raise Exception('两个列表的长度必须相同')
    n = len(x)
    avg_a = sum(x) / n #均值
    avg_b = sum(y) / n #均值
    mean_a = [round(i - avg_a, 5) for i in x] # 均差值,精确到小数点3位
    mean_b = [round(i - avg_b, 5) for i in y] # 均差值,精确到小数点3位
    return numpy.dot(mean_a, mean_b)/(n - 1) # dot 点乘，表示，对应索引的两个值相乘，最后在相加

# 标准差
def standard_deviation(x):
    avg = sum(x) / len(x)  # 样本期望值
    sumValue = sum([(avg - i) ** 2 for i in x])
    return math.sqrt(sumValue / (len(x) - 1))

# 相关系数
def correlation(x, y):
    stdev_x = standard_deviation(x)
    stdev_y = standard_deviation(y)
    if stdev_x > 0 and stdev_y > 0:
        return covariance(x, y) / stdev_x / stdev_y
    else:
        return 0 # 如果没有变动，相关系数为零

def get_row(A, i):
    return A[i] # A[i]是第i行

def get_column(A, j):
    return [A_i[j] # 第A_i行的第j个元素
        for A_i in A] # 对每个A_i行

def make_matrix(num_rows, num_cols, entry_fn):
    return [[entry_fn(i, j) # 根据i创建一个列表
        for j in range(num_cols)] # [entry_fn(i, 0), ... ]
        for i in range(num_rows)] # 为每一个i创建一个列表

def correlation_matrix(data):
    _, num_columns = numpy.shape(data)
    def matrix_entry(i, j):
        return correlation(get_column(data, i), get_column(data, j))
    return make_matrix(num_columns, num_columns, matrix_entry)

# 2. 一个更为直观的方法（如果维度不太多）是做散点图矩阵（图 10-4），以展示配对散点图。
# 通过命令 plt.subplots() 可以生成子图。我们给出了行数和列数，它返回一个 figure 对
# 象（我们不会用到它）和一个 axes 对象的二维数组（每个都会绘出）
data =[
    [
        random.randrange(10) for _ in range(4)
    ] for _ in range(10)
]

_, num_columns = numpy.shape(data)
fig, ax = plt.subplots(num_columns, num_columns)

for i in range(num_columns):
    for j in range(num_columns):
        if i != j:# x轴上column_j对y轴上column_i的散点
            ax[i][j].scatter(get_column(data, j), get_column(data, i))
        else:# 只有当 i == j时显示序列名
            ax[i][j].annotate("series " + str(i), (0.5, 0.5), xycoords='axes fraction', ha="center", va="center")
        # 除了图的左侧和下方之外，隐藏图的标记
        if i < num_columns - 1:
            ax[i][j].xaxis.set_visible(False)
        if j > 0:
            ax[i][j].yaxis.set_visible(False)
# 修复右下方和左上方的图标记
# 因为它们只有文本，是错误的
ax[-1][-1].set_xlim(ax[0][-1].get_xlim())
ax[0][0].set_ylim(ax[0][1].get_ylim())
plt.show()