Calculate an estimate of the var(theta)

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)) )