cardinal_pythonlib.rpm

Original code copyright (C) 2009-2022 Rudolf Cardinal (rudolf@pobox.com).

This file is part of cardinal_pythonlib.

Licensed under the Apache License, Version 2.0 (the “License”); you may not use this file except in compliance with the License. You may obtain a copy of the License at

Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an “AS IS” BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

Randomized probability matching (RPM).

As per:

  • Scott SL. A modern Bayesian look at the multi-armed bandit. Applied Stochastic Models in Business and Industry 26 (2010): 639–58. https://doi.org/10.1002/asmb.874.

An R version is in rpm.R within https://github.com/rudolfcardinal/rlib.

cardinal_pythonlib.rpm.bd0[source]

//github.com/atks/Rmath/blob/master/bd0.c.

Type:Per https
cardinal_pythonlib.rpm.beta_cdf_fast

Original license:

/*
 * zlib License
 *
 * Regularized Incomplete Beta Function
 *
 * Copyright (c) 2016, 2017 Lewis Van Winkle
 * https://CodePlea.com
 *
 * This software is provided 'as-is', without any express or implied
 * warranty. In no event will the authors be held liable for any damages
 * arising from the use of this software.
 *
 * Permission is granted to anyone to use this software for any purpose,
 * including commercial applications, and to alter it and redistribute it
 * freely, subject to the following restrictions:
 *
 * 1. The origin of this software must not be misrepresented; you must not
 *    claim that you wrote the original software. If you use this software
 *    in a product, an acknowledgement in the product documentation would be
 *    appreciated but is not required.
 * 2. Altered source versions must be plainly marked as such, and must not be
 *    misrepresented as being the original software.
 * 3. This notice may not be removed or altered from any source distribution.
 */
cardinal_pythonlib.rpm.beta_pdf_fast[source]

Beta probability distribution. Translated and adapted from https://en.wikipedia.org/wiki/Beta_distribution, but calculated in the log domain.

In R: plot(function(x) dbeta(x, shape1 = a, shape2 = b))

See https://github.com/SurajGupta/r-source/blob/master/src/nmath/dbeta.c. - For lower.tail = TRUE and log.p = FALSE (the defaults), R_DT_0 means 0. - For log.p = FALSE (the default), R_D_val(x) means x.

Original license:

/*
 *  AUTHOR
 *    Catherine Loader, catherine@research.bell-labs.com.
 *    October 23, 2000.
 *
 *  Merge in to R:
 *      Copyright (C) 2000, The R Core Team
 *  Changes to case a, b < 2, use logs to avoid underflow
 *      Copyright (C) 2006-2014 The R Core Team
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 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 General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, a copy is available at
 *  https://www.R-project.org/Licenses/
 *
 *
 *  DESCRIPTION
 *    Beta density,
 *                   (a+b-1)!     a-1       b-1
 *      p(x;a,b) = ------------ x     (1-x)
 *                 (a-1)!(b-1)!
 *
 *               = (a+b-1) dbinom(a-1; a+b-2,x)
 *
 *    The basic formula for the log density is thus
 *    (a-1) log x + (b-1) log (1-x) - lbeta(a, b)
 *    If either a or b <= 2 then 0 < lbeta(a, b) < 710 and so no
 *    term is large.  We use Loader's code only if both a and b > 2.
 */
cardinal_pythonlib.rpm.dbinom_raw_log[source]

Translated and adapted from https://github.com/atks/Rmath/blob/master/dbinom.c – the version where give_log is TRUE, for which:

  • R_D_exp(x) translates to x
  • R_D__0 translates to -inf
  • R_D__1 translates to 0

Original license:

/*
 * AUTHOR
 *   Catherine Loader, catherine@research.bell-labs.com.
 *   October 23, 2000.
 *
 *  Merge in to R and further tweaks :
 *      Copyright (C) 2000, The R Core Team
 *      Copyright (C) 2008, The R Foundation
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 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 General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, a copy is available at
 *  https://www.r-project.org/Licenses/
 *
 *
 * DESCRIPTION
 *
 *   To compute the binomial probability, call dbinom(x,n,p).
 *   This checks for argument validity, and calls dbinom_raw().
 *
 *   dbinom_raw() does the actual computation; note this is called by
 *   other functions in addition to dbinom().
 *     (1) dbinom_raw() has both p and q arguments, when one may be represented
 *         more accurately than the other (in particular, in df()).
 *     (2) dbinom_raw() does NOT check that inputs x and n are integers. This
 *         should be done in the calling function, where necessary.
 *         -- but is not the case at all when called e.g., from df() or dbeta() !
 *     (3) Also does not check for 0 <= p <= 1 and 0 <= q <= 1 or NaN's.
 *         Do this in the calling function.
 */
cardinal_pythonlib.rpm.dummy_jit_integrand_function_with_args(integrand_function: Callable[[numpy.ndarray], float]) → Callable[[...], float][source]

Dummy version of jit_integrand_function_with_args, for debugging. Use this instead if, for example, you want to be able to use the Python logger.

cardinal_pythonlib.rpm.incbeta[source]

Original license:

/*
 * zlib License
 *
 * Regularized Incomplete Beta Function
 *
 * Copyright (c) 2016, 2017 Lewis Van Winkle
 * https://CodePlea.com
 *
 * This software is provided 'as-is', without any express or implied
 * warranty. In no event will the authors be held liable for any damages
 * arising from the use of this software.
 *
 * Permission is granted to anyone to use this software for any purpose,
 * including commercial applications, and to alter it and redistribute it
 * freely, subject to the following restrictions:
 *
 * 1. The origin of this software must not be misrepresented; you must not
 *    claim that you wrote the original software. If you use this software
 *    in a product, an acknowledgement in the product documentation would be
 *    appreciated but is not required.
 * 2. Altered source versions must be plainly marked as such, and must not be
 *    misrepresented as being the original software.
 * 3. This notice may not be removed or altered from any source distribution.
 */
cardinal_pythonlib.rpm.jit_integrand_function_with_args(integrand_function: Callable[[numpy.ndarray], float]) → scipy._lib._ccallback.LowLevelCallable[source]

Decorator to wrap a function that will be integrated by Scipy. See https://stackoverflow.com/questions/51109429/.

carray:

cfunc:

CPointer:

  • Type class for pointers to other types.
  • Syntax: CPointer(dtype, addrspace=None).

Scipy’s quad() accepts either a Python function or a C callback wrapped in a ctypes callback object. This decorator converts the former to the latter.

Specifically (https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html), quad accepts a scipy.LowLevelCallable with one of these signatures:

double func(double x);
double func(double x, void *user_data);
double func(int n, double *xx);  // THIS ONE.
double func(int n, double *xx, void *user_data);

NOTE: “In the call forms with xx, n is the length of the xx array which contains xx[0] == x and the rest of the items are numbers contained in the args argument of quad.”

See specimen use below.

cardinal_pythonlib.rpm.rpm_probabilities_successes_failures(n_successes: Union[numpy.typing._array_like._SupportsArray[numpy.dtype][numpy.dtype], numpy.typing._nested_sequence._NestedSequence[numpy.typing._array_like._SupportsArray[numpy.dtype]][numpy.typing._array_like._SupportsArray[numpy.dtype][numpy.dtype]], bool, int, float, complex, str, bytes, numpy.typing._nested_sequence._NestedSequence[typing.Union[bool, int, float, complex, str, bytes]][Union[bool, int, float, complex, str, bytes]]], n_failures: Union[numpy.typing._array_like._SupportsArray[numpy.dtype][numpy.dtype], numpy.typing._nested_sequence._NestedSequence[numpy.typing._array_like._SupportsArray[numpy.dtype]][numpy.typing._array_like._SupportsArray[numpy.dtype][numpy.dtype]], bool, int, float, complex, str, bytes, numpy.typing._nested_sequence._NestedSequence[typing.Union[bool, int, float, complex, str, bytes]][Union[bool, int, float, complex, str, bytes]]]) → numpy.ndarray[source]

Calculate the optimal choice probabilities.

Note that Scott’s original R version, compute.probopt(), on Figure 3 (p648) has arguments y (number of successes) and n (number of trials, NOT the number of failures).

cardinal_pythonlib.rpm.rpm_probabilities_successes_failures_fast(n_successes: Union[numpy.typing._array_like._SupportsArray[numpy.dtype][numpy.dtype], numpy.typing._nested_sequence._NestedSequence[numpy.typing._array_like._SupportsArray[numpy.dtype]][numpy.typing._array_like._SupportsArray[numpy.dtype][numpy.dtype]], bool, int, float, complex, str, bytes, numpy.typing._nested_sequence._NestedSequence[typing.Union[bool, int, float, complex, str, bytes]][Union[bool, int, float, complex, str, bytes]]], n_failures: Union[numpy.typing._array_like._SupportsArray[numpy.dtype][numpy.dtype], numpy.typing._nested_sequence._NestedSequence[numpy.typing._array_like._SupportsArray[numpy.dtype]][numpy.typing._array_like._SupportsArray[numpy.dtype][numpy.dtype]], bool, int, float, complex, str, bytes, numpy.typing._nested_sequence._NestedSequence[typing.Union[bool, int, float, complex, str, bytes]][Union[bool, int, float, complex, str, bytes]]]) → numpy.ndarray[source]

Fast version of rpm_probabilities_successes_failures().

cardinal_pythonlib.rpm.rpm_probabilities_successes_failures_twochoice_fast(n_success_this: int, n_failure_this: int, n_success_other: int, n_failure_other: int) → float[source]

Calculate the optimal choice probability for the first of two options in two-choice RPM.

Curiously, copying rpm_probabilities_successes_failures from this library to user code, essentially unmodified, stopped a memory explosion (it looks like scipy is playing with docstrings, maybe?).

It was still very slow and that relates to quad().

Here we (a) only calculate one action (50% faster just from that!) and (b) use functions optimized using numba.jit.

Massively tedious optimization (translation from R’s C code to Python) but it works very well.

cardinal_pythonlib.rpm.rpm_probabilities_successes_totals(n_successes: Union[numpy.typing._array_like._SupportsArray[numpy.dtype][numpy.dtype], numpy.typing._nested_sequence._NestedSequence[numpy.typing._array_like._SupportsArray[numpy.dtype]][numpy.typing._array_like._SupportsArray[numpy.dtype][numpy.dtype]], bool, int, float, complex, str, bytes, numpy.typing._nested_sequence._NestedSequence[typing.Union[bool, int, float, complex, str, bytes]][Union[bool, int, float, complex, str, bytes]]], n_total: Union[numpy.typing._array_like._SupportsArray[numpy.dtype][numpy.dtype], numpy.typing._nested_sequence._NestedSequence[numpy.typing._array_like._SupportsArray[numpy.dtype]][numpy.typing._array_like._SupportsArray[numpy.dtype][numpy.dtype]], bool, int, float, complex, str, bytes, numpy.typing._nested_sequence._NestedSequence[typing.Union[bool, int, float, complex, str, bytes]][Union[bool, int, float, complex, str, bytes]]]) → numpy.ndarray[source]

Randomized probability matching (RPM).

Parameters:
  • n_successes – Number of successes (per option).
  • n_total – Total number of trials (per option).
Returns:

Optimal choice probabilities (per option) via RPM.
cardinal_pythonlib.rpm.stirlerr[source]

Stirling expansion error. Translated and adapted from https://github.com/atks/Rmath/blob/master/stirlerr.c

Original license:

/*
 *  AUTHOR
 *    Catherine Loader, catherine@research.bell-labs.com.
 *    October 23, 2000.
 *
 *  Merge in to R:
 *      Copyright (C) 2000, The R Core Team
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 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 General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, a copy is available at
 *  https://www.r-project.org/Licenses/
 *
 *
 *  DESCRIPTION
 *
 *    Computes the log of the error term in Stirling's formula.
 *      For n > 15, uses the series 1/12n - 1/360n^3 + ...
 *      For n <=15, integers or half-integers, uses stored values.
 *      For other n < 15, uses lgamma directly (don't use this to
 *        write lgamma!)
 *
 * Merge in to R:
 * Copyright (C) 2000, The R Core Team
 * R has lgammafn, and lgamma is not part of ISO C
 */