## Automatically adapted for numpy.oldnumeric Mar 26, 2007 by alter_code1.py
 
##
## Biskit, a toolkit for the manipulation of macromolecular structures
## Copyright (C) 2004 Michael Habeck
## Copyright (C) 2005 Raik Gruenberg
##
## This program is free software; you can redistribute it and/or
## modify it under the terms of the GNU General Public License as
## published by the Free Software Foundation; either version 3 of the
## License, or any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
## General Public License for more details.
##
## You find a copy of the GNU General Public License in the file
## license.txt along with this program; if not, write to the Free
## Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 
## Contributions: Michael Habeck, Raik Gruenberg
## $Revision: 987 $
## last $Date: 2011-09-03 10:09:26 -0400 (Sat, 03 Sep 2011) $
## last $Author: graik $
 
"""
lognormal distribution
"""
 
import numpy.oldnumeric as N
import numpy.oldnumeric.random_array as R
 
 
def rand_log_normal(alpha, beta, shape):
    return N.exp(R.normal(alpha, beta, shape))
 
 
def ln(r, alpha, beta):
    return N.exp(-0.5/beta**2 * (N.log(r) - alpha)**2 \
                 - 0.5*N.log(2*N.pi)-N.log(beta*r))
 
 
def erf(x):
    """
    Approximation to the erf-function with fractional error
    everywhere less than 1.2e-7
 
    @param x: value
    @type  x: float
 
    @return: value
    @rtype: float
    """
    if x > 10.: return 1.
    if x < -10.: return -1.
 
    z = abs(x)
    t = 1. / (1. + 0.5 * z)
 
    r = t * N.exp(-z * z - 1.26551223 + t * (1.00002368 + t * (0.37409196 + \
                                                               t * (0.09678418 + t * (-0.18628806 + t * (0.27886807 + t * \
                                                                                                         (-1.13520398 + t * (1.48851587 + t * (-0.82215223 + t * \
                                                                                                                                               0.17087277)))))))))
 
    if x >= 0.:
        return 1. - r
    else:
        return r - 1.
 
 
def logArea(x, alpha, beta):
    """
    Area of the smallest interval of a lognormal distribution that still
    includes x.
 
    @param x: border value
    @type  x: float
    @param alpha: mean of log-transformed distribution
    @type  alpha: float
    @param beta: standarddev of log-transformed distribution
    @type  beta: float
 
    @return: probability that x is NOT drawn from the given distribution
    @rtype: float
    """
    r_max = N.exp(alpha - beta**2)
 
    if x < r_max: x = r_max**2 / x
 
    upper = (N.log(x) - alpha) / beta 
 
    return 0.5 * (erf(upper / N.sqrt(2)) - erf(-(upper + 2*beta) / N.sqrt(2)))
 
 
def logMean( alpha, beta ):
    """
    @param alpha: mean of log-transformed distribution
    @type  alpha: float
    @param beta: standarddev of log-transformed distribution
    @type  beta: float
 
    @return: mean of the original lognormal distribution
    @rtype: float
    """
    return N.exp( alpha + (beta**2)/2. )
 
 
def logSigma( alpha, beta ):
    """
    @param alpha: mean of log-transformed distribution
    @type  alpha: float
    @param beta: standarddev of log-transformed distribution
    @type  beta: float
 
    @return: 'standard deviation' of the original lognormal distribution
    @rtype: float
    """
    return logMean( alpha, beta ) * N.sqrt( N.exp(beta**2) - 1.)
 
 
def logMedian( alpha, beta=None ):
    """
    @param alpha: mean of log-transformed distribution
    @type  alpha: float
    @param beta: not needed
    @type  beta: float
 
    @return: median of the original lognormal distribution
    @rtype: float
    """
    return N.exp( alpha )
 
 
def logConfidence( x, R, clip=0 ):
    """
    Estimate the probability of x NOT beeing a random observation from a
    lognormal distribution that is described by a set of random values.
 
    @param x: observed value
    @type  x: float
    @param R: sample of random values
    @type  R: [float]
    @param clip: clip zeros at this value  0->don't clip (default: 0)
    @type  clip: float
 
    @return: confidence that x is not random, median of random distr.
    @rtype: (float, float)
    """
    if clip and 0 in R:
        R = N.clip( R, clip, max( R ) )
    if clip and x == 0:
        x = clip
 
    ## remove 0 instead of clipping
    R = N.compress( R, R )
    if x == 0:
        return 0, 0
 
    ## get mean and stdv of log-transformed random sample
    alpha = N.average( N.log( R ) )
 
    n = len( R )
 
    beta = N.sqrt(N.sum(N.power(N.log( R ) - alpha, 2)) / (n - 1.))
 
    return logArea( x, alpha, beta ), logMedian( alpha )
 
 
 
#############
##  TESTING        
#############
import Biskit.test as BT
 
class Test(BT.BiskitTest):
    """
    Test class
    """
 
    def test_lognormal(self):
        """Statistics.lognormal test"""
        import random
        import Biskit.gnuplot as gnuplot
        import Biskit.hist as H
 
        cr = []
        for i in range( 10000 ):
            ## Some random values drawn from the same lognormal distribution 
 
            alpha = 1.5
            beta = .7
            x = 10.
 
            R = [ random.lognormvariate( alpha, beta ) for j in range( 10 ) ]
 
            cr += [ logConfidence( x, R )[0] ]
 
 
        ca = logArea( x, alpha, beta )
 
        if self.local:
            gnuplot.plot( H.density( N.array(cr) - ca, 100 ) )
 
            globals().update( locals() )
 
        self.assertAlmostEqual( ca,  0.86877651432955771, 7)
 
 
if __name__ == '__main__':
 
    BT.localTest()