<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Distribution Fitting
Goals
Step1:  2. Inspect the Distribution for a Feature
Step2:  Extract the values of the KDE curve for fitting purposes
Step3:  Let's plot it to make sure it looks like what we've seen above
Step4:  Good! It looks like it matches exactly.
3. Postulate a Fit Function
I'm guessing that this is four Gaussian functions added together.
Define a Gaussian function
Define our custom function made up of Gaussians
Make a guess at our parameters
Use scipy.optimize.curve_fit to optimize the parameters
$$ f_{G}(x) = a \cdot exp\left( -\frac{(x-\mu)^2}{2\sigma^2} \right) $$
Step5:  Here I make my ballpark guesses for the amplitudes, means, and deviations
Step6:  Take a look at what that gives us
Step7:  ... This looks not good, but it's a nice little nucleation point for an optimization routine.
4. Use scipy.stats to fit the function to the distribution
Step8:  Put in these optimized parameters and see what we get
Step9:  Compare to the actual KDE distribution
Step10:  I'd say that's pretty good for a first go
Calculate a $\chi^2$ value for this fit
With the hypothesis that these two distributions are the same, we calculate
Step11:  For 12 degrees of freedom (12 fit parameters), we can look at a $\chi^2$ table to find that we have a $p>0.995$ that this is a good fit to the distribution.
Or use scipy.stats to get the $\chi^2$ and p-value
Step12:  Overlay onto original data
<ASSISTANT_TASK:>
Python Code:
import pandas as pd
data_df = pd.read_csv('raw-data.csv', index_col='eventID')
data_df.head()
Explanation: Distribution Fitting
Goals:
Load up raw data from previous post
Inspect distribution for a feature
Postulate a fit function
Use scipy.stats to fit function to the distribution
1. Load the Raw Data
End of explanation
import matplotlib.pylab as plt
import seaborn as sns
# Show plots in notebook
%matplotlib inline
# Set some styling options
sns.set_style("darkgrid")
sns.set_context("paper", font_scale=1.4)
feature = "AfterInhMATRIX5"
sns.distplot(data_df.query("label == 1")[feature])
Explanation: 2. Inspect the Distribution for a Feature
End of explanation
from sklearn.neighbors import KernelDensity
import numpy as np
import matplotlib.pylab as plt
# Put the data you want to apply the KDE to into an array
data = data_df.query("label == 1")[feature].values[:, np.newaxis]
# Create a KDE object and then fit it to the data
kde = KernelDensity(kernel='gaussian', bandwidth=1400).fit(data)
Explanation: Extract the values of the KDE curve for fitting purposes
End of explanation
X_plot = np.linspace(10000, 100000, 1000)[:, np.newaxis]
# Get log density values for each point on the x-axis
log_dens = kde.score_samples(X_plot)
Y_plot = np.exp(log_dens)
# Plot the two against each other
sns.distplot(data_df.query("label == 1")[feature],
             hist=False, color='black',
             label='Seaborn KDE')
plt.plot(X_plot, Y_plot, '--', color='red', lw=2,
         label='scikit-learn KDE')
plt.legend(loc='best')
Explanation: Let's plot it to make sure it looks like what we've seen above
End of explanation
def gauss_dist(xdata, amp, mean, stddev):
    return (amp * np.exp( np.divide(-1 * np.square(xdata-mean),
                                    (2 * stddev**2))))
# Take four amplitudes, means, and standard deviations
# Compute sum of four Gaussians
def my_fit(xdata,
           a1, a2, a3, a4,
           m1, m2, m3, m4,
           s1, s2, s3, s4):
    exp1 = gauss_dist(xdata, a1, m1, s1)
    exp2 = gauss_dist(xdata, a2, m2, s2)
    exp3 = gauss_dist(xdata, a3, m3, s3)
    exp4 = gauss_dist(xdata, a4, m4, s4)
    return exp1 + exp2 + exp3 + exp4
Explanation: Good! It looks like it matches exactly.
3. Postulate a Fit Function
I'm guessing that this is four Gaussian functions added together.
Define a Gaussian function
Define our custom function made up of Gaussians
Make a guess at our parameters
Use scipy.optimize.curve_fit to optimize the parameters
$$ f_{G}(x) = a \cdot exp\left( -\frac{(x-\mu)^2}{2\sigma^2} \right) $$
End of explanation
p0 = [0.00001, 0.00002, 0.000061, 0.000005,
      31000, 51000, 66000, 83000,
      1000, 1500, 2000, 3000]
Explanation: Here I make my ballpark guesses for the amplitudes, means, and deviations
End of explanation
my_guess = my_fit(X_plot[:, 0], *p0)
plt.plot(X_plot, my_guess, '-')
Explanation: Take a look at what that gives us
End of explanation
from scipy.optimize import curve_fit
popt, pcov = curve_fit(my_fit, X_plot[:, 0], Y_plot, p0)
print popt
Explanation: ... This looks not good, but it's a nice little nucleation point for an optimization routine.
4. Use scipy.stats to fit the function to the distribution
End of explanation
optim_fit = my_fit(X_plot[:, 0], *popt)
# Plot the whole fit
plt.plot(X_plot, optim_fit, '-.', lw=3, color='black')
# Along with the consituent gaussians
for i in range(0,4):
    plt.plot(X_plot,
             gauss_dist(X_plot[:, 0], popt[i], popt[i+4], popt[i+8]),
             '-', lw=2)
Explanation: Put in these optimized parameters and see what we get
End of explanation
plt.plot(X_plot, Y_plot, '-', lw=3, color='black',
         label='scikit-learn KDE')
plt.plot(X_plot, optim_fit, '-', color='red',
         label=r'$\sum_{i=1}^4\ f_G(a_i, \mu_i, \sigma_i)$')
plt.legend(loc='best')
Explanation: Compare to the actual KDE distribution
End of explanation
def calc_chisq(obs, exp):
    chisq = 0.0
    for i in range(0,len(obs)):
        chisq += (obs[i] - exp[i])**2 / exp[i]
    return chisq
print calc_chisq(Y_plot, optim_fit)
Explanation: I'd say that's pretty good for a first go
Calculate a $\chi^2$ value for this fit
With the hypothesis that these two distributions are the same, we calculate:
$$ \chi^2 = \sum_i \frac{(O_i - E_i)^2}{E_i} $$
End of explanation
from scipy.stats import chisquare
chisq, pvalue = chisquare(Y_plot, optim_fit, ddof=12)
print ("Chi-squared: %.02f\np-value: %.02f" % (chisq, pvalue))
Explanation: For 12 degrees of freedom (12 fit parameters), we can look at a $\chi^2$ table to find that we have a $p>0.995$ that this is a good fit to the distribution.
Or use scipy.stats to get the $\chi^2$ and p-value
End of explanation
sns.distplot(data_df.query("label == 1")[feature],
             kde=False, hist=True, norm_hist=True)
for i in range(0,4):
    plt.plot(X_plot,
             gauss_dist(X_plot[:, 0], popt[i], popt[i+4], popt[i+8]),
             '-', lw=2)
Explanation: Overlay onto original data
End of explanation
<END_TASK>