I am trying to replicate a plot in Figure 8 of this "https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2012JD017979" article using precipitation data.
Basically, it shows the distribution of (in my case) precipitation where the y-axis (oftentimes in distribution plots, the density) is the normalized volume, allowing different precipitation products to be comparable in one plot.
I have tried the code below using Scipy's stats.binned_statistic function and numpy:
compute mean of values in each bin, bin edges, and bin number
#NB: sim_id = simulated, obs_id = observed
im_bin_means, sim_bin_edges, sim_binnumber = stats.binned_statistic(x = sim_id,values = sim_id,statistic='mean',bins=2000,range=(0.01,120))
sim_bin_cnt, sim_bin_edges, sim_binnumber = stats.binned_statistic(x = sim_id,values = sim_id,statistic='count',bins=2000,range=(0.01,120))
obs_bin_means, obs_bin_edges, obs_binnumber = stats.binned_statistic(x = obs_id,values = obs_id,statistic='mean',bins=2000,range=(0.01,120))
obs_bin_cnt, obs_bin_edges, obs_binnumber = stats.binned_statistic(x = obs_id,values = obs_id,statistic='count',bins=2000,range=(0.01,120))
compute a customized distribution as
sim_dist = (sim_bin_cnt*sim_bin_means)/np.nansum(sim_bin_cnt)
obs_dist = (obs_bin_cnt*obs_bin_means)/np.nansum(obs_bin_cnt)
the plot
f,ax = plt.subplots()
ax.plot(sim_bin_edges[:-1],sim_dist,c='r',label='SIM',lw=2,ls='-')
ax.plot(obs_bin_edges[:-1],obs_dist,c='k',label='OBS',lw=2,ls='--')
ax.set_xscale('log')
ax.legend(loc = 'best',frameon=False)
yielding the plot below
The problem is that the area under each curve should integrate to 1; i.e., np.nansum(sim_dist * np.diff(sim_bin_edges)) should be = 1. But this is not the case with my current method.
I would be glad if I could get help with an elegant method that leads to a more desirable answer. Thanks in advance