scipy-yli/tests/test_hdi.py

#   scipy-yli: Helpful SciPy utilities and recipes
#   Copyright © 2022  Lee Yingtong Li (RunasSudo)
#
#   This program is free software: you can redistribute it and/or modify
#   it under the terms of the GNU Affero General Public License as published by
#   the Free Software Foundation, either version 3 of the License, or
#   (at your option) any later version.
#
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU Affero General Public License for more details.
#
#   You should have received a copy of the GNU Affero General Public License
#   along with this program.  If not, see <https://www.gnu.org/licenses/>.

from pytest import approx

import numpy as np
from scipy import stats

import yli

def test_hdi_beta_vs_r():
	"""Compare yli.hdi for beta distribution with R binom.bayes"""
	
	# Inference on a beta-binomial posterior
	n, N = 12, 250
	prior_a, prior_b = 1, 1
	
	distribution = stats.beta(n + prior_a, N - n + prior_b)
	
	# Compute HDI
	hdi = yli.hdi(distribution, 0.95)
	
	#> library(binom)
	#> binom.bayes(12, 250, conf.level=0.95, type='highest', prior.shape1=1, prior.shape2=1)
	#  method  x   n shape1 shape2      mean      lower     upper  sig
	#1  bayes 12 250     13    239 0.0515873 0.02594348 0.0793006 0.05
	expected = np.array([0.02594348, 0.0793006])
	
	assert hdi == approx(expected)