Fitting lognormal distribution using Scipy vs Matlab [拟合的对数正态分布使用SciPy与MATLAB] - 问题-字节技术
Fitting lognormal distribution using Scipy vs Matlab
拟合的对数正态分布使用SciPy与MATLAB
问题 (Question)
I am trying to fit a lognormal distribution using Scipy. I've already done it using Matlab before but because of the need to extend the application beyond statistical analysis, I am in the process of trying to reproduce the fitted values in Scipy.
Below is the Matlab code I used to fit my data:
% Read input data (one value per line)
fid = fopen(file_path, 'r'); % reading is default action for fopen
disp('Reading network degree data...');
if fid == -1
disp('[ERROR] Unable to open data file.')
while ~feof(fid)
[x] = [x fscanf(fid, '%f', [1])];
c = fclose(fid);
disp('File closed successfully.');
disp('[ERROR] There was a problem with closing the file.');
[f,xx] = ecdf(x);
= lognfit(x); % MLE estimate
mu = parmhat(1);
sigma = parmhat(2);
And here's the fitted plot:
Now here's my Python code with the aim of achieving the same:
import math
from scipy import stats
from statsmodels.distributions.empirical_distribution import ECDF
# The same input is read as a list in Python
ecdf_func = ECDF(degrees)
x = ecdf_func.x
ccdf = 1-ecdf_func.y
# Fit data
shape, loc, scale = stats.lognorm.fit(degrees, floc=0)
# Parameters
sigma = shape # standard deviation
mu = math.log(scale) # meanlog of the distribution
fit_ccdf = stats.lognorm.sf(x, [sigma], floc=1, scale=scale)
Here's the fit using the Python code.
As you can see, both sets of code are capable of producing good fits, at least visually speaking.
Problem is that there is a huge difference in the estimated parameters mu and sigma.
From Matlab: mu = 1.62 sigma = 1.29.
From Python: mu = 2.78 sigma = 1.74.
Why is there such a difference?
Note: I have double checked that both sets of data fitted are exactly the same. Same number of points, same distribution.
Your help is much appreciated! Thanks in advance.
Other info:
import scipy
import numpy
import statsmodels
scipy.__version__
numpy.__version__
statsmodels.__version__
'0.5.0.dev-1bbd4ca'
Version of Matlab is R2011b.
As demonstrated in the answer below, the fault lies with Scipy 0.9. I am able to reproduce the mu and sigma results from Matlab using Scipy 11.0.
An easy way to update your Scipy is:
pip install --upgrade Scipy
If you don't have pip (you should!):
sudo apt-get install pip
我试图使用SciPy符合对数正态分布。我已经完成了使用MATLAB之前，但因为需要扩展应用程序超出了统计分析，我试图在SciPy再现过程拟合值。下面是我用我的数据拟合的MATLAB代码：% Read input data (one value per line)
fid = fopen(file_path, 'r'); % reading is default action for fopen
disp('Reading network degree data...');
if fid == -1
disp('[ERROR] Unable to open data file.')
while ~feof(fid)
[x] = [x fscanf(fid, '%f', [1])];
c = fclose(fid);
disp('File closed successfully.');
disp('[ERROR] There was a problem with closing the file.');
[f,xx] = ecdf(x);
= lognfit(x); % MLE estimate
mu = parmhat(1);
sigma = parmhat(2);
这里的装配图：现在这里是我的Python代码以实现同样的目的：import math
from scipy import stats
from statsmodels.distributions.empirical_distribution import ECDF
# The same input is read as a list in Python
ecdf_func = ECDF(degrees)
x = ecdf_func.x
ccdf = 1-ecdf_func.y
# Fit data
shape, loc, scale = stats.lognorm.fit(degrees, floc=0)
# Parameters
sigma = shape # standard deviation
mu = math.log(scale) # meanlog of the distribution
fit_ccdf = stats.lognorm.sf(x, [sigma], floc=1, scale=scale)
这是适合使用Python代码。正如你可以看到，代码集能够很好的配合，至少在视觉上来说。问题是，在估计的参数μ和σ的巨大差异。从MATLAB：μ= 1.62σ= 1.29。从Python：亩= 2.78σ=。为什么会有这样的差别呢？注：我有双重检查，数据拟合两套完全正确同样的。相同数目的点，相同的分布。非常感谢你的帮助！先谢谢了。其他信息：import scipy
import numpy
import statsmodels
scipy.__version__
numpy.__version__
statsmodels.__version__
'0.5.0.dev-1bbd4ca'
matlab版是r2011b。版：如下面的回答表明，错在SciPy 0.9。我能够重现使用SciPy 11μ和σ结果MATLAB。更新你的SciPy的简单方法是：pip install --upgrade Scipy
如果你没有画中画（你！）：sudo apt-get install pip
最佳答案 (Best Answer)
There is a bug in the fit method in scipy 0.9.0 that has been fixed in later versions of scipy.
The output of the script below should be:
Explicit formula:
mu = 4., sig = 0.
Fit log(x) to norm: mu = 4., sig = 0.
Fit x to lognorm:
mu = 4., sig = 0.
but with scipy 0.9.0, it is
Explicit formula:
mu = 4., sig = 0.
Fit log(x) to norm: mu = 4., sig = 0.
Fit x to lognorm:
mu = 4., sig = 1.
The following test script shows three ways to get the same results:
import numpy as np
from scipy import stats
def lognfit(x, ddof=0):
x = np.asarray(x)
logx = np.log(x)
mu = logx.mean()
sig = logx.std(ddof=ddof)
return mu, sig
# A simple data set for easy reproducibility
x = np.array([50., 50, 100, 200, 200, 300, 500])
# Explicit formula
my_mu, my_sig = lognfit(x)
# Fit a normal distribution to log(x)
norm_mu, norm_sig = stats.norm.fit(np.log(x))
# Fit the lognormal distribution
lognorm_sig, _, lognorm_expmu = stats.lognorm.fit(x, floc=0)
print "Explicit formula:
mu = %10.8f, sig = %10.8f" % (my_mu, my_sig)
print "Fit log(x) to norm: mu = %10.8f, sig = %10.8f" % (norm_mu, norm_sig)
print "Fit x to lognorm:
mu = %10.8f, sig = %10.8f" % (np.log(lognorm_expmu), lognorm_sig)
With the option ddof=1 in the std. dev. calculation to use the unbiased variance estimation:
In [104]: x
Out[104]: array([
In [105]: lognfit(x, ddof=1)
Out[105]: (4.66)
There is a note in matlab's
that says when censoring is not used, lognfit computes sigma using the square root of the unbiased estimator of the variance.
This corresponds to using ddof=1 in the above code.
有一个bugfit在SciPy 0.9.0方法已被固定在SciPy以后的版本。下面的脚本的输出应：Explicit formula:
mu = 4., sig = 0.
Fit log(x) to norm: mu = 4., sig = 0.
Fit x to lognorm:
mu = 4., sig = 0.
但SciPy 0.9.0，它是Explicit formula:
mu = 4., sig = 0.
Fit log(x) to norm: mu = 4., sig = 0.
Fit x to lognorm:
mu = 4., sig = 1.
下面的测试脚本显示得到相同的结果的三种方式：import numpy as np
from scipy import stats
def lognfit(x, ddof=0):
x = np.asarray(x)
logx = np.log(x)
mu = logx.mean()
sig = logx.std(ddof=ddof)
return mu, sig
# A simple data set for easy reproducibility
x = np.array([50., 50, 100, 200, 200, 300, 500])
# Explicit formula
my_mu, my_sig = lognfit(x)
# Fit a normal distribution to log(x)
norm_mu, norm_sig = stats.norm.fit(np.log(x))
# Fit the lognormal distribution
lognorm_sig, _, lognorm_expmu = stats.lognorm.fit(x, floc=0)
print "Explicit formula:
mu = %10.8f, sig = %10.8f" % (my_mu, my_sig)
print "Fit log(x) to norm: mu = %10.8f, sig = %10.8f" % (norm_mu, norm_sig)
print "Fit x to lognorm:
mu = %10.8f, sig = %10.8f" % (np.log(lognorm_expmu), lognorm_sig)
与选项ddof=1在标准开发使用无偏估计方差的计算：In [104]: x
Out[104]: array([
In [105]: lognfit(x, ddof=1)
Out[105]: (4.66)
在Matlab中有一个说审查时没有使用，计算西格玛使用lognfit方差的无偏估计量的平方根。这对应于使用自由度= 1在上面的代码。
本文翻译自StackoverFlow，英语好的童鞋可直接参考原文：Categories
January 2017
9101112131415
16171819202122
23242526272829
by pluskid, on , in &&&&
快三个月没有写日志了，大概是我开始认真写 blog 来第一次，也是因为发生了一些预料之外的事情，中断了许久，到后来又一直非常非常忙，不过我终于又爬上来冒个泡了，表明我还活着。
第二点要澄清的是，我这里并不是要讲“伪随机”、“真随机”这样的问题，而是关于如何生成服从某个概率分布的随机数（或者说 sample）的问题。比如，你想要从一个服从正态分布的随机变量得到 100 个样本，那么肯定抽到接近其均值的样本的概率要大许多，从而导致抽到的样本很多是集中在那附近的。当然，要解决这个问题，我们通常都假设我们已经有了一个生成 0 到 1 之间均匀分布的随机数的工具，就好像
给我们的结果那样，事实上许多时候我们也并不太关心它们是真随机数还是伪随机数，看起来差不多就行了。
现在再回到我们的问题，看起来似乎是很简单的，按照概率分布的话，只要在概率密度大的地方多抽一些样本不就行了吗？可是具体要怎么做呢？要真动起手来，似乎有不是那么直观了。实际上，这个问题曾经也是困扰了我很久，最近又被人问起，那我们不妨在这里一起来总结一下。为了避免一下子就陷入抽象的公式推导，那就还是从一个简单的具体例子出发好了，假设我们要抽样的概率分布其概率密度函数为
，并且被限制在区间
上，如右上图所示。
好了，假设现在我们要抽 100 个服从这个分布的随机数，直观上来讲，抽出来的接近 3 的数字肯定要比接近 0 的数字要多。那究竟要怎样抽才能得到这样的结果呢？由于我们实际上是不能控制最原始的随机数生成过程的，我们只能得到一组均匀分布的随机数，而这组随机数的生成过程对于我们完全是透明的，所以，我们能做的只有把这组均匀分布的随机数做一些变换让他符合我们的需求。找到下手的点了，可是究竟要怎样变换呢？有一个变换相信大家都是很熟悉的，假设我们有一组
之间的均匀分布的随机数
的话， 就是一组在
之间均匀分布的随机数了，不难想象， 等于某个数
的概率就是
的概率（“等于某个数的概率”这种说法对于连续型随机变量来说其实是不合适的，不过大概可以理解所表达的意思啦）。似乎有一种可以“逆转回去”的感觉了。
于是让我们来考虑更一般的变换。首先，我们知道
的概率密度函数是
，假设现在我们令
，不妨先假定
是严格单调递增的函数，这样我们可以求其逆函数
（也是严格单调递增的）。现在来看变换后的随机变量
会服从一个什么样的分布呢？
这里需要小心，因为这里都是连续型的随机变量，并不像离散型随机变量那样可以说成“等于某个值的概率”，因此我们需要转换为概率分布函数来处理，也就是求一个积分啦：
的概率分布函数为
小于或等于某个特定的值
这件事情是等价于
这件事情的。换句话说， 等于
。于是， 的概率分布函数就可以得到了：
再求导我们就能得到
的概率密度函数：
这样一来，我们就得到了对于一个随机变量进行一个映射
之后得到的随即变量的分布，那么，回到我们刚才的问题，我们想让这个结果分布就是我们所求的，然后再反推得
经过简单的化简就可以得到
。也就是说，把得到的随机数
带入到到函数
中所得到的结果，就是符合我们预期要求的随机数啦！
让我们来验证一下：
#!/usr/bin/python
import numpy as np
import matplotlib.pyplot as plot
X0 = np.random.rand(N)
= np.power(9*X1, 1.0/3)
t = np.arange(0.0, 3.0, 0.01)
plot.plot(t, y, 'r-', linewidth=1)
plot.hist(Y, bins=50, normed=1, facecolor='green', alpha=0.75)
plot.show()
这就没错啦，目的达成啦！让我们来总结一下。问题是这样的，我们有一个服从均匀分布的随机变量
，它的概率密度函数为一个常数
上的分布，那么常数
就直接等于 1 了。现在我们要得到一个随机变量
使其概率密度函数为
，做法就是构造出一个函数
满足（在这里加上了绝对值符号，这是因为
如果不是递增而是递减的话，推导的过程中有一处就需要反过来）
反推过来就是，对目标
的概率密度函数求一个积分（其实就是得到它的概率分布函数 CDF ，如果一开始就拿到的是 CDF 当然更好），然后求其反函数就可以得到需要的变换
了。实际上，这种方法有一个听起来稍微专业一点的名字： 。不过，虽然看起来很简单，但是实际操作起来却比较困难，因为对于许多函数来说，求逆是比较困难的，求积分就更困难了，如果写不出解析解，不得已只能用数值方法来逼近的话，计算效率就很让人担心了。可事实上也是如此，就连我们最常见的一维标准正态分布，也很难用这样的方法来抽样，因为它的概率密度函数
的不定积分没有一个解析形式。这可真是一点也不好玩，费了这么大劲，结果好像什么都干不了。看来这个看似简单的问题似乎还是比较复杂的，不过也不要灰心，至少对于高斯分布来说，我们还有一个叫做
的方法可以专门来做这个事情。因为高斯分布比较奇怪，虽然一维的时候概率分布函数无法写出解析式，但是二维的情况却可以通过一些技巧得出一个解析式来。
首先我们来考虑一个二维的且两个维度相互独立的高斯分布，它的概率密度函数为
这个分布是关于原点对称的，如果考虑使用极坐标
）的话，我们有 ， 这样的变换。这样，概率密度函数是写成：
注意到在给定
的情况下其概率密度是不依赖于
的，也就是说对于
来说是一个均匀分布，这和我们所了解的标准正态分布也是符合的：在一个圆上的点的概率是相等的。确定了
的分布，让我们再来看 ，用类似于前面的方法：
根据前面得出的结论，我现在得到了
的概率分布函数，是不是只要求一下逆就可以得到一个
现在只要把这一些线索串起来，假设我们有两个相互独立的平均分布在
上的随机变量
就可以得到
了（实际上，由于
实际上是相同的分布，所以通常直接写为 ）。再把极坐标换回笛卡尔坐标：
这样我们就能得到一个二维的正态分布的抽样了。可以直观地验证一下，二维不太好画，就画成 heatmap 了，看着比较热的区域就是概率比较大的，程序如下：
#!/usr/bin/python
import numpy as np
import matplotlib.pyplot as plot
T1 = np.random.rand(N)
T2 = np.random.rand(N)
r = np.sqrt(-2*np.log(T2))
theta = 2*np.pi*T1
X = r*np.cos(theta)
Y = r*np.sin(theta)
heatmap, xedges, yedges = np.histogram2d(X, Y, bins=80)
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
plot.imshow(heatmap, extent=extent)
plot.show()
画出来的图像这个样子：
不太好看，但是大概的形状是可以看出来的。其实有了二维的高斯分布，再注意到两个维度在我们这里是相互独立的，那么直接取其中任意一个维度，就是一个一维高斯分布了。如下：
即服从标准正态分布的话，则有
，也就是说，有了标准正态分布，其他所有的正态分布的抽样也都可以完成了。这下总算有点心满意足了。不过别急，还有最后一个问题：多元高斯分布。一般最常用不就是二元吗？二元不是我们一开始就推出来了吗？推出来了确实没错，不过我们考虑的是最简单的情形，当然同样可以通过
这样的方式来处理每一个维度，不过高维的情形还有一个需要考虑的就是各个维度之间的相关性——我们之前处理的都是两个维度相互独立的情况。对于一般的多维正态分布
，如果各个维度之间是相互独立的，就对应于协方差矩阵
是一个对角阵，但是如果
在非对角线的地方存在非零元素的话，就说明对应的两个维度之间存在相关性。
这个问题还是比较好解决的，高斯分布有这样的：类似于一维的情况，对于多维正态分布
，那么新的随机变量
所以，对于一个给定的高斯分布
来说，只要先生成一个对应维度的标准正态分布
即可，其中
的结果，即
结束之前让我们来看看 matlab 画个 3D 图来改善一下心情：
N = 50000;
T1 = rand(1, N);
T2 = rand(1, N);
r = sqrt(-2*log(T2));
theta = 2*pi*T1;
X = [r.*cos(theta); r.*sin(theta)];
mu = [1; 2];
Sigma = [5 2; 2 1];
L = chol(Sigma);
X1 = repmat(mu, 1, N) + L*X;
nbin = 30;
hist3(X1', [nbin nbin]);
set(gcf, 'renderer', 'opengl');
set(get(gca, 'child'), 'FaceColor', 'interp', 'CDataMode', 'auto');
[z c] = hist3(X1', [nbin nbin]);
[x y] = meshgrid(c{1}, c{2});
surfc(x,y,-z);
下面两幅图，哪幅好看一些（注意坐标比例不一样，所以看不出形状和旋转了）？似乎都不太好看，不过感觉还是比前面的 heatmap 要好一点啦！
然后，到这里为止，我们算是把高斯分布弄清楚了，不过这只是给一个介绍性的东西，里面的数学推导也并不严格，而 Box Muller 也并不是最高效的高斯采样的算法，不过，就算我们不打算再深入讨论高斯采样，采样这个问题本身也还有许多不尽人意的地方，我们推导出来的结论可以说只能用于一小部分简单的分布，连高斯分布都要通过 trick 来解决，另一些本身连概率密度函数都写不出来或者有各种奇怪数学特性的分布就更难处理了。所以本文的标题里也说了，这是上篇，如果什么时候有机会抽出时间来写下篇的话，我将会介绍一些更加通用和强大的方法，诸如 Rejection Sampling 、Gibbs Sampling 以及 Markov Chain Monte Carlo (MCMC) 等方法。如果你比较感兴趣，可以先自行 Google 一下解馋！
末了，还得废话两句，虽然是冒了一个泡，但是也是冒得很艰难，发了这篇日志也并不代表 blog 会恢复往日的发文频率，实际上这篇 blog 也差不多凑了半个多月的时间才折腾完，文章开头的“三个月”我想是不是也该更新为“四个月”了。
不过我想，如果有时间和有趣的 idea 的话，我还是会写的，毕竟这也是我的乐趣之一呢！
Leave a Reply
Except where otherwise noted, content on this site is licensed under a .matlab中normfit在正态分布中的使用技巧
我的图书馆
matlab中normfit在正态分布中的使用技巧
matlab中normfit在正态分布中的使用技巧如下:
函数 normfit
格式 [muhat,sigmahat,muci,sigmaci] = normfit(X)
[muhat,sigmahat,muci,sigmaci] = normfit(X,alpha)
说明 muhat,sigmahat分别为正态分布的参数μ和σ的估计值,muci,sigmaci分别为置信区间,其置信度为;alpha给出显著水平α,缺省时默认为0.05,即置信度为95%.
例4-62 有两组(每组100个元素)正态随机数据,其均值为10,均方差为2,求95%的置信区间和参数估计值.
解:&&r = normrnd (10,2,100,2); %产生两列正态随机数据
&&[mu,sigma,muci,sigmaci] = normfit(r)
10.7 %各列的均值的估计值
1.6 %各列的均方差的估计值
说明 muci,sigmaci中各列分别为原随机数据各列估计值的置信区间,置信度为95%.
例4-63 分别使用金球和铂球测定引力常数
(1)用金球测定观察值为:6.683 6.681 6.676 6.678 6.679 6.672
(2)用铂球测定观察值为:6.661 6.661 6.667 6.667 6.664
设测定值总体为,μ和σ为未知.对(1),(2)两种情况分别求μ和σ的置信度为0.9的置信区间.
解:建立M文件:LX0833.m
X=[6.683 6.681 6.676 6.678 6.679 6.672];
Y=[6.661 6.661 6.667 6.667 6.664];
[mu,sigma,muci,sigmaci]=normfit(X,0.1) %金球测定的估计
[MU,SIGMA,MUCI,SIGMACI]=normfit(Y,0.1) %铂球测定的估计
运行后结果显示如下:
由上可知,金球测定的μ估计值为6.6782,置信区间为[6.3];
σ的估计值为0.0039,置信区间为[0.1].
泊球测定的μ估计值为6.6640,置信区间为[6.9];
σ的估计值为0.0030,置信区间为[0.1].
发表评论：
TA的最新馆藏[转]&[转]&matlab 卡方检验用来 分布的拟合性 比如正太 对数 高斯 瑞利 ，韦伯
259万源代码下载-
&文件名称: matlab& & [
& & & & &&]
&&所属分类:
&&开发工具: matlab
&&文件大小: 2 KB
&&上传时间:
&&下载次数: 124
&&提 供 者:
&详细说明：卡方检验用来检验分布的拟合性 比如正太分布对数正太分布高斯分布瑞利分布，韦伯分布等里面包含了对数据的检测以及统计的原理和方法-The chi-square test used to test the The chi-square test used to test the distribution fitting, such as too the distribution of the lognormal distribution is too Gauss points Burui Li distribution, Weibull distribution, which contains the principles and methods of detection and statistical data
文件列表(点击判断是否您需要的文件，如果是垃圾请在下面评价投诉):
&&matlab.doc
&[]:纯粹是垃圾&[]:纯粹是垃圾&[]:纯粹是垃圾&[]:纯粹是垃圾
&近期下载过的用户:
&&&&&&&&&&&&[]
&相关搜索:
&输入关键字，在本站259万海量源码库中尽情搜索：
&[] - 利用Python检验一组数据是否为正态分布 用KS-test
&[] - 本matlab程序对t，f，卡方三大概率分布进行绘图分析，绘制了每种分布不同参数下的曲线，可以直观的看出不同参数对不同分布的影响。
&[] - 使用matlab模拟、检验和估计泊松过程、}