python - サンプルのサブセットのみにガンマ分布を当てはめる

Question

次のグラフに示す入力データ (黒) のヒストグラムがあります。

Gamma distributionデータ全体ではなく、ヒストグラムの最初の曲線 (最初のモード)に合わせようとしています。前のグラフの緑色のプロットは、を利用する次のコードをGamma distribution使用してすべてのサンプルに を当てはめた場合に対応します。pythonscipy.stats.gamma

img = IO.read(input_file)
data = img.flatten() + abs(np.min(img)) + 1

# calculate dB positive image
img_db = 10 * np.log10(img)
img_db_pos = img_db + abs(np.min(img_db))
data = img_db_pos.flatten() + 1

# data histogram
n, bins, patches = plt.hist(data, 1000, normed=True)

# slice histogram here

# estimation of the parameters of the gamma distribution
fit_alpha, fit_loc, fit_beta = gamma.fit(data, floc=0)
x = np.linspace(0, 100)
y = gamma.pdf(x, fit_alpha, fit_loc, fit_beta)
print '(alpha, beta): (%f, %f)' % (fit_alpha, fit_beta)

# plot estimated model
plt.plot(x, y, linewidth=2, color='g')
plt.show()

このデータの興味深いサブセットのみにフィッティングを制限するにはどうすればよいですか?

Update1 (スライス):

前のヒストグラムの最大値を下回る値のみを保持して入力データをスライスしましたが、結果はあまり説得力がありませんでした。

# slice histogram hereこれは、前のコードのコメントの下に次のコードを挿入することで実現されました。

max_data = bins[np.argmax(n)]
data = data[data < max_data]

Update2 (scipy.optimize.minimize):

以下のコードは、scipy.optimize.minimize()を使用してエネルギー関数を最小化して を見つける方法を示しています(alpha, beta)。

import matplotlib.pyplot as plt
import numpy as np
from geotiff.io import IO
from scipy.stats import gamma
from scipy.optimize import minimize


def truncated_gamma(x, max_data, alpha, beta):
    gammapdf = gamma.pdf(x, alpha, loc=0, scale=beta)
    norm = gamma.cdf(max_data, alpha, loc=0, scale=beta)
    return np.where(x < max_data, gammapdf / norm, 0)


# read image
img = IO.read(input_file)

# calculate dB positive image
img_db = 10 * np.log10(img)
img_db_pos = img_db + abs(np.min(img_db))
data = img_db_pos.flatten() + 1

# data histogram
n, bins = np.histogram(data, 100, normed=True)

# using minimize on a slice data below max of histogram
max_data = bins[np.argmax(n)]
data = data[data < max_data]

data = np.random.choice(data, 1000)
energy = lambda p: -np.sum(np.log(truncated_gamma(data, max_data, *p)))
initial_guess = [np.mean(data), 2.]
o = minimize(energy, initial_guess, method='SLSQP')
fit_alpha, fit_beta = o.x

# plot data histogram and model
x = np.linspace(0, 100)
y = gamma.pdf(x, fit_alpha, 0, fit_beta)
plt.hist(data, 30, normed=True)
plt.plot(x, y, linewidth=2, color='g')
plt.show()

上記のアルゴリズムは のサブセットに対して収束し、dataの出力oは次のとおりです。

x: array([ 16.66912781,   6.88105559])

しかし、下のスクリーンショットでわかるように、ガンマ プロットはヒストグラムに適合しません。

Answer 1

scipy.optimize.minimize目的の関数の切り詰められたバージョンを適合させる などの一般的な最適化ツールを使用して、うまく適合させることができます。

まず、変更された関数:

def truncated_gamma(x, alpha, beta):
    gammapdf = gamma.pdf(x, alpha, loc=0, scale=beta)
    norm = gamma.cdf(max_data, alpha, loc=0, scale=beta)
    return np.where(x<max_data, gammapdf/norm, 0)

これは、ガンマ分布から値を選択し、x < max_data他の場所ではゼロを選択します。とにかく、データはもっぱら左側にあるため、ここではそのnp.where部分は実際には重要ではありませmax_dataん。重要なのは正規化です。変化すると、元のガンマの切り捨てポイントの左側の領域が変更されるためalphaですbeta。

残りは最適化の技術です。

対数を扱うのが一般的であるため、「エネルギー」と呼ばれることもある、確率密度の逆数の対数を使用しました。

energy = lambda p: -np.sum(np.log(truncated_gamma(data, *p)))

最小化:

initial_guess = [np.mean(data), 2.]
o = minimize(energy, initial_guess, method='SLSQP')
fit_alpha, fit_beta = o.x

私の出力は(alpha, beta): (11.595208, 824.712481)です。元のように、これは最尤推定です。

収束率に満足できない場合は、

1. かなり大きなデータセットからサンプルを選択します。 data = np.random.choice(data, 10000)

2. methodキーワード引数を使用して別のアルゴリズムを試してください。

一部の最適化ルーチンは、不確実性の推定に役立つ逆ヘッセ行列の表現を出力します。パラメータの非負性を強制することも良い考えです。

切り捨てを行わない対数スケールのプロットは、分布全体を示しています。

Answer 2

これは、与えられたプロットと多かれ少なかれ一致する Excel で手動で作成されたデータセットを使用する別の可能なアプローチです。

生データ

概要

• Pandas データフレームにデータをインポートしました。
• 最大応答インデックスの後のインデックスをマスクします。
• 残りのデータのミラー イメージを作成します。
• 空のスペースのバッファーを残しながらミラー イメージを追加します。
• 変更されたデータに目的の分布を当てはめます。以下では、モーメント法で通常のフィットを行い、振幅と幅を調整します。

作業スクリプト

    # Import data to dataframe.
    df = pd.read_csv('sample.csv', header=0, index_col=0)
    # Mask indices after index at max Y.
    mask = df.index.values <= df.Y.argmax()
    df = df.loc[mask, :]
    scaled_y = 100*df.Y.values

    # Create new df with mirror image of Y appended.
    sep = 6
    app_zeroes = np.append(scaled_y, np.zeros(sep, dtype=np.float))
    mir_y = np.flipud(scaled_y)
    new_y = np.append(app_zeroes, mir_y)

    # Using Scipy-cookbook to fit a normal by method of moments.
    idxs = np.arange(new_y.size)  # idxs=[0, 1, 2,...,len(data)]
    mid_idxs = idxs.mean() # len(data)/2
    # idxs-mid_idxs is [-53.5, -52.5, ..., 52.5, len(data)/2]
    scaling_param = np.sqrt(np.abs(np.sum((idxs-mid_idxs)**2*new_y)/np.sum(new_y)))

    # adjust amplitude
    fmax = new_y.max()*1.2 # adjusted function max to 120% max y.
    # adjust width
    scaling_param = scaling_param*.7 # adjusted by 70%.
    # Fit normal.
    fit = lambda t: fmax*np.exp(-(t-mid_idxs)**2/(2*scaling_param**2))

    # Plot results.
    plt.plot(new_y, '.')
    plt.plot(fit(idxs), '--')
    plt.show()

結果

モーメント法を使用した法線のフィッティングの詳細については、scipy-cookbook フィッティング データページを参照してください。