Skip to content

Instantly share code, notes, and snippets.

@leandrolcampos

leandrolcampos/betaincinv_cephes32.ipynb Secret

Created August 3, 2022 00:15
Show Gist options
  • Save leandrolcampos/5dfbc66753310598ef531f59901a343b to your computer and use it in GitHub Desktop.
Save leandrolcampos/5dfbc66753310598ef531f59901a343b to your computer and use it in GitHub Desktop.
Download ZIP
betaincinv_cephes32.ipynb
Display the source blob
Display the rendered blob
Raw
betaincinv_cephes32.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "betaincinv_cephes32.ipynb",
"provenance": [],
"collapsed_sections": [
"rfxrWk_hRoyN",
"GC-0lYG0HYAp",
"d15s3qJNgkAc"
],
"authorship_tag": "ABX9TyN9w9yNL/ifaV18qUU29AuI",
"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/5dfbc66753310598ef531f59901a343b/betaincinv_cephes32.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",
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "c42caacd-4339-489c-e234-18d1e3e17102"
},
"execution_count": 1,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"\u001b[K |████████████████████████████████| 511.7 MB 6.5 kB/s \n",
"\u001b[K |████████████████████████████████| 438 kB 51.9 MB/s \n",
"\u001b[K |████████████████████████████████| 1.6 MB 44.6 MB/s \n",
"\u001b[K |████████████████████████████████| 5.8 MB 51.2 MB/s \n",
"\u001b[?25h"
]
}
]
},
{
"cell_type": "code",
"source": [
"!pip install -U -q tensorflow-probability"
],
"metadata": {
"id": "YgA81jppUk2T",
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "0617de78-528b-4bb4-dd90-118cb74ed9eb"
},
"execution_count": 2,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"\u001b[K |████████████████████████████████| 6.5 MB 31.7 MB/s \n",
"\u001b[?25h"
]
}
]
},
{
"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": [
"## Cephes implementation of `betaincinv` for `float32` [my port]"
],
"metadata": {
"id": "rfxrWk_hRoyN"
}
},
{
"cell_type": "code",
"source": [
"def cephes_bisection(a, b, y0, initial_candidate, use_bisection):\n",
" tiny = tf.constant(np.finfo(np.float32).tiny, dtype=tf.float32)\n",
" tolerance = 1e-3\n",
"\n",
" def bisection_step(should_stop, low, high, delta, error, candidate):\n",
" error_is_negative = (error < 0.)\n",
" new_low = tf.where(error_is_negative, candidate, low)\n",
" new_high = tf.where(~error_is_negative, candidate, high)\n",
" new_delta = tf.where(error_is_negative, 0.5, tf.math.square(delta))\n",
" new_delta = tf.where(tf.math.equal(new_delta, 0.), 0.5, new_delta)\n",
" new_candidate = new_delta * new_high + (1. - new_delta) * low\n",
" new_error = tf.math.betainc(a, b, new_candidate) - y0\n",
" yp = new_error / tf.math.maximum(y0, tiny)\n",
"\n",
" should_stop = should_stop | (tf.math.abs(yp) < tolerance)\n",
"\n",
" return should_stop, new_low, new_high, new_delta, new_error, new_candidate\n",
"\n",
" max_iterations = 20\n",
"\n",
" (_, _, _, _, _, result) = tf.while_loop(\n",
" cond=lambda stop, *_: tf.reduce_any(~stop),\n",
" body=bisection_step,\n",
" loop_vars=(\n",
" ~use_bisection,\n",
" tf.zeros_like(y0),\n",
" tf.ones_like(y0),\n",
" 0.5 * tf.ones_like(y0),\n",
" tf.math.betainc(a, b, initial_candidate) - y0,\n",
" initial_candidate),\n",
" maximum_iterations=max_iterations)\n",
"\n",
" return result"
],
"metadata": {
"id": "fXnplVzQcs8q"
},
"execution_count": 4,
"outputs": []
},
{
"cell_type": "code",
"source": [
"def cephes_newton(a, b, y0, initial_candidate, use_newton):\n",
" minlogf = -103.278929903431851103\n",
" machepf = 5.9604644775390625e-8\n",
" lgm = tf.math.lgamma(a + b) - tf.math.lgamma(a) - tf.math.lgamma(b)\n",
"\n",
" def newton_step(should_stop, i, candidate):\n",
" y = tf.math.betainc(a, b, candidate)\n",
" d = (a - 1.) * tf.math.log(candidate) + (b - 1.) * tf.math.log(\n",
" 1. - candidate) + lgm\n",
" d_lt_min_log = (d < minlogf)\n",
" new_candidate = tf.where(d_lt_min_log, 0., candidate)\n",
" should_stop = should_stop | d_lt_min_log\n",
" d = (y - y0) / tf.math.exp(d)\n",
" new_candidate = tf.where(should_stop, new_candidate, new_candidate - d)\n",
" new_candidate = tf.clip_by_value(new_candidate, 0., 1.)\n",
"\n",
" should_stop = should_stop | tf.math.equal(new_candidate, 0.) | tf.math.equal(\n",
" new_candidate, 1.)\n",
" should_stop = should_stop | ((i > 1) & (\n",
" tf.math.abs(d / new_candidate) < 256. * machepf))\n",
"\n",
" return should_stop, i + 1, new_candidate\n",
"\n",
" max_iterations = 10\n",
"\n",
" (_, _, result) = tf.while_loop(\n",
" cond=lambda stop, *_: tf.reduce_any(~stop),\n",
" body=newton_step,\n",
" loop_vars=(\n",
" ~use_newton,\n",
" 0,\n",
" initial_candidate),\n",
" maximum_iterations=max_iterations)\n",
"\n",
" return result"
],
"metadata": {
"id": "4kk9qbl3JFiB"
},
"execution_count": 5,
"outputs": []
},
{
"cell_type": "code",
"source": [
"def cephes_betaincinvf(a, b, y):\n",
" dtype = tf.float32\n",
" tiny = tf.constant(np.finfo(np.float32).tiny, dtype=dtype)\n",
" minlogf = -103.278929903431851103\n",
" machepf = 5.9604644775390625e-8\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",
" use_symmetry_relation = (y > 0.5)\n",
" a_orig = a\n",
" a = tf.where(use_symmetry_relation, b, a)\n",
" b = tf.where(use_symmetry_relation, a_orig, b)\n",
" y0 = tf.where(use_symmetry_relation, 1. - y, y)\n",
"\n",
" yp = -tf.math.ndtri(y)\n",
" lgm = (yp * yp - 3.) * 0.16666666666666667\n",
" x0 = 2. / (1. / (2. * a - 1.) + 1. / (2. * b - 1.))\n",
" y = (yp * tf.math.sqrt(x0 + lgm) / x0 \n",
" - ( 1. / (2. * b - 1.) - 1. / (2. * a - 1.))\n",
"\t * (lgm + 0.833333333333333333 - 2. / (3. * x0)))\n",
" y = tf.where((a < 1.) | (b < 1.), 0.5 * yp * yp, 2. * y)\n",
" y = tf.math.maximum(y, minlogf)\n",
"\n",
" initial_candidate = a / (a + b * tf.math.exp(y))\n",
" initial_error = tf.math.betainc(a, b, initial_candidate) - y0\n",
" yp = initial_error / y0\n",
"\n",
" use_bisection = tf.math.equal(y0, 0.) | tf.math.equal(y0, 1.) | (\n",
" tf.math.abs(yp) >= 0.1)\n",
" initial_candidate = cephes_bisection(\n",
" a, b, y0, initial_candidate, use_bisection)\n",
" \n",
" use_newton = tf.math.equal(y0, 0.) | tf.math.equal(y0, 1.) | tf.math.equal(\n",
" initial_candidate, 0.)\n",
" result = cephes_newton(a, b, y0, initial_candidate, use_newton)\n",
"\n",
" y = tf.where(use_symmetry_relation, 1. - y0, y0)\n",
" result = tf.where(use_symmetry_relation, 1. - result, result)\n",
" result = tf.clip_by_value(result, 0., 1.)\n",
"\n",
" result = tf.where(tf.equal(y, 0.) | tf.equal(y, 1.), y, result)\n",
"\n",
" return result"
],
"metadata": {
"id": "F7wGr55PjZoi"
},
"execution_count": 6,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## Sanity checks for `float32`"
],
"metadata": {
"id": "GC-0lYG0HYAp"
}
},
{
"cell_type": "code",
"source": [
"dtype = np.float32"
],
"metadata": {
"id": "Nzw2Z4EyYtcc"
},
"execution_count": 7,
"outputs": []
},
{
"cell_type": "code",
"source": [
"a = 1.2\n",
"b = 9.8\n",
"y = 0.05"
],
"metadata": {
"id": "QI5pr8VWYtcc"
},
"execution_count": 8,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "0c981ed2-2220-4404-9d0a-19fa73e51bc6",
"id": "1Jc84_MNYtcc"
},
"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": "e4b56136-2dd2-4db0-a35b-a8b31107421a",
"id": "XdZoldGlYtcd"
},
"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": "tsHJErcTYtcd"
},
"execution_count": 11,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = cephes_betaincinvf(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "784fca52-fe4f-45d0-ea47-8e70bd8f7c02",
"id": "XwDwhIzKYtcd"
},
"execution_count": 12,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=0.009365334>"
]
},
"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": "bd81f70f-f571-4994-a776-9799f0b20009",
"id": "CVLeA67cYtcd"
},
"execution_count": 13,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.009366301"
]
},
"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/"
},
"outputId": "f04a9916-3da0-4a39-c8f3-69eea8492354",
"id": "5_A4g3RfYtcd"
},
"execution_count": 14,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.00010321181"
]
},
"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/"
},
"outputId": "bc591db2-d1b4-47bd-ca63-80a420bb4540",
"id": "9LqoevwVYtcd"
},
"execution_count": 15,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.00011898577"
]
},
"metadata": {},
"execution_count": 15
}
]
},
{
"cell_type": "code",
"source": [
"a = 1.2\n",
"b = 9.8\n",
"y = 0.9"
],
"metadata": {
"id": "NK_PaL5SYtcd"
},
"execution_count": 16,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "8c65c36b-4c0b-44b3-f837-a52b632d0e22",
"id": "maOr4WwOYtcd"
},
"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": "06ebac7a-2c15-4208-cfe6-a9353a9229e4",
"id": "UKyvYHNZYtce"
},
"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": "zozfh1qTYtce"
},
"execution_count": 19,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = cephes_betaincinvf(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "7b535c78-f2cb-485a-e32e-79adb5cac5a6",
"id": "4r4vFk0_Ytce"
},
"execution_count": 20,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=0.2343002>"
]
},
"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": "388bcc70-c9a7-4f57-bebc-2b64637cdc39",
"id": "ZdqEdT4zYtce"
},
"execution_count": 21,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.23426864"
]
},
"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": "708352b5-2737-4e2f-aee9-afe49cb35f4a",
"id": "qOVrEZc4Ytce"
},
"execution_count": 22,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.00013471996"
]
},
"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": "aaef5f5b-f35f-40e7-8ae8-b5d4a3433c5c",
"id": "78_42wynYtce"
},
"execution_count": 23,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"4.2650438e-05"
]
},
"metadata": {},
"execution_count": 23
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.8\n",
"b = 1.7\n",
"y = 0.1"
],
"metadata": {
"id": "v0TD8biVYtce"
},
"execution_count": 24,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "8cb1c432-8c10-4615-d8b2-c9fd66603c0e",
"id": "ffumNbmFYtce"
},
"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": "d19930da-20c0-41ea-9e30-1a2f0c1ebcb4",
"id": "yJNdq2bLYtce"
},
"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": "0aKEOrPHYtcf"
},
"execution_count": 27,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = cephes_betaincinvf(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "092aa691-a943-4286-d1b0-e18bfd673823",
"id": "PYuvEtBwYtcf"
},
"execution_count": 28,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=0.032367595>"
]
},
"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": "87312bc5-47ba-4890-fb74-04548235efc8",
"id": "DX193sZiYtcf"
},
"execution_count": 29,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.032386925"
]
},
"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": "0d5e9125-1aee-4df0-ae15-fe1dc0b81667",
"id": "U-oruwvEYtcf"
},
"execution_count": 30,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.0005968622"
]
},
"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": "1fda4def-84aa-4a2d-b5f5-141b2a4e0403",
"id": "Lg31dLhhYtcf"
},
"execution_count": 31,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.00047154725"
]
},
"metadata": {},
"execution_count": 31
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.8\n",
"b = 1.7\n",
"y = 0.9"
],
"metadata": {
"id": "GwHlFs1vYtcf"
},
"execution_count": 32,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "0cfd3cc8-eaec-4cde-9315-6aa867a4397d",
"id": "tK6ycxpVYtcf"
},
"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": "138b35c7-90df-49d7-c213-4bad0733e1e1",
"id": "2gQEFO7VYtcf"
},
"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": "CRgWLhxfYtcf"
},
"execution_count": 35,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = cephes_betaincinvf(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "46f0aaa3-22c0-45ca-b961-36ea0a6b77f5",
"id": "8ESEzUH1Ytcf"
},
"execution_count": 36,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=0.69944507>"
]
},
"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": "aa88a1f7-f987-44ae-a1eb-455d8df24e7f",
"id": "PU5RIDmLYtcg"
},
"execution_count": 37,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.6994021"
]
},
"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": "ebc352a9-6a37-415f-f574-917cdb3a8cd6",
"id": "N_bgWT0kYtcg"
},
"execution_count": 38,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"6.1445266e-05"
]
},
"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": "4b706e84-4532-465f-a464-f911ce7ce490",
"id": "mfdoCSnNYtcg"
},
"execution_count": 39,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"2.7683047e-05"
]
},
"metadata": {},
"execution_count": 39
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.1\n",
"b = 0.7\n",
"y = 0.1"
],
"metadata": {
"id": "tW8uIcrdYtcg"
},
"execution_count": 40,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "197b6de8-9ece-4a35-c75b-73319554e76b",
"id": "AE6GIJbeYtcg"
},
"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": "8c9757d0-bf6f-42b5-c16d-b5dad33525cd",
"id": "xIXss844Ytcg"
},
"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": "kMWTtL4lYtcg"
},
"execution_count": 43,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = cephes_betaincinvf(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "9cdeb0d2-6f27-47f6-9457-fc18b1c5d1f4",
"id": "jbVgfZVuYtcg"
},
"execution_count": 44,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=1.7960505e-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": "d2ca2b13-c69b-4190-c8c1-1dcec06383cb",
"id": "TOahbeaVYtcg"
},
"execution_count": 45,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.8028194e-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": "2bf8ca41-416f-435d-90ba-ab8bfd3a772c",
"id": "7l3dqJzSYtcg"
},
"execution_count": 46,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.003754614"
]
},
"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": "57aa4fbd-8626-44c1-a249-5275c0090ef9",
"id": "FOqCOm6yYtch"
},
"execution_count": 47,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.0003761798"
]
},
"metadata": {},
"execution_count": 47
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.1\n",
"b = 0.7\n",
"y = 0.2"
],
"metadata": {
"id": "RAyFzGUBYtch"
},
"execution_count": 48,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "8811f91c-21ec-4587-9bad-38deefdbeb54",
"id": "iYts-njEYtch"
},
"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": "419e3617-114a-4e39-ab80-158e5713eff7",
"id": "TnTST3tdYtch"
},
"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": "LP06GvW2Ytch"
},
"execution_count": 51,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = cephes_betaincinvf(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "d5521e34-53ad-4791-8b19-c94454d867e2",
"id": "AojiasC4Ytci"
},
"execution_count": 52,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=1.860906e-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": "edb89d70-2a6d-4112-ce7e-76e9832fe953",
"id": "Vm424kynYtci"
},
"execution_count": 53,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.8460868e-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": "9994d162-3407-427b-faf5-147207ff6487",
"id": "Aqq2CIIwYtci"
},
"execution_count": 54,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.00802737"
]
},
"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": "62e1bfab-e23c-4999-f184-43987579c72a",
"id": "dp58WWKxYtci"
},
"execution_count": 55,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.0007998198"
]
},
"metadata": {},
"execution_count": 55
}
]
},
{
"cell_type": "code",
"source": [
"a = 0.1\n",
"b = 0.7\n",
"y = 0.9"
],
"metadata": {
"id": "0HIEFHsUYtci"
},
"execution_count": 56,
"outputs": []
},
{
"cell_type": "code",
"source": [
"mean = a / (a + b)\n",
"mean"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "eab636b3-e100-4785-a289-ad04964b58d0",
"id": "614IxAryYtci"
},
"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": "73e6eb47-346e-4e42-db77-5336fff5f70c",
"id": "XsYrSd_qYtcj"
},
"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": "GNIGIxOAYtcj"
},
"execution_count": 59,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_tf = cephes_betaincinvf(a, b, y)\n",
"x_tf"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "a52e3c40-adb7-47fc-9daf-b00353ccde9f",
"id": "2-m6ewESYtcj"
},
"execution_count": 60,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=0.5262923>"
]
},
"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": "d6e38a40-1d6b-453d-e2d0-a1644c71a450",
"id": "PS34OHCYYtcj"
},
"execution_count": 61,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.5259731"
]
},
"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": "461dd61b-58c1-41b9-b614-c16365185361",
"id": "dT6XkN0xYtcj"
},
"execution_count": 62,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0.0006069559"
]
},
"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": "be769ed2-3430-4c61-d779-f265ffb50b93",
"id": "Kha7NgpJYtcj"
},
"execution_count": 63,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"7.470449e-05"
]
},
"metadata": {},
"execution_count": 63
}
]
},
{
"cell_type": "code",
"source": [
"cephes_betaincinvf(a, b, 0.)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "fc5823d2-4643-4a01-9c97-5c7533aa120e",
"id": "fSYWwLuRYtcj"
},
"execution_count": 64,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=0.0>"
]
},
"metadata": {},
"execution_count": 64
}
]
},
{
"cell_type": "code",
"source": [
"cephes_betaincinvf(a, b, 1.)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "542a160c-7a4e-45dd-e838-9595148d384f",
"id": "gbXR4YuXYtcj"
},
"execution_count": 65,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<tf.Tensor: shape=(), dtype=float32, numpy=1.0>"
]
},
"metadata": {},
"execution_count": 65
}
]
},
{
"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": 66,
"outputs": []
},
{
"cell_type": "code",
"source": [
"len(SEEDS)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "2HR-JJbhXtAp",
"outputId": "69a8c935-1469-449f-8041-bcc30c97772e"
},
"execution_count": 67,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"10"
]
},
"metadata": {},
"execution_count": 67
}
]
},
{
"cell_type": "code",
"source": [
"SIZE_PER_SEED = 50_000\n",
"# SIZE = SIZE_PER_SEED * len(SEEDS)"
],
"metadata": {
"id": "d0UUVrpW7xIC"
},
"execution_count": 68,
"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": 69,
"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": 70,
"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": 71,
"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": 72,
"outputs": []
},
{
"cell_type": "code",
"source": [
"sp_betaincinv = sp_special.betaincinv\n",
"tf_betaincinv = tf.function(\n",
" cephes_betaincinvf, autograph=False, jit_compile=False)"
],
"metadata": {
"id": "fUEgyjLQbBmV"
},
"execution_count": 73,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## Test results for `float32`"
],
"metadata": {
"id": "nl4hOD0Lg7OJ"
}
},
{
"cell_type": "code",
"source": [
"DTYPE = np.float32"
],
"metadata": {
"id": "eXU_VXWlWDqG"
},
"execution_count": 74,
"outputs": []
},
{
"cell_type": "code",
"source": [
"TINY = np.finfo(DTYPE).tiny\n",
"TINY"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "c7341d07-2b6b-4fe3-c398-23b88a45bc1b",
"id": "wRbuQNJ7WDqH"
},
"execution_count": 75,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.1754944e-38"
]
},
"metadata": {},
"execution_count": 75
}
]
},
{
"cell_type": "code",
"source": [
"EPS = np.finfo(DTYPE).eps\n",
"EPS"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "uci_ntVBhBsn",
"outputId": "195222a7-10a6-413f-a81e-8d921ed1b0bb"
},
"execution_count": 76,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1.1920929e-07"
]
},
"metadata": {},
"execution_count": 76
}
]
},
{
"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": 77,
"outputs": []
},
{
"cell_type": "code",
"source": [
"sp_res = [\n",
" sp_betaincinv(a, b, y) for a, b, y in samples]"
],
"metadata": {
"id": "-luHQfz6WDqH"
},
"execution_count": 78,
"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": "93190e32-792a-416c-9217-26fe65b97519"
},
"execution_count": 79,
"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: Cephes implementation of betaincinv for float32 [my port]')\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": "iF7MHErXZI-e",
"outputId": "73ed979f-0c22-457c-9744-ab5149166235"
},
"execution_count": 80,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Method: Cephes implementation of betaincinv for float32 [my port]\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.4e-38 1.00 9.6e-04 \n",
"1 uniform 5 500000 0 0 1.7e-38 1.00 9.7e-04 \n",
"2 uniform 10 500000 0 0 0 1.00 9.7e-04 \n",
"3 uniform 30 500000 0 0 0 1.00 9.8e-04 \n",
"4 uniform 100 500000 0 0 0 1.00 9.8e-04 \n",
"5 uniform 1000 500000 0 0 0 1.00 9.8e-04 \n",
"6 uniform 10000 500000 0 0 7.4e-38 1.00 9.8e-04 \n",
"7 uniform 100000 500000 0 0 1.6e-29 1.00 9.8e-04 \n",
"\n",
" Mean Abs Error Max Rel Error Mean Rel Error \n",
"0 2.7e-05 9.6e+14 4.3e+10 \n",
"1 2.7e-05 3.9e+14 3.4e+09 \n",
"2 3.4e-05 9.3e+13 7.2e+08 \n",
"3 4.6e-05 1.2e+13 5.2e+07 \n",
"4 5.5e-05 5.0e+11 3.1e+06 \n",
"5 1.1e-04 51819.14 0.14 \n",
"6 1.7e-04 1044.32 4.8e-03 \n",
"7 1.8e-04 121.73 2.2e-03 "
],
"text/html": [
"\n",
" <div id=\"df-bd8d79fd-8768-4cc3-964b-8e3051ace498\">\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.4e-38</td>\n",
" <td>1.00</td>\n",
" <td>9.6e-04</td>\n",
" <td>2.7e-05</td>\n",
" <td>9.6e+14</td>\n",
" <td>4.3e+10</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.7e-38</td>\n",
" <td>1.00</td>\n",
" <td>9.7e-04</td>\n",
" <td>2.7e-05</td>\n",
" <td>3.9e+14</td>\n",
" <td>3.4e+09</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>0</td>\n",
" <td>1.00</td>\n",
" <td>9.7e-04</td>\n",
" <td>3.4e-05</td>\n",
" <td>9.3e+13</td>\n",
" <td>7.2e+08</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>0</td>\n",
" <td>1.00</td>\n",
" <td>9.8e-04</td>\n",
" <td>4.6e-05</td>\n",
" <td>1.2e+13</td>\n",
" <td>5.2e+07</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>0</td>\n",
" <td>1.00</td>\n",
" <td>9.8e-04</td>\n",
" <td>5.5e-05</td>\n",
" <td>5.0e+11</td>\n",
" <td>3.1e+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>0</td>\n",
" <td>1.00</td>\n",
" <td>9.8e-04</td>\n",
" <td>1.1e-04</td>\n",
" <td>51819.14</td>\n",
" <td>0.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>7.4e-38</td>\n",
" <td>1.00</td>\n",
" <td>9.8e-04</td>\n",
" <td>1.7e-04</td>\n",
" <td>1044.32</td>\n",
" <td>4.8e-03</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.6e-29</td>\n",
" <td>1.00</td>\n",
" <td>9.8e-04</td>\n",
" <td>1.8e-04</td>\n",
" <td>121.73</td>\n",
" <td>2.2e-03</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>\n",
" <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-bd8d79fd-8768-4cc3-964b-8e3051ace498')\"\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-bd8d79fd-8768-4cc3-964b-8e3051ace498 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-bd8d79fd-8768-4cc3-964b-8e3051ace498');\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": 80
}
]
},
{
"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": 81,
"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": "OzHeiy9bZd1v",
"outputId": "7ccda82c-fdb2-467f-f3dd-a8995c21961e"
},
"execution_count": 82,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"sample 0 | max_rtol: 2.7157089877466057e+33\n",
"sample 1 | max_rtol: 1.6581545302019176e+32\n",
"sample 2 | max_rtol: 1.063581020455313e+32\n",
"sample 3 | max_rtol: 2.4670244551569195e+31\n",
"sample 4 | max_rtol: 8.569336481572322e+30\n",
"sample 5 | max_rtol: 5.83373747059461e+34\n",
"sample 6 | max_rtol: 1044.321533203125\n",
"sample 7 | max_rtol: 121.73307037353516\n"
]
}
]
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment