Source Code for Module Biskit.Statistics.lognormal

  1  ## 
  2  ## Biskit, a toolkit for the manipulation of macromolecular structures 
  3  ## Copyright (C) 2004 Michael Habeck 
  4  ## Copyright (C) 2005 Raik Gruenberg 
  5  ## 
  6  ## This program is free software; you can redistribute it and/or 
  7  ## modify it under the terms of the GNU General Public License as 
  8  ## published by the Free Software Foundation; either version 2 of the 
  9  ## License, or any later version. 
 10  ## 
 11  ## This program is distributed in the hope that it will be useful, 
 12  ## but WITHOUT ANY WARRANTY; without even the implied warranty of 
 13  ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 
 14  ## General Public License for more details. 
 15  ## 
 16  ## You find a copy of the GNU General Public License in the file 
 17  ## license.txt along with this program; if not, write to the Free 
 18  ## Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 
 19   
 20  ## Contributions: Michael Habeck, Raik Gruenberg 
 21  ## $Revision: 2.4 $ 
 22  ## last $Date: 2006/09/12 06:34:36 $ 
 23  ## last $Author: leckner $ 
 24   
 25  """ 
 26  lognormal distribution 
 27  """ 
 28   
 29  import Numeric as N 
 30  import RandomArray as R 
 31   
 32   
 33 -def rand_log_normal(alpha, beta, shape): 
 34      return N.exp(R.normal(alpha, beta, shape)) 
 35   
 36   
 37 -def ln(r, alpha, beta): 
 38      return N.exp(-0.5/beta**2 * (N.log(r) - alpha)**2 \ 
 39                   - 0.5*N.log(2*pi)-N.log(beta*r)) 
 40   
 41   
 42 -def erf(x): 
 43      """ 
 44      Approximation to the erf-function with fractional error 
 45      everywhere less than 1.2e-7 
 46   
 47      @param x: value 
 48      @type  x: float 
 49       
 50      @return: value 
 51      @rtype: float 
 52      """ 
 53      if x > 10.: return 1. 
 54      if x < -10.: return -1. 
 55   
 56      z = abs(x) 
 57      t = 1. / (1. + 0.5 * z) 
 58   
 59      r = t * N.exp(-z * z - 1.26551223 + t * (1.00002368 + t * (0.37409196 + \ 
 60          t * (0.09678418 + t * (-0.18628806 + t * (0.27886807 + t * \ 
 61          (-1.13520398 + t * (1.48851587 + t * (-0.82215223 + t * \ 
 62          0.17087277))))))))) 
 63   
 64      if x >= 0.: 
 65          return 1. - r 
 66      else: 
 67          return r - 1. 
 68   
 69   
 70 -def logArea(x, alpha, beta): 
 71      """ 
 72      Area of the smallest interval of a lognormal distribution that still 
 73      includes x. 
 74       
 75      @param x: border value 
 76      @type  x: float 
 77      @param alpha: mean of log-transformed distribution 
 78      @type  alpha: float 
 79      @param beta: standarddev of log-transformed distribution 
 80      @type  beta: float 
 81       
 82      @return: probability that x is NOT drawn from the given distribution 
 83      @rtype: float 
 84      """ 
 85      r_max = N.exp(alpha - beta**2) 
 86   
 87      if x < r_max: x = r_max**2 / x 
 88   
 89      upper = (N.log(x) - alpha) / beta  
 90   
 91      return 0.5 * (erf(upper / N.sqrt(2)) - erf(-(upper + 2*beta) / N.sqrt(2))) 
 92   
 93   
 94 -def logMean( alpha, beta ): 
 95      """ 
 96      @param alpha: mean of log-transformed distribution 
 97      @type  alpha: float 
 98      @param beta: standarddev of log-transformed distribution 
 99      @type  beta: float 
100       
101      @return: mean of the original lognormal distribution 
102      @rtype: float 
103      """ 
104      return N.exp( alpha + (beta**2)/2. ) 
105   
106   
107 -def logSigma( alpha, beta ): 
108      """ 
109      @param alpha: mean of log-transformed distribution 
110      @type  alpha: float 
111      @param beta: standarddev of log-transformed distribution 
112      @type  beta: float 
113       
114      @return: 'standard deviation' of the original lognormal distribution 
115      @rtype: float 
116      """ 
117      return logMean( alpha, beta ) * N.sqrt( N.exp(beta**2) - 1.) 
118   
119   
120 -def logMedian( alpha, beta=None ): 
121      """ 
122      @param alpha: mean of log-transformed distribution 
123      @type  alpha: float 
124      @param beta: not needed 
125      @type  beta: float 
126       
127      @return: median of the original lognormal distribution 
128      @rtype: float 
129      """ 
130      return N.exp( alpha ) 
131   
132   
133 -def logConfidence( x, R, clip=0 ): 
134      """ 
135      Estimate the probability of x NOT beeing a random observation from a 
136      lognormal distribution that is described by a set of random values. 
137       
138      @param x: observed value 
139      @type  x: float 
140      @param R: sample of random values 
141      @type  R: [float] 
142      @param clip: clip zeros at this value  0->don't clip (default: 0) 
143      @type  clip: float 
144       
145      @return: confidence that x is not random, median of random distr. 
146      @rtype: (float, float) 
147      """ 
148      if clip and 0 in R: 
149          R = N.clip( R, clip, max( R ) ) 
150      if clip and x == 0: 
151          x = clip 
152   
153      ## remove 0 instead of clipping 
154      R = N.compress( R, R ) 
155      if x == 0: 
156          return 0, 0 
157   
158      ## get mean and stdv of log-transformed random sample 
159      alpha = N.average( N.log( R ) ) 
160   
161      n = len( R ) 
162   
163      beta = N.sqrt(N.sum(N.power(N.log( R ) - alpha, 2)) / (n - 1.)) 
164   
165      return logArea( x, alpha, beta ), logMedian( alpha ) 
166   
167   
168   
169  ############# 
170  ##  TESTING         
171  ############# 
172           
173 -class Test: 
174      """ 
175      Test class 
176      """ 
177       
178 -    def run( self, local=0 ): 
179          """ 
180          run function test 
181           
182          @param local: transfer local variables to global and perform 
183                        other tasks only when run locally 
184          @type  local: 1|0 
185           
186          @return: 1 
187          @rtype: int 
188          """ 
189          import random 
190          import Biskit.gnuplot as gnuplot 
191          import Biskit.hist as H 
192   
193          cr = [] 
194          for i in range( 10000 ): 
195              ## Some random values drawn from the same lognormal distribution  
196   
197              alpha = 1.5 
198              beta = .7 
199              x = 10. 
200   
201              R = [ random.lognormvariate( alpha, beta ) for i in range( 10 ) ] 
202   
203              cr += [ logConfidence( x, R )[0] ] 
204   
205   
206          ca = logArea( x, alpha, beta ) 
207   
208          if local: 
209              gnuplot.plot( H.density( N.array(cr) - ca, 100 ) ) 
210               
211              globals().update( locals() ) 
212                                 
213          return ca 
214   
215   
216   
217 -    def expected_result( self ): 
218          """ 
219          Precalculated result to check for consistent performance. 
220   
221          @return: 1 
222          @rtype:  int 
223          """ 
224          return 0.86877651432955771 
225       
226   
227  if __name__ == '__main__': 
228   
229      test = Test() 
230   
231      assert abs( test.run( local=1 ) - test.expected_result() ) < 1e-8 
232