#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Fri Mar 31 15:41:58 2023 @author: Konstantinos Paschos, Dimitris Stavropoulos """ import numpy as np import pandas as pd import matplotlib.pyplot as plt from sklearn.utils import shuffle from sklearn.model_selection import train_test_split from sklearn.ensemble import RandomForestClassifier, HistGradientBoostingClassifier, AdaBoostClassifier from sklearn import metrics from sklearn.metrics import ConfusionMatrixDisplay from sklearn.metrics import roc_curve, auc, precision_recall_curve from sklearn.model_selection import RandomizedSearchCV from sklearn.feature_selection import SelectKBest, f_classif #from xgboost import XGBClassifier from scipy.stats import chisquare, ks_2samp import seaborn as sns features = pd.read_csv("features_Em2_pEnDist.csv") """ features dataset shuffling and splitting """ print("---> shuffling dataset...") # shuffle dataset features_shuffled = shuffle(features, random_state=9) print(features_shuffled) # feature extraction #feature_values = features_shuffled.values X = features_shuffled.iloc[:, :-5] quant = features_shuffled.iloc[:, -5:-1] y = features_shuffled.iloc[:, -1] print("---> ",sum(y=='signal')) ''' # feature importance print("--->Estimating feature importance...") test = SelectKBest(score_func=f_classif, k=10) fit = test.fit(X, Y) # summarize scores plt.figure(1) plt.bar(features.columns.to_numpy()[:-1], fit.scores_) plt.show() ''' ''' # tuning a bdt to find the best parameters bdt_scan = HistGradientBoostingClassifier() tuned_parameters = [{'max_depth': [3, 4, 5, 6, 7, 8], 'max_iter': [50, 100, 150, 200, 250, 300, 350, 400, 450, 500], "learning_rate": [0.05, 0.1, 0.2, 0.3, 0.4, 0.5]}] # Use random search to find the best hyperparameters rand_search = RandomizedSearchCV(bdt_scan, param_distributions=tuned_parameters, n_iter=50, cv=5) rand_search.fit(X_train, y_train) means = rand_search.cv_results_["mean_test_score"] stds = rand_search.cv_results_["std_test_score"] for mean, std, params in zip(means, stds, rand_search.cv_results_['params']): print(f"{mean:0.3f} (+/-{std * 2:0.03f}) for {params}") ''' # create train and test datasets X_train, X_test, y_train, y_test, quant_train, quant_test = train_test_split(X, y, quant, test_size=0.2 , random_state=42 ) # Create and fit the model bdt = HistGradientBoostingClassifier(max_iter=150, max_depth=3, learning_rate=0.2, random_state=42) #HistGradientBoostingClassifier(max_iter=500, max_depth=5, learning_rate=0.05) #bdt = AdaBoostClassifier(n_estimators=40, learning_rate=0.2, random_state=42) bdt.fit(X_train, y_train) y_pred = bdt.predict(X_test) print('The accuracy of the Decision Tree is', "{:.5f}".format(metrics.accuracy_score(y_pred, y_test))) # Calculate probabilities for test and train samples scores_train = bdt.predict_proba(X_train) scores_test = bdt.predict_proba(X_test) # Create dataframes with the information of train and test samples test_concat = pd.concat( [pd.DataFrame(quant_test), pd.DataFrame(y_test)] , axis=1) test_concat['scores'] = scores_test[:,1] train_set = pd.DataFrame(y_train) train_set['scores'] = scores_train[:,1] # Get test samples for signal and noise train_signal = train_set[train_set['label']=='signal'] train_noise = train_set[train_set['label']=='noise'] test_signal = test_concat[ test_concat['label']=='signal'] test_noise = test_concat[ test_concat['label']=='noise'] # Binnings for visualisation score_binning = np.arange( 0.0 , 1.01 , 0.02 ) energy_binning = np.arange( 10.0 , 11.01 , 0.1 ) # BDT score distribution for test smple plt.hist(test_signal.loc[:,'scores'], bins=score_binning, alpha=1.0, label="signal", histtype='step', stacked=False, linewidth=2) plt.hist(test_noise.loc[:,'scores'], bins=score_binning, alpha=1.0, label="noise", histtype='step', stacked=False, linewidth=2) plt.title('BDT score for test samples') plt.ylim( 0.0 , 1.0 ) plt.ylim( 0.5 , 50000 ) plt.yscale("log") plt.grid() plt.legend(loc='upper left') plt.show() # ROC curve fpr, tpr, thresholds = roc_curve(test_concat['label'], test_concat['scores'] , pos_label="signal" ) roc_auc = auc(fpr, tpr) fig, ax = plt.subplots() ax.plot(fpr, tpr, label='%s ROC (area = %0.2f)' % ("BDT", roc_auc)) ax.plot([0, 1], [0, 1], 'k--') ax.plot([0, 0], [0, 1], 'r--') ax.plot([0, 1], [1, 1], 'r--') ax.set_xlabel('False Positive Rate', fontsize=16) ax.set_ylabel('True Positive Rate', fontsize=16) ax.legend(loc=0, fontsize=20) plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.show() # =================================================== # Overtraining check # =================================================== train_signal_entries, junk = np.histogram( train_signal.loc[:,'scores'] , bins=score_binning ) train_noise_entries, junk = np.histogram( train_noise.loc[:,'scores'] , bins=score_binning ) test_signal_entries, junk = np.histogram( test_signal.loc[:,'scores'] , bins=score_binning ) test_noise_entries, junk = np.histogram( test_noise.loc[:,'scores'] , bins=score_binning ) ks_stat, p_value = ks_2samp(train_signal.loc[:,'scores'],test_signal.loc[:,'scores']) print("Signal:") print("KS statistic:", ks_stat) print("KS Test p-value:", p_value) print("==============") print("xxxxxxxxxxxxxx") print("==============") ks_stat, p_value = ks_2samp(train_noise.loc[:,'scores'],test_noise.loc[:,'scores']) print("Noise:") print("KS statistic:", ks_stat) print("KS Test p-value:", p_value) print("==============") print("xxxxxxxxxxxxxx") print("==============") print("==============") print("train_signal_entries: ",sum(train_signal_entries)) print("train_noise_entries: ",sum(train_noise_entries)) print("test_signal_entries: ",sum(test_signal_entries)) print("test_noise_entries: ",sum(test_noise_entries)) print("==============") train_signal_entries = train_signal_entries/sum(train_signal_entries) train_noise_entries = train_noise_entries/sum(train_noise_entries) test_signal_entries = test_signal_entries/sum(test_signal_entries) test_noise_entries = test_noise_entries/sum(test_noise_entries) # KS test to check overtraining ks_stat, p_value = ks_2samp(train_signal_entries, test_signal_entries) print("KS statistic:", ks_stat) print("KS Test p-value:", p_value) # X^2 crashes due to zero bin entries #chisq, p_value = chisquare(f_obs=train_signal_entries, f_exp=test_signal_entries) #print("Chi-Square Test p-value:", p_value) plt.hist( train_signal.loc[:,'scores'], bins=score_binning , alpha=1.0,label="train-signal", histtype='step', stacked=False, linewidth=1, density=True ) plt.hist( train_noise.loc[:,'scores'], bins=score_binning , alpha=1.0,label="train-noise", histtype='step', stacked=False, linewidth=1, density=True , color='blue') plt.hist( test_signal.loc[:,'scores'], bins=score_binning , alpha=1.0,label="test-signal", histtype='step', stacked=False, linewidth=1, density=True ) plt.hist( test_noise.loc[:,'scores'], bins=score_binning , alpha=1.0,label="test-noise", histtype='step', stacked=False, linewidth=1, density=True ) plt.legend(loc='upper left') plt.show() # =================================================== # F1 score to decide optimal cut # =================================================== precisions, recalls, thresholds = precision_recall_curve( test_concat['label'], test_concat['scores'] , pos_label="signal") f1_scores = 2 * (precisions[:-1] * recalls[:-1]) / (precisions[:-1] + recalls[:-1]) optimal_idx = f1_scores.argmax() optimal_threshold = thresholds[optimal_idx] print('F1-score') print("For cut at ",optimal_threshold,", Efficinecy= ", recalls[optimal_idx], " purity= ", precisions[optimal_idx]) # Plot F1-scores as a function of thresholds plt.plot(thresholds, precisions[:-1], label='precision (purity)') plt.plot(thresholds, recalls[:-1], label='recall (efficiency)') plt.plot(thresholds, f1_scores,label='F1-Score') plt.ylim(0.45 , 1.02 ) x_thres = [ optimal_threshold, optimal_threshold] y_thres = [ 0 ,1 ] plt.plot( x_thres , y_thres, label='optimal score cut', linestyle='dashed',color='black') plt.xlabel('BDT score cut') #plt.ylabel('F1-Score') plt.title('F1-Score vs. BDT score cut value') plt.legend() plt.grid() plt.show() # =================================================== # Efficiency = f(energy), for several distances # =================================================== test_signal = pd.DataFrame(test_signal.copy()) # calculate distances test_signal.loc[:,'dist'] = np.sqrt( test_signal.loc[:,'x']**2 + test_signal.loc[:,'y']**2 + test_signal.loc[:,'z']**2 ) # filter dataframe for different distances and BDT score cut pulses_50m = test_signal.query('dist <= 75.0') pulses_100m = test_signal.query('75.0 <= dist < 125.0') pulses_150m = test_signal.query('125.0 <= dist') pulses_200m = test_signal.query('175.0 <= dist < 225.0') pulses_250m = test_signal.query('225.0 <= dist') # plot, divided with total number/bin plt.hist( pulses_50m.loc[:,'energy'] , bins=energy_binning , histtype='step', label='all') plt.hist( pulses_50m[pulses_50m.loc[:,'scores']>optimal_threshold].loc[:,'energy'] , bins=energy_binning , histtype='step', label='selected') plt.legend() plt.title('Pulses @ 50 m') plt.show() plt.hist( pulses_100m.loc[:,'energy'] , bins=energy_binning , histtype='step', label='all') plt.hist( pulses_100m[pulses_100m.loc[:,'scores']>optimal_threshold].loc[:,'energy'] , bins=energy_binning , histtype='step', label='selected') plt.legend() plt.title('Pulses @ 100 m') plt.show() plt.hist( pulses_150m.loc[:,'energy'] , bins=energy_binning , histtype='step', label='all') plt.hist( pulses_150m[pulses_150m.loc[:,'scores']>optimal_threshold].loc[:,'energy'] , bins=energy_binning , histtype='step', label='selected') plt.legend() plt.title('Pulses @ 150 m') plt.show() plt.hist( pulses_200m.loc[:,'energy'] , bins=energy_binning , histtype='step', label='all') plt.hist( pulses_200m[pulses_200m.loc[:,'scores']>optimal_threshold].loc[:,'energy'] , bins=energy_binning , histtype='step', label='selected') plt.legend() plt.title('Pulses @ 200 m') plt.show() plt.hist( pulses_250m.loc[:,'energy'] , bins=energy_binning , histtype='step', label='all') plt.hist( pulses_250m[pulses_250m.loc[:,'scores']>optimal_threshold].loc[:,'energy'] , bins=energy_binning , histtype='step', label='selected') plt.legend() plt.title('Pulses @ 250 m') plt.show() bin_centers = (energy_binning[:-1] + energy_binning[1:]) / 2 # Centers of the bins # Compute efficiency per bin efficiency_all = [] for i in range(len(energy_binning) - 1): # Events in the current energy bin in_bin = (test_signal.loc[:,'energy'] >= energy_binning[i]) & (test_signal.loc[:,'energy'] < energy_binning[i + 1]) total_in_bin = np.sum( in_bin ) # Total events in bin passed_in_bin = np.sum( (test_signal.loc[:,'scores'] >= optimal_threshold) & in_bin ) # Events passing the BDT cut in bin efficiency_all.append(passed_in_bin/total_in_bin) efficiency_50m = [] for i in range(len(energy_binning) - 1): # Events in the current energy bin in_bin = (pulses_50m.loc[:,'energy'] >= energy_binning[i]) & (pulses_50m.loc[:,'energy'] < energy_binning[i + 1]) total_in_bin = np.sum( in_bin ) # Total events in bin passed_in_bin = np.sum( (pulses_50m.loc[:,'scores'] >= optimal_threshold) & in_bin ) # Events passing the BDT cut in bin efficiency_50m.append(passed_in_bin/total_in_bin) efficiency_100m = [] for i in range(len(energy_binning) - 1): # Events in the current energy bin in_bin = (pulses_100m.loc[:,'energy'] >= energy_binning[i]) & (pulses_100m.loc[:,'energy'] < energy_binning[i + 1]) total_in_bin = np.sum( in_bin ) # Total events in bin passed_in_bin = np.sum( (pulses_100m.loc[:,'scores'] >= optimal_threshold) & in_bin ) # Events passing the BDT cut in bin efficiency_100m.append(passed_in_bin/total_in_bin) efficiency_150m = [] for i in range(len(energy_binning) - 1): # Events in the current energy bin in_bin = (pulses_150m.loc[:,'energy'] >= energy_binning[i]) & (pulses_150m.loc[:,'energy'] < energy_binning[i + 1]) total_in_bin = np.sum( in_bin ) # Total events in bin passed_in_bin = np.sum( (pulses_150m.loc[:,'scores'] >= optimal_threshold) & in_bin ) # Events passing the BDT cut in bin efficiency_150m.append(passed_in_bin/total_in_bin) efficiency_200m = [] for i in range(len(energy_binning) - 1): # Events in the current energy bin in_bin = (pulses_200m.loc[:,'energy'] >= energy_binning[i]) & (pulses_200m.loc[:,'energy'] < energy_binning[i + 1]) total_in_bin = np.sum( in_bin ) # Total events in bin passed_in_bin = np.sum( (pulses_200m.loc[:,'scores'] >= optimal_threshold) & in_bin ) # Events passing the BDT cut in bin efficiency_200m.append(passed_in_bin/total_in_bin) efficiency_250m = [] for i in range(len(energy_binning) - 1): # Events in the current energy bin in_bin = (pulses_250m.loc[:,'energy'] >= energy_binning[i]) & (pulses_250m.loc[:,'energy'] < energy_binning[i + 1]) total_in_bin = np.sum( in_bin ) # Total events in bin passed_in_bin = np.sum( (pulses_250m.loc[:,'scores'] >= optimal_threshold) & in_bin ) # Events passing the BDT cut in bin efficiency_250m.append(passed_in_bin/total_in_bin) plt.plot(bin_centers, efficiency_50m, marker='o', linestyle='-', label='50 m') plt.plot(bin_centers, efficiency_100m, marker='o', linestyle='-', label='100 m') plt.plot(bin_centers, efficiency_150m, marker='o', linestyle='-', label='>100 m') '''plt.plot(bin_centers, efficiency_200m, marker='o', linestyle='-', label='200 m') plt.plot(bin_centers, efficiency_250m, marker='o', linestyle='-', label='250 m')''' plt.legend() plt.ylim( 0.0, 1.05) plt.grid() plt.show() plt.plot(bin_centers, efficiency_all, marker='o', linestyle='-', label='all pulses') plt.legend() plt.ylim( 0.0, 1.05) plt.grid() plt.show() ''' # =================================================== # ROC curve and Jouden's index to decide optimal cut # =================================================== # not confident about the following fpr, tpr, thresholds = roc_curve(test_concat['label'], test_concat['scores'] , pos_label="signal" ) # Compute Youden's Index youden_index = tpr - fpr - 1.0 optimal_idx = youden_index.argmax() optimal_threshold = thresholds[optimal_idx] plt.plot(thresholds,youden_index) plt.show() print("Optimal Threshold J-index:", optimal_threshold, " Efficinecy= ", recalls[optimal_idx], " purity= ", precisions[optimal_idx]) ''' ''' confusion matrix ''' ''' #disp = ConfusionMatrixDisplay.from_predictions(y_test, y_pred, normalize="true", display_labels=["background", "signal"]) #plt.show() noise = sum(y_pred[y_test == "noise"] == "noise") / sum(y_test == "noise") signal = sum(y_pred[y_test == "signal"] == "signal") / sum(y_test == "signal") heat = pd.DataFrame( {"Predicted": ["signal", "background"], "background": [1 - noise, noise], "signal": [signal, 1 - signal]}) heat = heat.set_index("Predicted") heat = heat.rename_axis('True', axis=1) sns.set(font_scale=1.3) sns.heatmap(heat, annot=True, cmap="Blues") plt.show() ''' ''' confusion matrix (george) ''' '''heat = pd.DataFrame({"Predicted": ["signal", "background"], "background": [0.047, 0.953], "signal": [0.759, 0.241]}) heat = heat.set_index("Predicted") heat = heat.rename_axis('True', axis=1) sns.set(font_scale=1.3) sns.heatmap(heat, annot=True, cmap="Blues") bdt.probability = True probs = bdt.predict_proba(X_test) fpr, tpr, thresholds = roc_curve(y_test, probs[:, 1], pos_label="signal") roc_auc = auc(fpr, tpr) fig, ax = plt.subplots() ax.plot(fpr, tpr, label='%s ROC (area = %0.2f)' % ("BDT", roc_auc)) ax.plot([0, 1], [0, 1], 'k--') ax.plot([0, 0], [0, 1], 'r--') ax.plot([0, 1], [1, 1], 'r--') ax.set_xlabel('False Positive Rate', fontsize=16) ax.set_ylabel('True Positive Rate', fontsize=16) ax.legend(loc=0, fontsize=20) plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.show()''' ''' DStavropoulos # Plot efficiency, purity and significance by hand test_signal_entries, junk = np.histogram( test_signal.loc[:,'scores'] , bins=score_binning ) test_noise_entries, junk = np.histogram( test_noise.loc[:,'scores'] , bins=score_binning ) test_signal_entries_cum = np.cumsum(test_signal_entries[::-1]) test_signal_entries_cum_norm = test_signal_entries_cum/sum(test_signal_entries) test_signal_entries_cum_norm = np.append( test_signal_entries_cum_norm , 1.0 ) # to match number of elements with bin edges test_signal_entries_invcum_norm = test_signal_entries_cum_norm[::-1] test_noise_entries_cum = np.cumsum(test_noise_entries[::-1]) test_noise_entries_cum_norm = test_noise_entries_cum/sum(test_noise_entries) test_noise_entries_cum_norm = np.append( test_noise_entries_cum_norm , 1.0 ) # to match number of elements with bin edges test_noise_entries_invcum_norm = test_noise_entries_cum_norm[::-1] test_signal_entries_cum = np.append( test_signal_entries_cum , test_signal_entries_cum[-1] ) test_noise_entries_cum = np.append( test_noise_entries_cum , test_noise_entries_cum[-1] ) significance = test_signal_entries_cum/np.sqrt((test_signal_entries_cum**2)+(test_noise_entries_cum**2)) significance = significance[::-1] #significance = test_signal_entries_cum/np.sqrt(test_noise_entries_cum) plt.plot( score_binning, test_signal_entries_invcum_norm ) plt.plot( score_binning, test_noise_entries_invcum_norm ) plt.plot( score_binning, significance ) plt.show() ''' ''' # tuning a bdt to find the best parameters bdt_scan = HistGradientBoostingClassifier() tuned_parameters = [{'max_depth': [3, 4, 5, 6, 7, 8], 'max_iter': [50, 100, 150, 200, 250, 300, 350, 400, 450, 500], "learning_rate": [0.05, 0.1, 0.2, 0.3, 0.4, 0.5]}] # Use random search to find the best hyperparameters rand_search = RandomizedSearchCV(bdt_scan, param_distributions=tuned_parameters, n_iter=50, cv=5) rand_search.fit(X_train, y_train) means = rand_search.cv_results_["mean_test_score"] stds = rand_search.cv_results_["std_test_score"] for mean, std, params in zip(means, stds, rand_search.cv_results_['params']): print(f"{mean:0.3f} (+/-{std * 2:0.03f}) for {params}") ''' ''' # random forest classifier generation and training rft = RandomForestClassifier(n_estimators=1000, oob_score=True, max_depth=3) rft.fit(X_train, y_train) # test the tree y_pred = rft.predict(X_test) print('The accuracy of the Decision Tree is', "{:.3f}".format(metrics.accuracy_score(y_pred, y_test))) # importance of each feature rft.feature_importances_ plt.figure(2) disp = ConfusionMatrixDisplay.from_predictions(y_test, y_pred, normalize="true") # predict probabilities: the predict_proba method returns the probabilities for each data point to belong in a class rft.probability = True probs = rft.predict_proba(X_test) fpr, tpr, thresholds = roc_curve(y_test, probs[:, 1], pos_label="signal") roc_auc = auc(fpr, tpr) fig, ax = plt.subplots() ax.plot(fpr, tpr, label='%s ROC (area = %0.2f)' % ("RandomForest", roc_auc)) ax.plot([0, 1], [0, 1], 'k--') ax.set_xlabel('False Positive Rate') ax.set_ylabel('True Positive Rate') ax.legend(loc=0, fontsize='small') # classify the evaluation dataset X_eval, y_eval = evaluate.iloc[:, :-1], evaluate.iloc[:, -1] y_eval_pred = rft.predict(X_eval) disp_eval = ConfusionMatrixDisplay.from_predictions(y_eval, y_eval_pred, normalize="true") # tuning a rft to find the best parameters rf = RandomForestClassifier(n_estimators=500) tuned_parameters = [{'max_depth': [3, 4, 5], 'min_samples_split': [2, 3, 4]}] # Use random search to find the best hyperparameters rand_search = RandomizedSearchCV(rf, param_distributions=tuned_parameters, n_iter=9, cv=5) rand_search.fit(X_train, y_train) means = rand_search.cv_results_["mean_test_score"] stds = rand_search.cv_results_["std_test_score"] for mean, std, params in zip(means, stds, rand_search.cv_results_['params']): print(f"{mean:0.3f} (+/-{std * 2:0.03f}) for {params}") # generate preferred tree rft_final = RandomForestClassifier(n_estimators=1000, oob_score=True, max_depth=6, min_samples_split=5) rft_final.fit(X_train, y_train) y_pred_final = rft_final.predict(X_test) print('The accuracy of the Decision Tree is', "{:.3f}".format(metrics.accuracy_score(y_pred_final, y_test))) disp_final = ConfusionMatrixDisplay.from_predictions(y_test, y_pred_final, normalize="true") # classify the evaluation dataset X_eval, y_eval = evaluate.iloc[:, :-1], evaluate.iloc[:, -1] y_eval_pred = rft_final.predict(X_eval) disp_eval = ConfusionMatrixDisplay.from_predictions(y_eval, y_eval_pred, normalize="true") ''' ''' bins = np.linspace(0, 100, 2101) fig, ax = plt.subplots() # histogramms for the noise and the signal chunks has = ax.hist(maxmean_signal, bins) no = ax.hist(maxmean_noise, bins) # cumulative curves: efficiency, contamination, purity eff = [] cont = [] eff_noise = [] signal = features[features['label'] == "signal"]["col_sum"] noise = features[features['label'] == "noise"]["col_sum"] bins = np.linspace(0, 100, 2101) fig, ax = plt.subplots() has = ax.hist(signal, bins) no = ax.hist(noise, bins) for i in range(len(has[0][:])): eff.append(sum(has[0][i:])) eff_noise.append(sum(no[0][i:])) cont.append(sum(no[0][i:]) / sum(has[0][i:])) eff = np.array(eff) / eff[0] eff_noise = np.array(eff_noise) / eff_noise[0] pur = 1 - np.array(cont) eff70 = eff[np.logical_and(eff < 0.72, eff > 0.68)] index = np.where(np.isin(eff, eff70)) observ = np.mean(pur[index]) # plots of p2p heights of pulses identified as noise by the bdt y_pred_df = pd.Series(data=y_pred, index=y_test.index) ind = np.array(y_pred_df[y_test == "signal"][y_pred_df[y_test == "signal"] == "noise"].index) ind_pulses = list(ind // 2) ind2 = np.array(y_pred_df[y_test == "signal"][y_pred_df[y_test == "signal"] == "signal"].index) ind_pulses2 = list(ind2 // 2) header = ['t_' + str(i + 1) for i in range(147)] + \ ['x', 'y', 'z', 'E', 'width'] pulses = pd.read_csv('./Simulations/pulse_production.csv', names=header, index_col=None) for index, row in pulses.iterrows(): pulses.loc[index, "width_resampled"] = np.max(row[:147]) - np.min(row[:147]) pulses_test = pulses.iloc[np.array(y_test[y_test == "signal"].index) // 2] fig, ax = plt.subplots() color = 'tab:blue' bins = np.linspace(0, 1, 500) detected = pulses.iloc[ind_pulses2, -1] not_detected = pulses.iloc[ind_pulses, -1] det, bins, patches_det = ax.hist(detected, bins=bins, alpha=0.8, color='tab:blue', label="detected") not_det, bins, patches_not_det = ax.hist(not_detected, bins=bins, alpha=0.8, color='tab:orange', label="not detected") ax.set_xticks(np.linspace(0, 1, 11)) detected_cum = [] not_detected_cum = [] for i in range(len(det[:])): detected_cum.append(sum(det[:i])) for i in range(len(not_det)): not_detected_cum.append(sum(not_det[:i])) detected_cum = np.array(detected_cum) / detected_cum[-1] not_detected_cum = np.array(not_detected_cum) / not_detected_cum[-1] ax.grid() ax.set_xlabel("peak to peak width [Pa]") ax.set_ylabel("count [-]", color=color) ax.tick_params(axis='y', labelcolor=color) ax2 = ax.twinx() # instantiate a second axes that shares the same x-axis color = 'tab:red' det_cum, = ax2.plot(bins[:-1], detected_cum, color='tab:blue', label="detected") not_det_cum, = ax2.plot(bins[:-1], not_detected_cum, color='tab:orange', label="detected") leg = fig.legend([(patches_det, det_cum), (patches_not_det, not_det_cum)], ["detected", "not detected"], 'right') leg.set_draggable(state=True) ax2.set_ylabel('cumulative curves [-]', color=color) # we already handled the x-label with ax1 ax2.tick_params(axis='y', labelcolor=color) plt.show() p2p_noise = [] for i in ind: i_true = i while np.isin(i_true, fish): i_true = i_true + i_end chunk = y[space * i_true: space * (i_true + 1)] chunk = chunk / 32768.0 # normalization of raw file chunk = chunk - dc_offset yMicroPascals = chunk / sensitivityLinear yPascals = yMicroPascals / 1000000 Filtered_5 = filtfilt(a, b, yPascals) p2p_noise.append(max(Filtered_5) - min(Filtered_5)) p2p_noise2 = [] for i in ind2: i_true = i while np.isin(i_true, fish): i_true = i_true + i_end chunk = y[space * i_true: space * (i_true + 1)] chunk = chunk / 32768.0 # normalization of raw file chunk = chunk - dc_offset yMicroPascals = chunk / sensitivityLinear yPascals = yMicroPascals / 1000000 Filtered_5 = filtfilt(a, b, yPascals) p2p_noise2.append(max(Filtered_5) - min(Filtered_5)) sub_tree_42 = bdt.estimators_[42, 0] '''