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()
