[Numpy-discussion] Upper and lower envelopes
Bevan
bevan07@gmail....
Tue Sep 30 17:35:21 CDT 2008
Anne Archibald <peridot.faceted <at> gmail.com> writes:
>
> 2008/9/30 bevan <bevan07 <at> gmail.com>:
> > Hello,
> >
> > I have some XY data. I would like to generate the equations for an upper
and
> > lower envelope that excludes a percentage of the data points.
> >
> > I would like to define the slope of the envelope line (say 3) and then
have my
> > code find the intercept that fits my requirements (say 5% of data below the
> > lower envelope). This would then give me the equation and I could plot the
> > upper and lower envelopes.
> >
> >
> > I hope this makes sense. Thanks for any help.
>
> For this particular problem - where you know the slope - it's not too
> hard. If the slope is b, and your points are x and y, compute y-b*x,
> then sort that array, and choose the 5th and 95th percentile values.
>
> Anne
>
David and Anne,
Thanks for your help. I first tried your suggestion David but did not get it
working by the time Anne's post arrived. I plan on trying to get the
optimize.fmin to work as i can see some cool uses for it in the future.
I managed to get Anne's suggestion to work relatively quickly (only slowed by
the dullard operating the machine...)
Apologies for the lack of signature in my first post and thanks again.
My code (for what its worth) is as follows:
import numpy as np
import pylab
from scipy import stats,polyval
#If the slope is b, and your points are x and y, compute y-b*x,
#then sort that array, and choose the 5th and 95th percentile values.
def Envelope(x,y,slpe,percntExclude):
ans = y-slpe*x
ans.sort()
intercpt =stats.scoreatpercentile(ans,percntExclude)
return intercpt
slpeUpper = 1.0
slpeLower = 3.0
percntUpper = 95
percntLower = 5
x = np.random.rand(100)
y = x*3+1.5+ np.random.rand(100)
#Linear regression using stats.linregress
#Returns: slope, intercept, r, two-tailed prob, stderr-of-the-estimate
slpe,intercpt,r,tt,stderr=stats.linregress(x,y)
yRegress=polyval([slpe,intercpt],x)
intercptLower = Envelope(x,y,slpeLower,percntLower)
intercptUpper = Envelope(x,y,slpeUpper,percntUpper)
yLower = polyval([slpeLower,intercptLower],x)
yUpper = polyval([slpeUpper,intercptUpper],x)
pylab.figure()
pylab.plot(x,y,'k.')
pylab.plot(x,yRegress,'b-')
pylab.plot(x,yLower,'r--')
pylab.plot(x,yUpper,'r--')
pylab.show()
Bevan Jenkins
More information about the Numpy-discussion
mailing list