仅将伽马分布拟合到样本的子集
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
的左侧。关键是 归一化 ,因为改变 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)
。和原来一样,是最大似然估计。
如果您对收敛速度不满意,您可能想要
Select 来自您相当大的数据集的样本:
data = np.random.choice(data, 10000)
使用 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 页面了解更多关于拟合正常使用矩的方法。
我的输入数据(黑色)的直方图如下图所示:
我正在尝试拟合 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
的左侧。关键是 归一化 ,因为改变 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)
。和原来一样,是最大似然估计。
如果您对收敛速度不满意,您可能想要
Select 来自您相当大的数据集的样本:
data = np.random.choice(data, 10000)
使用
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 页面了解更多关于拟合正常使用矩的方法。