python 计算概率密度、累计分布、逆函数的例子
计算概率分布的相关参数时,一般使用scipy包,常用的函数包括以下几个:
pdf:连续随机分布的概率密度函数
pmf:离散随机分布的概率密度函数
cdf:累计分布函数
百分位函数(累计分布函数的逆函数)
生存函数的逆函数(1-cdf的逆函数)
函数里面不仅能跟一个数据,还能跟一个数组。下面用正态分布举例说明:
>>>importscipy.statsasst >>>st.norm.cdf(0)#标准正态分布在0处的累计分布概率值 0.5 >>>st.norm.cdf([-1,0,1])#标准正态分布分别在-1,0,1处的累计分布概率值 array([0.15865525,0.5,0.84134475]) >>>st.norm.pdf(0)#标准正态分布在0处的概率密度值 0.3989422804014327 >>>st.norm.ppf(0.975)#标准正态分布在0.975处的逆函数值 1.959963984540054 >>>st.norm.lsf(0.975)#标准正态分布在0.025处的生存函数的逆函数值 1.959963984540054
对于非标准正态分布,通过更改参数loc与scale来改变均值与标准差:
>>>st.norm.cdf(0,loc=2,scale=1)#均值为2,标准差为1的正态分布在0处的累计分布概率值 0.022750131948179195
对于其他随机分布,可能更改的参数不一样,具体需要查官方文档。下面我们举一些常用分布的例子:
>>>st.binom.pmf(4,n=100,p=0.05)#参数值n=100,p=0.05的二项分布在4处的概率密度值 0.17814264156968956 >>>st.geom.pmf(4,p=0.05)#参数值p=0.05的几何分布在4处的概率密度值 0.04286875 >>>st.poisson.pmf(2,mu=3)#参数值mu=3的泊松分布在2处的概率密度值 0.22404180765538775 >>>st.chi2.ppf(0.95,df=10)#自由度为10的卡方分布在0.95处的逆函数值 18.307038053275146 >>>st.t.ppf(0.975,df=10)#自由度为10的t分布在0.975处的逆函数值 2.2281388519649385 >>>st.f.ppf(0.95,dfn=2,dfd=12)#自由度为2,12的F分布在0.95处的逆函数值 3.8852938346523933
补充拓展:给定概率密度,生成随机数python实现
实现的方法可以不止一种:
rejectionsampling
invertthecdf
MetropolisAlgorithm(MCMC)
本篇介绍根据累积概率分布函数的逆函数(2:inverttheCDF)生成的方法。
自己的理解不一定正确,有错误望指正。
目标:
已知y=pdf(x),现想由给定的pdf,生成对应分布的x
PDF是概率分布函数,对其积分或者求和可以得到CDF(累积概率分布函数),PDF积分或求和的结果始终为1
步骤(具体解释后面会说):
1、根据pdf得到cdf
2、由cdf得到inverseofthecdf
3、对于给定的均匀分布[0,1),带入inversecdf,得到的结果即是我们需要的x
求cdf逆函数的具体方法:
对于上面的第二步,可以分成两类:
1、当CDF的逆函数好求时,直接根据公式求取,
2、反之当CDF的逆函数不好求时,用数值模拟方法
自己的理解:为什么需要根据cdf的逆去获得x?
原因一:
因为cdf是单调函数因此一定存在逆函数(cdf是s型函数,而pdf则不一定,例如正态分布,不单调,对于给定的y,可能存在两个对应的x,就不可逆)
原因二:
这仅是我自己的直观理解,根据下图所示(左上为pdf,右上为cdf)
由步骤3可知,我们首先生成[0,1)的均匀随机数,此随机数作为cdf的y,去映射到cdf的x(若用cdf的逆函数表示则是由x映射到y),可以参考上图的右上,既然cdf的y是均匀随机的,那么对于cdf中同样范围的x,斜率大的部分将会有更大的机会被映射,因为对应的y范围更大(而y是随即均匀分布的),那么,cdf的斜率也就等同于pdf的值,这正好符合若x的pdf较大,那么有更大的概率出现(即重复很多次后,该x会出现的次数最多)
代码实现——方法一,公式法
importnumpyasnp importmath importrandom importmatplotlib.pyplotasplt importcollections count_dict=dict() bin_count=20 definverseCDF(): """ returnthexvalueinPDF """ uniform_random=random.random() returninverse_cdf(uniform_random) defpdf(x): return2*x #cdf=x^2,其逆函数很好求,因此直接用公式法 definverse_cdf(x): returnmath.sqrt(x) defdraw_pdf(D): globalbin_count D=collections.OrderedDict(sorted(D.items())) plt.bar(range(len(D)),list(D.values()),align='center') #因为映射bin的时候采用的floor操作,因此加上0.5 value_list=[(key+0.5)/bin_countforkeyinD.keys()] plt.xticks(range(len(D)),value_list) plt.xlabel('x',fontsize=5) plt.ylabel('counts',fontsize=5) plt.title('countingbits') plt.show() foriinrange(90000): x=inverseCDF() #用bin去映射,否则不好操作 bin=math.floor(x*bin_count)#type(bin):int count_dict[bin]=count_dict.get(bin,0)+1 draw_pdf(count_dict)
结果:
代码实现——方法二,数值法
数值模拟cdf的关键是创建lookuptable,
table的size越大则结果越真实(即区间划分的个数)
importnumpyasnp importmath importrandom importmatplotlib.pyplotasplt importcollections lookup_table_size=40 CDFlookup_table=np.zeros((lookup_table_size)) count_dict=dict() bin_count=20 definverse_cdf_numerically(y): globallookup_table_size globalCDFlookup_table value=0.0 foriinrange(lookup_table_size): x=i*1.0/(lookup_table_size-1) value+=pdf2(x) CDFlookup_table[i]=value CDFlookup_table/=value#normalizethecdf ify=y: index=j break #linearinterpolation t=(y-CDFlookup_table[index-1])/\ (CDFlookup_table[index]-CDFlookup_table[index-1]) fractional_index=index+t#因为index从0开始,所以不是(index-1)+t returnfractional_index/lookup_table_size definverseCDF(): """ returnthexvalueinPDF """ uniform_random=random.random() returninverse_cdf_numerically(uniform_random) defpdf2(x): return(x*x*x-10.0*x*x+5.0*x+11.0)/(10.417) defdraw_pdf(D): globalbin_count D=collections.OrderedDict(sorted(D.items())) plt.bar(range(len(D)),list(D.values()),align='center') value_list=[(key+0.5)/bin_countforkeyinD.keys()] plt.xticks(range(len(D)),value_list) plt.xlabel('x',fontsize=5) plt.ylabel('counts',fontsize=5) plt.title('countingbits') plt.show() foriinrange(90000): x=inverseCDF() bin=math.floor(x*bin_count)#type(bin):int count_dict[bin]=count_dict.get(bin,0)+1 draw_pdf(count_dict)
真实函数与模拟结果
扩展:生成伯努利、正太分布
importnumpyasnp importmatplotlib.pyplotasplt """ reference: https://blog.demofox.org/2017/07/25/counting-bits-the-normal-distribution/ """ defplot_bar_x(): #thisisforplottingpurpose index=np.arange(counting.shape[0]) plt.bar(index,counting) plt.xlabel('x',fontsize=5) plt.ylabel('counts',fontsize=5) plt.title('countingbits') plt.show() #ifdice_side=2,isbinomialdistribution #ifdice_side>2,ismultinomialdistribution dice_side=2 #ifNbecomeslarger,thenmultinomialdistributionwillmorelikenormaldistribution N=100 counting=np.zeros(((dice_side-1)*N+1)) foriinrange(30000): sum=0 forjinrange(N): dice_result=np.random.randint(0,dice_side) sum+=dice_result counting[sum]+=1 #normalization counting/=np.sum(counting) plot_bar_x()
以上这篇python计算概率密度、累计分布、逆函数的例子就是小编分享给大家的全部内容了,希望能给大家一个参考,也希望大家多多支持毛票票。
声明:本文内容来源于网络,版权归原作者所有,内容由互联网用户自发贡献自行上传,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任。如果您发现有涉嫌版权的内容,欢迎发送邮件至:czq8825#qq.com(发邮件时,请将#更换为@)进行举报,并提供相关证据,一经查实,本站将立刻删除涉嫌侵权内容。