diff options
| author | tlatorre <tlatorre@uchicago.edu> | 2020-04-13 12:53:08 -0500 |
|---|---|---|
| committer | tlatorre <tlatorre@uchicago.edu> | 2020-04-13 12:53:08 -0500 |
| commit | ba6c298b213d87153d56e12494675ce01207a440 (patch) | |
| tree | f35d0eb87146e402fa0c510282e93bd6ddc6f2e6 | |
| parent | aa2699809e9cd6cd48b57a011af278660e30d9ab (diff) | |
| download | sddm-ba6c298b213d87153d56e12494675ce01207a440.tar.gz sddm-ba6c298b213d87153d56e12494675ce01207a440.tar.bz2 sddm-ba6c298b213d87153d56e12494675ce01207a440.zip | |
fix a few bugs in plot-energy
| -rwxr-xr-x | utils/plot-energy | 35 |
1 files changed, 18 insertions, 17 deletions
diff --git a/utils/plot-energy b/utils/plot-energy index 5079c15..4ec77f8 100755 --- a/utils/plot-energy +++ b/utils/plot-energy @@ -435,7 +435,7 @@ def plot_corner_plot(ev, title, save=None): ev = ev.dropna(subset=variables) - f = plt.figure(figsize=(FIGSIZE[0],FIGSIZE[0])) + fig = plt.figure(figsize=(FIGSIZE[0],FIGSIZE[0])) despine(fig,trim=True) for i in range(len(variables)): for j in range(len(variables)): @@ -446,25 +446,26 @@ def plot_corner_plot(ev, title, save=None): plt.hist(ev[variables[i]],bins=np.linspace(limits[i][0],limits[i][1],100),histtype='step') plt.gca().set_xlim(limits[i]) else: - p_i_lo = np.count_nonzero(ev[variables[i]] < cuts[i])/len(ev) - p_j_lo = np.count_nonzero(ev[variables[j]] < cuts[j])/len(ev) - p_lolo = p_i_lo*p_j_lo - p_lohi = p_i_lo*(1-p_j_lo) - p_hilo = (1-p_i_lo)*p_j_lo - p_hihi = (1-p_i_lo)*(1-p_j_lo) - n_lolo = np.count_nonzero((ev[variables[i]] < cuts[i]) & (ev[variables[j]] < cuts[j])) - n_lohi = np.count_nonzero((ev[variables[i]] < cuts[i]) & (ev[variables[j]] >= cuts[j])) - n_hilo = np.count_nonzero((ev[variables[i]] >= cuts[i]) & (ev[variables[j]] < cuts[j])) - n_hihi = np.count_nonzero((ev[variables[i]] >= cuts[i]) & (ev[variables[j]] >= cuts[j])) - n = len(ev) - observed = np.array([n_lolo,n_lohi,n_hilo,n_hihi]) - expected = n*np.array([p_lolo,p_lohi,p_hilo,p_hihi]) - psi = -poisson.logpmf(observed,expected).sum() + poisson.logpmf(observed,observed).sum() - psi /= np.std(-poisson.logpmf(np.random.poisson(observed,size=(10000,4)),observed).sum(axis=1) + poisson.logpmf(observed,observed).sum()) plt.scatter(ev[variables[j]],ev[variables[i]],s=0.5) plt.gca().set_xlim(limits[j]) plt.gca().set_ylim(limits[i]) - plt.title(r"$\psi = %.1f$" % psi) + n = len(ev) + if n: + p_i_lo = np.count_nonzero(ev[variables[i]] < cuts[i])/n + p_j_lo = np.count_nonzero(ev[variables[j]] < cuts[j])/n + p_lolo = p_i_lo*p_j_lo + p_lohi = p_i_lo*(1-p_j_lo) + p_hilo = (1-p_i_lo)*p_j_lo + p_hihi = (1-p_i_lo)*(1-p_j_lo) + n_lolo = np.count_nonzero((ev[variables[i]] < cuts[i]) & (ev[variables[j]] < cuts[j])) + n_lohi = np.count_nonzero((ev[variables[i]] < cuts[i]) & (ev[variables[j]] >= cuts[j])) + n_hilo = np.count_nonzero((ev[variables[i]] >= cuts[i]) & (ev[variables[j]] < cuts[j])) + n_hihi = np.count_nonzero((ev[variables[i]] >= cuts[i]) & (ev[variables[j]] >= cuts[j])) + observed = np.array([n_lolo,n_lohi,n_hilo,n_hihi]) + expected = n*np.array([p_lolo,p_lohi,p_hilo,p_hihi]) + psi = -poisson.logpmf(observed,expected).sum() + poisson.logpmf(observed,observed).sum() + psi /= np.std(-poisson.logpmf(np.random.poisson(observed,size=(10000,4)),observed).sum(axis=1) + poisson.logpmf(observed,observed).sum()) + plt.title(r"$\psi = %.1f$" % psi) if i == len(variables) - 1: plt.xlabel(labels[j]) else: |