How can I fit multiple Gaussian curved to mass spectrometry data in Python? -
i fit multiple gaussian curves mass spectrometry data in python. right i'm fitting data 1 gaussian @ time -- literally 1 range @ time.
is there more streamlined way this? there way can run data through loop plot gaussian @ each peak? i'm guessing there's gotta better way, i've combed through internet.
my graph 2 gaussians shown below.
my example data can found at: http://txt.do/dooxv
and here's current code:
import numpy np import matplotlib.pyplot plt import scipy.optimize opt scipy.interpolate import interp1d rgadata = np.loadtxt("/users/ilenemitchell/desktop/rgascan.txt", skiprows=14) rgadata=rgadata.transpose() x=rgadata[0] y=rgadata[1] # graph labels plt.ylabel('ion current') plt.xlabel('mass/charge ratio') plt.xticks(np.arange(min(rgadata[0]), max(rgadata[0])+2, 2.0)) plt.ylim([10**-12.5, 10**-9]) plt.title('rga data jul 25, 2017') plt.semilogy(x, y,'b') #fitting guassian peak def gauss(x, a, mu, sig): return a*np.exp(-(x-mu)**2/(2*sig**2)) fitx=x[(x>40)*(x<43)] fity=y[(x>40)*(x<43)] mu=np.sum(fitx*fity)/np.sum(fity) sig=np.sqrt(np.sum(fity*(fitx-mu)**2)/np.sum(fity)) print (mu, sig, max(fity)) popt, pcov = opt.curve_fit(gauss, fitx, fity, p0=[max(fity),mu, sig]) plt.semilogy(x, gauss(x, popt[0],popt[1],popt[2]), 'r-', label='fit') #second guassian fitx2=x[(x>26)*(x<31)] fity2=y[(x>26)*(x<31)] mu=np.sum(fitx2*fity2)/np.sum(fity2) sig=np.sqrt(np.sum(fity2*(fitx2-mu)**2)/np.sum(fity2)) print (mu, sig, max(fity2)) popt2, pcov2 = opt.curve_fit(gauss, fitx2, fity2, p0=[max(fity2),mu, sig]) plt.semilogy(x, gauss(x, popt2[0],popt2[1],popt2[2]), 'm', label='fit2') plt.show()
in addition alex f's answer, need identify peaks , analyze surroundings identify xmin , xmax values.
if have done that, can use refactored code , loop within plot relevant data
import numpy np import matplotlib.pyplot plt import scipy.optimize opt scipy.interpolate import interp1d def _gauss(x, a, mu, sig): return a*np.exp(-(x-mu)**2/(2*sig**2)) def gauss(x, y, xmin, xmax): fitx = x[(x>xmin)*(x<xmax)] fity = y[(x>xmin)*(x<xmax)] mu = np.sum(fitx*fity)/np.sum(fity) sig = np.sqrt(np.sum(fity*(fitx-mu)**2)/np.sum(fity)) print (mu, sig, max(fity)) popt, pcov = opt.curve_fit(_gauss, fitx, fity, p0=[max(fity), mu, sig]) return _gauss(x, popt[0], popt[1], popt[2]) # load data , define x - y rgadata = np.loadtxt("/users/ilenemitchell/desktop/rgascan.txt", skiprows=14) x, y = rgadata.t # create plot fig, ax = plt.subplots() ax.semilogy(x, y, 'b') # plot gaussian's between xmin , xmax xmin, xmax in [(40, 43), (26, 31)]: yg = gauss(x, y, xmin, xmax) ax.semilogy(x, yg) # prettify graph ax.set_xlabel("mass/charge ratio") ax.set_ylabel("ion current") ax.set_xticks(np.arange(min(x), max(x)+2, 2.0)) ax.set_ylim([10**-12.5, 10**-9]) ax.set_title("rga data jul 25, 2017") plt.show() 
Comments
Post a Comment