Back to index

python-biopython  1.60
Functions
Bio.Statistics.lowess Namespace Reference

Functions

def lowess
def _test

Function Documentation

def Bio.Statistics.lowess._test ( ) [private]
Run the Bio.Statistics.lowess module's doctests.

Definition at line 94 of file lowess.py.

00094 
00095 def _test():
00096     """Run the Bio.Statistics.lowess module's doctests."""
00097     print "Running doctests..."
00098     import doctest
00099     doctest.testmod()
00100     print "Done"

def Bio.Statistics.lowess.lowess (   x,
  y,
  f = 2./3.,
  iter = 3 
)
lowess(x, y, f=2./3., iter=3) -> yest

Lowess smoother: Robust locally weighted regression.
The lowess function fits a nonparametric regression curve to a scatterplot.
The arrays x and y contain an equal number of elements; each pair
(x[i], y[i]) defines a data point in the scatterplot. The function returns
the estimated (smooth) values of y.

The smoothing span is given by f. A larger value for f will result in a
smoother curve. The number of robustifying iterations is given by iter. The
function will run faster with a smaller number of iterations.

x and y should be numpy float arrays of equal length.  The return value is
also a numpy float array of that length.

e.g.
>>> import numpy
>>> x = numpy.array([4,  4,  7,  7,  8,  9, 10, 10, 10, 11, 11, 12, 12, 12,
...                 12, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 16, 16,
...                 17, 17, 17, 18, 18, 18, 18, 19, 19, 19, 20, 20, 20, 20,
...                 20, 22, 23, 24, 24, 24, 24, 25], numpy.float)
>>> y = numpy.array([2, 10,  4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24,
...                 28, 26, 34, 34, 46, 26, 36, 60, 80, 20, 26, 54, 32, 40,
...                 32, 40, 50, 42, 56, 76, 84, 36, 46, 68, 32, 48, 52, 56,
...                 64, 66, 54, 70, 92, 93, 120, 85], numpy.float)
>>> result = lowess(x, y)
>>> len(result)
50
>>> print "[%0.2f, ..., %0.2f]" % (result[0], result[-1])
[4.85, ..., 84.98]

Definition at line 33 of file lowess.py.

00033 
00034 def lowess(x, y, f=2./3., iter=3):
00035     """lowess(x, y, f=2./3., iter=3) -> yest
00036 
00037     Lowess smoother: Robust locally weighted regression.
00038     The lowess function fits a nonparametric regression curve to a scatterplot.
00039     The arrays x and y contain an equal number of elements; each pair
00040     (x[i], y[i]) defines a data point in the scatterplot. The function returns
00041     the estimated (smooth) values of y.
00042 
00043     The smoothing span is given by f. A larger value for f will result in a
00044     smoother curve. The number of robustifying iterations is given by iter. The
00045     function will run faster with a smaller number of iterations.
00046 
00047     x and y should be numpy float arrays of equal length.  The return value is
00048     also a numpy float array of that length.
00049 
00050     e.g.
00051     >>> import numpy
00052     >>> x = numpy.array([4,  4,  7,  7,  8,  9, 10, 10, 10, 11, 11, 12, 12, 12,
00053     ...                 12, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 16, 16,
00054     ...                 17, 17, 17, 18, 18, 18, 18, 19, 19, 19, 20, 20, 20, 20,
00055     ...                 20, 22, 23, 24, 24, 24, 24, 25], numpy.float)
00056     >>> y = numpy.array([2, 10,  4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24,
00057     ...                 28, 26, 34, 34, 46, 26, 36, 60, 80, 20, 26, 54, 32, 40,
00058     ...                 32, 40, 50, 42, 56, 76, 84, 36, 46, 68, 32, 48, 52, 56,
00059     ...                 64, 66, 54, 70, 92, 93, 120, 85], numpy.float)
00060     >>> result = lowess(x, y)
00061     >>> len(result)
00062     50
00063     >>> print "[%0.2f, ..., %0.2f]" % (result[0], result[-1])
00064     [4.85, ..., 84.98]
00065     """
00066     n = len(x)
00067     r = int(numpy.ceil(f*n))
00068     h = [numpy.sort(abs(x-x[i]))[r] for i in range(n)]
00069     w = numpy.clip(abs(([x]-numpy.transpose([x]))/h),0.0,1.0)
00070     w = 1-w*w*w
00071     w = w*w*w
00072     yest = numpy.zeros(n)
00073     delta = numpy.ones(n)
00074     for iteration in range(iter):
00075         for i in xrange(n):
00076             weights = delta * w[:,i]
00077             weights_mul_x = weights * x
00078             b1 = numpy.dot(weights,y)
00079             b2 = numpy.dot(weights_mul_x,y)
00080             A11 = sum(weights)
00081             A12 = sum(weights_mul_x)
00082             A21 = A12
00083             A22 = numpy.dot(weights_mul_x,x)
00084             determinant = A11*A22 - A12*A21
00085             beta1 = (A22*b1-A12*b2) / determinant
00086             beta2 = (A11*b2-A21*b1) / determinant
00087             yest[i] = beta1 + beta2*x[i]
00088         residuals = y-yest
00089         s = median(abs(residuals))
00090         delta[:] = numpy.clip(residuals/(6*s),-1,1)
00091         delta[:] = 1-delta*delta
00092         delta[:] = delta*delta
00093     return yest

Here is the call graph for this function:

Here is the caller graph for this function: