Gradiente di dati rumorosi, python

Question

Ho uno spettro di energia da un rivelatore di raggi cosmici. Lo spettro segue una curva esponenziale ma avrà grumi ampi (e forse molto lievi). I dati, ovviamente, contengono un elemento di rumore.

Sto cercando di appianare i dati e quindi tracciare il gradiente. Finora ho usato la funzione scipy per lisciarlo e poi np.gradient().

Come si può vedere dall'immagine, il metodo della funzione gradiente è trovare le differenze tra ogni punto e non mostra i grumi in modo molto chiaro.

Fondamentalmente ho bisogno di un grafico gradiente uniforme. Qualsiasi aiuto sarebbe fantastico!

che ho provato metodi 2 spline:

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 

Edit: Provato differenziazione numerica:

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 

Tuttavia, appena ritorna di nuovo nello stesso grafico!

Answer 1

C'è un interessante metodo pubblicato su questo: Numerical Differentiation of Noisy Data. Dovrebbe darti una bella soluzione al tuo problema. Ulteriori dettagli sono forniti in un altro, accompanying paper. L'autore dà anche Matlab code that implements it; è disponibile anche un'alternativa implementation in Python.

Se si vuole perseguire la interpolazione con spline metodo, vorrei suggerire di regolare il fattore di livellamento s di scipy.interpolate.UnivariateSpline().

Un'altra soluzione sarebbe semplificare la funzione tramite la convoluzione (ad esempio con un gaussiano).

Il documento che ho allegato afferma di impedire alcuni degli artefatti che si presentano con l'approccio della convoluzione (l'approccio spline potrebbe soffrire di difficoltà simili).

Answer 2

Non cercherò la validità matematica di questo; sembra che la carta di LANL che EOL ha citato varrebbe la pena esaminare. Ad ogni modo, ho ottenuto risultati decenti usando la differenziazione integrata di SciPy's splines quando si utilizza splev.

%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()