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
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
