npcompare package documentation

class npcompare.Compare(f1, f2, psamples1, psamples2, lower=0, upper=1, weights1=None, weights2=None, metric=None)

Bases: object

Compare two samples.

Parameters:

f1 : function

Density function for the first population.

f2 : function

Density function for the second population.

psamples1 : array like

Must be either a one dimensional numpy.array (in which case, each element will be passed one at a time to f1) or a two dimensional numpy.array (in which case, each row will be passed one at a time to f1).

psamples2 : array like

Analog to psamples1.

weights1 : array like

Give weights to posterior psamples1. Set to None if each posterior sample has the same weight (the usual case for MCMC methods).

weights2 : array like

Analog to weights2.

metric : function

Metric function to be used.

Defaults to:

def metric(f1, f2, param1, param2): return quad(lambda x:
(f1(x, param1) - f2(x, param2))**2, a, b)[0]

Can be set to a user-defined function of the same signature.

lower : float

Lower integration limit passed to default metric function.

upper : float

Upper integration limit passed to default metric function.

classmethod frombfs(bfsobj1, bfsobj2, transformation=True, metric=None)

Create a class Compare from two EstimateBFS objects.

Parameters:

bfsobj1 : EstimateBFS object

bfsobj2 : EstimateBFS object

transformation :

If set to False, metric will be evaluated in [0, 1] without transformation (EstimateBFS class documentation). Otherwise, transformation will be applied.

Ignored if bfsobj1 has no transformation (which implies sample space of [0, 1]).

metric : function

Metric function to be used, defaults to class’s default.

Returns:

New instance of class Compare

plot(ax=None, **kwargs)

Plot empirical CDF of metric samples.

Parameters:

ax : matplotlib axes

Axis to plot, defaults to axes of a new figure.

show : bool

If True, calls matplotlib.pyplot plt.show() at end.

**kwargs :

Aditional named arguments passed to matplotlib.axes.Axes.step.

Returns:

matplotlib.axes.Axes object

sampleposterior(niter=1000, refresh=100)

Compare two samples.

Parameters:

niter : integer

Number of simulations to be draw.

refresh : integer

Interval of samples to print the amount of samples obtained so far. Set to 0 to disable printing.

Returns:

None

class npcompare.EstimateBFS(obs=None, nmaxcomp=10, hpp=1, hpgamma=0, transformation=None, mixture=True, **kwargs)

Bases: object

Estimate univariate density using Bayesian Fourier Series. This method only works with data the lives in [0, 1], however, the class implements methods to automatically transform user inputted data to [0, 1]. See parameter transform below.

Parameters:

obs : array like

List/tuple or numpy.array 1D array of observations. If unset, the constructor will not call method fit, and you will need to call it manually later to configure your model data (this can be usefull for using it with scikit-learns GridSearchCV). nmaxcomp : maximum number of components of the Fourier series expansion.

nmaxcomp : integer

Maximum number of components of the Fourier series expansion.

hpp : float

Hyperparameter p (defaults to “conservative” value of 1).

hpgamma : float

Hyperparameter p (defaults to “conservative” value of 0).

transform :

Transformation function to use. Can be either:

String “logit” to use logit transformation. Usefull if the sample space is the real line. However the obs must actually assume only low values like, for example, between [-5, 5], this is due to the fact that inverse logit starts getting (computationally) very close to 1 (or 0) after 10 (or -10).

Dictionary {“transf”: “fixed”, “vmin”: vmin, “vmax”: vmax} for fixed value transformation, where vmin and vmax are the minimun and maximum values of the sample space (for obs), respectivelly.

A user defined transformation function with a dict with elements transf, itransf and laditransf, where transf is a function that transforms [0, 1] to sample space, itransf is its inverse, and laditransf is the log absolute derivative of itransf. These 3 functions must accept and return numpy 1D arrays.

mixture : bool

If True, will work with a mixture of Fourier series models with up to nmaxcomponent components (that is, a prior will be set on the number of components of the Fourier Series). If False, will work with a single Fourier series model with exactly nmaxcomponent components.

**kwargs:

Additional keyword arguments passed to method fit.

compilestanmodel()

Compile Stan model necessary for method sample. This method is called automatically by obj.sample().

Returns:self
evalgrid(gridsize=1000)

Calculates posterior estimated mean density value at grid points so they can be used later by method plot or directly by the user.

Parameters:gridsize : size of the grid.
Returns:self

Notes

The results are store in an instance variable called egresults which is a dictionary with the following objects:

If object was initialized with mixture=True, you will have:

Log densities for full mixture of components stored as logdensitymixmean.

Densities for full mixture of components stored as densitymixmean.

Log densities for individual components stored as logdensityindivmean.

Densities for individual components stored as densityindivmean.

If object was initialized with mixture=False, you will have:

Log densities stored as logdensitymean.

Densities stored as densitymean.

evaluate(points, logdensity=False, transformed=True, component=None)

Predict posterior density for estimated model.

Parameters:

points : array like or scalar

Scalar or 1D numpy array of points where density will be evaluated.

transformed : bool

True (default) if points live in the sample space of your inputted data. False if your points live in the sample space of Bayesian Fourier Series sample space ([0, 1]).

component : NoneType or integer

Which individual component of the mixture to use. Set to None (default) to full mixture with posterior with sieve prior (strongly recomended). Ignored if object was initialized with mixture=False.

logdensity : bool

If True, will return the logdensity instead of the density.

Returns:

Numpy 1D array with predicted mean density.

fit(obs, nmaxcomp=None, hpp=None, hpgamma=None, transformation=None, niter=5000, **kwargs)

Configure object model data and automatically calls method sampleposterior (to obtain MCMC samples from the posterior).

Parameters:

obs : array like

List/tuple or numpy.array 1D array of observations.

nmaxcomp : integer

Maximum number of components of the Fourier series expansion. Defaults to the one already set in object construction.

hpp : float

Hyperparameter p (defaults to “conservative” value of 1). Defaults to the one already set in object construction.

hpgamma : float

Hyperparameter p (defaults to “conservative” value of 0). Defaults to the one already set in object construction.

transform :

Transformation function to use. Can be either:

String “logit” to use logit transformation. Usefull if the sample space is the real line. However the obs must actually assume only low values like, for example, between [-5, 5], this is due to the fact that inverse logit starts getting (computationally) very close to 1 (or 0) after 10 (or -10).

Dictionary {“transf”: “fixed”, “vmin”: vmin, “vmax”: vmax} for fixed value transformation, where vmin and vmax are the minimun and maximum values of the sample space (for obs), respectivelly.

A user defined transformation function with a dict with elements transf, itransf and laditransf, where transf is a function that transforms [0, 1] to sample space, itransf is its inverse, and laditransf is the log absolute derivative of itransf. These 3 functions must accept and return numpy 1D arrays.

Defaults to the one already set in object construction.

niter : integer

Number of iterations to sample. If set to 0, then will not call method sample.

**kwargs:

Additional keyword arguments passed to method sample.

Returns:

self

plot(ax=None, component=None, **kwargs)

Plot samples.

Parameters:

ax : matplotlib axes

Axis to plot, defaults to axes of a new figure.

component : NoneType or integer

Which individual component of the mixture to use. Set to None (default) to full mixture with posterior with sieve prior (strongly recomended). Ignored if object was initialized with mixture=False.

**kwargs :

Aditional named arguments passed to matplotlib.axes.Axes.step.

Returns:

matplotlib.axes.Axes object

sampleposterior(niter=5000, nchains=4, njobs=-1, tolrhat=0.02, **kwargs)

Samples from posterior.

Parameters:

niter :

Number of simulations to be draw.

nchains :

Number of MCMC chains.

njobs :

Number of CPUs to be used in parallel. If -1 (default), all CPUs will be used.

tolrhat :

Maximum tolerable distance an Rhat of any sampled parameter can have from 1 (we will resample model approximatedly 10% more iterations until this convergence criteria is met).

**kwargs :

Additional keyword arguments passed to pystan (e.g.: refresh parameter to configure sampler printing status).

Returns:

self

score(points, transformed=True, component=None)

Return sum of average logdensity of points.

This is equivalent to calling:

obj.evaluate(points, logdensity=True, transformed,
component).sum()[()]
Parameters:

points : array like or scalar

Scalar or 1D numpy array of points where density will be evaluated.

transformed : bool

True (default) if points live in the sample space of your inputted data. False if your points live in the sample space of Bayesian Fourier Series sample space ([0, 1]).

component : NoneType or integer

Which individual component of the mixture to use. Set to None (default) to full mixture with posterior with sieve prior (strongly recomended). Ignored if object was initialized with mixture=False.

Returns:

Numpy float with sum of predicted mean log densities.

class npcompare.EstimateLindleyBFS(obs0=None, obs1=None, nmaxcomp=10, hpp=1, hpgamma=0, hplindley=0.5, transformation=None, mixture=True, **kwargs)

Bases: object

Estimate univariate density using Bayesian Fourier Series. This method only works with data the lives in [0, 1], however, the class implements methods to automatically transform user inputted data to [0, 1]. See parameter transform below.

Parameters:

obs0 : array like

List/tuple or numpy.array 1D array of observations for model 1.

obs1 : array like

List/tuple or numpy.array 1D array of observations for model 2.

hplindley : float

A priori probability that obs0 and obs1 came from different populations.

**kwargs:

Additional keyword arguments passed to method fit.

Notes

For other parameters, see documentation for EstimateBFS.

compilestanmodel()

Compile Stan model necessary for method sample. This method is called automatically by obj.sample().

Returns:self
evalgrid(gridsize=1000)

Call method evalgrid for each EstimateBFS object that the self holds, that is, this is equivalent to calling:

obj.bfs0.evalgrid(gridsize)
obj.bfs1.evalgrid(gridsize)
obj.bfsconcat.evalgrid(gridsize)
Parameters:gridsize : size of the grid.
Returns:self
fit(obs0, obs1, nmaxcomp=None, hpp=None, hpgamma=None, hplindley=None, transformation=None, niter=5000, **kwargs)

Configure object model data and automatically calls method sampleposterior (to obtain MCMC samples from the posterior).

Parameters:

obs0 : array like

List/tuple or numpy.array 1D array of observations for model 1.

obs1 : array like

List/tuple or numpy.array 1D array of observations for model 2.

hplindley : float

A priori probability that obs0 and obs1 came from different populations.

**kwargs:

Additional keyword arguments passed to method sample.

Returns:

self

Notes

For other parameters, see documentation for method fit of EstimateBFS.

sampleposterior(niter=5000, nchains=4, njobs=-1, tolrhat=0.02, **kwargs)

Samples from posterior.

Parameters:

niter :

Number of simulations to be draw.

nchains :

Number of MCMC chains.

njobs :

Number of CPUs to be used in parallel. If -1 (default), all CPUs will be used.

tolrhat :

Maximum tolerable distance an Rhat of any sampled parameter can have from 1 (we will resample model approximatedly 10% more iterations until this convergence criteria is met).

**kwargs :

Additional keyword arguments passed to pystan (e.g.: refresh parameter to configure sampler printing status).

Returns:

self

npcompare.fourierseries(x, ncomponents)

Calculate Fourier Series Expansion.

Parameters:

x : 1D numpy.array, list or tuple of numbers to calculate

fourier series expansion

ncomponents : int

number of components of the series

Returns:

2D numpy.array where each line is the Fourier series expansion of

each component of x.