ノイズの多いデータの勾配、python

Question

私は宇宙線検出器からのエネルギースペクトルを持っています。スペクトラムは指数関数的な曲線に従いますが、それには幅広い（そしておそらくごくわずかな）塊があります。明らかに、データにはノイズの要素が含まれています。

私はデータを滑らかにし、グラデーションをプロットしようとしています。 これまでは、scipy sline関数を使って平滑化してから、np.gradient（）を使っていました。

画像からわかるように、グラデーション関数の方法は、各点の違いを見つけることです。塊を非常にはっきりと見せません。

基本的に滑らかな勾配グラフが必要です。どんな助けも素晴らしいだろう！

私が試した2スプライン方法：

def smooth_data(y,x,factor): 
    print "smoothing data by interpolation..." 
    xnew=np.linspace(min(x),max(x),factor*len(x)) 
    smoothy=spline(x,y,xnew) 
    return smoothy,xnew 

def smooth2_data(y,x,factor): 
    xnew=np.linspace(min(x),max(x),factor*len(x)) 
    f=interpolate.UnivariateSpline(x,y) 
    g=interpolate.interp1d(x,y) 
    return g(xnew),xnew 

編集：試みた数値微分：

def smooth_data(y,x,factor): 
    print "smoothing data by interpolation..." 
    xnew=np.linspace(min(x),max(x),factor*len(x)) 
    smoothy=spline(x,y,xnew) 
    return smoothy,xnew 

def minim(u,f,k): 
    """"functional to be minimised to find optimum u. f is original, u is approx""" 
    integral1=abs(np.gradient(u)) 
    part1=simps(integral1) 
    part2=simps(u) 
    integral2=abs(part2-f)**2. 
    part3=simps(integral2) 
    F=k*part1+part3 
    return F 


def fit(data_x,data_y,denoising,smooth_fac): 
    smy,xnew=smooth_data(data_y,data_x,smooth_fac) 
    y0,xnnew=smooth_data(smy,xnew,1./smooth_fac) 
    y0=list(y0) 
    data_y=list(data_y) 
    data_fit=fmin(minim, y0, args=(data_y,denoising), maxiter=1000, maxfun=1000) 
    return data_fit 

はしかし、それだけで同じグラフをもう一度返します！

Answer 1

は、この上で公開され興味深い方法あり：Numerical Differentiation of Noisy Data。それはあなたにあなたの問題に対する素晴らしい解決策を与えるはずです。詳細は別のaccompanying paperに記載されています。著者はまたMatlab code that implements itを与える;代わりのimplementation in Pythonも利用可能です。

あなたはスプライン方法で補間を追求したい場合は、私はscipy.interpolate.UnivariateSpline()の平滑化係数sを調整することをお勧め。

もう1つの解決策は、畳み込み（ガウス関数を使用して）を使用して関数を円滑にすることです。

私は畳み込みアプローチを思い付く成果物の一部を防ぐために、クレームにリンクされている紙（スプラインアプローチは、同様の困難に苦しむかもしれません）。

Answer 2

私はこれの数学的妥当性を保証しません。 LANLからEOLが引用した論文は、調べる価値があるようです。とにかく、splevを使用しているとき、SciPyのスプラインの組み込み差別化を使用して、まともな結果を得ました。あなたのために理にかなってスムージングのどのレベル

%matplotlib inline 
from matplotlib import pyplot as plt 
import numpy as np 
from scipy.interpolate import splrep, splev 

x = np.arange(0,2,0.008) 
data = np.polynomial.polynomial.polyval(x,[0,2,1,-2,-3,2.6,-0.4]) 
noise = np.random.normal(0,0.1,250) 
noisy_data = data + noise 

f = splrep(x,noisy_data,k=5,s=3) 
#plt.plot(x, data, label="raw data") 
#plt.plot(x, noise, label="noise") 
plt.plot(x, noisy_data, label="noisy data") 
plt.plot(x, splev(x,f), label="fitted") 
plt.plot(x, splev(x,f,der=1)/10, label="1st derivative") 
#plt.plot(x, splev(x,f,der=2)/100, label="2nd derivative") 
plt.hlines(0,0,2) 
plt.legend(loc=0) 
plt.show()