diff --git a/swp2.py b/swp2.py index 346362f..8c28299 100644 --- a/swp2.py +++ b/swp2.py @@ -262,14 +262,6 @@ def plot_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc # divide by N0 y[i] = y[i]/N0 x[i] = x[i]/N0 - sum_theta_i = 0 - print(epoch, x, y) - for i in range(2, len(y)-1): - sum_theta_i=y[i] / (i-1) - prop = [] - for k in range(2, len(y)-1): - prop.append(y[k+1] / (k - 1) / sum_theta_i) - #print(epoch, prop) plt.plot(x, y, 'o', linestyle = "-", alpha=0.75, lw=2, label = str(epoch)+' BrkPt | Lik='+greatest_likelihood) if theta_scale: plt.xlabel("Coal. time") @@ -294,15 +286,15 @@ def plot_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc # number of monomorphic sites L = L_stored S0 = L-S - print("SFS", SFS_stored) - print("S", S, "L", L, "S0=", S0) + # print("SFS", SFS_stored) + # print("S", S, "L", L, "S0=", S0) # compute Ln Ln = log_facto(S+S0) - log_facto(S0) + np.log(float(S0)/(S+S0)) * S0 for xi in range(0, len(SFS_stored)): p_i = SFS_stored[xi] / float(S+S0) Ln += np.log(p_i) * SFS_stored[xi] - log_facto(SFS_stored[xi]) res = Ln - print(res) + # print(res) # basic plot likelihood plt.figure(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi) plt.rcParams['font.size'] = '18' @@ -362,6 +354,14 @@ def plot_test_theta(folder_path, mu, tgen, title = "Title", theta_scale = True, N0 = y[0] for i in range(len(y)): y[i] = y[i]/N0 + # compute the proportion of information used at each bin of the SFS + sum_theta_i = 0 + for i in range(2, len(y)-1): + sum_theta_i+=y[i] / (i-1) + prop = [] + for k in range(2, len(y)-1): + prop.append(y[k] / (k - 1) / sum_theta_i) + # plot plt.plot(x, y, 'o', linestyle="dotted", alpha=0.75, lw=2, label = str(epoch)+' brks') plt.xlabel("# breaks") plt.ylabel("theta")