仅将伽马分布拟合到样本的子集

Fit the gamma distribution only to a subset of the samples

我的输入数据(黑色)的直方图如下图所示:

我正在尝试拟合 Gamma distribution 但不是针对整个数据,而是针对直方图的第一条曲线(第一种模式)。上图中的绿色图对应于我使用以下使用 scipy.stats.gamma:

python 代码在所有样本上安装 Gamma distribution
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]

更新 2 (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])

但正如下面的屏幕截图所示,伽马图不符合直方图:

您可以使用 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 的左侧。关键是 归一化 ,因为改变 alphabeta 会改变原始伽玛中截断点左侧的区域。

剩下的只是优化技术。

使用对数是常见的做法,所以我使用了有时称为“能量”的东西,或概率密度倒数的对数。

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. Select 来自您相当大的数据集的样本: data = np.random.choice(data, 10000)

  2. 使用 method 关键字参数尝试不同的算法。

一些优化例程输出逆海森矩阵的表示,这对不确定性估计很有用。对参数强制执行非负性也可能是个好主意。

没有截断的对数标度图显示了整个分布:

这是另一种可能的方法,使用 excel 中手动创建的数据集,或多或少与给定的图相匹配。

原始数据

大纲

  • 将数据导入 Pandas 数据框。
  • 屏蔽后的索引 最大反应指数
  • 创建剩余数据的镜像。
  • 附加镜像,同时留下空缓冲区 space。
  • 将所需的分布拟合到修改后的数据。下面我用矩量的方法做了一个正常的拟合,调整了幅值和宽度。

工作脚本

    # 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 fitting data 页面了解更多关于拟合正常使用矩的方法。