0
$\begingroup$

EDIT: I have optimized values of sigma and mu for a lognormal distribution. I want to represent that data as a Gaussian curve, for which the x-axis will be logarithmic. How do I find the normalization constant for my Gaussian and for my lognormal distribution? For the Gaussian distribution (pictured below), I used cnorm = 1 / ( sigma * (2 * pi)^.5 ).

TL;DR:

I am trying to perfect a routine in python on a random-number-generated data sample so that I can apply the "perfected" routine to my actual dataset. I was able to make the routine work smoothly for a Gaussian distribution but not for a lognormal distribution.

The code uses an initial guess value of mu and sigma to calculate an initial Chi-Square value. The parameters (mu and sigma) that produce the minimized Chi-Square are optimized. I was able to optimize these parameters for a lognormal distribution and a Gaussian distribution. By converting the x-axis of a lognormal distribution to a log(x)-axis, the lognormal distribution should look more Gaussian. However, I cannot figure out how to find the proper normalization constant for this new Gaussian representation of the data.

I have 2 sample plots. The one on the left (in the picture) is for the case of a random-Gaussian-number-generator; this one works properly. The one on the right (in the picture) is from running the entire code. The histogram data is normalized but the normalization constant for my lognormal Gaussian fit is clearly wrong. How can I find the proper normalization constant when displaying a normal representation of a lognormal distribution?

That's the TL;DR. The full problem and code is explained in-depth below.

""" IMPORTS """

## GENERATE DATA SAMPLE (( TEST CASE ))
from numpy.random import lognormal as randomlgn
## MATH SYNTAX
from math import exp
from math import log
from math import pi
## INTEGRATATION ROUTINE
from scipy.integrate import quad
## OPTIMIZATION ROUTINES
from scipy.optimize import minimize
from scipy.optimize import basinhopping
## CALCULATE CHISQ ++ PROBABILITY
from scipy.stats import chisquare
## NORMALIZE DISTRIBUTION
from scipy.stats import norm
## PLOT
import numpy as np
import matplotlib.pyplot as plt

""" GENERATE DATA """

mu_true, sigma_true = 48, 7
randlog = randomlgn( mu_true , sigma_true , size = 25000)
datatrue = sorted(randlog)
data = [log(val) for val in datatrue]

""" DEFINE x-range """

small = 1 # min(data) ## smallest value in data
big = 100 # max(data) ## largest value in data
domain = np.linspace(small,big,10000) ## 10000 steps

I then defined functions to help me count the number of observations per bin and to help me tick the midpoints of the bins along the x-axis.

""" DEFINE BIN FUNCTIONS ++ ORGANIZE BINS """

numbins = 40 ## number of bins

def binbounder( small , big , numbins ):
    ## generates list of bound bins
    ## used for histogram ++ bincount
    binwide = ( big - small ) / numbins
    binleft = []
    for index in range( numbins ):
        binleft.append( small + index * binwide )
    binbound = [val for val in binleft]
    binbound.append( big )
    return binbound

def countperbin( binbound , values ):
    ## calculates multiplicity of observed values per bin
    values = sorted( values )
    bincount = []
    for jndex in range( len( binbound ) ):
        if jndex != len( binbound ) - 1:
            summ = 0
            for val in values:
                if val > binbound[ jndex ] and val <= binbound[ jndex + 1 ]:
                    summ += 1
            bincount.append( summ )
        if jndex == len( binbound ) - 1:
            pass
    return bincount

binborders = binbounder( small , big , numbins )
obsperbin = countperbin( binborders , data )
numobs  = sum( obsperbin )
maxbin = max(obsperbin)

To calculate the expectation values per bin, one must find the area under the curve, per bin, from x1 = left bin edge to x2 = right bin edge. The area under the curve is the same as the integral of the distribution function from the left bin edge to the right bin edge. Luckily, the bin edges are defined above to simplify the integrations per bin.

""" DEFINE DISTRIBUTION ++ EXPECTATION VALUE CALCULATOR """

def distribution( x , mu , sigma ): ## LOGNORMAL
    sigma = abs(sigma)
    cnorm = 1 / ( sigma * ( 2*pi )**(1/2) )
    return cnorm / x * exp( (-1) * ( log(x) - mu )**2 / ( 2 * (sigma **2) ) )

def expectperbin( args ):
    ## args[0] = mu
    ## args[1] = sigma
    ## calculates expectation values per bin
    ## expectation value of single bin is equal to area under Gaussian
    ## from left bin-edge to right bin-edge multiplied by the sample size
    xpct = []
    for i in range(len(binborders)-1): # ith i does not exist for rightmost boundary
        xpct.append( quad( distribution , binborders[ i ] , binborders[ i + 1 ],
                         args = ( args[0] , args[1] ))[0] * numobs )
    return xpct

Using the expected values per bin and the observed multiplicities per bin, one can now calculate the Chi-Squared value of the fit. The goal is to minimize Chi-Squared. Two minimization methods are posted, though "Basin Hopping" takes much longer and is unnecessary. The minimization methods require an initial guess of Chi-Squared, which in turn depends on an initial guess of mu and sigma. The code should spit out the optimized parameters (mu and sigma) that yield the minimized value of Chi-Squared.

""" DEFINE CHI SQUARE ++ MINIMIZATION ROUTINE """

def chisq( args ):
    ## first subscript [0] gives chi-squared value, [1] gives 0 = p-value = 1
    return chisquare( obsperbin , expectperbin( args ))[0]

def miniz( chisq , parameterguess, ):
    ## MINIMIZATION ROUTINES    --    https://docs.scipy.org/doc/scipy-0.14.0/reference/optimize.html
    ##      via SCIPY

    """ ROUTINE 1:    BASIN HOP """
    # globmin = basinhopping( chisq , parameterguess , niter = 200 )
    # while globmin.fun <= 0:
    #     ## self-correcting mechanism if ChiSq is negative
    #     ## GLOBAL MIN ==> DOES NOT WORK WITH 'SUCCESS TEST' of SCIPY MODULE
    #
    #     try:
    #         globmin = basinhopping( chisq , parameterguess , niter = 200 )
    #         print("ERROR:   BASINHOP LOOPING AGAIN")
    #         break
    #     except globmin.fun > 0:
    #         print("TA DAA")
    #         break
    #
    # return globmin

    """ ROUTINE 2:    MINIMIZE (SCALAR MULTIVARIABLE FUNCTION) """
    globmin = minimize( chisq , parameterguess)
    while globmin.success == False:
        ## self-correcting mechanism if 'success test' fails

        try:
            globmin = minimize( chisq , parameterguess)
            print("ERROR:   MINIMIZE LOOPING AGAIN")
            break
        except globmin.success == True:
            print("TA DAA")
            break

    return globmin

""" PICK INITIAL PARAMETERS TO PICK GUESS VALUE """

## INITIAL GUESS OF mu ++ sigma (ex:    100, 20)
## TRUE VALUES: mu_true = 48, sigma_true = 7
## PICK WRONG GUESSES TO TEST STRENGTH OF MINIMIZATION METHOD (limits = ??)
initial_mu, initial_sigma = 100,20

## INITIAL GUESS OF CHI SQ MINIMUM (DEPENDS ON GUESS VALUES OF mu ++ sigma)
chisqguess = chisquare( obsperbin , expectperbin( [initial_mu , initial_sigma] ))[0]

""" OPTIMIZATION """

## CALCULATION
res = miniz( chisq, [initial_mu , initial_sigma] ) ## FULL OPTIMIZED RESULT
mulog = res.x[0] ## OPTIMIZED MU FOR LGN DISTRIBUTION
sigmalog = res.x[1] ## OPTIMIZED SIGMA FOR LGN DISTRIBUTION

The optimized values of mu and sigma must be recalibrated from the lognormal distribution to fit the normal distribution.

## SOURCE >> https://www.mathworks.com/help/stats/lognstat.html
mu_opt = exp( mulog + (sigmalog**2 /2) )  ## OPTIMIZED MU
sigma_opt = ( exp(2 * mulog + sigmalog**2) * (exp( sigmalog**2 ) - 1) )**(1/2) ## OPTIMIZED SIGMA
cnorm = 1 / ( sigma_opt * ( 2*pi )**(1/2) ) ## OPTIMIZED NORMALIZATION CONSTANT

chisqstats = chisquare( obsperbin, expectperbin( [mu_opt , sigma_opt] ) ) ## RECALCULATE CHISQ via (mu_opt , sigma_opt)
chisqmin = chisqstats[0]
pval = chisqstats[1]

numconstraints = 3 ## mu, sigma, numobs
dof = numbins - numconstraints ## degrees of freedom

## CHECK RESULTS
print("""""")
chisqmin = chisquare( obsperbin, expectperbin( [mu_opt , sigma_opt] ) )[0]
chisqred = chisqmin/dof

Below is a bunch of print statements that made/make debugging easier. Everything should be easy to read from the terminal window.

print("                 CHI SQUARE MINIMIZATION")
print("")
print("OPTIMIZED PARAMETERS:")
print("     guess values          --      mu =", initial_mu,",","sigma =", initial_sigma)
print("     actual values         --      mu =", mu_true,",","sigma =", sigma_true)
print("     optimized values      --      mu =", mu_opt,",","sigma =", sigma_opt)
print("")
print("CHISQ")
print("     guess value           --      chisq =", chisqguess)
print("     minimized value       --      chisq =", chisqmin)
print("     probability           --      p-value =", pval)
print("")
print("     number of bins        =", numbins)
print("     number of constraints =", numconstraints, "   ","(mu, sigma, number of observed values)")
print("     degrees of freedom    =", dof)
print("     reduced value         --      chisq_red =", chisqred)
# print("     probability         --      p-value =", predval)


print("")
print("MINIMIZATION SPECS")
print("""
""")

I then want to plot a histogram of the data and my fit to the data that minimized Chi-Squared.

""" DEFINE FITTED DISTRIBUTIONS TO PLOT """

binwide = (big - small) / numbins ## binwidth
numdata = len(data) ## number of observations
cscale = binwide * numdata * cnorm ## rescale from normalization

def distribfitted(x):
    ## PLOT DISTRIBUTION (( NORMALIZED ))
    cnorm = 1 / ( sigma_opt * ( 2*pi )**(1/2) )
    return [(( cnorm * exp( (-1) * (x[index] - mu_opt)**2 / ( 2 * (sigma_opt **2) ) ) )) for index in range(len(x))]

def distribscaled(x):
    ## PLOT DISTRIBUTION (( NOT NORMALIZED ))
    cnorm = 1 / ( sigma_opt * ( 2*pi )**(1/2) )
    cscale = binwide * numdata * cnorm
    return [(( cscale * exp( (-1) * (x[index] - mu_opt)**2 / ( 2 * (sigma_opt **2) ) ) )) for index in range(len(x))]

""" PLOT PROBABILITY HISTOGRAM ++ FITTED OVERLAY """

## >> NORMALIZED

## DISPLAY VALUES IN LEGEND
plt.plot(np.NaN, np.NaN, '-', color='none', label = '$\mu_{opt}=%.3f$, $\mu=%.3f$' %(mu_opt,mu_true))
plt.plot(np.NaN, np.NaN, '-', color='none', label = '$\sigma_{opt}=%.3f$, $\sigma=%.3f$' %(sigma_opt,sigma_true))
plt.plot(np.NaN, np.NaN, '-', color='none', label = '$\chi^2_{opt}=%.3f$' %chisqmin)

## PLOT HISTOGRAM ++ OVERLAY FIT
plt.hist(data, bins = binborders, alpha = 0.25, label = 'Binned Data', normed = True)
plt.plot(domain,distribfitted(domain),'r-',label = 'Gaussian Fit')

## FORMATTING
plt.title('Histogram: Normalized Probability Distribution of Data with Fitted Overlay') ## put equation f(x)=N*exp(-...)
plt.ylabel('Probabilities')
plt.xlabel('ln (x)')
plt.legend(loc='best', fontsize = 'medium')
plt.tight_layout()
plt.grid(True)

plt.show()
$\endgroup$
  • $\begingroup$ While there is a subtle difference in the way your data fits, I'm not sure why you want to need better fitting data or why you are trying to use the $\chi^2$ to create data when there are much simpler and straightforward ways to create a perfect representation of the data. For example, using the percentile point function and inverse distributions as described in this question yield perfect results with much less hassle: stats.stackexchange.com/questions/259654/… $\endgroup$ – Tavrock Mar 6 '17 at 6:51
  • $\begingroup$ I am trying to replicate the plots in this research paper. This particular one is the second plot (on page 2). The researchers calculated Chi Square for their data fit. $\endgroup$ – MPath Mar 6 '17 at 6:56
  • $\begingroup$ Also to clarify, I am not familiar at all with percentile point functions or inverse distributions. $\endgroup$ – MPath Mar 6 '17 at 21:58
  • $\begingroup$ Inverse distributions allow you to create the pdf from the percent of the distribution and the distribution parameters. Percentile point functions allow you to create a uniform percentile that works well with inverse distributions. I'm not sure what the process is for Python to use inverse distributions, but I am pretty sure they have them. Also, looking back at the paper referenced, it seems that they went with an inverse distribution that was "good enough" for their data. As a result, your results may be as good as theirs were. $\endgroup$ – Tavrock Mar 6 '17 at 22:36
  • $\begingroup$ I don't understand the method just yet. But it appears that pdf (via scipy module for python) uses this to calculate the curve that best approximates the probability density function. In any case, I realized that my mu and sigma for the lognormal distribution are not the same mu and sigma as the normal distribution. Pretty obvious in hindsight, but I think it's possible to convert one to the other. $\endgroup$ – MPath Mar 7 '17 at 11:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.