[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