python - ピークの半値全幅を見つける

Question

青色のピークの半値全幅 (FWHM) を把握しようとしています (画像を参照)。緑のピークとマゼンタのピークを合わせて、青のピークを作ります。緑とマゼンタのピークの FWHM を求めるために、次の式を使用しています。fwhm = 2*np.sqrt(2*(math.log(2)))*sdここで、sd = 標準偏差。緑とマゼンタのピークを作成し、標準偏差を知っているので、その方程式を使用できます。

次のコードを使用して、緑とマゼンタのピークを作成しました。

def make_norm_dist(self, x, mean, sd):
    import numpy as np

    norm = []
    for i in range(x.size):
        norm += [1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x[i] - mean)**2/(2*sd**2))]
    return np.array(norm) 

青色のピークが 2 つのピークで構成されていることを知らず、データに青色のピークしかない場合、どのように FWHM を見つけることができますか?

私はこのコードを使用してピークトップを見つけました：

peak_top = 0.0e-1000
for i in x_axis:
    if i > peak_top:
        peak_top = i

を 2 で割ってpeak_top高さの半分を見つけ、高さの半分に対応する y 値を見つけようとすることもできますが、高さの半分に正確に一致する x 値がない場合、問題が発生します。

私が試しているものには、よりエレガントな解決策があると確信しています。

Answer 1

スプラインを使用して [青い曲線 - ピーク/2] に合わせてから、根を見つけることができます。

import numpy as np
from scipy.interpolate import UnivariateSpline

def make_norm_dist(x, mean, sd):
    return 1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x - mean)**2/(2*sd**2))

x = np.linspace(10, 110, 1000)
green = make_norm_dist(x, 50, 10)
pink = make_norm_dist(x, 60, 10)

blue = green + pink   

# create a spline of x and blue-np.max(blue)/2 
spline = UnivariateSpline(x, blue-np.max(blue)/2, s=0)
r1, r2 = spline.roots() # find the roots

import pylab as pl
pl.plot(x, blue)
pl.axvspan(r1, r2, facecolor='g', alpha=0.5)
pl.show()

結果は次のとおりです。

Answer 2

これはiPythonでうまくいきました（素早く汚い、3行に減らすことができます）：

def FWHM(X,Y):
    half_max = max(Y) / 2.
    #find when function crosses line half_max (when sign of diff flips)
    #take the 'derivative' of signum(half_max - Y[])
    d = sign(half_max - array(Y[0:-1])) - sign(half_max - array(Y[1:]))
    #plot(X[0:len(d)],d) #if you are interested
    #find the left and right most indexes
    left_idx = find(d > 0)[0]
    right_idx = find(d < 0)[-1]
    return X[right_idx] - X[left_idx] #return the difference (full width)

解像度をより正確にするためにいくつかの追加を行うことができますが、X 軸に沿って多くのサンプルがあり、データにノイズが多すぎないという制限では、これはうまく機能します。

データがガウスではなく、少しノイズが多い場合でも、うまくいきました（最初と最後の半分の最大値がデータを横切るだけです）。

Answer 3

データにノイズがある場合 (現実の世界では常にノイズがある場合)、より堅牢なソリューションは、ガウス分布をデータに当てはめ、そこから FWHM を抽出することです。

import numpy as np
import scipy.optimize as opt

def gauss(x, p): # p[0]==mean, p[1]==stdev
    return 1.0/(p[1]*np.sqrt(2*np.pi))*np.exp(-(x-p[0])**2/(2*p[1]**2))

# Create some sample data
known_param = np.array([2.0, .7])
xmin,xmax = -1.0, 5.0
N = 1000
X = np.linspace(xmin,xmax,N)
Y = gauss(X, known_param)

# Add some noise
Y += .10*np.random.random(N)

# Renormalize to a proper PDF
Y /= ((xmax-xmin)/N)*Y.sum()

# Fit a guassian
p0 = [0,1] # Inital guess is a normal distribution
errfunc = lambda p, x, y: gauss(x, p) - y # Distance to the target function
p1, success = opt.leastsq(errfunc, p0[:], args=(X, Y))

fit_mu, fit_stdev = p1

FWHM = 2*np.sqrt(2*np.log(2))*fit_stdev
print "FWHM", FWHM

プロットされた画像は、次の方法で生成できます。

from pylab import *
plot(X,Y)
plot(X, gauss(X,p1),lw=3,alpha=.5, color='r')
axvspan(fit_mu-FWHM/2, fit_mu+FWHM/2, facecolor='g', alpha=0.5)
show()

さらに優れた近似では、適合前に特定のしきい値を下回るノイズの多いデータを除外します。

Answer 4

これは、スプライン アプローチを使用した素敵な小さな関数です。

from scipy.interpolate import splrep, sproot, splev

class MultiplePeaks(Exception): pass
class NoPeaksFound(Exception): pass

def fwhm(x, y, k=10):
    """
    Determine full-with-half-maximum of a peaked set of points, x and y.

    Assumes that there is only one peak present in the datasset.  The function
    uses a spline interpolation of order k.
    """

    half_max = amax(y)/2.0
    s = splrep(x, y - half_max, k=k)
    roots = sproot(s)

    if len(roots) > 2:
        raise MultiplePeaks("The dataset appears to have multiple peaks, and "
                "thus the FWHM can't be determined.")
    elif len(roots) < 2:
        raise NoPeaksFound("No proper peaks were found in the data set; likely "
                "the dataset is flat (e.g. all zeros).")
    else:
        return abs(roots[1] - roots[0])

Answer 5

私は、ノイズが多くガウス的ではないデータに対してかなりうまく機能する経験的な解決策を実装しましたhaggis.math.full_width_half_max。使用方法は非常に簡単です。

fwhm = full_width_half_max(x, y)

この関数は堅牢です。要求された補間スキームを使用して、データの最大値と「半分下」のしきい値を超える最も近い点を見つけるだけです。

他の回答のデータを使用した例をいくつか示します。

@HYRYの滑らかなデータ

def make_norm_dist(x, mean, sd):
    return 1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x - mean)**2/(2*sd**2))

x = np.linspace(10, 110, 1000)
green = make_norm_dist(x, 50, 10)
pink = make_norm_dist(x, 60, 10)

blue = green + pink   

# create a spline of x and blue-np.max(blue)/2 
spline = UnivariateSpline(x, blue-np.max(blue)/2, s=0)
r1, r2 = spline.roots() # find the roots

# Compute using my function
fwhm, (x1, y1), (x2, y2) = full_width_half_max(x, blue, return_points=True)

# Print comparison
print('HYRY:', r2 - r1, 'MP:', fwhm)

plt.plot(x, blue)
plt.axvspan(r1, r2, facecolor='g', alpha=0.5)
plt.plot(x1, y1, 'r.')
plt.plot(x2, y2, 'r.')

滑らかなデータの場合、結果はかなり正確です。

HYRY: 26.891157007233254 MP: 26.891193606203814

@Hooked のノイズの多いデータ

def gauss(x, p): # p[0]==mean, p[1]==stdev
    return 1.0/(p[1]*np.sqrt(2*np.pi))*np.exp(-(x-p[0])**2/(2*p[1]**2))

# Create some sample data
known_param = np.array([2.0, .7])
xmin,xmax = -1.0, 5.0
N = 1000
X = np.linspace(xmin,xmax,N)
Y = gauss(X, known_param)

# Add some noise
Y += .10*np.random.random(N)

# Renormalize to a proper PDF
Y /= ((xmax-xmin)/N)*Y.sum()

# Fit a guassian
p0 = [0,1] # Inital guess is a normal distribution
errfunc = lambda p, x, y: gauss(x, p) - y # Distance to the target function
p1, success = opt.leastsq(errfunc, p0[:], args=(X, Y))

fit_mu, fit_stdev = p1

FWHM = 2*np.sqrt(2*np.log(2))*fit_stdev

# Compute using my function
fwhm, (x1, y1), (x2, y2) = full_width_half_max(X, Y, return_points=True)

# Print comparison
print('Hooked:', FWHM, 'MP:', fwhm)

plt.plot(X, Y)
plt.plot(X, gauss(X, p1), lw=3, alpha=.5, color='r')
plt.axvspan(fit_mu - FWHM / 2, fit_mu + FWHM / 2, facecolor='g', alpha=0.5)
plt.plot(x1, y1, 'r.')
plt.plot(x2, y2, 'r.')

ノイズの多いデータ (ベースラインが偏っている) の場合、結果はそれほど一貫していません。

Hooked: 1.9903193212254346 MP: 1.5039676990530118

一方では、ガウス近似はデータに対してあまり最適ではありませんが、他方では、半値しきい値と交差する最も近い点を選択する戦略も最適ではない可能性があります。