[frames] | no frames]

# Source Code for Module Bio.Statistics.lowess

```  1  # Copyright 2004-2008 by M de Hoon.
3  # This code is part of the Biopython distribution and governed by its
5  # as part of this package.
6  """Implements the Lowess function for nonparametric regression.
7
8  Functions:
9  lowess        Fit a smooth nonparametric regression curve to a scatterplot.
10
12
13  William S. Cleveland: "Robust locally weighted regression and smoothing
14  scatterplots", Journal of the American Statistical Association, December 1979,
15  volume 74, number 368, pp. 829-836.
16
17  William S. Cleveland and Susan J. Devlin: "Locally weighted regression: An
18  approach to regression analysis by local fitting", Journal of the American
19  Statistical Association, September 1988, volume 83, number 403, pp. 596-610.
20  """
21
22  from __future__ import print_function
23
24  from Bio._py3k import range
25
26  import numpy
27
28  try:
29      from Bio.Cluster import median
30      # The function median in Bio.Cluster is faster than the function median
31      # in NumPy, as it does not require a full sort.
32  except ImportError as x:
33      # Use the median function in NumPy if Bio.Cluster is not available
34      from numpy import median
35
36
37 -def lowess(x, y, f=2. / 3., iter=3):
38      """lowess(x, y, f=2./3., iter=3) -> yest
39
40      Lowess smoother: Robust locally weighted regression.
41      The lowess function fits a nonparametric regression curve to a scatterplot.
42      The arrays x and y contain an equal number of elements; each pair
43      (x[i], y[i]) defines a data point in the scatterplot. The function returns
44      the estimated (smooth) values of y.
45
46      The smoothing span is given by f. A larger value for f will result in a
47      smoother curve. The number of robustifying iterations is given by iter. The
48      function will run faster with a smaller number of iterations.
49
50      x and y should be numpy float arrays of equal length.  The return value is
51      also a numpy float array of that length.
52
53      e.g.
54      >>> import numpy
55      >>> x = numpy.array([4,  4,  7,  7,  8,  9, 10, 10, 10, 11, 11, 12, 12, 12,
56      ...                 12, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 16, 16,
57      ...                 17, 17, 17, 18, 18, 18, 18, 19, 19, 19, 20, 20, 20, 20,
58      ...                 20, 22, 23, 24, 24, 24, 24, 25], numpy.float)
59      >>> y = numpy.array([2, 10,  4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24,
60      ...                 28, 26, 34, 34, 46, 26, 36, 60, 80, 20, 26, 54, 32, 40,
61      ...                 32, 40, 50, 42, 56, 76, 84, 36, 46, 68, 32, 48, 52, 56,
62      ...                 64, 66, 54, 70, 92, 93, 120, 85], numpy.float)
63      >>> result = lowess(x, y)
64      >>> len(result)
65      50
66      >>> print("[%0.2f, ..., %0.2f]" % (result[0], result[-1]))
67      [4.85, ..., 84.98]
68      """
69      n = len(x)
70      r = int(numpy.ceil(f * n))
71      h = [numpy.sort(abs(x - x[i]))[r] for i in range(n)]
72      w = numpy.clip(abs(([x] - numpy.transpose([x])) / h), 0.0, 1.0)
73      w = 1 - w * w * w
74      w = w * w * w
75      yest = numpy.zeros(n)
76      delta = numpy.ones(n)
77      for iteration in range(iter):
78          for i in range(n):
79              weights = delta * w[:, i]
80              weights_mul_x = weights * x
81              b1 = numpy.dot(weights, y)
82              b2 = numpy.dot(weights_mul_x, y)
83              A11 = sum(weights)
84              A12 = sum(weights_mul_x)
85              A21 = A12
86              A22 = numpy.dot(weights_mul_x, x)
87              determinant = A11 * A22 - A12 * A21
88              beta1 = (A22 * b1 - A12 * b2) / determinant
89              beta2 = (A11 * b2 - A21 * b1) / determinant
90              yest[i] = beta1 + beta2 * x[i]
91          residuals = y - yest
92          s = median(abs(residuals))
93          delta[:] = numpy.clip(residuals / (6 * s), -1, 1)
94          delta[:] = 1 - delta * delta
95          delta[:] = delta * delta
96      return yest
97
98
99 -def _test():
100      """Run the Bio.Statistics.lowess module's doctests."""
101      print("Running doctests...")
102      import doctest
103      doctest.testmod()
104      print("Done")
105
106  if __name__ == "__main__":
107      _test()
108
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Wed Dec 17 16:12:03 2014 http://epydoc.sourceforge.net