docs: def plot(self, agebins=50, p=(2.5, 97.5), ax=None):
          """Age-depth plot"""
          if ax is None:
              ax = plt.gca()
          ax.hist2d(np.repeat(self.depth, self.age_ensemble.shape[1]), self.age_ensemble.flatten(),
                     (len(self.depth), agebins), cmin=1)
          ax.step(self.depth, self.age_median(), where='mid', color='red')
          ax.step(self.depth, self.age_percentile(p[0]), where='mid', color='red', linestyle=':')
          ax.step(self.depth, self.age_percentile(p[1]), where='mid', color='red', linestyle=':')
          ax.set_ylabel('Age (cal yr BP)')
          ax.set_xlabel('Depth (cm)')
          ax.grid(True)
          return ax