Polya-Gamma

https://github.com/zoj613/polyagamma
Admin
polyagamma is a Python package for efficiently generating samples from the Polya-Gamma distribution...
Tokens:4,795
Snippets:22
Trust Score:6.1
License:BSD-3-Clause
Update:3 weeks ago
Show doc for...
Tokens:
Raw
### Install polyagamma from source

Source: https://github.com/zoj613/polyagamma/blob/main/README.md

Install polyagamma from source. Requires cloning the repository, installing Cython 3.0.x, running cythonize on the .pyx files, and then installing the package.

```shell
$ git clone https://github.com/zoj613/polyagamma.git
$ cd polyagamma/
$ pip install cython==3.0.*
$ cythonize polyagamma/*.pyx
$ pip install .
```

--------------------------------

### Install polyagamma with pip

Source: https://github.com/zoj613/polyagamma/blob/main/README.md

Install the latest version of the polyagamma package from PyPI using pip. The --pre flag allows installation of pre-release versions, and -U upgrades to the latest version.

```shell
$ pip install --pre -U polyagamma
```

--------------------------------

### Generate benchmark plots locally

Source: https://github.com/zoj613/polyagamma/blob/main/README.md

Commands to generate the benchmark plots locally. Requires installing the scripts requirements and running the benchmark script with a size and z value. Runtimes may differ from reported ones depending on the machine.

```shell
$ pip install -r scripts/requirements.txt
$ python scripts/benchmark.py --size=<some size> --z=<z value>
```

--------------------------------

### Install polyagamma with conda

Source: https://github.com/zoj613/polyagamma/blob/main/README.md

Install the polyagamma package from the conda-forge channel using conda.

```shell
$ conda install -c conda-forge polyagamma
```

--------------------------------

### Benchmark random_polyagamma with explicit random_state

Source: https://github.com/zoj613/polyagamma/blob/main/README.md

IPython timing example showing the performance benefit of passing a random generator instance (np.random.SFC64) to random_polyagamma. Passing the generator reduces runtime from 90 µs to 1.69 µs per loop, mitigating overhead from bitgenerator creation, thread lock acquisition, and parameter validation.

```shell
In [3]: rng = np.random.SFC64(1)

In [4]: %timeit random_polyagamma()
90 µs ± 1.65 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [5]: %timeit random_polyagamma(random_state=rng)
1.69 µs ± 6.96 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
```

--------------------------------

### Use low-level Cython function random_polyagamma_fill

Source: https://github.com/zoj613/polyagamma/blob/main/README.md

Example of using the low-level Cython functions provided by the package, specifically `random_polyagamma_fill`. It shows how to obtain a pointer to the underlying bit generator struct from a `numpy.random.Generator` instance, set parameters, and fill a memory view of an array while holding the generator lock and releasing the GIL. The `...` indicates omitted code.

```cython
from cpython.pycapsule cimport PyCapsule_GetPointer
from polyagamma cimport random_polyagamma_fill, DEVROYE
from numpy.random cimport bitgen_t
import numpy as np

# assuming there exists an instance of the Generator class called `rng`.
bitgenerator = rng._bit_generator
# get pointer to the underlying bitgenerator struct
cdef bitgen_t* bitgen = <bitgen_t*>PyCapsule_GetPointer(bitgenerator.capsule, "BitGenerator")
# set distribution parameters
cdef double h = 1, z = 0
# get a memory view of the array to store samples in
cdef double[:] out = np.empty(300)
with bitgenerator.lock, nogil:
    random_polyagamma_fill(bitgen, h, z, DEVROYE, <size_t>out.shape[0], &out[0])
print(out.base)
...
```

--------------------------------

### Compute Polya-Gamma PDF and CDF with polyagamma_pdf and polyagamma_cdf

Source: https://github.com/zoj613/polyagamma/blob/main/README.md

Shows how to compute the probability density function (PDF) and cumulative distribution function (CDF) of the Polya-Gamma distribution. Inputs support broadcasting, and the `return_log` parameter returns log values. The displayed outputs are commented out in the source.

```python
from polyagamma import polyagamma_pdf, polyagamma_cdf

>>> polyagamma_pdf(0.1)
# 3.613955566329298
>>> polyagamma_cdf([1, 2], h=2, z=1)
# array([0.95637847, 0.99963397])
>>> polyagamma_pdf([2, 0.1], h=[[1, 2], [3, 4]], return_log=True)
# array([[   -8.03172733,  -489.17101125]
#        [   -3.82023942, -1987.09156971]])
>>> polyagamma_cdf(4, z=[-100, 0, 2], return_log=True)
# array([ 3.72007598e-44, -3.40628215e-09, -1.25463528e-12])
```

--------------------------------

### Generate Polya-Gamma samples with random_polyagamma

Source: https://github.com/zoj613/polyagamma/blob/main/README.md

Demonstrates the main Python API for generating Polya-Gamma random variates. Shows scalar, array, and broadcasting inputs, the `out` parameter for output arrays, the `method` parameter (e.g., "saddle"), using a `numpy.random.Generator` or its bit generator via `random_state`, and the `disable_checks` option for performance. The default method is a hybrid sampler that chooses the most efficient algorithm based on input values.

```python
import array
import numpy as np
from polyagamma import random_polyagamma

# generate a PG(1, 0) sample
o = random_polyagamma()

# Get a 5 by 1 array of PG(1, 2) variates.
o = random_polyagamma(z=2, size=5)

# We can pass sequences as input. Numpy's broadcasting rules apply here.
# Get a 10 by 2 array where column 1 is PG(2, -10) and column 2 is PG(1, 10)
o = random_polyagamma([2, 1], [-10, 10], size=(10, 2))
z = [[1.5, 2, -0.75, 4, 5],
     [9.5, -8, 7, 6, -0.9]]
o = random_polyagamma(1, z)

# We can pass an output array using the `out` parameter. It does not have to be
# a numpy array. it can be any object that implements the array or buffer protocols.
# As long as its type is 64bit float, contiguous in memory and aligned (e.g. Python's array object).
numpy_out = np.empty(5)
array_out = array.array('d', [0] * 5)
random_polyagamma(out=numpy_out)
print(numpy_out)
random_polyagamma(out=array_out)
print(array_out)

# one can choose a sampling method from {devroye, alternate, gamma, saddle}.
# If not given, the default behaviour is a hybrid sampler that picks the most
# efficient method based on the input values.
o = random_polyagamma(method="saddle")

# one can also use an existing instance of `numpy.random.Generator` as a parameter.
# This is useful to reproduce samples generated via a given seed.
rng = np.random.default_rng(12345)
o = random_polyagamma(random_state=rng)

# If one is using a `numpy.random.RandomState` instance instead of the `Generator`
# class, the object's underlying bitgenerator can be passed as the value of random_state
bit_gen = np.random.RandomState(12345)._bit_generator
o = random_polyagamma(random_state=bit_gen)

# When passing a large input array for the shape parameter `h`, parameter value
# validation checks can be disabled if the values are guaranteed to be positive
# to avoid some overhead, which may boost performance.
large_h = np.ones(1000000)
o = random_polyagamma(large_h, disable_checks=True)
```
Tokens:
Raw