{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Bernoulli and binomial distribution\n", "- References: [Bernoulli distribution on Wikipedia](https://en.wikipedia.org/wiki/Bernoulli_distribution) and [Binomial distribution on Wikipedia](https://en.wikipedia.org/wiki/Binomial_distribution#Generating_binomial_random_variates).\n", "\n", "The Bernoulli distribution of mean $p\\in[0,1]$ is defined as the distribution on $\\{0,1\\}$ such that $\\mathbb{P}(X=1) = p$ and $\\mathbb{P}(X=0) = 1-p$.\n", "\n", "If $X$ follows a Binomial distribution of mean $p\\in[0,1]$ and $n$ samples, $X$ is defined as the sum of $n$ independent and identically distributed (iid) samples from a Bernoulli distribution of mean $p$, that is $X\\in\\{0,\\dots,n\\}$ ($X\\in\\mathbb{N}$) and $\\forall k\\in\\{0,\\dots,n\\}, \\mathbb{P}(X=k) = {n \\choose k} p^k (1-p)^{n-k}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Requirements\n", "Let's import the modules required for this notebook." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "%load_ext cython" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Lilian Besson (Naereen) 2019-02-28T14:09:13+01:00\n", "\n", "CPython 3.6.7\n", "IPython 7.2.0\n", "\n", "numpy 1.15.4\n", "matplotlib 3.0.2\n", "cython 0.29.2\n" ] } ], "source": [ "%load_ext watermark\n", "%watermark -a \"Lilian Besson (Naereen)\" -i -v -p numpy,matplotlib,cython" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A naive generator\n", "Using the pseudo-random generator of (`float`) random numbers in $[0,1]$ from the `random` or `numpy.random` module, we can easily generate a sample from a Bernoulli distribution." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "import random\n", "\n", "def uniform_01() -> float:\n", " return random.random()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0.02240955612218365,\n", " 0.7770015834088284,\n", " 0.8358045673716832,\n", " 0.4095374878779069,\n", " 0.987872833694527]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[ uniform_01() for _ in range(5) ]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's very quick now:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "def bernoulli(p: float) -> int:\n", " return 1 if uniform_01() <= p else 0" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0, 0, 0, 0, 0]" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[ bernoulli(0) for _ in range(5) ]" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1, 0, 0, 0, 0]" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[ bernoulli(0.12345) for _ in range(5) ]" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1, 1, 1, 1, 1]" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[ bernoulli(1) for _ in range(5) ]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So we can naively generate samples from a Binomial distribution by summing iid samples generated using this `bernoulli` function." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "def naive_binomial(n: int, p: float) -> int:\n", " result = 0\n", " for k in range(n): # sum of n iid samples from Bernoulli(p)\n", " result += bernoulli(p)\n", " return result" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example :" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[0, 1, 2, 0, 1]" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[ naive_binomial(10, 0.1) for _ in range(5) ]" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[8, 4, 6, 5, 8]" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[ naive_binomial(10, 0.5) for _ in range(5) ]" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[10, 9, 8, 10, 10]" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[ naive_binomial(10, 0.9) for _ in range(5) ]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can quickly illustrate the generated distribution, to check it has the correct \"shape\":" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Generated by Cython 0.29.2
\n", "\n",
" Yellow lines hint at Python interaction.
\n",
" Click on a line that starts with a \"+\" to see the C code that Cython generated for it.\n",
"
01: \n",
"+02: import random\n", "
__pyx_t_1 = __Pyx_Import(__pyx_n_s_random, 0, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_1);\n", " if (PyDict_SetItem(__pyx_d, __pyx_n_s_random, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)\n", " __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;\n", "
03: \n",
"+04: def cython_inversion_binomial(int n, double p) -> int:\n", "
/* Python wrapper */\n",
"static PyObject *__pyx_pw_46_cython_magic_beb04c003b39b7afc4e8a1cfd3d3afad_1cython_inversion_binomial(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/\n",
"static PyMethodDef __pyx_mdef_46_cython_magic_beb04c003b39b7afc4e8a1cfd3d3afad_1cython_inversion_binomial = {\"cython_inversion_binomial\", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_46_cython_magic_beb04c003b39b7afc4e8a1cfd3d3afad_1cython_inversion_binomial, METH_VARARGS|METH_KEYWORDS, 0};\n",
"static PyObject *__pyx_pw_46_cython_magic_beb04c003b39b7afc4e8a1cfd3d3afad_1cython_inversion_binomial(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {\n",
" int __pyx_v_n;\n",
" double __pyx_v_p;\n",
" PyObject *__pyx_r = 0;\n",
" __Pyx_RefNannyDeclarations\n",
" __Pyx_RefNannySetupContext(\"cython_inversion_binomial (wrapper)\", 0);\n",
" {\n",
" static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_n,&__pyx_n_s_p,0};\n",
" PyObject* values[2] = {0,0};\n",
" if (unlikely(__pyx_kwds)) {\n",
" Py_ssize_t kw_args;\n",
" const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);\n",
" switch (pos_args) {\n",
" case 2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);\n",
" CYTHON_FALLTHROUGH;\n",
" case 1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);\n",
" CYTHON_FALLTHROUGH;\n",
" case 0: break;\n",
" default: goto __pyx_L5_argtuple_error;\n",
" }\n",
" kw_args = PyDict_Size(__pyx_kwds);\n",
" switch (pos_args) {\n",
" case 0:\n",
" if (likely((values[0] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_n)) != 0)) kw_args--;\n",
" else goto __pyx_L5_argtuple_error;\n",
" CYTHON_FALLTHROUGH;\n",
" case 1:\n",
" if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_p)) != 0)) kw_args--;\n",
" else {\n",
" __Pyx_RaiseArgtupleInvalid(\"cython_inversion_binomial\", 1, 2, 2, 1); __PYX_ERR(0, 4, __pyx_L3_error)\n",
" }\n",
" }\n",
" if (unlikely(kw_args > 0)) {\n",
" if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, \"cython_inversion_binomial\") < 0)) __PYX_ERR(0, 4, __pyx_L3_error)\n",
" }\n",
" } else if (PyTuple_GET_SIZE(__pyx_args) != 2) {\n",
" goto __pyx_L5_argtuple_error;\n",
" } else {\n",
" values[0] = PyTuple_GET_ITEM(__pyx_args, 0);\n",
" values[1] = PyTuple_GET_ITEM(__pyx_args, 1);\n",
" }\n",
" __pyx_v_n = __Pyx_PyInt_As_int(values[0]); if (unlikely((__pyx_v_n == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 4, __pyx_L3_error)\n",
" __pyx_v_p = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_p == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 4, __pyx_L3_error)\n",
" }\n",
" goto __pyx_L4_argument_unpacking_done;\n",
" __pyx_L5_argtuple_error:;\n",
" __Pyx_RaiseArgtupleInvalid(\"cython_inversion_binomial\", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 4, __pyx_L3_error)\n",
" __pyx_L3_error:;\n",
" __Pyx_AddTraceback(\"_cython_magic_beb04c003b39b7afc4e8a1cfd3d3afad.cython_inversion_binomial\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
" __Pyx_RefNannyFinishContext();\n",
" return NULL;\n",
" __pyx_L4_argument_unpacking_done:;\n",
" __pyx_r = __pyx_pf_46_cython_magic_beb04c003b39b7afc4e8a1cfd3d3afad_cython_inversion_binomial(__pyx_self, __pyx_v_n, __pyx_v_p);\n",
"\n",
" /* function exit code */\n",
" __Pyx_RefNannyFinishContext();\n",
" return __pyx_r;\n",
"}\n",
"\n",
"static PyObject *__pyx_pf_46_cython_magic_beb04c003b39b7afc4e8a1cfd3d3afad_cython_inversion_binomial(CYTHON_UNUSED PyObject *__pyx_self, int __pyx_v_n, double __pyx_v_p) {\n",
" int __pyx_v_result;\n",
" double __pyx_v_q;\n",
" double __pyx_v_current_proba;\n",
" double __pyx_v_cum_proba;\n",
" double __pyx_v_one_uniform_sample;\n",
" PyObject *__pyx_r = NULL;\n",
" __Pyx_RefNannyDeclarations\n",
" __Pyx_RefNannySetupContext(\"cython_inversion_binomial\", 0);\n",
"/* … */\n",
" /* function exit code */\n",
" __pyx_L1_error:;\n",
" __Pyx_XDECREF(__pyx_t_2);\n",
" __Pyx_XDECREF(__pyx_t_3);\n",
" __Pyx_XDECREF(__pyx_t_4);\n",
" __Pyx_XDECREF(__pyx_t_5);\n",
" __Pyx_XDECREF(__pyx_t_6);\n",
" __Pyx_XDECREF(__pyx_t_7);\n",
" __Pyx_XDECREF(__pyx_t_9);\n",
" __Pyx_AddTraceback(\"_cython_magic_beb04c003b39b7afc4e8a1cfd3d3afad.cython_inversion_binomial\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
" __pyx_r = NULL;\n",
" __pyx_L0:;\n",
" __Pyx_XGIVEREF(__pyx_r);\n",
" __Pyx_RefNannyFinishContext();\n",
" return __pyx_r;\n",
"}\n",
"/* … */\n",
" __pyx_tuple_ = PyTuple_Pack(7, __pyx_n_s_n, __pyx_n_s_p, __pyx_n_s_result, __pyx_n_s_q, __pyx_n_s_current_proba, __pyx_n_s_cum_proba, __pyx_n_s_one_uniform_sample); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 4, __pyx_L1_error)\n",
" __Pyx_GOTREF(__pyx_tuple_);\n",
" __Pyx_GIVEREF(__pyx_tuple_);\n",
"/* … */\n",
" __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_beb04c003b39b7afc4e8a1cfd3d3afad_1cython_inversion_binomial, NULL, __pyx_n_s_cython_magic_beb04c003b39b7afc4); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 4, __pyx_L1_error)\n",
" __Pyx_GOTREF(__pyx_t_1);\n",
" if (PyDict_SetItem(__pyx_d, __pyx_n_s_cython_inversion_binomial, __pyx_t_1) < 0) __PYX_ERR(0, 4, __pyx_L1_error)\n",
" __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;\n",
"+05: if p <= 1e-9:\n", "
__pyx_t_1 = ((__pyx_v_p <= 1e-9) != 0);\n",
" if (__pyx_t_1) {\n",
"/* … */\n",
" }\n",
"+06: return 0\n", "
__Pyx_XDECREF(__pyx_r);\n", " __Pyx_INCREF(__pyx_int_0);\n", " __pyx_r = __pyx_int_0;\n", " goto __pyx_L0;\n", "
+07: if p >= 1 - 1e-9:\n", "
__pyx_t_1 = ((__pyx_v_p >= (1.0 - 1e-9)) != 0);\n",
" if (__pyx_t_1) {\n",
"/* … */\n",
" }\n",
"+08: return n\n", "
__Pyx_XDECREF(__pyx_r);\n", " __pyx_t_2 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 8, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_2);\n", " __pyx_r = __pyx_t_2;\n", " __pyx_t_2 = 0;\n", " goto __pyx_L0;\n", "
+09: if p > 0.5: # speed up by computing for q and then substracting\n", "
__pyx_t_1 = ((__pyx_v_p > 0.5) != 0);\n",
" if (__pyx_t_1) {\n",
"/* … */\n",
" }\n",
"+10: return n - cython_inversion_binomial(n, 1.0 - p)\n", "
__Pyx_XDECREF(__pyx_r);\n", " __pyx_t_2 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 10, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_2);\n", " __Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_n_s_cython_inversion_binomial); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 10, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_4);\n", " __pyx_t_5 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 10, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_5);\n", " __pyx_t_6 = PyFloat_FromDouble((1.0 - __pyx_v_p)); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 10, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_6);\n", " __pyx_t_7 = NULL;\n", " __pyx_t_8 = 0;\n", " if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {\n", " __pyx_t_7 = PyMethod_GET_SELF(__pyx_t_4);\n", " if (likely(__pyx_t_7)) {\n", " PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);\n", " __Pyx_INCREF(__pyx_t_7);\n", " __Pyx_INCREF(function);\n", " __Pyx_DECREF_SET(__pyx_t_4, function);\n", " __pyx_t_8 = 1;\n", " }\n", " }\n", " #if CYTHON_FAST_PYCALL\n", " if (PyFunction_Check(__pyx_t_4)) {\n", " PyObject *__pyx_temp[3] = {__pyx_t_7, __pyx_t_5, __pyx_t_6};\n", " __pyx_t_3 = __Pyx_PyFunction_FastCall(__pyx_t_4, __pyx_temp+1-__pyx_t_8, 2+__pyx_t_8); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 10, __pyx_L1_error)\n", " __Pyx_XDECREF(__pyx_t_7); __pyx_t_7 = 0;\n", " __Pyx_GOTREF(__pyx_t_3);\n", " __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;\n", " __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;\n", " } else\n", " #endif\n", " #if CYTHON_FAST_PYCCALL\n", " if (__Pyx_PyFastCFunction_Check(__pyx_t_4)) {\n", " PyObject *__pyx_temp[3] = {__pyx_t_7, __pyx_t_5, __pyx_t_6};\n", " __pyx_t_3 = __Pyx_PyCFunction_FastCall(__pyx_t_4, __pyx_temp+1-__pyx_t_8, 2+__pyx_t_8); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 10, __pyx_L1_error)\n", " __Pyx_XDECREF(__pyx_t_7); __pyx_t_7 = 0;\n", " __Pyx_GOTREF(__pyx_t_3);\n", " __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;\n", " __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;\n", " } else\n", " #endif\n", " {\n", " __pyx_t_9 = PyTuple_New(2+__pyx_t_8); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 10, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_9);\n", " if (__pyx_t_7) {\n", " __Pyx_GIVEREF(__pyx_t_7); PyTuple_SET_ITEM(__pyx_t_9, 0, __pyx_t_7); __pyx_t_7 = NULL;\n", " }\n", " __Pyx_GIVEREF(__pyx_t_5);\n", " PyTuple_SET_ITEM(__pyx_t_9, 0+__pyx_t_8, __pyx_t_5);\n", " __Pyx_GIVEREF(__pyx_t_6);\n", " PyTuple_SET_ITEM(__pyx_t_9, 1+__pyx_t_8, __pyx_t_6);\n", " __pyx_t_5 = 0;\n", " __pyx_t_6 = 0;\n", " __pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_9, NULL); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 10, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_3);\n", " __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0;\n", " }\n", " __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;\n", " __pyx_t_4 = PyNumber_Subtract(__pyx_t_2, __pyx_t_3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 10, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_4);\n", " __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;\n", " __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;\n", " __pyx_r = __pyx_t_4;\n", " __pyx_t_4 = 0;\n", " goto __pyx_L0;\n", "
+11: cdef int result = 0\n", "
__pyx_v_result = 0;\n",
"+12: cdef double q = 1.0 - p\n", "
__pyx_v_q = (1.0 - __pyx_v_p);\n",
"+13: cdef double current_proba = q**n\n", "
__pyx_v_current_proba = pow(__pyx_v_q, ((double)__pyx_v_n));\n",
"+14: cdef double cum_proba = current_proba\n", "
__pyx_v_cum_proba = __pyx_v_current_proba;\n",
"+15: cdef double one_uniform_sample = random.random()\n", "
__Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_random); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 15, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_3);\n", " __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_random); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 15, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_2);\n", " __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;\n", " __pyx_t_3 = NULL;\n", " if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_2))) {\n", " __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_2);\n", " if (likely(__pyx_t_3)) {\n", " PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_2);\n", " __Pyx_INCREF(__pyx_t_3);\n", " __Pyx_INCREF(function);\n", " __Pyx_DECREF_SET(__pyx_t_2, function);\n", " }\n", " }\n", " __pyx_t_4 = (__pyx_t_3) ? __Pyx_PyObject_CallOneArg(__pyx_t_2, __pyx_t_3) : __Pyx_PyObject_CallNoArg(__pyx_t_2);\n", " __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;\n", " if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 15, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_4);\n", " __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;\n", " __pyx_t_10 = __pyx_PyFloat_AsDouble(__pyx_t_4); if (unlikely((__pyx_t_10 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 15, __pyx_L1_error)\n", " __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;\n", " __pyx_v_one_uniform_sample = __pyx_t_10;\n", "
+16: while cum_proba < one_uniform_sample:\n", "
while (1) {\n",
" __pyx_t_1 = ((__pyx_v_cum_proba < __pyx_v_one_uniform_sample) != 0);\n",
" if (!__pyx_t_1) break;\n",
"+17: current_proba *= (p * (n - result)) / (q * (result + 1))\n", "
__pyx_t_10 = (__pyx_v_p * (__pyx_v_n - __pyx_v_result));\n",
" __pyx_t_11 = (__pyx_v_q * (__pyx_v_result + 1));\n",
" if (unlikely(__pyx_t_11 == 0)) {\n",
" PyErr_SetString(PyExc_ZeroDivisionError, \"float division\");\n",
" __PYX_ERR(0, 17, __pyx_L1_error)\n",
" }\n",
" __pyx_v_current_proba = (__pyx_v_current_proba * (__pyx_t_10 / __pyx_t_11));\n",
"+18: cum_proba += current_proba\n", "
__pyx_v_cum_proba = (__pyx_v_cum_proba + __pyx_v_current_proba);\n",
"+19: result += 1\n", "
__pyx_v_result = (__pyx_v_result + 1);\n",
" }\n",
"+20: return result\n", "
__Pyx_XDECREF(__pyx_r);\n", " __pyx_t_4 = __Pyx_PyInt_From_int(__pyx_v_result); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 20, __pyx_L1_error)\n", " __Pyx_GOTREF(__pyx_t_4);\n", " __pyx_r = __pyx_t_4;\n", " __pyx_t_4 = 0;\n", " goto __pyx_L0;\n", "