python二维正态分布

#  -*- coding:utf-8 -*-
import numpy as np
from scipy import stats
import math
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

def calc_statistics(x):
    n = x.shape[0] #样本个数

    m = 0#期望
    m2 = 0#平方的期望
    m3 = 0#三次方的期望
    m4 = 0#四次方的期望
    for t in x:
        #向量的加法
        m += t
        m2 += t*t
        m3 += t**3
        m4 += t**4
    m /= n
    m2 /= n
    m3 /= n
    m4 /= n

    #标准差 = E((X - E(X))^2) = E(X^2) - E(X)^2
    sigma = np.sqrt(m2-m*m)
    #求偏度  = E((X-E(X))^3) = (m3 - 3*m*m2 + 2*m**3) / sigma**3
    skew = (m3 - 3*m*sigma**2 - m**3) / sigma**3
    #求峰度
    kurtosis = m4 / sigma**4 - 3
    print('手动计算均值、标准差、偏度、峰度：', m, sigma, skew, kurtosis)

    #使用系统函数验证
    mu = np.mean(x,axis=0)
    sigma = np.std(x,axis=0)
    skew = stats.skew(x)
    kurtosis = stats.kurtosis(x)
    return mu,sigma,skew,kurtosis



if __name__ == '__main__':
    # d = np.random.randn(100000)
    # print(d)
    # mu, sigma, skew, kurtosis = calc_statistics(d)
    # print('函数计算均值、标准差、偏度、峰度：', mu, sigma, skew, kurtosis)
    # # 一维直方图
    # mpl.rcParams[u'font.sans-serif'] = 'SimHei'
    # mpl.rcParams[u'axes.unicode_minus'] = False
    # #画出统计直方图
    # #bins直方图的条数
    # #density=True 画出趋势图
    # #y1:x1每个中每个值出现的次数的度量
    # #x1:d的值的范围
    # y1, x1, dummy = plt.hist(d, bins=50,density=True, color='g', alpha=0.75)
    # t = np.arange(x1.min(), x1.max(), 0.05)
    # #绘制标准正态分布的曲线
    # y = np.exp(-t ** 2 / 2) / math.sqrt(2 * math.pi)
    # plt.plot(t, y, 'r-', lw=2)
    # plt.title(u'高斯分布，样本个数：%d' % d.shape[0])
    # plt.grid(True)
    # plt.show()

    #二维
    print("二维正态分布")
    d = np.random.randn(100000, 2)
    mu, sigma, skew, kurtosis = calc_statistics(d)
    print('函数计算均值、标准差、偏度、峰度：', mu, sigma, skew, kurtosis)
    # 二维图像
    N = 50
    #density：edges中每个值出现的次数的度量
    density, edges = np.histogramdd(d, bins=[N, N])
    print('样本总数：', np.sum(density))
    density /= density.max()
    x = y = np.arange(N)
    t = np.meshgrid(x, y)

    fig = plt.figure(facecolor='gray')
    ax = fig.add_subplot(111, projection='3d')
    #x,y,z
    ax.scatter(t[0], t[1], density, c='r', s=15 * density, marker='o', depthshade=True)
    # ax.plot_surface(t[0], t[1], density, cmap=cm.Accent, rstride=2, cstride=2, alpha=0.9, lw=0.75)
    # ax.set_xlabel(u'X')
    # ax.set_ylabel(u'Y')
    # ax.set_zlabel(u'Z')
    # plt.title(u'二元高斯分布，样本个数：%d' % d.shape[0], fontsize=20)
    # plt.tight_layout(0.1)
    plt.show()

#  -*- coding:utf-8 -*-
import numpy as np
from scipy import stats
import math
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

def calc_statistics(x):
    n = x.shape[0] #样本个数

    m = 0#期望
    m2 = 0#平方的期望
    m3 = 0#三次方的期望
    m4 = 0#四次方的期望
    for t in x:
        #向量的加法
        m += t
        m2 += t*t
        m3 += t**3
        m4 += t**4
    m /= n
    m2 /= n
    m3 /= n
    m4 /= n

    #标准差 = E((X - E(X))^2) = E(X^2) - E(X)^2
    sigma = np.sqrt(m2-m*m)
    #求偏度  = E((X-E(X))^3) = (m3 - 3*m*m2 + 2*m**3) / sigma**3
    skew = (m3 - 3*m*sigma**2 - m**3) / sigma**3
    #求峰度
    kurtosis = m4 / sigma**4 - 3
    print('手动计算均值、标准差、偏度、峰度：', m, sigma, skew, kurtosis)

    #使用系统函数验证
    mu = np.mean(x,axis=0)
    sigma = np.std(x,axis=0)
    skew = stats.skew(x)
    kurtosis = stats.kurtosis(x)
    return mu,sigma,skew,kurtosis



if __name__ == '__main__':
    # d = np.random.randn(100000)
    # print(d)
    # mu, sigma, skew, kurtosis = calc_statistics(d)
    # print('函数计算均值、标准差、偏度、峰度：', mu, sigma, skew, kurtosis)
    # # 一维直方图
    # mpl.rcParams[u'font.sans-serif'] = 'SimHei'
    # mpl.rcParams[u'axes.unicode_minus'] = False
    # #画出统计直方图
    # #bins直方图的条数
    # #density=True 画出趋势图
    # #y1:x1每个中每个值出现的次数的度量
    # #x1:d的值的范围
    # y1, x1, dummy = plt.hist(d, bins=50,density=True, color='g', alpha=0.75)
    # t = np.arange(x1.min(), x1.max(), 0.05)
    # #绘制标准正态分布的曲线
    # y = np.exp(-t ** 2 / 2) / math.sqrt(2 * math.pi)
    # plt.plot(t, y, 'r-', lw=2)
    # plt.title(u'高斯分布，样本个数：%d' % d.shape[0])
    # plt.grid(True)
    # plt.show()

    #二维
    print("二维正态分布")
    d = np.random.randn(100000, 2)
    mu, sigma, skew, kurtosis = calc_statistics(d)
    print('函数计算均值、标准差、偏度、峰度：', mu, sigma, skew, kurtosis)
    # 二维图像
    N = 50
    #density：edges中每个值出现的次数的度量
    density, edges = np.histogramdd(d, bins=[N, N])
    print('样本总数：', np.sum(density))
    density /= density.max()
    x = y = np.arange(N)
    t = np.meshgrid(x, y)

    fig = plt.figure(facecolor='gray')
    ax = fig.add_subplot(111, projection='3d')
    #x,y,z
    ax.scatter(t[0], t[1], density, c='r', s=15 * density, marker='o', depthshade=True)
    # ax.plot_surface(t[0], t[1], density, cmap=cm.Accent, rstride=2, cstride=2, alpha=0.9, lw=0.75)
    # ax.set_xlabel(u'X')
    # ax.set_ylabel(u'Y')
    # ax.set_zlabel(u'Z')
    # plt.title(u'二元高斯分布，样本个数：%d' % d.shape[0], fontsize=20)
    # plt.tight_layout(0.1)
    plt.show()