[SciPy-User] FW: curve fitting by a sum of gaussian with scipy

[Charles R Harris]

On Thu, Apr 18, 2013 at 9:41 AM, Charles R Harris <charlesr.harris at gmail.com
> wrote:

>
>
> On Thu, Apr 18, 2013 at 9:18 AM, Stéphanie haaaaaaaa <
> flower_des_iles at hotmail.com> wrote:
>
>> Hi charles,
>> Yes, i left out all the distance data points that had zero matches. So, I
>> correct this, and i obtain the graph attached.
>> I've done an interpolation with the following code using python and scipy
>> :
>>
>> #!/usr/bin/env python
>> # -*- coding:Utf-8 -*-
>>
>> import numpy as np
>> import matplotlib.pyplot as plt
>> from matplotlib import pylab
>> from scipy import interpolate
>> from scipy import stats
>>
>> data = np.dtype( [ ('distance',np.int), ('effectif',np.int) ] )
>> enreg = np.loadtxt('tmp.txt', dtype=data, skiprows=True)
>>
>> data=zip(enreg['distance'],enreg['effectif'])
>>
>>
>> xmin = -60
>> xmax=60
>> ymin=0
>> ymax=max(enreg['effectif'])+5
>>
>> ax1=plt.subplot(2,1,1)
>> ax1.set_autoscaley_on(False)
>> ax1.set_ylim([ymin,ymax])
>> ax1.set_xlim([xmin,xmax])
>>
>> plt.plot(enreg['distance'],enreg['effectif'],color='#347C2C',lw=2)
>> plt.fill_between(enreg['distance'],enreg['effectif'],0,color='#4CC417')
>>
>> #tck= interpolate.splrep(enreg['distance'],enreg['effectif'])
>> xnew= np.linspace(enreg['distance'][0],enreg['distance'][-1], 1000)
>> #ynew = interpolate.splev(xnew,tck,der=0)
>> #plt.plot(xnew, ynew, 'r')
>>
>> #kernel
>> #kde1 = stats.gaussian_kde(enreg['distance'])
>> #kdepdf = kde1.evaluate(xnew)
>>
>>
>> #plt.plot(xnew, kdepdf, label='kde', color="r")
>>
>> tck= interpolate.splrep(enreg['distance'],enreg['effectif'])
>> xnew= np.linspace(enreg['distance'][0],enreg['distance'][-1], 1000)
>> ynew = interpolate.splev(xnew,tck,der=0)
>> plt.plot(xnew, ynew, 'r')
>>
>>
>> plt.subplot(2,1,2)
>> ax2 = plt.subplot(212, sharex=ax1)
>>
>> for dist, eff in data: # distance dans mon histogramme + effectif dans
>> chaque classe
>>     for e in np.arange(eff): # np.arrange(3) = array([0, 1, 2])
>>         #print "%d,%d : %d" %(dist,eff,e)
>>         emod=e+1
>>         ax2.scatter(dist,-emod,marker='s',color='#B048B5') # pour qu'il
>> n'y ait pas d'overlappe
>>
>> #plt.xlim(0,100)
>>
>> ax1.xaxis.tick_top()
>> # Set the X Axis label.
>> ax1.set_xlabel('Position from significant crosslink site',fontsize=12)
>> ax1.xaxis.set_label_position('top')
>> ax1.axvline(0, color='black', lw=2)
>> # Set the Y Axis label.
>> ax1.set_ylabel('# of piRNAs',fontsize=12)
>>
>>
>> ax2.axes.get_xaxis().set_visible(False)
>> ax2.axes.get_yaxis().set_visible(False)
>> ax2.axvline(0, color='black', lw=2)
>>
>>
>>
>>
>> plt.setp( ax2.get_xticklabels(), visible=False)
>> plt.setp( ax2.get_yticklabels(), visible=False)
>> plt.show()
>>
>> What do you think about it ?
>>
>>
> Certainly looks better ;) Not my field so I can't say more.
>
> For gaussian kernel density estimation there is scipy.stats.gaussian_kdethat might be worth fooling with, although the notes imply that multimodal
> distributions tend to be over smoothed.
>
> <snip>
>
>
And you might try
pchip<http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.pchip.html#scipy.interpolate.pchip>to
cut down on the ringing.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20130418/8429ea22/attachment.html>