草庐IT

正态分布,二维正态分布,卡方分布,学生t分布——概率分布学习 python

nature1949 2023-08-14 原文

目录

基本概念

概率密度函数(PDF: Probability Density Function)

累积分布函数(CDF: Cumulative Distribution Function)

核密度估计((kernel density estimation)

1.正态分布

概率密度函数(pdf)

正态分布累积分布函数(CDF)

正态分布核密度估计(kde)

正态分布四则运算

二维正态分布(逐渐补充)

马氏距离

2.卡方分布

概率密度函数(pdf):

 卡方分布表:

卡方分布相关计算

生成卡方分布随机数

3.学生t分布

概率密度函数(pdf):


基本概念

概率密度函数(PDF: Probability Density Function)

连续随机变量的概率分布特性。

累积分布函数(CDF: Cumulative Distribution Function)

在x点左侧事件发生的总和。

CDF特性:

①因为累计分布函数是计算x点左侧的点的数量,所以累计分布函数CDF是单调递增的。

②所有的CDF中,在x趋近-∞时,CDF趋近于0,当x趋近+∞时,CDF趋近于1。

③对于给定的数据集,CDF是唯一的

核密度估计((kernel density estimation)

核密度估计(kernel density estimation,KDE)是在概率论中用来估计未知的密度函数,属于非参数检验方法之一,通过核密度估计图可以比较直观的看出数据样本本身的分布特征。

scipy中的stats.gaussian_kde可以计算高斯核函数的密度函数,而且提供了直接计算区间的累计密度函数,integrate_box_1d(low=-np.Inf, high=x)。

1.正态分布

表示为:,其中期望为μ,方差为

概率密度函数(pdf)

python画图效果及代码(包含随机数生成):

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import matplotlib.cm as cm
import math
import scipy.stats as stats
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

################################                正态分布             ###########################
# 根据均值、标准差,求指定范围的正态分布概率值
def normfun(x, mu, sigma):
    pdf = np.exp(-((x - mu)**2)/(2*sigma**2)) / (sigma * np.sqrt(2*np.pi))
    return pdf                 
np.random.seed(0)   ##  定义一个随机数种子
result = np.random.normal(loc=10, scale=16, size=1000) # 均值为10,标准差为16
##           !!!强调,以上参数中scale为标准差(方差的根号),不是方差,
# 设定 x,y 轴,载入刚才的正态分布函数
x = np.arange(min(result), max(result), 0.1)
y = normfun(x, result.mean(), result.std())
plt.plot(x, y) # 这里画出理论的正态分布概率曲线
plt.hist(result, bins=20, rwidth=0.8, density=True)     ##  柱状图
plt.title('distribution')
plt.xlabel('temperature')
plt.ylabel('probability')
plt.show()

正态分布累积分布函数(CDF)

################################                累积分布函数cdf             ###########################
#计算正态概率密度函数在x处的值
def norm_dist_prob(theta):
    y = stats.norm.pdf(theta, loc=np.mean(data), scale=np.std(data))
    return y

#计算正态分布累积概率值
def norm_dist_cdf(theta):
    y = stats.norm.cdf(theta,loc=np.mean(data), scale=np.std(data))
    return y
##  数据生成
data = np.random.normal(loc=0.0, scale=10, size=1000)

x = np.linspace(stats.norm.ppf(0.01,loc=np.mean(data), scale=np.std(data)),
                stats.norm.ppf(0.99,loc=np.mean(data), scale=np.std(data)), len(data))  #linspace() 函数返回指定间隔内均匀间隔数字的 ndarray。

y1=norm_dist_prob(x)
y2=norm_dist_cdf(x)
plt.plot(x, y1,'g', label='pdf')
plt.plot(x, y2,'r', label='cdf1')
#或
sns.kdeplot(data,cumulative=True, label='cdf2')
plt.legend()

正态分布核密度估计(kde)

################################                核密度估计             ###########################
##  数据生成
data = np.random.normal(loc=0.0, scale=10, size=1000)
##  本程序是根据数据进行概率密度估计
density = stats.gaussian_kde(data)   #, bw_method=None, weights=[i[4] for i in data1]
density.covariance_factor = lambda : .25    #   lambda : .25
density._compute_covariance()
density.set_bandwidth(bw_method='silverman')        ##  调用set_bandwidth 后计算的新带宽用于估计密度的后续评估。可选‘scott’, ‘silverman’
xs = np.linspace(min(data), max(data), 200)
fig, ax = plt.subplots()
ax.plot(xs, density(xs), 'r')
ax.fill_between(xs, density(xs), color="r", alpha=0.1)
ax.hist(data, bins=30, rwidth=0.96, density =True, alpha=0.6,color = 'steelblue', edgecolor = 'w', label = 'dimensional histogram statistic ')

##  或者用seaborn
fig, ax = plt.subplots()
sns.distplot(data, hist=True, kde=True, rug=True, bins=20, ax=ax)
#   通过hist和kde参数调节是否显示直方图及核密度估计(默认hist,kde均为True)
#   bins:int或list,控制直方图的划分
#   rug:控制是否生成观测数值的小细条
#   ax = sns.distplot(x, rug=True, rug_kws={"color": "g"},
#        ...                   kde_kws={"color": "k", "lw": 3, "label": "KDE"},
#        ...                   hist_kws={"histtype": "step", "linewidth": 3,
#        ...                             "alpha": 1, "color": "g"})fig, ax = plt.subplots()

正态分布四则运算

 两个相互独立的正态分布分别满足

则:

二维正态分布(逐渐补充)

其生成及协方差椭圆的python实现如下:

################################                二维正态分布             ###########################
from matplotlib.patches import Ellipse
def get_error_ellipse_parameters(cov, confidence=None, sigma=None):
    """Returns parameters of an ellipse which contains a specified
    amount of normally-distributed 2D data, where the data is
    characterised by its covariance matrix.
    
    Parameters
    ----------
    cov : array_like
        Input covariance matrix of shape (2,2)
    confidence : float
        Fraction of data points within ellipse. 0 < confidence < 1.
        If confidence is not given, it is calculated according to sigma.
    sigma : float
        Length of axes of the ellipse in standard deviations. If 
        confidence is also given, sigma is ignored.
    
    Returns
    -------
    semi_major : float
        Length of major semiaxis of ellipse.
    semi_minor : float
        Length of minor semiaxis of ellipse.
    angle : float
        Rotation angle of ellipse in radian.
    confidence : float
        Fraction of data expected to lie within the ellipse.
    sigma : float
        Length of major and minor semiaxes in standard deviations.
    """
    cov = np.array(cov)
    if(cov.shape != (2,2)):
        raise ValueError("The covariance matrix needs to be of shape (2,2)")
    if(confidence == None and sigma == None):
        raise RuntimeError("One of confidence and sigma is needed as input argument")
    if(confidence and sigma):
        print("Argument sigma is ignored as confidence is also provided!")
    
    if(confidence == None):
        if(sigma < 0):
            raise ValueError("Sigma needs to be positive")
        #scaling = np.square(sigma)
        scaling = sigma
        confidence = stats.chi2.cdf(scaling, 2)
    if(sigma == None):
        if(confidence > 1 or confidence < 0):
            raise ValueError("Ensure that confidence lies between 0 and 1")
        scaling = stats.chi2.ppf(confidence, 2)
        #sigma = np.sqrt(scaling)
        sigma = scaling
    eigenvalues, eigenvectors = np.linalg.eig(cov)
    maxindex = np.argmax(eigenvalues)
    vx, vy = eigenvectors[:, maxindex]
    angle = np.arctan2(vy, vx)
    semi_minor, semi_major = np.sqrt(np.sort(eigenvalues) * scaling)
    print("With sigma = {:.3f}, {:.1f}% of data points lie within ellipse.".format(sigma, confidence * 100))
    return semi_major, semi_minor, angle, confidence, sigma

mu = [1,2]
cov = [[50,30],[30,50]] #sigma
#   随机数生成
z = stats.multivariate_normal(mu, cov)
data_points = z.rvs(size = 5000)

fig, ax = plt.subplots()
plt.scatter(data_points[:,0], data_points[:,1], alpha = .5)

#   画置信度椭圆
confidence = 0.95
semi_major, semi_minor, angle, confidence, sigma = get_error_ellipse_parameters(cov, confidence = confidence)
ax.add_patch(Ellipse(mu, 2*semi_major, 2*semi_minor, 180*angle/np.pi, facecolor = 'none', edgecolor = 'red', label = 'Confidence = {:.0f}% (sigma = {:.2f})'.format(confidence * 100, sigma)))
sigma = 1
semi_major, semi_minor, angle, confidence, sigma, = get_error_ellipse_parameters(cov, sigma = sigma)
ax.add_patch(Ellipse(mu, 2*semi_major, 2*semi_minor, 180*angle/np.pi, facecolor = 'none', edgecolor = 'yellow', label = 'Sigma = {:.0f} (confidence = {:.1f}%)'.format(sigma, confidence * 100)))
plt.legend()
plt.show()

马氏距离

计算马氏距离(Mahalanobis Distance)。一维马氏距离定义为:

iv = [[1, 0.5, 0.5], [0.5, 1, 0.5], [0.5, 0.5, 1]]
md = distance.mahalanobis([1, 0, 0], [0, 1, 0], iv)
print(md)
#   或
p = np.array([1,1])
distr = np.array([2,2])
cov = [[1,0.2],
        [0.2,1]]
dis = distance.mahalanobis(p, distr, cov)
# p: 一个点    
# distr : 一个分布    
# 计算分布的协方差矩阵    
#cov = np.cov(distr, rowvar=False)    
# 选取分布中各维度均值所在点    
#avg_distri = np.average(distr, axis=0)    
print(dis)

2.卡方分布

卡方分布,也写作:分布。服从自由度为n的卡方分布,记作,其均值为 n,方差为2n。

若n个相互独立的随机变量ξ₁,ξ₂,...,ξn ,均服从标准正态分布N(0,1),则这n个服从标准正态分布的随机变量的平方和构成一新的随机变量,其分布规律称为卡方分布(chi-square distribution)。

 直观说:如果 X1,X2,X3...X„是 n个具有标准正态分布的独立变量,那么其平方和,满足具有n个自由度的分布。

概率密度函数(pdf):

其中,是Gamma函数,n为自由度,一般情况

################################                卡方分布             ###########################
for PDF in range(1,8):
    plt.plot(np.linspace(0,15,100),stats.chi2.pdf(np.linspace(0,15,100),df=PDF),label='k='+str(PDF))
plt.tick_params(axis="both",which="major",labelsize=18)
plt.axhline(y=0,color="black",linewidth=1.3,alpha=.7)
plt.legend()

 卡方分布表:

卡方分布相关计算

##  卡方分布相关计算
#   累积分布函数
x = stats.chi2.cdf(5.991, df=2)
#   百分比点函数(与cdf—百分位数相反)
a = stats.chi2.ppf(0.95, df=2)  
print(x,a)

生成卡方分布随机数

#生成随机数
r = stats.chi2.rvs(df=df, size=1000)

3.学生t分布

Student's t-distribution,简称为t分布。

假设随机变量Z服从标准正态分布N(0,1),另一随机变量V服从m自由度的分布,进一步假设Z和 V 彼此独立,则下列的数量t服从自由度为m的学生t分布:

概率密度函数(pdf):

################################                t分布             ###########################
x = np.linspace( -3, 3, 100)
plt.plot(x, stats.t.pdf(x,1), label='df=1')
plt.plot(x, stats.t.pdf(x,2), label='df=20')
plt.plot(x, stats.t.pdf(x,100), label = 'df=100')
plt.plot( x[::5], stats.norm.pdf(x[::5]),'kx',  label='normal')
##  累积分布函数cdf
y = stats.t.cdf(x,df=100, loc=0, scale=1)
plt.plot(x,y, label='cdf')
plt.legend()

有关正态分布,二维正态分布,卡方分布,学生t分布——概率分布学习 python的更多相关文章

  1. python - 如何使用 Ruby 或 Python 创建一系列高音调和低音调的蜂鸣声? - 2

    关闭。这个问题是opinion-based.它目前不接受答案。想要改进这个问题?更新问题,以便editingthispost可以用事实和引用来回答它.关闭4年前。Improvethisquestion我想在固定时间创建一系列低音和高音调的哔哔声。例如:在150毫秒时发出高音调的蜂鸣声在151毫秒时发出低音调的蜂鸣声200毫秒时发出低音调的蜂鸣声250毫秒的高音调蜂鸣声有没有办法在Ruby或Python中做到这一点?我真的不在乎输出编码是什么(.wav、.mp3、.ogg等等),但我确实想创建一个输出文件。

  2. ruby - 分布式事务和队列,ruby,erlang,scala - 2

    我有一个涉及多台机器、消息队列和事务的问题。因此,例如用户点击网页,点击将消息发送到另一台机器,该机器将付款添加到用户的帐户。每秒可能有数千次点击。事务的所有方面都应该是容错的。我以前从未遇到过这样的事情,但一些阅读表明这是一个众所周知的问题。所以我的问题。我假设安全的方法是使用两阶段提交,但协议(protocol)是阻塞的,所以我不会获得所需的性能,我是否正确?我通常写Ruby,但似乎Redis之类的数据库和Rescue、RabbitMQ等消息队列系统对我的帮助不大——即使我实现某种两阶段提交,如果Redis崩溃,数据也会丢失,因为它本质上只是内存。所有这些让我开始关注erlang和

  3. Python 相当于 Perl/Ruby ||= - 2

    这个问题在这里已经有了答案:关闭10年前。PossibleDuplicate:Pythonconditionalassignmentoperator对于这样一个简单的问题表示歉意,但是谷歌搜索||=并不是很有帮助;)Python中是否有与Ruby和Perl中的||=语句等效的语句?例如:foo="hey"foo||="what"#assignfooifit'sundefined#fooisstill"hey"bar||="yeah"#baris"yeah"另外,类似这样的东西的通用术语是什么?条件分配是我的第一个猜测,但Wikipediapage跟我想的不太一样。

  4. java - 什么相当于 ruby​​ 的 rack 或 python 的 Java wsgi? - 2

    什么是ruby​​的rack或python的Java的wsgi?还有一个路由库。 最佳答案 来自Python标准PEP333:Bycontrast,althoughJavahasjustasmanywebapplicationframeworksavailable,Java's"servlet"APImakesitpossibleforapplicationswrittenwithanyJavawebapplicationframeworktoruninanywebserverthatsupportstheservletAPI.ht

  5. 华为OD机试用Python实现 -【明明的随机数】 2023Q1A - 2

    华为OD机试题本篇题目:明明的随机数题目输入描述输出描述:示例1输入输出说明代码编写思路最近更新的博客华为od2023|什么是华为od,od薪资待遇,od机试题清单华为OD机试真题大全,用Python解华为机试题|机试宝典【华为OD机试】全流程解析+经验分享,题型分享,防作弊指南华为o

  6. python - 如何读取 MIDI 文件、更改其乐器并将其写回? - 2

    我想解析一个已经存在的.mid文件,改变它的乐器,例如从“acousticgrandpiano”到“violin”,然后将它保存回去或作为另一个.mid文件。根据我在文档中看到的内容,该乐器通过program_change或patch_change指令进行了更改,但我找不到任何在已经存在的MIDI文件中执行此操作的库.他们似乎都只支持从头开始创建的MIDI文件。 最佳答案 MIDIpackage会为您完成此操作,但具体方法取决于midi文件的原始内容。一个MIDI文件由一个或多个音轨组成,每个音轨是十六个channel中任何一个上的

  7. 「Python|Selenium|场景案例」如何定位iframe中的元素? - 2

    本文主要介绍在使用Selenium进行自动化测试或者任务时,对于使用了iframe的页面,如何定位iframe中的元素文章目录场景描述解决方案具体代码场景描述当我们在使用Selenium进行自动化测试的时候,可能会遇到一些界面或者窗体是使用HTML的iframe标签进行承载的。对于iframe中的标签,如果直接查找是无法找到的,会抛出没有找到元素的异常。比如近在咫尺的例子就是,CSDN的登录窗体就是使用的iframe,大家可以尝试通过F12开发者模式查看到的tag_name,class_name,id或者xpath来定位中的页面元素,会抛出NoSuchElementException异常。解决

  8. 计算机毕业设计ssm+vue基本微信小程序的小学生兴趣延时班预约小程序 - 2

    项目介绍随着我国经济迅速发展,人们对手机的需求越来越大,各种手机软件也都在被广泛应用,但是对于手机进行数据信息管理,对于手机的各种软件也是备受用户的喜爱小学生兴趣延时班预约小程序的设计与开发被用户普遍使用,为方便用户能够可以随时进行小学生兴趣延时班预约小程序的设计与开发的数据信息管理,特开发了小程序的设计与开发的管理系统。小学生兴趣延时班预约小程序的设计与开发的开发利用现有的成熟技术参考,以源代码为模板,分析功能调整与小学生兴趣延时班预约小程序的设计与开发的实际需求相结合,讨论了小学生兴趣延时班预约小程序的设计与开发的使用。开发环境开发说明:前端使用微信微信小程序开发工具:后端使用ssm:VU

  9. python ffmpeg 使用 pyav 转换 一组图像 到 视频 - 2

    2022/8/4更新支持加入水印水印必须包含透明图像,并且水印图像大小要等于原图像的大小pythonconvert_image_to_video.py-f30-mwatermark.pngim_dirout.mkv2022/6/21更新让命令行参数更加易用新的命令行使用方法pythonconvert_image_to_video.py-f30im_dirout.mkvFFMPEG命令行转换一组JPG图像到视频时,是将这组图像视为MJPG流。我需要转换一组PNG图像到视频,FFMPEG就不认了。pyav内置了ffmpeg库,不需要系统带有ffmpeg工具因此我使用ffmpeg的python包装p

  10. Python 刷Leetcode题库,顺带学英语单词(31) - 2

    ValidPalindromeGivenastring,determineifitisapalindrome,consideringonlyalphanumericcharactersandignoringcases. [#125]Example:"Aman,aplan,acanal:Panama"isapalindrome."raceacar"isnotapalindrome.Haveyouconsiderthatthestringmightbeempty?Thisisagoodquestiontoaskduringaninterview.Forthepurposeofthisproblem

随机推荐