Skip to content

Instantly share code, notes, and snippets.

@leandrolcampos

leandrolcampos/betaincinv.ipynb Secret

Created August 2, 2022 23:24
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save leandrolcampos/ccb74333a424613d9337349666dbec33 to your computer and use it in GitHub Desktop.
Save leandrolcampos/ccb74333a424613d9337349666dbec33 to your computer and use it in GitHub Desktop.
Download ZIP
betaincinv.ipynb
Display the source blob
Display the rendered blob
Raw
betaincinv.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "betaincinv.ipynb",
"provenance": [],
"collapsed_sections": [
"qCxzCA-zR7Bq",
"3LrNhs84Xkyd",
"GC-0lYG0HYAp",
"d15s3qJNgkAc",
"Jyts6uECw7Cp"
],
"authorship_tag": "ABX9TyNxK4CZLODS5aW7aRGOV0oR",
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
},
"language_info": {
"name": "python"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/leandrolcampos/ccb74333a424613d9337349666dbec33/betaincinv.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "code",
"source": [
"!pip install -U -q tensorflow"
],
"metadata": {
"id": "5xAGAPPG8wcL"
},
"execution_count": 1,
"outputs": []
},
{
"cell_type": "code",
"source": [
"!pip install -U -q tensorflow-probability"
],
"metadata": {
"id": "YgA81jppUk2T"
},
"execution_count": 2,
"outputs": []
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"id": "KA-hrzh_7myi"
},
"outputs": [],
"source": [
"import functools\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"import scipy.special as sp_special\n",
"import tensorflow as tf\n",
"import tensorflow_probability as tfp\n",
"\n",
"from tensorflow_probability.python import math as tfp_math\n",
"from tensorflow_probability.python.internal import dtype_util\n",
"from tensorflow_probability.python.internal import prefer_static as ps"
]
},
{
"cell_type": "markdown",
"source": [
"## My implementation of `betaincinv`"
],
"metadata": {
"id": "qCxzCA-zR7Bq"
}
},
{
"cell_type": "code",
"source": [
"# The implementation of the inverse of the regularized incomplete beta function\n",
"# is based on ideas and equations available in the following references:\n",
"# [1] Milton Abramowitz and Irene A. Stegun\n",
"# Handbook of Mathematical Functions with Formulas, Graphs, and\n",
"# Mathematical Tables.\n",
"# US Government Printing Office, 1964 (reprinted 1972).\n",
"# https://archive.org/details/AandS-mono600\n",
"# [2] William Press, Saul Teukolsky, William Vetterling and Brian Flannery\n",
"# Numerical Recipes: The Art of Scientific Computing\n",
"# Cambridge University Press, 2007 (Third Edition)\n",
"# http://numerical.recipes/book/book.html\n",
"# [3] John Maddock, Paul A. Bristow, et al.\n",
"# The Incomplete Beta Function Inverses\n",
"# https://www.boost.org/doc/libs/1_79_0/libs/math/doc/html/special.html\n",
"# [4] Stephen L. Moshier\n",
"# Cephes Mathematical Library\n",
"# https://netlib.org/cephes/"
],
"metadata": {
"id": "U9zirLQRbQB2"
},
"execution_count": 4,
"outputs": []
},
{
"cell_type": "code",
"source": [
"def _betaincinv_initial_approx(a, b, y, dtype):\n",
" \"\"\"Computes an initial approximation for `betaincinv(a, b, y)`.\"\"\"\n",
" numpy_dtype = dtype_util.as_numpy_dtype(dtype)\n",
" tiny = tf.constant(np.finfo(numpy_dtype).tiny, dtype=dtype)\n",
" eps = tf.constant(np.finfo(numpy_dtype).eps, dtype=dtype)\n",
" one = tf.constant(1., dtype=dtype)\n",
" two = tf.constant(2., dtype=dtype)\n",
" three = tf.constant(3., dtype=dtype)\n",
" five = tf.constant(5., dtype=dtype)\n",
" six = tf.constant(6., dtype=dtype)\n",
" min_log = tf.math.log(\n",
" tf.constant(2. ** np.finfo(numpy_dtype).minexp, dtype))\n",
" max_log = tf.math.log(\n",
" tf.constant(2. ** (np.finfo(numpy_dtype).maxexp - 1.), dtype))\n",
"\n",
" # When min(a, b) >= 1, we use the approximation proposed by [1].\n",
"\n",
" # Equation 26.5.22 [1, page 945].\n",
" yp = -tf.math.ndtri(y)\n",
" inv_2a_minus_one = tf.math.reciprocal(two * a - one)\n",
" inv_2b_minus_one = tf.math.reciprocal(two * b - one)\n",
" lmb = (tf.math.square(yp) - three) / six\n",
" h = two * tf.math.reciprocal(inv_2a_minus_one + inv_2b_minus_one)\n",
" w = (yp * tf.math.sqrt(h + lmb) / h - \n",
" (inv_2b_minus_one - inv_2a_minus_one) * \n",
" (lmb + five / six - two / (three * h)))\n",
" result_for_large_a_and_b = a / (a + b * tf.math.exp(two * w))\n",
"\n",
" # When min(a, b) < 1 and max(a, b) >= 1, we use the approximation proposed by\n",
" # [2]. This approximation depends on the following approximation for betainc:\n",
" # betainc(a, b, x) ~=\n",
" # x ** a / (integral_approx * a) , when x <= mean ,\n",
" # (1 - x) ** b / (integral_approx * b) , when x > mean ,\n",
" # where:\n",
" # integral_approx = (mean ** a) / a + (mean_complement ** b) / b ,\n",
" # mean = a / (a + b) ,\n",
" # mean_complement = 1 - mean = b / (a + b) .\n",
" # We invert betainc(a, b, x) with respect to x in the proper regime.\n",
"\n",
" # Equation 6.4.7 [2, page 271].\n",
" a_plus_b = a + b\n",
" mean = a / a_plus_b\n",
" mean_complement = b / a_plus_b\n",
" integral_approx_part_a = tf.math.exp(tf.math.xlogy(a, mean) - tf.math.log(a))\n",
" integral_approx_part_b = tf.math.exp(tf.math.xlogy(b, mean_complement) -\n",
" tf.math.log(b))\n",
" integral_approx = integral_approx_part_a + integral_approx_part_b\n",
"\n",
" # Solve Equation 6.4.8 [2, page 271] for x in the respective regimes.\n",
" inv_a = tf.math.reciprocal(a)\n",
" inv_b = tf.math.reciprocal(b)\n",
" result_for_small_a_or_b = tf.where(\n",
" tf.math.less_equal(y, (integral_approx_part_a / integral_approx)),\n",
" tf.math.exp(tf.math.xlogy(inv_a, y) + tf.math.xlogy(inv_a, a) +\n",
" tf.math.xlogy(inv_a, integral_approx)),\n",
" one - tf.math.exp(tf.math.xlog1py(inv_b, -y) + tf.math.xlogy(inv_b, b) +\n",
" tf.math.xlogy(inv_b, integral_approx)))\n",
"\n",
" # And when max(a, b) < 1, we use the approximation proposed by [3] for the\n",
" # same domain:\n",
" # betaincinv(a, b, y) ~= xg / (1 + xg) ,\n",
" # where:\n",
" # xg = (a * y * Beta(a, b)) ** (1 / a) .\n",
" log_xg = tf.math.xlogy(inv_a, a) + tf.math.xlogy(inv_a, y) + (\n",
" inv_a * tfp_math.lbeta(a, b))\n",
" xg = tf.math.exp(tf.clip_by_value(log_xg, min_log, max_log))\n",
" result_for_small_a_and_b = xg / (one + xg)\n",
"\n",
" # Return the appropriate result for parameters a and b.\n",
" result = tf.where(\n",
" tf.math.greater_equal(tf.math.minimum(a, b), one),\n",
" result_for_large_a_and_b,\n",
" tf.where(\n",
" tf.math.maximum(a, b) < one,\n",
" result_for_small_a_and_b,\n",
" result_for_small_a_or_b))\n",
" \n",
" return tf.clip_by_value(result, tiny, one - eps)"
],
"metadata": {
"id": "SZeloSmIz52v"
},
"execution_count": 5,
"outputs": []
},
{
"cell_type": "code",
"source": [
"def _betaincinv_computation(a, b, y):\n",
" \"\"\"Returns the inverse of `betainc(a, b, x)` with respect to `x`.\"\"\"\n",
" dtype = dtype_util.common_dtype([a, b, y], tf.float32)\n",
" numpy_dtype = dtype_util.as_numpy_dtype(dtype)\n",
" zero = tf.constant(0., dtype=dtype)\n",
" tiny = tf.constant(np.finfo(numpy_dtype).tiny, dtype=dtype)\n",
" eps = tf.constant(np.finfo(numpy_dtype).eps, dtype=dtype)\n",
" half = tf.constant(0.5, dtype=dtype)\n",
" one = tf.constant(1., dtype=dtype)\n",
" two = tf.constant(2., dtype=dtype)\n",
" halley_correction_min = tf.constant(0.5, dtype=dtype)\n",
" halley_correction_max = tf.constant(1.5, dtype=dtype)\n",
"\n",
" a = tf.convert_to_tensor(a, dtype=dtype)\n",
" b = tf.convert_to_tensor(b, dtype=dtype)\n",
" y = tf.convert_to_tensor(y, dtype=dtype)\n",
"\n",
" broadcast_shape = functools.reduce(\n",
" ps.broadcast_shape, [ps.shape(a), ps.shape(b), ps.shape(y)])\n",
"\n",
" a = tf.broadcast_to(a, broadcast_shape)\n",
" b = tf.broadcast_to(b, broadcast_shape)\n",
" y = tf.broadcast_to(y, broadcast_shape)\n",
"\n",
" # When tfp_math.betainc(a, b, 0.5) < y, we apply the symmetry relation given\n",
" # here: https://dlmf.nist.gov/8.17.E4\n",
" # betainc(a, b, x) = 1 - betainc(b, a, 1 - x) .\n",
" # If dtype != float64, we have additional conditions to apply this relation:\n",
" # (a < 1) & (b < 1) & (tfp_math.betainc(a, b, a / (a + b)) < y) .\n",
" error_at_half = tfp_math.betainc(a, b, half) - y\n",
" if numpy_dtype == np.float64:\n",
" use_symmetry_relation = (error_at_half < zero)\n",
" else:\n",
" a_and_b_are_small = (a < one) & (b < one)\n",
" error_at_mean = tfp_math.betainc(a, b, a / (a + b)) - y\n",
" use_symmetry_relation = (error_at_half < zero) & a_and_b_are_small & (\n",
" error_at_mean < zero)\n",
"\n",
" a_orig = a\n",
" a = tf.where(use_symmetry_relation, b, a)\n",
" b = tf.where(use_symmetry_relation, a_orig, b)\n",
" y = tf.where(use_symmetry_relation, one - y, y)\n",
"\n",
" a_minus_1 = a - one\n",
" b_minus_1 = b - one\n",
" one_minus_eps = one - eps\n",
"\n",
" # tolerance was set by experimentation and max_iterations was taken from [4].\n",
" if numpy_dtype == np.float64:\n",
" tolerance = 1e-12\n",
" max_iterations = 8\n",
" else:\n",
" tolerance = 1e-6\n",
" max_iterations = 10\n",
"\n",
" def root_finding_iteration(should_stop, low, high, candidate):\n",
" error = tfp_math.betainc(a, b, candidate) - y\n",
" error_over_der = error / tf.math.exp(\n",
" tf.math.xlogy(a_minus_1, candidate) +\n",
" tf.math.xlog1py(b_minus_1, -candidate) -\n",
" tfp_math.lbeta(a, b))\n",
" second_der_over_der = a_minus_1 / candidate - b_minus_1 / (one - candidate)\n",
" # Following [2, section 9.4.2, page 463], we limit the influence of the\n",
" # Halley's correction to the Newton's method, since this correction can\n",
" # reduce the Newton's region of convergence. We set minimum and maximum\n",
" # values for this correction by experimentation.\n",
" halley_correction = tf.clip_by_value(\n",
" one - half * error_over_der * second_der_over_der,\n",
" halley_correction_min,\n",
" halley_correction_max)\n",
" halley_delta = error_over_der / halley_correction\n",
" halley_candidate = tf.where(\n",
" should_stop, candidate, candidate - halley_delta)\n",
"\n",
" # Fall back to bisection if the current step would take the new candidate\n",
" # out of bounds.\n",
" new_candidate = tf.where(\n",
" tf.less_equal(halley_candidate, low),\n",
" half * (candidate + low),\n",
" tf.where(\n",
" tf.greater_equal(halley_candidate, high),\n",
" half * (candidate + high),\n",
" halley_candidate))\n",
"\n",
" new_delta = candidate - new_candidate\n",
" new_delta_is_negative = (new_delta < zero)\n",
" new_low = tf.where(new_delta_is_negative, candidate, low)\n",
" new_high = tf.where(new_delta_is_negative, high, candidate)\n",
"\n",
" adjusted_tolerance = tf.maximum(tolerance * new_candidate, two * tiny)\n",
" should_stop = (should_stop | (tf.math.abs(new_delta) < adjusted_tolerance) |\n",
" tf.math.equal(new_low, new_high))\n",
"\n",
" return should_stop, new_low, new_high, new_candidate\n",
"\n",
" initial_candidate = _betaincinv_initial_approx(a, b, y, dtype)\n",
" # Bracket the solution with the interval (low, high).\n",
" initial_low = tf.zeros_like(y)\n",
" if numpy_dtype == np.float64:\n",
" initial_high = tf.ones_like(y) * half\n",
" else:\n",
" initial_high = tf.ones_like(y) * tf.where(\n",
" a_and_b_are_small & (error_at_mean < zero), half, one)\n",
"\n",
" (_, _, _, result) = tf.while_loop(\n",
" cond=lambda stop, *_: tf.reduce_any(~stop),\n",
" body=root_finding_iteration,\n",
" loop_vars=(\n",
" tf.equal(y, initial_low) | tf.equal(y, initial_high),\n",
" initial_low,\n",
" initial_high,\n",
" initial_candidate),\n",
" maximum_iterations=max_iterations)\n",
"\n",
" # If we are taking advantage of the symmetry relation, we have to adjust the\n",
" # input y and the solution.\n",
" y = tf.where(use_symmetry_relation, one - y, y)\n",
" result = tf.where(use_symmetry_relation, one - result, result)\n",
" result = tf.clip_by_value(result, tiny, one_minus_eps)\n",
"\n",
" # Handle trivial cases.\n",
" result = tf.where(tf.equal(y, zero) | tf.equal(y, one), y, result)\n",
"\n",
" # Determine if the inputs are out of range (should return NaN output).\n",
" result_is_nan = (a <= zero) | (b <= zero) | (y < zero) | (y > one)\n",
" result = tf.where(result_is_nan, numpy_dtype(np.nan), result)\n",
"\n",
" return result"
],
"metadata": {
"id": "b4-v9FBS1w_z"
},
"execution_count": 6,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## Sanity checks for `float64`"
],
"metadata": {
"id": "3LrNhs84Xkyd"
}
},
{
"cell_type": "code",
"source": [
"dtype = np.float64"
],
"metadata": {
"id": "7IP4aGZTy2Fd"
},
"execution_count": 7,
"outputs": []
},
{
"cell_type": "code",
"source": [
"a = 1.2\n",
"b = 9.8\n",
"y = 0.05"
],
"metadata": {
"id": "2E1wtw0Wy2Fk"
},
"execution_count": 8,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "c4fd8157-345d-4d3f-e3da-a513845418dd",
"id": "g45UbCsby2Fl"
},
"execution_count": 9,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.10909090909090909"
]
},
"metadata": {},
"execution_count": 9
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "87e4ebd0-79cf-4325-cde0-c0204cfdb2b3",
"id": "t6lfHBXoy2Fl"
},
"execution_count": 10,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"True"
]
},
"metadata": {},
"execution_count": 10
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "1YCHfznoy2Fl"
},
"execution_count": 11,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "71a2ad32-006e-43f8-c5ec-075f135cc759",
"id": "W0-qLGZ6y2Fl"
},
"execution_count": 12,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.009366299700584217>"
]
},
"metadata": {},
"execution_count": 12
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "4882ab3d-ab6e-4656-e0cc-ba22f8c9be6c",
"id": "BQb-h3gyy2Fm"
},
"execution_count": 13,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.009366299700584236"
]
},
"metadata": {},
"execution_count": 13
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "8v0cTM7RVpv2",
"outputId": "bbeb5674-8766-499e-8b60-09eda2623198"
},
"execution_count": 14,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"2.037299557535471e-15"
]
},
"metadata": {},
"execution_count": 14
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "cdh5nynJVyCc",
"outputId": "c62fa193-c8cb-4784-c3e2-acc6e8d99448"
},
"execution_count": 15,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.3877787807814457e-16"
]
},
"metadata": {},
"execution_count": 15
}
]
},
{
"cell_type": "code",
"source": [
"a = 1.2\n",
"b = 9.8\n",
"y = 0.9"
],
"metadata": {
"id": "ObMxZkOby2Fm"
},
"execution_count": 16,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "a8968d6f-611c-4a24-b358-2626d6f8f4a6",
"id": "HaIJB67ey2Fm"
},
"execution_count": 17,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.10909090909090909"
]
},
"metadata": {},
"execution_count": 17
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "aab76b84-1bcd-44df-fe7b-d3a2e6b4478a",
"id": "265FgAWgy2Fn"
},
"execution_count": 18,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"False"
]
},
"metadata": {},
"execution_count": 18
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "PfxulgB6y2Fn"
},
"execution_count": 19,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "63e11b82-643e-4197-ce5b-e9760721559e",
"id": "cNDG2cnzWae0"
},
"execution_count": 20,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.23426865827003515>"
]
},
"metadata": {},
"execution_count": 20
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "e37d867b-8dc5-4d1b-8b81-0959a058e5ff",
"id": "EHHIcm4LWae1"
},
"execution_count": 21,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.234268658270035"
]
},
"metadata": {},
"execution_count": 21
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "56d8febc-5396-4f9c-b411-bd37aa775cd5",
"id": "AjkIs724Wae1"
},
"execution_count": 22,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"5.92387727419257e-16"
]
},
"metadata": {},
"execution_count": 22
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "2bbf195c-526f-44d9-ff7e-893d1cbbc777",
"id": "_SCum_kaWae1"
},
"execution_count": 23,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.0"
]
},
"metadata": {},
"execution_count": 23
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.8\n",
"b = 1.7\n",
"y = 0.1"
],
"metadata": {
"id": "sJTSQzgUy2Fn"
},
"execution_count": 24,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "a851ee77-12eb-4354-d332-a16d86bc8fdb",
"id": "t8U-mXpZy2Fo"
},
"execution_count": 25,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.32"
]
},
"metadata": {},
"execution_count": 25
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "a6d3576d-e377-4b8c-9940-d42224730aca",
"id": "uSBO7tTmy2Fo"
},
"execution_count": 26,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"True"
]
},
"metadata": {},
"execution_count": 26
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "mAOvYqbjy2Fo"
},
"execution_count": 27,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "0bc4edd1-b69f-4a36-fb71-f9089dece728",
"id": "Uisnjop7WkAR"
},
"execution_count": 28,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.03238692564298561>"
]
},
"metadata": {},
"execution_count": 28
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "4b5b4238-cff7-4450-a04e-1bd810539248",
"id": "Ug4FaNTBWkAR"
},
"execution_count": 29,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.03238692564298564"
]
},
"metadata": {},
"execution_count": 29
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "4322307f-9e16-4c6a-ba5c-4d26623e66c4",
"id": "p-AM6KFgWkAR"
},
"execution_count": 30,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"8.569993929522673e-16"
]
},
"metadata": {},
"execution_count": 30
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "4bf2ab3a-2cb8-4763-a811-e2772166ad2e",
"id": "9dlAvTOVWkAR"
},
"execution_count": 31,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"2.7755575615628914e-16"
]
},
"metadata": {},
"execution_count": 31
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.8\n",
"b = 1.7\n",
"y = 0.9"
],
"metadata": {
"id": "I5AQ2YfJy2Fp"
},
"execution_count": 32,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "bdb38003-c4bb-4ad2-9c31-d912ddb52a2b",
"id": "XT08_3aRy2Fp"
},
"execution_count": 33,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.32"
]
},
"metadata": {},
"execution_count": 33
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "f011bfb8-2f25-4bbe-e8ec-5d59123b374e",
"id": "nyNlT5KDy2Fp"
},
"execution_count": 34,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"False"
]
},
"metadata": {},
"execution_count": 34
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "8a7S0b_fy2Fp"
},
"execution_count": 35,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "50019658-dbca-432c-8d09-90f5f80e3499",
"id": "ooj438JZWptP"
},
"execution_count": 36,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.6994021675222115>"
]
},
"metadata": {},
"execution_count": 36
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "44807b41-9242-4626-b0c3-0f10dc45f238",
"id": "q3SCpje4WptQ"
},
"execution_count": 37,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.6994021675222115"
]
},
"metadata": {},
"execution_count": 37
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "33ee3f32-aa88-4423-9dd9-c4dc4f64bc2a",
"id": "2msnvY5gWptQ"
},
"execution_count": 38,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.0"
]
},
"metadata": {},
"execution_count": 38
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "34e312e1-a280-47d4-8559-b2e4fc04abb8",
"id": "QCwij7x3WptQ"
},
"execution_count": 39,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.0"
]
},
"metadata": {},
"execution_count": 39
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.1\n",
"b = 0.7\n",
"y = 0.1"
],
"metadata": {
"id": "s1-PQr3iBYSr"
},
"execution_count": 40,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "d990d89c-4085-4964-f6b2-f81e19d0c475",
"id": "3RZacerABYSs"
},
"execution_count": 41,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.12500000000000003"
]
},
"metadata": {},
"execution_count": 41
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "40233af5-051b-42e2-f0e9-05c0d3dbabe6",
"id": "2RJS-6JQBYSs"
},
"execution_count": 42,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"True"
]
},
"metadata": {},
"execution_count": 42
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "nw-sTDuvBYSt"
},
"execution_count": 43,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "9ac8b058-61d4-4417-d1e7-0a14dc44a621",
"id": "fmd4SfYZWui7"
},
"execution_count": 44,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=1.8028184911194362e-10>"
]
},
"metadata": {},
"execution_count": 44
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "513470ce-9fd4-4d1e-940f-4aa0f27f9289",
"id": "2BkRTJn0Wui7"
},
"execution_count": 45,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.8028184911194486e-10"
]
},
"metadata": {},
"execution_count": 45
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "04a44bca-9354-4fdc-ee2f-7914fa1a96d0",
"id": "lVTS9l2BWui7"
},
"execution_count": 46,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"6.882395121536016e-15"
]
},
"metadata": {},
"execution_count": 46
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "4bf4dca6-0e31-46b0-bdd2-623f6b7777e6",
"id": "kbdiMK6CWui7"
},
"execution_count": 47,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.1102230246251565e-15"
]
},
"metadata": {},
"execution_count": 47
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.1\n",
"b = 0.7\n",
"y = 0.2"
],
"metadata": {
"id": "eccvxwbaBYSt"
},
"execution_count": 48,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "08da3c91-cca7-4084-89ad-b5a5f5e356be",
"id": "IZQ9ydE_BYSu"
},
"execution_count": 49,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.12500000000000003"
]
},
"metadata": {},
"execution_count": 49
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "f9abaf26-1d2a-43a9-8198-7d324b9166ee",
"id": "DtCVX5BABYSu"
},
"execution_count": 50,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"False"
]
},
"metadata": {},
"execution_count": 50
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "Pn95B27SBYSu"
},
"execution_count": 51,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "6013d606-0498-41e3-c5db-95cc9795c012",
"id": "jxrsUCGcWzCq"
},
"execution_count": 52,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=1.846086042050682e-07>"
]
},
"metadata": {},
"execution_count": 52
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "2a56287c-27f2-4f42-c9a4-c5a7b9fb4198",
"id": "QkY6eWoAWzCq"
},
"execution_count": 53,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.8460860420506996e-07"
]
},
"metadata": {},
"execution_count": 53
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "50817412-ca54-46ed-abf3-e13f895f0509",
"id": "UwuT_SrgWzCq"
},
"execution_count": 54,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"9.606677007013447e-15"
]
},
"metadata": {},
"execution_count": 54
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "796f5e70-8695-406e-ae23-5cb26fe74dcc",
"id": "6qbuEvJNWzCq"
},
"execution_count": 55,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"2.7755575615628914e-16"
]
},
"metadata": {},
"execution_count": 55
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.1\n",
"b = 0.7\n",
"y = 0.9"
],
"metadata": {
"id": "0zsIbmFQK4lG"
},
"execution_count": 56,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "115f48bf-9c59-4173-d09d-686b62078b82",
"id": "HoHUslCSK4lH"
},
"execution_count": 57,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.12500000000000003"
]
},
"metadata": {},
"execution_count": 57
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "a977cd0d-414f-49db-a618-e4076c352351",
"id": "l7bgK-w6K4lI"
},
"execution_count": 58,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"False"
]
},
"metadata": {},
"execution_count": 58
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "UjnltRCEK4lI"
},
"execution_count": 59,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "d60fd9e4-2a27-47cb-de4b-25b086bff038",
"id": "34MqhVdGW39K"
},
"execution_count": 60,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.5259731898104099>"
]
},
"metadata": {},
"execution_count": 60
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "a1b3ba30-75af-4a8a-b5bb-3f68d9ac2eaa",
"id": "6tYQ35bGW39L"
},
"execution_count": 61,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.5259731898104131"
]
},
"metadata": {},
"execution_count": 61
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "740af2fb-70f3-4d35-df71-b8f9203849f8",
"id": "CI2LuvXEW39L"
},
"execution_count": 62,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"6.121313469558163e-15"
]
},
"metadata": {},
"execution_count": 62
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "539cdc2b-12bd-4c08-8528-8a6c72358f62",
"id": "SfeRHj0ZW39M"
},
"execution_count": 63,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"8.635067969306773e-16"
]
},
"metadata": {},
"execution_count": 63
}
]
},
{
"cell_type": "code",
"source": [
"_betaincinv_computation(a, b, 0.)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "_61gtRwByfT7",
"outputId": "40d3f3f0-6ecb-40ae-87aa-210c98679726"
},
"execution_count": 64,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.0>"
]
},
"metadata": {},
"execution_count": 64
}
]
},
{
"cell_type": "code",
"source": [
"_betaincinv_computation(a, b, 1.)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "6LRqC0Coeoh9",
"outputId": "743d3c9d-e1a5-4c84-ece1-02c21f21fa1e"
},
"execution_count": 65,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=1.0>"
]
},
"metadata": {},
"execution_count": 65
}
]
},
{
"cell_type": "markdown",
"source": [
"## Sanity checks for `float32`"
],
"metadata": {
"id": "GC-0lYG0HYAp"
}
},
{
"cell_type": "code",
"source": [
"dtype = np.float32"
],
"metadata": {
"id": "Nzw2Z4EyYtcc"
},
"execution_count": 66,
"outputs": []
},
{
"cell_type": "code",
"source": [
"a = 1.2\n",
"b = 9.8\n",
"y = 0.05"
],
"metadata": {
"id": "QI5pr8VWYtcc"
},
"execution_count": 67,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "ce7154bc-b117-4000-bc35-1a0d3e9462d6",
"id": "1Jc84_MNYtcc"
},
"execution_count": 68,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.10909090909090909"
]
},
"metadata": {},
"execution_count": 68
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "93f01f36-1c6c-4077-d7d5-06555b9d929e",
"id": "XdZoldGlYtcd"
},
"execution_count": 69,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"True"
]
},
"metadata": {},
"execution_count": 69
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "tsHJErcTYtcd"
},
"execution_count": 70,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "c4985d6d-2fa1-4d20-ab44-69770ff5a64a",
"id": "XwDwhIzKYtcd"
},
"execution_count": 71,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=0.0093663065>"
]
},
"metadata": {},
"execution_count": 71
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "65387ddb-2707-40b2-b22e-ef36e7080729",
"id": "CVLeA67cYtcd"
},
"execution_count": 72,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.009366301"
]
},
"metadata": {},
"execution_count": 72
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "2ebb7990-9b72-45c3-88b8-d197f2f59e1d",
"id": "5_A4g3RfYtcd"
},
"execution_count": 73,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"5.9660005e-07"
]
},
"metadata": {},
"execution_count": 73
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "fdbde3c9-e6a9-42f8-9b8a-e18f68902893",
"id": "9LqoevwVYtcd"
},
"execution_count": 74,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"7.4505806e-08"
]
},
"metadata": {},
"execution_count": 74
}
]
},
{
"cell_type": "code",
"source": [
"a = 1.2\n",
"b = 9.8\n",
"y = 0.9"
],
"metadata": {
"id": "NK_PaL5SYtcd"
},
"execution_count": 75,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "5f512f6c-7ad8-44e0-f6b6-c9c3a7050f19",
"id": "maOr4WwOYtcd"
},
"execution_count": 76,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.10909090909090909"
]
},
"metadata": {},
"execution_count": 76
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "df471c81-f776-4781-9c7c-2a45312ad397",
"id": "UKyvYHNZYtce"
},
"execution_count": 77,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"False"
]
},
"metadata": {},
"execution_count": 77
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "zozfh1qTYtce"
},
"execution_count": 78,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "5a565124-5c2e-438c-9c74-582d62236603",
"id": "4r4vFk0_Ytce"
},
"execution_count": 79,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=0.23426865>"
]
},
"metadata": {},
"execution_count": 79
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "b54ba212-bb37-4d4a-b8a9-8b114736f318",
"id": "ZdqEdT4zYtce"
},
"execution_count": 80,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.23426864"
]
},
"metadata": {},
"execution_count": 80
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "010ec521-c1af-431c-f4f5-6c447040bf4b",
"id": "qOVrEZc4Ytce"
},
"execution_count": 81,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"6.360715e-08"
]
},
"metadata": {},
"execution_count": 81
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "92034978-f986-45db-8840-dfd85cd8140c",
"id": "78_42wynYtce"
},
"execution_count": 82,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"6.6227386e-08"
]
},
"metadata": {},
"execution_count": 82
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.8\n",
"b = 1.7\n",
"y = 0.1"
],
"metadata": {
"id": "v0TD8biVYtce"
},
"execution_count": 83,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "60ba4652-2cee-4b8c-c44e-5925f19b1038",
"id": "ffumNbmFYtce"
},
"execution_count": 84,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.32"
]
},
"metadata": {},
"execution_count": 84
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "84d3667f-3885-49c5-a061-bc6fd9342db9",
"id": "yJNdq2bLYtce"
},
"execution_count": 85,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"True"
]
},
"metadata": {},
"execution_count": 85
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "0aKEOrPHYtcf"
},
"execution_count": 86,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "ea6b49c4-e81b-49aa-8afd-d9654069538f",
"id": "PYuvEtBwYtcf"
},
"execution_count": 87,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=0.032386925>"
]
},
"metadata": {},
"execution_count": 87
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "6d1744a7-9f8d-4375-c964-95f7700092ce",
"id": "DX193sZiYtcf"
},
"execution_count": 88,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.032386925"
]
},
"metadata": {},
"execution_count": 88
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "5efc6cb2-d224-463a-940e-80e7cae95f0a",
"id": "U-oruwvEYtcf"
},
"execution_count": 89,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.0"
]
},
"metadata": {},
"execution_count": 89
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "8213c5a4-28ad-4a97-cea0-f4fe972c0a24",
"id": "Lg31dLhhYtcf"
},
"execution_count": 90,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.4901161e-07"
]
},
"metadata": {},
"execution_count": 90
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.8\n",
"b = 1.7\n",
"y = 0.9"
],
"metadata": {
"id": "GwHlFs1vYtcf"
},
"execution_count": 91,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "1ab1078b-85ea-4e66-d6f2-b17459a15667",
"id": "tK6ycxpVYtcf"
},
"execution_count": 92,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.32"
]
},
"metadata": {},
"execution_count": 92
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "bf1f8743-258e-4986-b387-3002c562e8eb",
"id": "2gQEFO7VYtcf"
},
"execution_count": 93,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"False"
]
},
"metadata": {},
"execution_count": 93
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "CRgWLhxfYtcf"
},
"execution_count": 94,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "3b9b7c36-ccb5-4727-9c13-e430490375b9",
"id": "8ESEzUH1Ytcf"
},
"execution_count": 95,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=0.6994022>"
]
},
"metadata": {},
"execution_count": 95
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "fa86cea6-bb98-4e62-e017-5d772a0dee28",
"id": "PU5RIDmLYtcg"
},
"execution_count": 96,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.6994021"
]
},
"metadata": {},
"execution_count": 96
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "f1aec1c5-864d-4530-8656-295dbd54ab93",
"id": "N_bgWT0kYtcg"
},
"execution_count": 97,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.7044457e-07"
]
},
"metadata": {},
"execution_count": 97
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "5a460c82-76dd-4823-816b-504ac5917f9d",
"id": "mfdoCSnNYtcg"
},
"execution_count": 98,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.0"
]
},
"metadata": {},
"execution_count": 98
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.1\n",
"b = 0.7\n",
"y = 0.1"
],
"metadata": {
"id": "tW8uIcrdYtcg"
},
"execution_count": 99,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "65486c9d-c5ad-4454-e01f-f031e437d00a",
"id": "AE6GIJbeYtcg"
},
"execution_count": 100,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.12500000000000003"
]
},
"metadata": {},
"execution_count": 100
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "8d5d6d88-eff9-4795-b6a5-28991cc8dae1",
"id": "xIXss844Ytcg"
},
"execution_count": 101,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"True"
]
},
"metadata": {},
"execution_count": 101
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "kMWTtL4lYtcg"
},
"execution_count": 102,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "453172b9-0f59-443e-f150-02cfeca3b759",
"id": "jbVgfZVuYtcg"
},
"execution_count": 103,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=1.8028251e-10>"
]
},
"metadata": {},
"execution_count": 103
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "700f3bf5-fe12-4e96-e9fa-73c8f816f5b5",
"id": "TOahbeaVYtcg"
},
"execution_count": 104,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.8028194e-10"
]
},
"metadata": {},
"execution_count": 104
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "436a3dc2-d4eb-47f8-95ef-9330daa00bfd",
"id": "7l3dqJzSYtcg"
},
"execution_count": 105,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"3.156108e-06"
]
},
"metadata": {},
"execution_count": 105
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "66ecfffa-602e-4229-baee-9a5e607630b0",
"id": "FOqCOm6yYtch"
},
"execution_count": 106,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"7.4505806e-08"
]
},
"metadata": {},
"execution_count": 106
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.1\n",
"b = 0.7\n",
"y = 0.2"
],
"metadata": {
"id": "RAyFzGUBYtch"
},
"execution_count": 107,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "aaae2e83-1cd6-4844-c508-e8a6069db7e8",
"id": "iYts-njEYtch"
},
"execution_count": 108,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.12500000000000003"
]
},
"metadata": {},
"execution_count": 108
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "2c7dc8c7-e765-4342-a4bc-db40d1f7febe",
"id": "TnTST3tdYtch"
},
"execution_count": 109,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"False"
]
},
"metadata": {},
"execution_count": 109
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "LP06GvW2Ytch"
},
"execution_count": 110,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "692077df-01db-41b0-9e09-765d4ca919bd",
"id": "AojiasC4Ytci"
},
"execution_count": 111,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=1.846088e-07>"
]
},
"metadata": {},
"execution_count": 111
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "d058bcee-37cf-4ccd-af06-400e83946a6e",
"id": "Vm424kynYtci"
},
"execution_count": 112,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.8460868e-07"
]
},
"metadata": {},
"execution_count": 112
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "2cfcf892-9681-4735-952c-e544148635f4",
"id": "Aqq2CIIwYtci"
},
"execution_count": 113,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"6.158261e-07"
]
},
"metadata": {},
"execution_count": 113
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "06d73b9e-8269-42d6-fdcb-42548d0f6b49",
"id": "dp58WWKxYtci"
},
"execution_count": 114,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"7.4505806e-08"
]
},
"metadata": {},
"execution_count": 114
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.1\n",
"b = 0.7\n",
"y = 0.9"
],
"metadata": {
"id": "0HIEFHsUYtci"
},
"execution_count": 115,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "e8f30d99-e578-42e9-f547-3e0d8cb6be6d",
"id": "614IxAryYtci"
},
"execution_count": 116,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.12500000000000003"
]
},
"metadata": {},
"execution_count": 116
}
]
},
{
"cell_type": "code",
"source": [
"y <= mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "11688901-35a7-49d4-9ee1-c4103836f9e5",
"id": "XsYrSd_qYtcj"
},
"execution_count": 117,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"False"
]
},
"metadata": {},
"execution_count": 117
}
]
},
{
"cell_type": "code",
"source": [
"a, b, y = [dtype(z) for z in [a, b, y]]"
],
"metadata": {
"id": "GNIGIxOAYtcj"
},
"execution_count": 118,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = _betaincinv_computation(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "b16ca6eb-f190-40c4-f74c-391b949129cc",
"id": "2-m6ewESYtcj"
},
"execution_count": 119,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=0.5259731>"
]
},
"metadata": {},
"execution_count": 119
}
]
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(a, b, y)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "564f6e3e-7bbc-4ffa-9212-84d88e4dec33",
"id": "PS34OHCYYtcj"
},
"execution_count": 120,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.5259731"
]
},
"metadata": {},
"execution_count": 120
}
]
},
{
"cell_type": "code",
"source": [
"rerr_x_tf = np.abs(x_tf - x_sp) / x_sp\n",
"rerr_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "d003f25b-a1f4-46fc-a1f9-8d4b9a765426",
"id": "dT6XkN0xYtcj"
},
"execution_count": 121,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.0"
]
},
"metadata": {},
"execution_count": 121
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_tf = np.abs(tfp_math.betainc(a, b, x_tf) - y) / y\n",
"rerr_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "96d203e4-5c9f-44f6-f5c3-d2134fafa4a1",
"id": "Kha7NgpJYtcj"
},
"execution_count": 122,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"6.6227386e-08"
]
},
"metadata": {},
"execution_count": 122
}
]
},
{
"cell_type": "code",
"source": [
"_betaincinv_computation(a, b, 0.)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "a8041602-14e5-4977-ee8e-20689e93f776",
"id": "fSYWwLuRYtcj"
},
"execution_count": 123,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=0.0>"
]
},
"metadata": {},
"execution_count": 123
}
]
},
{
"cell_type": "code",
"source": [
"_betaincinv_computation(a, b, 1.)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "53516f2f-9373-4e84-b810-8e0498eb29ee",
"id": "gbXR4YuXYtcj"
},
"execution_count": 124,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=1.0>"
]
},
"metadata": {},
"execution_count": 124
}
]
},
{
"cell_type": "markdown",
"source": [
"## Test settings and auxiliary functions"
],
"metadata": {
"id": "d15s3qJNgkAc"
}
},
{
"cell_type": "code",
"source": [
"SEEDS = [1, 17, 42, 51, 184, 301, 346, 448, 733, 985]"
],
"metadata": {
"id": "bFBYwIZhtV-G"
},
"execution_count": 125,
"outputs": []
},
{
"cell_type": "code",
"source": [
"len(SEEDS)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "2HR-JJbhXtAp",
"outputId": "8bce3929-106f-4fda-e7f8-a9509d755ba2"
},
"execution_count": 126,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"10"
]
},
"metadata": {},
"execution_count": 126
}
]
},
{
"cell_type": "code",
"source": [
"SIZE_PER_SEED = 50_000\n",
"# SIZE = SIZE_PER_SEED * len(SEEDS)"
],
"metadata": {
"id": "d0UUVrpW7xIC"
},
"execution_count": 127,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"`sample_specs` specifies how to simulate values for parameters `a` and `b`. We call each specification a domain.\n",
"\n",
"Values for the parameter `y` are simulated from a Uniform distribution on the interval (0, 1)."
],
"metadata": {
"id": "O9ynlPQzA2sR"
}
},
{
"cell_type": "code",
"source": [
"sample_specs = {\n",
" np.float64: [\n",
" ('uniform', 1.),\n",
" ('uniform', 5.),\n",
" ('uniform', 10.),\n",
" ('uniform', 30.),\n",
" ('uniform', 100.),\n",
" ('uniform', 1000.),\n",
" ('uniform', 10000.),\n",
" ('uniform', 100000.)],\n",
" np.float32: [\n",
" ('uniform', 1.),\n",
" ('uniform', 5.),\n",
" ('uniform', 10.),\n",
" ('uniform', 30.),\n",
" ('uniform', 100.),\n",
" ('uniform', 1000.),\n",
" ('uniform', 10000.),\n",
" ('uniform', 100000.)]}"
],
"metadata": {
"id": "GEB-pRUz7yHH"
},
"execution_count": 128,
"outputs": []
},
{
"cell_type": "code",
"source": [
"def make_samples(sample_specs, size_per_seed, seeds, dtype):\n",
" tiny = tf.constant(np.finfo(dtype).tiny, dtype)\n",
" one = tf.constant(1., dtype)\n",
" eps = tf.constant(np.finfo(dtype).eps, dtype)\n",
"\n",
" # This function supports uniform and half_normal. \n",
" def sample(rng, name, scale):\n",
" if name == 'uniform':\n",
" return np.clip(\n",
" scale * rng.uniform(size=size_per_seed), a_min=tiny, a_max=None)\n",
" return np.clip(\n",
" scale * np.abs(rng.randn(size_per_seed)), a_min=tiny, a_max=None)\n",
"\n",
" samples = []\n",
" for spec in sample_specs:\n",
" list_a = []\n",
" list_b = []\n",
" list_y = []\n",
"\n",
" for seed in seeds:\n",
" rng = np.random.RandomState(seed)\n",
" a = sample(rng, *spec).astype(dtype)\n",
" list_a.append(a)\n",
" b = sample(rng, *spec).astype(dtype)\n",
" list_b.append(b)\n",
" y_min = tfp_math.betainc(a, b, tiny)\n",
" y_max = tfp_math.betainc(a, b, one - eps)\n",
" # Simulates values for y such that tiny < betaincinv(a, b, y) < 1 - eps\n",
" y = rng.uniform(size=size_per_seed, low=y_min, high=y_max).astype(dtype)\n",
" list_y.append(y)\n",
" \n",
" a = np.concatenate(list_a)\n",
" b = np.concatenate(list_b)\n",
" y = np.concatenate(list_y)\n",
" samples.append([a, b, y])\n",
"\n",
" return samples"
],
"metadata": {
"id": "0EfoCKxF70Sw"
},
"execution_count": 129,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"`format_number` and `get_metrics_dataframe` functions are used to print test results."
],
"metadata": {
"id": "vsIEqUoPQtaE"
}
},
{
"cell_type": "code",
"source": [
"def format_number(num):\n",
" if num == 0:\n",
" return '0'\n",
" elif -2 < np.log10(num) < 6: \n",
" if num == int(num):\n",
" return f'{num:.0f}'\n",
" else:\n",
" return f'{num:.2f}'\n",
" return f'{num:.1e}'"
],
"metadata": {
"id": "Biw3tiGVucgS"
},
"execution_count": 130,
"outputs": []
},
{
"cell_type": "code",
"source": [
"def get_metrics_dataframe(\n",
" sample_specs,\n",
" sp_results,\n",
" tf_results,\n",
"):\n",
" records = []\n",
"\n",
" idx_cmp = 0\n",
" for idx, (distr, domain) in enumerate(sample_specs[DTYPE]):\n",
" size = sp_results[idx].shape[0]\n",
" record = {'Distribution': distr}\n",
" record['Domain'] = format_number(domain)\n",
" record['#Trials'] = format_number(size)\n",
"\n",
" perc_nan = np.sum(np.isnan(tf_results[idx])) / size * 100\n",
" perc_inf = np.sum(np.isinf(tf_results[idx])) / size * 100\n",
" record['%NaN'] = format_number(perc_nan)\n",
" record['%Inf'] = format_number(perc_inf)\n",
"\n",
" if (perc_nan + perc_inf) == 100:\n",
" record['Min Result'] = 'NaN'\n",
" record['Max Result'] = 'NaN'\n",
" else:\n",
" res = tf_results[idx]\n",
" min_res = np.min(res, initial=np.inf, where=np.isfinite(res))\n",
" max_res = np.max(res, initial=-np.inf, where=np.isfinite(res))\n",
" record['Min Result'] = format_number(min_res)\n",
" record['Max Result'] = format_number(max_res)\n",
"\n",
" if (perc_nan + perc_inf) == 100:\n",
" record['Max Abs Error'] = 'NaN'\n",
" record['Mean Abs Error'] = 'NaN'\n",
" record['Max Rel Error'] = 'NaN'\n",
" record['Mean Rel Error'] = 'NaN'\n",
" else:\n",
" aerr = np.abs(sp_results[idx] - tf_results[idx])\n",
" aerr_valid = np.isfinite(aerr)\n",
" max_aerr = np.max(aerr, initial=0., where=aerr_valid)\n",
" mean_aerr = np.mean(aerr, where=aerr_valid)\n",
" record['Max Abs Error'] = format_number(max_aerr)\n",
" record['Mean Abs Error'] = format_number(mean_aerr)\n",
"\n",
" div_by_zero = (np.abs(sp_results[idx]) < np.sqrt(TINY))\n",
" rerr = aerr / np.where(div_by_zero, 1., np.abs(sp_results[idx]))\n",
" rerr_valid = np.isfinite(rerr)\n",
" max_rerr = np.max(rerr, initial=0., where=rerr_valid)\n",
" mean_rerr = np.mean(rerr, where=rerr_valid)\n",
" record['Max Rel Error'] = format_number(max_rerr)\n",
" record['Mean Rel Error'] = format_number(mean_rerr)\n",
"\n",
" records.append(record)\n",
"\n",
" return pd.DataFrame.from_records(records)"
],
"metadata": {
"id": "XRvmj1-9Kukh"
},
"execution_count": 131,
"outputs": []
},
{
"cell_type": "code",
"source": [
"sp_betaincinv = sp_special.betaincinv\n",
"tf_betaincinv = tf.function(\n",
" _betaincinv_computation, autograph=False, jit_compile=False)"
],
"metadata": {
"id": "fUEgyjLQbBmV"
},
"execution_count": 132,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## Test results for `float64`"
],
"metadata": {
"id": "KebKjtrOgwgx"
}
},
{
"cell_type": "code",
"source": [
"DTYPE = np.float64"
],
"metadata": {
"id": "KSS4xfUr7odg"
},
"execution_count": 133,
"outputs": []
},
{
"cell_type": "code",
"source": [
"TINY = np.finfo(DTYPE).tiny\n",
"TINY"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "ShOeQ13V7pqG",
"outputId": "017ec1d9-ac0d-4f78-8a95-e843b76cbee0"
},
"execution_count": 134,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"2.2250738585072014e-308"
]
},
"metadata": {},
"execution_count": 134
}
]
},
{
"cell_type": "code",
"source": [
"EPS = np.finfo(DTYPE).eps\n",
"EPS"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "3Vws5VIjhinP",
"outputId": "85c51f90-4034-4c6b-832b-f60e5fd9e504"
},
"execution_count": 135,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"2.220446049250313e-16"
]
},
"metadata": {},
"execution_count": 135
}
]
},
{
"cell_type": "code",
"source": [
"samples = make_samples(\n",
" sample_specs[DTYPE], size_per_seed=SIZE_PER_SEED, seeds=SEEDS, dtype=DTYPE)"
],
"metadata": {
"id": "umRhJjde8Cp7"
},
"execution_count": 136,
"outputs": []
},
{
"cell_type": "code",
"source": [
"sp_res = [\n",
" sp_betaincinv(a, b, y) for a, b, y in samples]"
],
"metadata": {
"id": "JUw64yHs-dj8"
},
"execution_count": 137,
"outputs": []
},
{
"cell_type": "code",
"source": [
"tf_res = []\n",
"\n",
"for i, (a, b, y) in enumerate(samples):\n",
" result = tf_betaincinv(a, b, y)\n",
" tf_res.append(result)\n",
" print(f'sample {i} is done.')"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Bpohh_aLJZpG",
"outputId": "65aa7ee0-09d6-45b8-9d8e-816f894a9c3e"
},
"execution_count": 138,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"sample 0 is done.\n",
"sample 1 is done.\n",
"sample 2 is done.\n",
"sample 3 is done.\n",
"sample 4 is done.\n",
"sample 5 is done.\n",
"sample 6 is done.\n",
"sample 7 is done.\n"
]
}
]
},
{
"cell_type": "code",
"source": [
"print(f'Method: My implementation of betaincinv')\n",
"print(f'Benchmark: Cephes implementation of betaincinv for float64 [SciPy]')\n",
"print(f'Dtype: {np.dtype(DTYPE).name}')\n",
"get_metrics_dataframe(sample_specs, sp_res, tf_res)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 354
},
"id": "i1j6PxtaUEnB",
"outputId": "9f4a0043-3662-4776-d472-aff537b65743"
},
"execution_count": 139,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Method: My implementation of betaincinv\n",
"Benchmark: Cephes implementation of betaincinv for float64 [SciPy]\n",
"Dtype: float64\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": [
" Distribution Domain #Trials %NaN %Inf Min Result Max Result Max Abs Error \\\n",
"0 uniform 1 500000 0 0 3.2e-308 1.00 7.6e-12 \n",
"1 uniform 5 500000 0 0 2.2e-308 1.00 6.7e-13 \n",
"2 uniform 10 500000 0 0 1.9e-307 1.00 3.6e-13 \n",
"3 uniform 30 500000 0 0 1.9e-307 1.00 2.6e-13 \n",
"4 uniform 100 500000 0 0 2.1e-307 1.00 2.0e-13 \n",
"5 uniform 1000 500000 0 0 5.1e-307 1.00 9.5e-13 \n",
"6 uniform 10000 500000 0 0 6.1e-235 1.00 4.3e-13 \n",
"7 uniform 100000 500000 0 0 1.5e-29 1.00 1.9e-12 \n",
"\n",
" Mean Abs Error Max Rel Error Mean Rel Error \n",
"0 2.0e-16 2.8e-10 5.4e-15 \n",
"1 9.2e-17 4.4e-11 1.1e-15 \n",
"2 1.7e-16 1.3e-10 1.3e-15 \n",
"3 4.2e-16 6.8e-11 2.2e-15 \n",
"4 1.3e-15 7.9e-11 5.4e-15 \n",
"5 6.4e-15 2.5e-10 2.1e-14 \n",
"6 3.6e-14 4.1e-10 1.1e-13 \n",
"7 1.3e-13 2.4e-10 3.9e-13 "
],
"text/html": [
"\n",
" <div id=\"df-a54e86a7-f18d-49b5-ae4e-d58db5bc19cf\">\n",
" <div class=\"colab-df-container\">\n",
" <div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Distribution</th>\n",
" <th>Domain</th>\n",
" <th>#Trials</th>\n",
" <th>%NaN</th>\n",
" <th>%Inf</th>\n",
" <th>Min Result</th>\n",
" <th>Max Result</th>\n",
" <th>Max Abs Error</th>\n",
" <th>Mean Abs Error</th>\n",
" <th>Max Rel Error</th>\n",
" <th>Mean Rel Error</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>uniform</td>\n",
" <td>1</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>3.2e-308</td>\n",
" <td>1.00</td>\n",
" <td>7.6e-12</td>\n",
" <td>2.0e-16</td>\n",
" <td>2.8e-10</td>\n",
" <td>5.4e-15</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>uniform</td>\n",
" <td>5</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>2.2e-308</td>\n",
" <td>1.00</td>\n",
" <td>6.7e-13</td>\n",
" <td>9.2e-17</td>\n",
" <td>4.4e-11</td>\n",
" <td>1.1e-15</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>uniform</td>\n",
" <td>10</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1.9e-307</td>\n",
" <td>1.00</td>\n",
" <td>3.6e-13</td>\n",
" <td>1.7e-16</td>\n",
" <td>1.3e-10</td>\n",
" <td>1.3e-15</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>uniform</td>\n",
" <td>30</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1.9e-307</td>\n",
" <td>1.00</td>\n",
" <td>2.6e-13</td>\n",
" <td>4.2e-16</td>\n",
" <td>6.8e-11</td>\n",
" <td>2.2e-15</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>uniform</td>\n",
" <td>100</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>2.1e-307</td>\n",
" <td>1.00</td>\n",
" <td>2.0e-13</td>\n",
" <td>1.3e-15</td>\n",
" <td>7.9e-11</td>\n",
" <td>5.4e-15</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>uniform</td>\n",
" <td>1000</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>5.1e-307</td>\n",
" <td>1.00</td>\n",
" <td>9.5e-13</td>\n",
" <td>6.4e-15</td>\n",
" <td>2.5e-10</td>\n",
" <td>2.1e-14</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>uniform</td>\n",
" <td>10000</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>6.1e-235</td>\n",
" <td>1.00</td>\n",
" <td>4.3e-13</td>\n",
" <td>3.6e-14</td>\n",
" <td>4.1e-10</td>\n",
" <td>1.1e-13</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>uniform</td>\n",
" <td>100000</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1.5e-29</td>\n",
" <td>1.00</td>\n",
" <td>1.9e-12</td>\n",
" <td>1.3e-13</td>\n",
" <td>2.4e-10</td>\n",
" <td>3.9e-13</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>\n",
" <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-a54e86a7-f18d-49b5-ae4e-d58db5bc19cf')\"\n",
" title=\"Convert this dataframe to an interactive table.\"\n",
" style=\"display:none;\">\n",
" \n",
" <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
" width=\"24px\">\n",
" <path d=\"M0 0h24v24H0V0z\" fill=\"none\"/>\n",
" <path d=\"M18.56 5.44l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94zm-11 1L8.5 8.5l.94-2.06 2.06-.94-2.06-.94L8.5 2.5l-.94 2.06-2.06.94zm10 10l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94z\"/><path d=\"M17.41 7.96l-1.37-1.37c-.4-.4-.92-.59-1.43-.59-.52 0-1.04.2-1.43.59L10.3 9.45l-7.72 7.72c-.78.78-.78 2.05 0 2.83L4 21.41c.39.39.9.59 1.41.59.51 0 1.02-.2 1.41-.59l7.78-7.78 2.81-2.81c.8-.78.8-2.07 0-2.86zM5.41 20L4 18.59l7.72-7.72 1.47 1.35L5.41 20z\"/>\n",
" </svg>\n",
" </button>\n",
" \n",
" <style>\n",
" .colab-df-container {\n",
" display:flex;\n",
" flex-wrap:wrap;\n",
" gap: 12px;\n",
" }\n",
"\n",
" .colab-df-convert {\n",
" background-color: #E8F0FE;\n",
" border: none;\n",
" border-radius: 50%;\n",
" cursor: pointer;\n",
" display: none;\n",
" fill: #1967D2;\n",
" height: 32px;\n",
" padding: 0 0 0 0;\n",
" width: 32px;\n",
" }\n",
"\n",
" .colab-df-convert:hover {\n",
" background-color: #E2EBFA;\n",
" box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
" fill: #174EA6;\n",
" }\n",
"\n",
" [theme=dark] .colab-df-convert {\n",
" background-color: #3B4455;\n",
" fill: #D2E3FC;\n",
" }\n",
"\n",
" [theme=dark] .colab-df-convert:hover {\n",
" background-color: #434B5C;\n",
" box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
" filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
" fill: #FFFFFF;\n",
" }\n",
" </style>\n",
"\n",
" <script>\n",
" const buttonEl =\n",
" document.querySelector('#df-a54e86a7-f18d-49b5-ae4e-d58db5bc19cf button.colab-df-convert');\n",
" buttonEl.style.display =\n",
" google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
"\n",
" async function convertToInteractive(key) {\n",
" const element = document.querySelector('#df-a54e86a7-f18d-49b5-ae4e-d58db5bc19cf');\n",
" const dataTable =\n",
" await google.colab.kernel.invokeFunction('convertToInteractive',\n",
" [key], {});\n",
" if (!dataTable) return;\n",
"\n",
" const docLinkHtml = 'Like what you see? Visit the ' +\n",
" '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
" + ' to learn more about interactive tables.';\n",
" element.innerHTML = '';\n",
" dataTable['output_type'] = 'display_data';\n",
" await google.colab.output.renderOutput(dataTable, element);\n",
" const docLink = document.createElement('div');\n",
" docLink.innerHTML = docLinkHtml;\n",
" element.appendChild(docLink);\n",
" }\n",
" </script>\n",
" </div>\n",
" </div>\n",
" "
]
},
"metadata": {},
"execution_count": 139
}
]
},
{
"cell_type": "code",
"source": [
"cases = {}\n",
"\n",
"for sample_id in range(8):\n",
" a, b, y = samples[sample_id]\n",
" x_tf = tf_res[sample_id]\n",
" x_sp = sp_res[sample_id]\n",
" data = np.column_stack([a, b, y, x_tf, x_sp])\n",
" df = pd.DataFrame(data=data, columns=['a', 'b', 'y', 'x_tf', 'x_sp'])\n",
" df['aerr'] = np.abs(df.x_tf - df.x_sp)\n",
" df['rerr'] = df.aerr / df.x_sp\n",
" df['aerr_y_tf'] = np.abs(tfp_math.betainc(df.a, df.b, df.x_tf) - df.y)\n",
" df['rerr_y_tf'] = df.aerr_y_tf / df.y\n",
" df['aerr_y_sp'] = np.abs(sp_special.betainc(df.a, df.b, df.x_sp) - df.y)\n",
" df['rerr_y_sp'] = df.aerr_y_sp / df.y\n",
" df['rtol'] = np.maximum(df.aerr - 1e-12, 0) / df.x_sp\n",
" df.sort_values(by='aerr', ascending=False, inplace=True)\n",
" cases[sample_id] = df"
],
"metadata": {
"id": "8-cqYX0YZpYg"
},
"execution_count": 140,
"outputs": []
},
{
"cell_type": "code",
"source": [
"for sample_id in range(8):\n",
" max_rtol = cases[sample_id].rtol.max()\n",
" print(f'sample {sample_id} | max_rtol: {max_rtol}')"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Z11dHJXOZtd5",
"outputId": "0bcb01ed-bc64-4d02-8154-4150e609c337"
},
"execution_count": 141,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"sample 0 | max_rtol: 7.665715989998892e-12\n",
"sample 1 | max_rtol: 0.0\n",
"sample 2 | max_rtol: 0.0\n",
"sample 3 | max_rtol: 0.0\n",
"sample 4 | max_rtol: 0.0\n",
"sample 5 | max_rtol: 0.0\n",
"sample 6 | max_rtol: 0.0\n",
"sample 7 | max_rtol: 1.69155856600542e-12\n"
]
}
]
},
{
"cell_type": "markdown",
"source": [
"### Some comments"
],
"metadata": {
"id": "Jyts6uECw7Cp"
}
},
{
"cell_type": "markdown",
"source": [
"When parameters `a` and `b` take very small values, the functions `sp_special.betainc` and `tfp_math.betainc` go from 0 to 1 extremely fast. Thus, a small approximation error between these functions can result in a large approximation error between the `sp_special.betaincinv` and `tfp_math.betaincinv` functions.\n",
"\n",
"Furthermore, there may be a nondegenerate interval `I` in `(0, 1)` such that `x1` and `x2` belong to `I`, `x1 < x2`, but `sp_special.betainc(a, b, x1)` is equal to `sp_special.betainc(a , b, x2)`. The same goes for `tfp_math.betainc`.\n",
"\n",
"See Example 1."
],
"metadata": {
"id": "p6VPLcIAL5pf"
}
},
{
"cell_type": "markdown",
"source": [
"**Example 1**"
],
"metadata": {
"id": "Aj1NTKrewcxC"
}
},
{
"cell_type": "code",
"source": [
"aa = 1.3323706817470157e-18\n",
"bb = 8.748581246342964e-13\n",
"yy = 0.9999984770460288\n",
"\n",
"aa, bb, yy = [DTYPE(z) for z in (aa, bb, yy)]"
],
"metadata": {
"id": "Y9AQ1pNQKic0"
},
"execution_count": 142,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_left = DTYPE(0.01)\n",
"x_middle = DTYPE(0.5)\n",
"x_right = DTYPE(0.99)"
],
"metadata": {
"id": "3d9frtxge6ru"
},
"execution_count": 143,
"outputs": []
},
{
"cell_type": "code",
"source": [
"sp_special.betainc(aa, bb, x_middle) - sp_special.betainc(aa, bb, x_left)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Htakm1RjfF53",
"outputId": "48a5f0ee-2f68-4069-fa05-70e282597dfc"
},
"execution_count": 144,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.0"
]
},
"metadata": {},
"execution_count": 144
}
]
},
{
"cell_type": "code",
"source": [
"sp_special.betainc(aa, bb, x_right) - sp_special.betainc(aa, bb, x_left)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "lIZA5fJKoeKx",
"outputId": "6a4cf1d3-5bf6-4cb0-f93c-fabae005c59c"
},
"execution_count": 145,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"2.220446049250313e-16"
]
},
"metadata": {},
"execution_count": 145
}
]
},
{
"cell_type": "code",
"source": [
"tfp_math.betainc(aa, bb, x_middle) - tfp_math.betainc(aa, bb, x_left)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "PttEG8VQpV-7",
"outputId": "589bf17a-49ec-4fe7-c17d-c2722f3a700d"
},
"execution_count": 146,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.0>"
]
},
"metadata": {},
"execution_count": 146
}
]
},
{
"cell_type": "code",
"source": [
"tfp_math.betainc(aa, bb, x_right) - tfp_math.betainc(aa, bb, x_left)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "BZ8E42C3olL_",
"outputId": "2deda2ea-ab80-418c-eb57-ddad8121c3a3"
},
"execution_count": 147,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=2.886579864025407e-15>"
]
},
"metadata": {},
"execution_count": 147
}
]
},
{
"cell_type": "code",
"source": [
"tfp_math.betainc(aa, bb, x_right) - tfp_math.betainc(aa, bb, x_left)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "_4CL_yMPfM0v",
"outputId": "203a7dca-ac1f-4ab5-d420-eed6989fe97b"
},
"execution_count": 148,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=2.886579864025407e-15>"
]
},
"metadata": {},
"execution_count": 148
}
]
},
{
"cell_type": "markdown",
"source": [
"When parameters `a` and `b` take very large values, the function `tfp_math.betainc` presents numerical instability. In some cases, this function doesn't preserve the monotonicity of the theoretical function. In fact, there are values `x0` and `x1` in the interval `(0, 1)` such that `x0 < x1` and `tfp_math.betainc(a, b, x0) > tfp_math.betainc(a, b, x1)`. See Example 2. In other cases, the function `tfp_math.betainc` may return a value out of the interval [0, 1]. See Example 3.\n",
"\n",
"In either case, the approximation error of `tfp_math.betaincinv` can be large.\n",
"\n",
"The function `sp_special.betainc` also presents numerical instability when parameters `a` and `b` take very large values. See Example 4."
],
"metadata": {
"id": "m44eqCg5R1M7"
}
},
{
"cell_type": "markdown",
"source": [
"**Example 2**"
],
"metadata": {
"id": "-WoH_4xRnHVD"
}
},
{
"cell_type": "code",
"source": [
"aa = 134038921201746.05\n",
"bb = 2025417281289407.0\n",
"yy = 0.9999727386535785\n",
"\n",
"aa, bb, yy = [DTYPE(z) for z in (aa, bb, yy)]"
],
"metadata": {
"id": "F7WbTKht9KpQ"
},
"execution_count": 149,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_low = 0.06207070385156131\n",
"x = 0.062070703856789776\n",
"x_high = 0.06207070386201824\n",
"\n",
"x_low, x, x_high = [DTYPE(z) for z in (x_low, x, x_high)]"
],
"metadata": {
"id": "uGbzj_fH9wQa"
},
"execution_count": 150,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_low < x < x_high"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "CPGAGEKL9wF7",
"outputId": "9bbc4c1c-9804-4d20-99d1-2941759d93ea"
},
"execution_count": 151,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"True"
]
},
"metadata": {},
"execution_count": 151
}
]
},
{
"cell_type": "code",
"source": [
"y_low = tfp_math.betainc(aa, bb, x_low)\n",
"y = tfp_math.betainc(aa, bb, x)\n",
"y_high = tfp_math.betainc(aa, bb, x_high)"
],
"metadata": {
"id": "PbGxubp_m-cw"
},
"execution_count": 152,
"outputs": []
},
{
"cell_type": "code",
"source": [
"y_low < y < y_high"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "STk9SUF9nL39",
"outputId": "b7fde213-7335-4bdf-e9cf-9a65107cbb14"
},
"execution_count": 153,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=bool, numpy=False>"
]
},
"metadata": {},
"execution_count": 153
}
]
},
{
"cell_type": "code",
"source": [
"y_low"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "qTAXjg5B9eFJ",
"outputId": "a2916129-27e3-4f80-bfc6-baaf2f70e6af"
},
"execution_count": 154,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.9984144385935695>"
]
},
"metadata": {},
"execution_count": 154
}
]
},
{
"cell_type": "code",
"source": [
"y"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "D5EY54Lg-jNI",
"outputId": "b8686daf-6c11-430c-a007-0284020b8aad"
},
"execution_count": 155,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.6784307401997556>"
]
},
"metadata": {},
"execution_count": 155
}
]
},
{
"cell_type": "code",
"source": [
"y_high"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "iBP_u3DG-n0n",
"outputId": "7ef60766-ba0e-407c-a8f9-db7ef189da91"
},
"execution_count": 156,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.999992660705851>"
]
},
"metadata": {},
"execution_count": 156
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_low = np.abs(yy - y_low) / yy\n",
"rerr_y_low"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "4_uiF8f8_9wt",
"outputId": "5c17762e-eb13-47e8-faf5-01499893467e"
},
"execution_count": 157,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.0015583425425249053"
]
},
"metadata": {},
"execution_count": 157
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y = np.abs(yy - y) / yy\n",
"rerr_y"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "VXebtkxC_xu-",
"outputId": "a8f12a91-07cb-4eb3-ad8c-129e54bef040"
},
"execution_count": 158,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.32155076436060226"
]
},
"metadata": {},
"execution_count": 158
}
]
},
{
"cell_type": "code",
"source": [
"rerr_y_high = np.abs(yy - y_high) / yy\n",
"rerr_y_high"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "2uL53p0X_9fI",
"outputId": "215c3a13-1dd5-4ae4-f5a8-739b3c4d11d1"
},
"execution_count": 159,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.9922595389299715e-05"
]
},
"metadata": {},
"execution_count": 159
}
]
},
{
"cell_type": "code",
"source": [
"initial_x_tf = _betaincinv_initial_approx(aa, bb, yy, DTYPE)\n",
"initial_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "wBUYqy6n_9Yq",
"outputId": "a17d45bd-45e4-4b79-e530-b2faeaef502a"
},
"execution_count": 160,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.06207070386201824>"
]
},
"metadata": {},
"execution_count": 160
}
]
},
{
"cell_type": "code",
"source": [
"initial_x_tf == x_high"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "kKCstSRRnqr8",
"outputId": "ac0bc375-4a2d-40b1-9262-9f34a1a26ce4"
},
"execution_count": 161,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=bool, numpy=True>"
]
},
"metadata": {},
"execution_count": 161
}
]
},
{
"cell_type": "code",
"source": [
"final_x_tf = _betaincinv_computation(aa, bb, yy)\n",
"final_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "IkmIFSU-nvWx",
"outputId": "5445bfd9-32bd-4643-82e2-32d98eef7c7a"
},
"execution_count": 162,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.062070703856789776>"
]
},
"metadata": {},
"execution_count": 162
}
]
},
{
"cell_type": "code",
"source": [
"final_x_tf == x"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Dlyo0reOoqCu",
"outputId": "e5e82397-00c3-4190-877e-6f147ccb2af9"
},
"execution_count": 163,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=bool, numpy=True>"
]
},
"metadata": {},
"execution_count": 163
}
]
},
{
"cell_type": "markdown",
"source": [
"**Example 3**"
],
"metadata": {
"id": "E-EVDnEdqLU4"
}
},
{
"cell_type": "code",
"source": [
"aa = 2287949019054372.5\n",
"bb = 2416753116337517.0\n",
"yy = 0.1549079401918569\n",
"\n",
"aa, bb, yy = [DTYPE(z) for z in (aa, bb, yy)]"
],
"metadata": {
"id": "tR9q6YYo7urZ"
},
"execution_count": 164,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(aa, bb, yy)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "bZ3W_bt275wL",
"outputId": "2ba95887-c7f5-441e-c3c5-0494c3d0ddb8"
},
"execution_count": 165,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.48631113107520124"
]
},
"metadata": {},
"execution_count": 165
}
]
},
{
"cell_type": "code",
"source": [
"initial_x_tf = _betaincinv_initial_approx(aa, bb, yy, DTYPE)\n",
"initial_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "ahtLtKLnLXUP",
"outputId": "b4eaeaab-403e-4ca6-b5e1-53720f78b492"
},
"execution_count": 166,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.48631112414646716>"
]
},
"metadata": {},
"execution_count": 166
}
]
},
{
"cell_type": "code",
"source": [
"rerr_initial_x_tf = np.abs(x_sp - initial_x_tf) / x_sp\n",
"rerr_initial_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "JJvdenlQSo9i",
"outputId": "9d5548ed-3302-4c79-e381-91a48126a871"
},
"execution_count": 167,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.4247533407587163e-08"
]
},
"metadata": {},
"execution_count": 167
}
]
},
{
"cell_type": "code",
"source": [
"initial_y_tf = tfp_math.betainc(aa, bb, initial_x_tf)\n",
"initial_y_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "UAMepEwIrVN5",
"outputId": "ad1c8c1c-6ef9-4554-9351-436ca5a141ad"
},
"execution_count": 168,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=2.055420439727784e+19>"
]
},
"metadata": {},
"execution_count": 168
}
]
},
{
"cell_type": "code",
"source": [
"0. < initial_y_tf < 1."
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "19QIBx2prrwQ",
"outputId": "4d23e303-6ce5-4fc1-dfad-3ec19ce9bee7"
},
"execution_count": 169,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=bool, numpy=False>"
]
},
"metadata": {},
"execution_count": 169
}
]
},
{
"cell_type": "code",
"source": [
"final_x_tf = _betaincinv_computation(aa, bb, yy)\n",
"final_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "yAWV6ViFOgYT",
"outputId": "b7b69556-e107-4b8e-d84d-d3846a8196b5"
},
"execution_count": 170,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float64, numpy=0.48441147131777>"
]
},
"metadata": {},
"execution_count": 170
}
]
},
{
"cell_type": "code",
"source": [
"rerr_final_x_tf = np.abs(x_sp - final_x_tf) / x_sp\n",
"rerr_final_x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "wPqLi-lcsA1O",
"outputId": "8d8b3b0c-79a9-484c-de40-5164b57e422b"
},
"execution_count": 171,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.003906264191878986"
]
},
"metadata": {},
"execution_count": 171
}
]
},
{
"cell_type": "markdown",
"source": [
"**Example 4**"
],
"metadata": {
"id": "ELhr49p6uXq7"
}
},
{
"cell_type": "code",
"source": [
"aa = 1561191471113026.5\n",
"bb = 4251271024489255.0\n",
"yy = 0.6428796345575062\n",
"\n",
"aa, bb, yy = [DTYPE(z) for z in (aa, bb, yy)]"
],
"metadata": {
"id": "rcsrscDruqF5"
},
"execution_count": 172,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_sp = sp_special.betaincinv(aa, bb, yy)\n",
"x_sp"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "h15zkVaf0bI7",
"outputId": "b376c3fa-fcd4-4d56-8e58-d39f920fe9f3"
},
"execution_count": 173,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.26859381624699485"
]
},
"metadata": {},
"execution_count": 173
}
]
},
{
"cell_type": "code",
"source": [
"sp_special.betainc(aa, bb, x_sp)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "ZDKqM67E0mAJ",
"outputId": "c8dabf39-b8ff-49ef-bce3-c58ffdf01308"
},
"execution_count": 174,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"-2.5194361402037033e+18"
]
},
"metadata": {},
"execution_count": 174
}
]
},
{
"cell_type": "markdown",
"source": [
"## Test results for `float32`"
],
"metadata": {
"id": "nl4hOD0Lg7OJ"
}
},
{
"cell_type": "code",
"source": [
"DTYPE = np.float32"
],
"metadata": {
"id": "eXU_VXWlWDqG"
},
"execution_count": 175,
"outputs": []
},
{
"cell_type": "code",
"source": [
"TINY = np.finfo(DTYPE).tiny\n",
"TINY"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "065d2ff9-fcee-4dfd-8d22-b0808abfbed4",
"id": "wRbuQNJ7WDqH"
},
"execution_count": 176,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.1754944e-38"
]
},
"metadata": {},
"execution_count": 176
}
]
},
{
"cell_type": "code",
"source": [
"EPS = np.finfo(DTYPE).eps\n",
"EPS"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "uci_ntVBhBsn",
"outputId": "160d7dff-6e99-438e-8ab3-82335e9ac9db"
},
"execution_count": 177,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.1920929e-07"
]
},
"metadata": {},
"execution_count": 177
}
]
},
{
"cell_type": "code",
"source": [
"samples = make_samples(\n",
" sample_specs[DTYPE], size_per_seed=SIZE_PER_SEED, seeds=SEEDS, dtype=DTYPE)"
],
"metadata": {
"id": "TocNiKaXWDqH"
},
"execution_count": 178,
"outputs": []
},
{
"cell_type": "code",
"source": [
"sp_res = [\n",
" sp_betaincinv(a, b, y) for a, b, y in samples]"
],
"metadata": {
"id": "-luHQfz6WDqH"
},
"execution_count": 179,
"outputs": []
},
{
"cell_type": "code",
"source": [
"tf_res = []\n",
"\n",
"for i, (a, b, y) in enumerate(samples):\n",
" result = tf_betaincinv(a, b, y)\n",
" tf_res.append(result)\n",
" print(f'sample {i} is done.')"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "rcaLGwOJKQWs",
"outputId": "9f0aed83-3d5f-446c-f6b3-0e208c2ce184"
},
"execution_count": 180,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"sample 0 is done.\n",
"sample 1 is done.\n",
"sample 2 is done.\n",
"sample 3 is done.\n",
"sample 4 is done.\n",
"sample 5 is done.\n",
"sample 6 is done.\n",
"sample 7 is done.\n"
]
}
]
},
{
"cell_type": "code",
"source": [
"print(f'Method: My implementation of betaincinv')\n",
"print(f'Benchmark: Cephes implementation of betaincinv for float64 [SciPy]')\n",
"print(f'Dtype: {np.dtype(DTYPE).name}')\n",
"get_metrics_dataframe(sample_specs, sp_res, tf_res)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 354
},
"id": "x2AKnR2EVFem",
"outputId": "adfa38f3-546b-4365-a5fa-fd150acf3763"
},
"execution_count": 181,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Method: My implementation of betaincinv\n",
"Benchmark: Cephes implementation of betaincinv for float64 [SciPy]\n",
"Dtype: float32\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": [
" Distribution Domain #Trials %NaN %Inf Min Result Max Result Max Abs Error \\\n",
"0 uniform 1 500000 0 0 1.2e-38 1.00 1.0e-04 \n",
"1 uniform 5 500000 0 0 1.2e-38 1.00 8.3e-04 \n",
"2 uniform 10 500000 0 0 1.2e-38 1.00 6.8e-04 \n",
"3 uniform 30 500000 0 0 1.2e-38 1.00 2.1e-04 \n",
"4 uniform 100 500000 0 0 1.2e-38 1.00 2.9e-05 \n",
"5 uniform 1000 500000 0 0 1.2e-38 1.00 3.1e-05 \n",
"6 uniform 10000 500000 0 0 1.2e-38 1.00 1.6e-04 \n",
"7 uniform 100000 500000 0 0 2.2e-29 1.00 4.8e-04 \n",
"\n",
" Mean Abs Error Max Rel Error Mean Rel Error \n",
"0 2.5e-08 0.01 4.6e-07 \n",
"1 5.5e-08 0.02 3.1e-07 \n",
"2 1.1e-07 0.04 5.9e-07 \n",
"3 2.7e-07 0.02 1.2e-06 \n",
"4 7.5e-07 0.09 2.9e-06 \n",
"5 3.1e-06 0.09 1.0e-05 \n",
"6 1.6e-05 0.23 4.8e-05 \n",
"7 6.1e-05 0.21 1.8e-04 "
],
"text/html": [
"\n",
" <div id=\"df-9d579f69-9612-4eaa-87d7-fea3558c2d6d\">\n",
" <div class=\"colab-df-container\">\n",
" <div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Distribution</th>\n",
" <th>Domain</th>\n",
" <th>#Trials</th>\n",
" <th>%NaN</th>\n",
" <th>%Inf</th>\n",
" <th>Min Result</th>\n",
" <th>Max Result</th>\n",
" <th>Max Abs Error</th>\n",
" <th>Mean Abs Error</th>\n",
" <th>Max Rel Error</th>\n",
" <th>Mean Rel Error</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>uniform</td>\n",
" <td>1</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1.2e-38</td>\n",
" <td>1.00</td>\n",
" <td>1.0e-04</td>\n",
" <td>2.5e-08</td>\n",
" <td>0.01</td>\n",
" <td>4.6e-07</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>uniform</td>\n",
" <td>5</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1.2e-38</td>\n",
" <td>1.00</td>\n",
" <td>8.3e-04</td>\n",
" <td>5.5e-08</td>\n",
" <td>0.02</td>\n",
" <td>3.1e-07</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>uniform</td>\n",
" <td>10</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1.2e-38</td>\n",
" <td>1.00</td>\n",
" <td>6.8e-04</td>\n",
" <td>1.1e-07</td>\n",
" <td>0.04</td>\n",
" <td>5.9e-07</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>uniform</td>\n",
" <td>30</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1.2e-38</td>\n",
" <td>1.00</td>\n",
" <td>2.1e-04</td>\n",
" <td>2.7e-07</td>\n",
" <td>0.02</td>\n",
" <td>1.2e-06</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>uniform</td>\n",
" <td>100</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1.2e-38</td>\n",
" <td>1.00</td>\n",
" <td>2.9e-05</td>\n",
" <td>7.5e-07</td>\n",
" <td>0.09</td>\n",
" <td>2.9e-06</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>uniform</td>\n",
" <td>1000</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1.2e-38</td>\n",
" <td>1.00</td>\n",
" <td>3.1e-05</td>\n",
" <td>3.1e-06</td>\n",
" <td>0.09</td>\n",
" <td>1.0e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>uniform</td>\n",
" <td>10000</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1.2e-38</td>\n",
" <td>1.00</td>\n",
" <td>1.6e-04</td>\n",
" <td>1.6e-05</td>\n",
" <td>0.23</td>\n",
" <td>4.8e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>uniform</td>\n",
" <td>100000</td>\n",
" <td>500000</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>2.2e-29</td>\n",
" <td>1.00</td>\n",
" <td>4.8e-04</td>\n",
" <td>6.1e-05</td>\n",
" <td>0.21</td>\n",
" <td>1.8e-04</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>\n",
" <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-9d579f69-9612-4eaa-87d7-fea3558c2d6d')\"\n",
" title=\"Convert this dataframe to an interactive table.\"\n",
" style=\"display:none;\">\n",
" \n",
" <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
" width=\"24px\">\n",
" <path d=\"M0 0h24v24H0V0z\" fill=\"none\"/>\n",
" <path d=\"M18.56 5.44l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94zm-11 1L8.5 8.5l.94-2.06 2.06-.94-2.06-.94L8.5 2.5l-.94 2.06-2.06.94zm10 10l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94z\"/><path d=\"M17.41 7.96l-1.37-1.37c-.4-.4-.92-.59-1.43-.59-.52 0-1.04.2-1.43.59L10.3 9.45l-7.72 7.72c-.78.78-.78 2.05 0 2.83L4 21.41c.39.39.9.59 1.41.59.51 0 1.02-.2 1.41-.59l7.78-7.78 2.81-2.81c.8-.78.8-2.07 0-2.86zM5.41 20L4 18.59l7.72-7.72 1.47 1.35L5.41 20z\"/>\n",
" </svg>\n",
" </button>\n",
" \n",
" <style>\n",
" .colab-df-container {\n",
" display:flex;\n",
" flex-wrap:wrap;\n",
" gap: 12px;\n",
" }\n",
"\n",
" .colab-df-convert {\n",
" background-color: #E8F0FE;\n",
" border: none;\n",
" border-radius: 50%;\n",
" cursor: pointer;\n",
" display: none;\n",
" fill: #1967D2;\n",
" height: 32px;\n",
" padding: 0 0 0 0;\n",
" width: 32px;\n",
" }\n",
"\n",
" .colab-df-convert:hover {\n",
" background-color: #E2EBFA;\n",
" box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
" fill: #174EA6;\n",
" }\n",
"\n",
" [theme=dark] .colab-df-convert {\n",
" background-color: #3B4455;\n",
" fill: #D2E3FC;\n",
" }\n",
"\n",
" [theme=dark] .colab-df-convert:hover {\n",
" background-color: #434B5C;\n",
" box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
" filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
" fill: #FFFFFF;\n",
" }\n",
" </style>\n",
"\n",
" <script>\n",
" const buttonEl =\n",
" document.querySelector('#df-9d579f69-9612-4eaa-87d7-fea3558c2d6d button.colab-df-convert');\n",
" buttonEl.style.display =\n",
" google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
"\n",
" async function convertToInteractive(key) {\n",
" const element = document.querySelector('#df-9d579f69-9612-4eaa-87d7-fea3558c2d6d');\n",
" const dataTable =\n",
" await google.colab.kernel.invokeFunction('convertToInteractive',\n",
" [key], {});\n",
" if (!dataTable) return;\n",
"\n",
" const docLinkHtml = 'Like what you see? Visit the ' +\n",
" '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
" + ' to learn more about interactive tables.';\n",
" element.innerHTML = '';\n",
" dataTable['output_type'] = 'display_data';\n",
" await google.colab.output.renderOutput(dataTable, element);\n",
" const docLink = document.createElement('div');\n",
" docLink.innerHTML = docLinkHtml;\n",
" element.appendChild(docLink);\n",
" }\n",
" </script>\n",
" </div>\n",
" </div>\n",
" "
]
},
"metadata": {},
"execution_count": 181
}
]
},
{
"cell_type": "code",
"source": [
"cases = {}\n",
"\n",
"for sample_id in range(8):\n",
" a, b, y = samples[sample_id]\n",
" x_tf = tf_res[sample_id]\n",
" x_sp = sp_res[sample_id]\n",
" data = np.column_stack([a, b, y, x_tf, x_sp])\n",
" df = pd.DataFrame(data=data, columns=['a', 'b', 'y', 'x_tf', 'x_sp'])\n",
" df['aerr'] = np.abs(df.x_tf - df.x_sp)\n",
" df['rerr'] = df.aerr / df.x_sp\n",
" df['aerr_y_tf'] = np.abs(tfp_math.betainc(df.a, df.b, df.x_tf) - df.y)\n",
" df['rerr_y_tf'] = df.aerr_y_tf / df.y\n",
" df['aerr_y_sp'] = np.abs(sp_special.betainc(df.a, df.b, df.x_sp) - df.y)\n",
" df['rerr_y_sp'] = df.aerr_y_sp / df.y\n",
" df['rtol'] = np.maximum(df.aerr - 1e-12, 0) / df.x_sp\n",
" df.sort_values(by='aerr', ascending=False, inplace=True)\n",
" cases[sample_id] = df"
],
"metadata": {
"id": "Ra2HBhqCZypb"
},
"execution_count": 182,
"outputs": []
},
{
"cell_type": "code",
"source": [
"for sample_id in range(8):\n",
" max_rtol = cases[sample_id].rtol.max()\n",
" print(f'sample {sample_id} | max_rtol: {max_rtol}')"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "AwynJH8wZ05b",
"outputId": "4526df08-bc8a-4bcf-f650-bf199bfc67bd"
},
"execution_count": 183,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"sample 0 | max_rtol: 0.009682435542345047\n",
"sample 1 | max_rtol: 0.017569657415151596\n",
"sample 2 | max_rtol: 0.03620864078402519\n",
"sample 3 | max_rtol: 0.024280576035380363\n",
"sample 4 | max_rtol: 0.08601529151201248\n",
"sample 5 | max_rtol: 0.08581182360649109\n",
"sample 6 | max_rtol: 0.23033356666564941\n",
"sample 7 | max_rtol: 0.21363873779773712\n"
]
}
]
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment