# Package: ci_rvm¶

## Installation¶

The package can be installed via pip. To install the package, you can use

```pip install ci_rvm
```

## Usage¶

The most convenient way to use the algorithm is to use the method find_CI. An example is below.

```# Example for finding profile likelihood confidence intervals for a
# negative binomial model

# We import some packages for convenience
import numpy as np               # for numerical operations
from scipy import stats          # for stats functions
from scipy import optimize as op # to maximize the likelihood

import numdifftools as nd        # to compute gradient and Hessian numerically;
# the package can be found on pypi.
# Another good package for that purpose
# (using automatic differentiation) is autograd

from ci_rvm import find_CI       # to determine confidence intervals

# Define the size of the data set
n = 100

# Define the true parameters
k, p = 5, 0.1

# Generate the data set
data = np.random.negative_binomial(k, p, size=n)

# Because the parameters are constrained to the positive range and the
# interval (0, 1), respectively, we work on a transformed parameter space
# with unbounded domain.
def transform_parameters(params):
k, p = params
return np.exp(k), 1/(1+np.exp(-p))

# Log-Likelihood function for a negative binomial model
def logL(params):
k, p = transform_parameters(params)
return stats.nbinom.logpmf(data, k, p).sum()

# negative log-Likelihood function for optimization (because we use
# minimization algorithms instead of maximization algorithms)
negLogL = lambda params: -logL(params)

# Initial guess
x0 = [0, 0]

# Maximize the likelihood
result = op.minimize(negLogL, x0)

# Print the result (we need to transform the parameters to the original
# parameter space to make them interpretable)
print("The estimate is: k={:5.3f}, p={:5.3f}".format(*transform_parameters(result.x)))

hess = nd.Hessian(logL)

# Find confidence intervals for all parameters.
# Note: For complicated problems, it is worthwile doing this in parallel.
#       However, then we would need to encapsulate the procedure in a
#       method and define the likelihood function, gradient, and Hessian
#       on the top level of the module.
CIs = find_CI(result.x, logL, jac, hess,
disp=True) # the disp argument lets the algorithm print
# status messages.

# CIs is a 2D numpy array with CIs[i, 0] containing the lower bound of the
# confidence interval for the i-th parameter and CIs[i, 1] containing the
# respective upper bound.

# Print the confidence intervals. Note: we need to transform the parameters
# back to the original parameter space.
original_lower = transform_parameters(CIs[:,0])
original_upper = transform_parameters(CIs[:,1])
print("Confidence interval for k: [{:5.3f}, {:5.3f}]".format(
original_lower, original_upper))
print("Confidence interval for n: [{:5.3f}, {:5.3f}]".format(
original_lower, original_upper))
```

## Scientific Publication¶

The theory behind the algorithm implemented in this package is explained in the paper “A robust and efficient algorithm to find profile likelihood confidence intervals” (preprint; accepted in “Statistics and Computing”). Please cite this publication if you have used the package in your own research.