1. 程式人生 > 程式設計 >使用python模擬高斯分佈例子

使用python模擬高斯分佈例子

正態分佈(Normal distribution),也稱“常態分佈”,又名高斯分佈(Gaussian distribution)

正態曲線呈鍾型,兩頭低,中間高,左右對稱因其曲線呈鐘形,因此人們又經常稱之為鐘形曲線。

若隨機變數X服從一個數學期望為μ、方差為σ^2的正態分佈。其概率密度函式為正態分佈的期望值μ決定了其位置,其標準差σ決定了分佈的幅度。當μ = 0,σ = 1時的正態分佈是標準正態分佈。

用python 模擬

#!/usr/bin/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
import seaborn


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

  mu = m
  sigma = np.sqrt(m2 - mu*mu)
  skew = (m3 - 3*mu*m2 + 2*mu**3) / sigma**3
  kurtosis = (m4 - 4*mu*m3 + 6*mu*mu*m2 - 4*mu**3*mu + mu**4) / sigma**4 - 3
  print('手動計算均值、標準差、偏度、峰度:',mu,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,kurtosis


if __name__ == '__main__':
  d = np.random.randn(10000)
  print(d)
  print(d.shape)
  mu,kurtosis = calc_statistics(d)
  print('函式庫計算均值、標準差、偏度、峰度:',kurtosis)
  # 一維直方圖
  mpl.rcParams['font.sans-serif'] = 'SimHei'
  mpl.rcParams['axes.unicode_minus'] = False
  plt.figure(num=1,facecolor='w')
  y1,x1,dummy = plt.hist(d,bins=30,normed=True,color='g',alpha=0.75,edgecolor='k',lw=0.5)
  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('高斯分佈,樣本個數:%d' % d.shape[0])
  plt.grid(b=True,ls=':',color='#404040')
  # plt.show()

  d = np.random.randn(100000,2)
  mu,kurtosis)

  # 二維影象
  N = 30
  density,edges = np.histogramdd(d,bins=[N,N])
  print('樣本總數:',np.sum(density))
  density /= density.max()
  x = y = np.arange(N)
  print('x = ',x)
  print('y = ',y)
  t = np.meshgrid(x,y)
  print(t)
  fig = plt.figure(facecolor='w')
  ax = fig.add_subplot(111,projection='3d')
  # ax.scatter(t[0],t[1],density,c='r',s=50*density,marker='o',depthshade=True,edgecolor='k')
  ax.plot_surface(t[0],cmap=cm.Accent,rstride=1,cstride=1,alpha=0.9,lw=0.75,edgecolor='k')
  ax.set_xlabel('X')
  ax.set_ylabel('Y')
  ax.set_zlabel('Z')
  plt.title('二元高斯分佈,樣本個數:%d' % d.shape[0],fontsize=15)
  plt.tight_layout(0.1)
  plt.show()

來個6的

二元高斯分佈方差比較

#!/usr/bin/python
# -*- coding:utf-8 -*-

import numpy as np
from scipy import stats
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm


if __name__ == '__main__':
  x1,x2 = np.mgrid[-5:5:51j,-5:5:51j]
  x = np.stack((x1,x2),axis=2)
  print('x1 = \n',x1)
  print('x2 = \n',x2)
  print('x = \n',x)

  mpl.rcParams['axes.unicode_minus'] = False
  mpl.rcParams['font.sans-serif'] = 'SimHei'
  plt.figure(figsize=(9,8),facecolor='w')
  sigma = (np.identity(2),np.diag((3,3)),np.diag((2,5)),np.array(((2,1),(1,5))))
  for i in np.arange(4):
    ax = plt.subplot(2,2,i+1,projection='3d')
    norm = stats.multivariate_normal((0,0),sigma[i])
    y = norm.pdf(x)
    ax.plot_surface(x1,x2,lw=0.3,edgecolor='#303030')
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
  plt.suptitle('二元高斯分佈方差比較',fontsize=18)
  plt.tight_layout(1.5)
  plt.show()

影象好看嗎?

以上這篇使用python模擬高斯分佈例子就是小編分享給大家的全部內容了,希望能給大家一個參考,也希望大家多多支援我們。