3

Scipyを使用して対数正規分布を近似しようとしています。以前にMatlabを使用してこれを実行しましたが、統計分析を超えてアプリケーションを拡張する必要があるため、Scipyで近似値を再現しようとしています。

以下は、データを適合させるために使用したMatlabコードです。

% Read input data (one value per line)
x = [];
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.')
else
    while ~feof(fid)
        [x] = [x fscanf(fid, '%f', [1])];

    end
    c = fclose(fid);
    if c == 0
         disp('File closed successfully.');
    else
        disp('[ERROR] There was a problem with closing the file.');
    end
end

[f,xx] = ecdf(x);
y = 1-f;

parmhat  = 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コードを使用した適合は次のとおりです。

ここに画像の説明を入力してください

ご覧のとおり、両方のコードセットは、少なくとも視覚的には、適切に適合させることができます。

問題は、推定されたパラメータmuとsigmaに大きな違いがあることです。

Matlabから：mu = 1.62 sigma=1.29。Pythonから：mu = 2.78 sigma=1.74。

なぜそのような違いがあるのですか？

注：適合した両方のデータセットが完全に同じであることを再確認しました。同じ数のポイント、同じ分布。

あなたの助けは大歓迎です！前もって感謝します。

他の情報：

import scipy
import numpy
import statsmodels

scipy.__version__
'0.9.0'

numpy.__version__
'1.6.1'

statsmodels.__version__
'0.5.0.dev-1bbd4ca'

MatlabのバージョンはR2011bです。

版：

以下の回答に示されているように、障害はScipy0.9にあります。Scipy 11.0を使用して、Matlabからのmuとsigmaの結果を再現できます。

Scipyを更新する簡単な方法は次のとおりです。

pip install --upgrade Scipy

あなたがピップを持っていない場合（あなたはすべきです！）：

sudo apt-get install pip
4

1 に答える 1

6

scipy 0.9.0のメソッドにバグがありfit、それ以降のバージョンのscipyで修正されています。

以下のスクリプトの出力は次のようになります。

Explicit formula:   mu = 4.99203450, sig = 0.81691086
Fit log(x) to norm: mu = 4.99203450, sig = 0.81691086
Fit x to lognorm:   mu = 4.99203468, sig = 0.81691081

しかし、scipy 0.9.0では、

Explicit formula:   mu = 4.99203450, sig = 0.81691086
Fit log(x) to norm: mu = 4.99203450, sig = 0.81691086
Fit x to lognorm:   mu = 4.23197270, sig = 1.11581240

次のテストスクリプトは、同じ結果を得る3つの方法を示しています。

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=1stdのオプション付き。開発者 不偏分散推定を使用するための計算：

In [104]: x
Out[104]: array([  50.,   50.,  100.,  200.,  200.,  300.,  500.])

In [105]: lognfit(x, ddof=1)
Out[105]: (4.9920345004312647, 0.88236457185021866)

matlabのlognfitドキュメントには、打ち切りが使用されていない場合、lognfitは分散の偏りのない推定量の平方根を使用してシグマを計算するという注記があります。これは、上記のコードでddof=1を使用することに対応します。

于 2013-03-26T08:59:35.280 に答える