Docstrings make the post look longer than it actually is. Also, my question is regarding a function at the top of a chain of function calls.
The part that works:
I would like to make a contour plot of chi squared values for a specified distribution. I understand the basics of how to make a contour plot work, but am unable to apply the techniques outside of basic examples. The problem may lie in vectorizing my functions. As a sample, consider a sample Gaussian dataset of 1000 points, the average and spread of which are 48 and 7, respectively.
# imports: import numpy as np, import random, from math import pi, from scipy.integrate import quad, from scipy.stats import chisquare, from scipy.optimize import minimize
dataset_gauss = [random.gauss(48, 7) for index in range(1000)]
My functions and variable names are the way they are because my full code takes multiple distributions (gaussian, lognormal)
def equation_gauss(x, a, b):
"""
This function returns the equation for the Gaussian distribution.
"""
cnorm = 1 / (b* (2*pi)**(1/2))
return cnorm * np.exp((-1) * (x - a)**2 / (2* b**2))
Using maximum log-likelihood, my script (not relevant to question, so code not shown) returns params_gauss = [47.972906400237889, 7.0241339595841286].
In order to calculate chi square, one must first make a list of bin bounds. Then one can equate each expectation value to the integral of the distribution equation from a bins left side to its right side for each bin. The observed values of each bin are the number of observed values within that bin. One can calculate chi square by summing the quotients of the square differences of the expected and observed values per bin divided by the expected values.
def get_bins(distribution, num_bins=50):
"""
This function returns a specified number of equally sized bins over
the domain of the distribution.
"""
if distribution == 'gauss':
dataset = dataset_gauss
return np.linspace(min(dataset), max(dataset), num_bins)
def get_binned_expectations(distribution, args):
"""
This function returns the expectation values per bin for a dataset
given by the specified distribution.
"""
if distribution == 'gauss':
dataset = dataset_gauss
func = equation_gauss
num_obs = len(dataset)
bins = get_bins(distribution)
res = []
for idx in range(len(bins)):
if idx != len(bins)-1:
res.append(quad(func, bins[idx] , bins[idx+1], args = (args[0] , args[1]))[0] * num_obs)
return res
def get_binned_observations(distribution):
"""
This function returns the observation values per bin for a dataset
given by the specified distribution.
"""
if distribution == 'gauss':
dataset = dataset_gauss
bins = get_bins(distribution)
bin_count = []
for idx in range(len(bins)):
if idx != len(bins)-1:
summ = 0
for datum in dataset:
if datum > bins[idx] and datum <= bins[idx+1]:
summ += 1
bin_count.append(summ)
if idx == len(bins)-1:
pass
return bin_count
def get_chi_square(distribution, params):
"""
This function returns the chi square value for a specified
distribution.
EX:
distribution : 'gauss', 'lognormal'
params : [a, b] for parameters a and b
'opt' (for optimized parameters)
"""
values_observation = get_binned_observations(distribution)
if params == 'opt':
if distribution == 'gauss':
params = params_gauss
values_expectation = get_binned_expectations(distribution, params)
return chisquare(values_observation, values_expectation)
As a check, let's try:
res = get_chi_square('gauss', params='opt')
print(res)
new_params = [40, 10]
new_res = get_chi_square('gauss', params=new_params)
print(new_res)
>> Power_divergenceResult(statistic=55.465132812431413, pvalue=0.21391356257718666)
>> Power_divergenceResult(statistic=14950.604250041084, pvalue=0.0)
The first value statistic is the chi squared value obtained with the corresponding parameters, while the second value pvalue is the probability of the parameters fit. For my purposes, it's best to call only the first element as print(new_res[0]). (The probabilities aren't very accurate as the degrees of freedom have not been specified).
In order to make a contour plot, my understanding is that I need to generate a grid space via dim-2 arrays. First, I write a function to return a list of numbers for each parameter. This is the function that returns x, y such that X, Y are its meshgrid.
def get_axis_data(param, frac, size):
"""
This function returns a specified number of elements in a range
centered around the value of the inputted parameter. The extrema
of this range are specified as:
param ± param * frac
"""
update = frac * param
return np.linspace(param - update, param + update, size)
My problem:
I know that I can use plt.contourf(X, Y, Z, cmap). But, I do not know how to format get_chi_square to receive the meshgrid-ed parameters as inputs since it calls the scipy module to (efficiently) calculate chi square via a list of optimizable parameters. I've commented out the things that I've tried unsuccessfully.
def get_grid_data(distribution, frac=1/4, size=9, func=get_chi_square, cmap='plasma'):
"""
This function returns the grid values for a contour plot of the
error metric as a function of the parameters of a specified
distribution.
EX:
func: 'chi square', 'maximum log-likelihood' (error metric)
"""
if distribution == 'gauss':
opt_params = params_gauss
a_vals = get_axis_data(opt_params[0], frac, size)
b_vals = get_axis_data(opt_params[1], frac, size)
X, Y = np.meshgrid(a_vals, b_vals)
# func = np.vectorize(func)
# Z = func(distribution, [X, Y])[0]
return X, Y#, Z
X, Y = get_grid_data('gauss')
print("X")
print(X)
print("")
print("Y")
print(Y)
Running the above gives:
X
[[ 35.9796798 38.97798645 41.9762931 44.97459975 47.9729064
50.97121305 53.9695197 56.96782635 59.966133 ]
[ 35.9796798 38.97798645 41.9762931 44.97459975 47.9729064
50.97121305 53.9695197 56.96782635 59.966133 ]
[ 35.9796798 38.97798645 41.9762931 44.97459975 47.9729064
50.97121305 53.9695197 56.96782635 59.966133 ]
[ 35.9796798 38.97798645 41.9762931 44.97459975 47.9729064
50.97121305 53.9695197 56.96782635 59.966133 ]
[ 35.9796798 38.97798645 41.9762931 44.97459975 47.9729064
50.97121305 53.9695197 56.96782635 59.966133 ]
[ 35.9796798 38.97798645 41.9762931 44.97459975 47.9729064
50.97121305 53.9695197 56.96782635 59.966133 ]
[ 35.9796798 38.97798645 41.9762931 44.97459975 47.9729064
50.97121305 53.9695197 56.96782635 59.966133 ]
[ 35.9796798 38.97798645 41.9762931 44.97459975 47.9729064
50.97121305 53.9695197 56.96782635 59.966133 ]
[ 35.9796798 38.97798645 41.9762931 44.97459975 47.9729064
50.97121305 53.9695197 56.96782635 59.966133 ]]
Y
[[ 5.26810047 5.26810047 5.26810047 5.26810047 5.26810047 5.26810047
5.26810047 5.26810047 5.26810047]
[ 5.70710884 5.70710884 5.70710884 5.70710884 5.70710884 5.70710884
5.70710884 5.70710884 5.70710884]
[ 6.14611721 6.14611721 6.14611721 6.14611721 6.14611721 6.14611721
6.14611721 6.14611721 6.14611721]
[ 6.58512559 6.58512559 6.58512559 6.58512559 6.58512559 6.58512559
6.58512559 6.58512559 6.58512559]
[ 7.02413396 7.02413396 7.02413396 7.02413396 7.02413396 7.02413396
7.02413396 7.02413396 7.02413396]
[ 7.46314233 7.46314233 7.46314233 7.46314233 7.46314233 7.46314233
7.46314233 7.46314233 7.46314233]
[ 7.9021507 7.9021507 7.9021507 7.9021507 7.9021507 7.9021507
7.9021507 7.9021507 7.9021507 ]
[ 8.34115908 8.34115908 8.34115908 8.34115908 8.34115908 8.34115908
8.34115908 8.34115908 8.34115908]
[ 8.78016745 8.78016745 8.78016745 8.78016745 8.78016745 8.78016745
8.78016745 8.78016745 8.78016745]]
I would like to print Z in a format identical to that of X or Y from the code above. How can I obtain the chi square function value in this way?
EDIT:
If I change the function get_grid_data to get_grid_params and redefine get_grid_data as is done below, I can generate 81 values of chi square. I think this is a step forward, but I'm not sure about the order of the array elements in res (aka Z above) necessary for the contour plot.
def get_grid_params(distribution, frac, size):
"""
This function returns the grid values for a contour plot of the
error metric as a function of the parameters of a specified
distribution.
EX:
func: 'chi square', 'maximum log-likelihood' (error metric)
"""
if distribution == 'gauss':
opt_params = params_gauss
a_vals = get_axis_data(opt_params[0], frac, size)
b_vals = get_axis_data(opt_params[1], frac, size)
X, Y = np.meshgrid(a_vals, b_vals)
# func = np.vectorize(func)
# Z = func(distribution, [X, Y])
return X, Y
def get_grid_data(distribution, frac=1/4, size=9, func=get_chi_square):
"""
"""
X, Y = get_grid_params(distribution, frac, size)
res = []
for idx in range(len(X)):
for jdx in range(len(Y)):
res.append(func(distribution, [X[idx][jdx], Y[idx][jdx]])[0])
print(res)
get_grid_data('gauss')
This prints
# 81 elements ==> 9x9 grid
[4208765217.1232886, 79756867.433148235, 2102012.2187297232, 77845.812346977109, 4299.2223157168837, 2529.7286507333743, 20486.858965000847, 257923.37090704756, 4854102.2912357552, 93281349.868633255, 3214630.1060019895, 149308.23999474355, 9526.0996064385563, 892.28204593366377, 1078.7222202890009, 6755.3095776326609, 53291.09528539874, 588864.18413363863, 4691132.998034155, 266721.46912966535, 20459.717521392733, 2093.3255539124393, 279.78284725132187, 577.3737260040574, 3111.9705345888774, 17462.38755758019, 125880.4188491786, 450519.22715869371, 40667.241172187212, 5020.7992346344054, 744.8798302729781, 116.9962855442742, 364.63898596547921, 1791.3456214870084, 7916.7426067634342, 40972.313769493878, 76104.092836489493, 10798.249475713539, 2013.1185415524558, 381.52353083113587, 66.126519584745949, 264.93942984225561, 1200.5798834763946, 4482.867919608283, 18107.837200860213, 21572.225934943446, 4551.094178016996, 1136.7099239043926, 253.51850353558262, 54.455759914884304, 218.13425049819415, 897.03841272531849, 2952.9334085022683, 9936.4277408736034, 9337.1516297669732, 2622.2698023608255, 789.26686546629082, 202.78664001629076, 60.365012999827258, 199.40257099587109, 726.84333101567586, 2159.6632005396755, 6339.5377293121628, 5372.7483380962221, 1815.8139713332946, 620.16531689499118, 184.61780691354744, 75.563465535153725, 196.96163816097214, 626.64757117448494, 1701.8233311097256, 4494.3117008380068, 3664.4699687203392, 1400.0096023072927, 527.65588603959168, 182.94718825996048, 96.20249715692033, 204.59025315045054, 566.75361531867895, 1416.8609878368447, 3434.8994517014899]
# reshape as 9x9 shows the order of params is wrong.
