sciapy.regress.statistics.waic_loo()

Complexity

Conditions

15

Size

Total Lines	157
Code Lines	74

Duplication

Lines	157
Ratio	100 %

Code Coverage

Tests	1
CRAP Score	229.7642

Importance

Changes

0

# -*- coding: utf-8 -*-
def waic_loo(model, times, data, errs,

This code seems to be duplicated in your project.

		samples,
		method="likelihood",
		train_data=None,
		noisy_targets=True,
		nthreads=1,
		use_dask=False,
		dask_cluster=None,
		):
	"""Watanabe-Akaike information criterion (WAIC) and LOO IC of the (GP) model

	Calculates the WAIC and leave-one-out (LOO) cross validation scores and
	information criteria (IC) from the MCMC samples of the posterior parameter
	distributions. Uses the posterior point-wise (per data point)
	probabilities and the formulae from [1]_ and [2]_.

	.. [1] Vehtari, Gelman, and Gabry, Stat Comput (2017) 27:1413–1432,
		doi: 10.1007/s11222-016-9696-4

	.. [2] Vehtari and Gelman, (unpublished)
		http://www.stat.columbia.edu/~gelman/research/unpublished/waic_stan.pdf
		http://www.stat.columbia.edu/~gelman/research/unpublished/loo_stan.pdf

	Parameters
	----------
	model : `celerite.GP`, `george.GP` or `CeleriteModelSet` instance
		The model instance whose parameter distribution was drawn.
	times : (M,) array_like
		The test coordinates to predict or evaluate the model on.
	data : (M,) array_like
		The test data to test the model against.
	errs : (M,) array_like
		The errors (variances) of the test data.
	samples : (K, L) array_like
		The `K` MCMC samples of the `L` parameter distributions.
	method : str ("likelihood" or "predict"), optional
		The method to "predict" the data, the default uses the (log)likelihood
		in the same way as is done when fitting (training) the model.
		"predict" uses the actual GP prediction, might be useful if the IC
		should be estimated for actual test data that was not used to train
		the model.
	train_data : (N,) array_like, optional
		The data on which the model was trained, needed if method="predict" is
		used, otherwise None is the default and the likelihood is used.
	noisy_targets : bool, optional
		Include the given errors when calculating the predictive probability.
	nthreads : int, optional
		Number of threads to distribute the point-wise probability
		calculations to (default: 1).
	use_dask : boot, optional
		Use `dask.distributed` to distribute the point-wise probability
		calculations to `nthreads` workers. The default is to use
		`multiprocessing.pool.Pool()`.
	dask_cluster: str, or `dask.distributed.Cluster` instance, optional
		Will be passed to `dask.distributed.Client()`
		This can be the address of a Scheduler server like a string
		'127.0.0.1:8786' or a cluster object like `dask.distributed.LocalCluster()`.

	Returns
	-------
	waic, waic_se, p_waic, loo_ic, loo_se, p_loo : tuple
		The WAIC and its standard error as well as the
		estimated effective number of parameters, p_waic.
		The LOO IC, its standard error, and the estimated
		effective number of parameters, p_loo.
	"""
	from functools import partial
	from multiprocessing import pool
	from scipy.special import logsumexp
	try:
		from tqdm.autonotebook import tqdm
	except ImportError:
		tqdm = None
	try:
		from dask.distributed import Client, LocalCluster, progress
	except ImportError:
		use_dask = False

	# the predictive covariance should include the data variance
	# set to a small value if we don't want to account for them
	if not noisy_targets or errs is None:
		errs = 1.123e-12

	# point-wise posterior/predictive probabilities
	_log_p_pt = partial(_log_lh_pt, model, times, data, errs)
	if method == "predict" and train_data is not None:
		_log_p_pt = partial(_log_pred_pt, model, train_data, times, data, errs)

	# calculate the point-wise probabilities and stack them together
	if nthreads > 1:
		if use_dask:
			if dask_cluster is None:
				# start local dask cluster
				_cl = LocalCluster(n_workers=nthreads, threads_per_worker=1)
			else:
				# use provided dask cluster
				_cl = dask_cluster
			_c = Client(_cl)
			_log_pred = _c.map(_log_p_pt, samples)
			progress(_log_pred)
			log_pred = np.stack(_c.gather(_log_pred))
			_c.close()
			if dask_cluster is None:
				_cl.close()
		else:
			# multiprocessing.pool
			_p = pool.Pool(processes=nthreads)
			_mapped = _p.imap_unordered(_log_p_pt, samples)
			if tqdm is not None:
				_mapped = tqdm(_mapped, total=len(samples))
			log_pred = np.stack(list(_mapped))
			_p.close()
			_p.join()
	else:
		if tqdm is not None:
			samples = tqdm(samples, total=len(samples))
		log_pred = np.stack(list(map(_log_p_pt, samples)))

	lppd_i = logsumexp(log_pred, b=1. / log_pred.shape[0], axis=0)
	p_waic_i = np.nanvar(log_pred, ddof=1, axis=0)
	if np.any(p_waic_i > 0.4):
		logging.warn("""For one or more samples the posterior variance of the
		log predictive densities exceeds 0.4. This could be indication of
		WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details
		""")
	elpd_i = lppd_i - p_waic_i
	waic_i = -2. * elpd_i
	waic_se = np.sqrt(len(waic_i) * np.nanvar(waic_i, ddof=1))
	waic = np.nansum(waic_i)
	p_waic = np.nansum(p_waic_i)
	if 2. * p_waic > len(waic_i):
		logging.warn("""p_waic > n / 2,
		the WAIC approximation is unreliable.
		""")
	logging.info("WAIC: %s, waic_se: %s, p_w: %s", waic, waic_se, p_waic)

	# LOO
	loo_ws = 1. / np.exp(log_pred - np.nanmax(log_pred, axis=0))
	loo_ws_n = loo_ws / np.nanmean(loo_ws, axis=0)
	loo_ws_r = np.clip(loo_ws_n, None, np.sqrt(log_pred.shape[0]))
	elpd_loo_i = logsumexp(log_pred,
			b=loo_ws_r / np.nansum(loo_ws_r, axis=0),
			axis=0)
	p_loo_i = lppd_i - elpd_loo_i
	loo_ic_i = -2 * elpd_loo_i
	loo_ic_se = np.sqrt(len(loo_ic_i) * np.nanvar(loo_ic_i))
	loo_ic = np.nansum(loo_ic_i)
	p_loo = np.nansum(p_loo_i)
	logging.info("loo IC: %s, se: %s, p_loo: %s", loo_ic, loo_ic_se, p_loo)

	# van der Linde, 2005, Statistica Neerlandica, 2005
	# https://doi.org/10.1111/j.1467-9574.2005.00278.x
	hy1 = -np.nanmean(lppd_i)
	hy2 = -np.nanmedian(lppd_i)
	logging.info("H(Y): mean %s, median: %s", hy1, hy2)

	return waic, waic_se, p_waic, loo_ic, loo_ic_se, p_loo