docs: def _calc_var_theta(self):
          """Calculate an estimate of the var(theta)"""
          if self.decaying_prior:
              n_sampled = np.clip(self.alpha_ + self.beta_, 1, np.inf)
              prior_weight = 1/n_sampled
              alpha = self.alpha_ + prior_weight * self.alpha_0
              beta = self.beta_ + prior_weight * self.beta_0
          else:
              alpha = self.alpha_ + self.alpha_0
              beta = self.beta_ + self.beta_0
          # Variance of Beta-distributed rv
          self.var_theta_ = ( alpha * beta /
                              ((alpha + beta)**2 * (alpha + beta + 1)) )