{
"cells": [
{
"cell_type": "code",
"execution_count": 37,
"source": [
"import Pkg; Pkg.add(\"SpecialFunctions\")\n",
"using SpecialFunctions, Plots\n",
"default(size=(480,200))"
],
"outputs": [
{
"output_type": "stream",
"name": "stderr",
"text": [
"\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
"\u001b[32m\u001b[1m No Changes\u001b[22m\u001b[39m to `~/.julia/environments/v1.6/Project.toml`\n",
"\u001b[32m\u001b[1m No Changes\u001b[22m\u001b[39m to `~/.julia/environments/v1.6/Manifest.toml`\n"
]
}
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Let's take a [standard](https://erikbern.com/2018/03/27/waiting-time-load-factor-and-queueing-theory.html) [result](https://www.sigarch.org/three-other-models-of-computer-system-performance-part-2) from queueing theory and say that the latency is a function of the load factor, $f$, defined as the ratio of arrival rate to service rate.\n",
"\n",
"$$Latency \\propto \\frac{1}{1-f}$$\n",
"\n",
"Latencies fly off to infinity when the load factor approaches 1. Makes sense because if the requests arrive faster than they can be serviced then the queue will grow until something else gives.\n",
"\n",
"Now we might assume there is some dispersal in the load factor due to the natural burstiness of incoming traffic. We don't want to model this noise with a normal distribution because our load factor is strictly between 0 and 1. A better choice is the [beta distribution](https://en.wikipedia.org/wiki/Beta_distribution) which has exactly that property and just so happens to nicely approximate the normal distribution if we [choose the right parameters](https://www.johndcook.com/blog/normal_approx_to_beta/).\n",
"\n",
"Now what happens to the distribution of latencies? Well, luck would have it that the transformation of the beta distribution by the formula above gives us the so called [beta-prime distribution](https://en.wikipedia.org/wiki/Beta_prime_distribution) (up to a constant translation representing the minimum latency). And that approximates a LogNormal distribution surprisingly well, especially around the tail."
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": 38,
"source": [
"p=10; q=10\n",
"betad(x) = (x^(p-1) * (1 - x)^(q-1))/beta(p, q)\n",
"\n",
"plot(0:0.01:1, betad, label=\"Beta($p,$q)\")"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"image/svg+xml": "\n\n"
},
"metadata": {},
"execution_count": 38
}
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": 39,
"source": [
"betaprime(x) = (x^(p - 1) * (x + 1)^(-p - q))/beta(p, q)\n",
"m=p/(q-1)\n",
"v=(p*(p + q - 1))/((q - 2)*(q - 1)^2)\n",
"\n",
"plot(0:0.01:4, ylim=1.0, betaprime, label=\"Beta Prime\")\n",
"\n",
"a = log(m^2/sqrt(v+m^2))\n",
"s = sqrt(log((v+m^2)/m^2))\n",
"lognormal(x) = (x*s*√(2*π))^-1 * exp(-(log(x) - a)^2/(2s^2))\n",
"\n",
"plot!(0:0.01:4, lognormal, label=\"LogNormal\")\n",
"\n",
"k=m^2/v; l=m/v\n",
"gammad(x) = (l^k)*x^(k-1)*exp(-l*x)/gamma(k)\n",
"\n",
"plot!(0:0.01:4, ylim=1.0, gammad, label=\"Gamma\")\n",
"\n",
"\n",
"vline!([m], label=\"Mean\")\n"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"image/svg+xml": "\n\n"
},
"metadata": {},
"execution_count": 39
}
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"The LogNormal is a good approximation to the BetaPrime for this situation."
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": 40,
"source": [
"plot(0:0.01:5, x -> lognormal(x)-betaprime(x), label=\"Absolute error (BetaPrime)\")\n",
"plot!(0:0.01:5, x -> lognormal(x)-gammad(x), label=\"Absolute error (Gamma)\")"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"image/svg+xml": "\n\n"
},
"metadata": {},
"execution_count": 40
}
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": 41,
"source": [
"p=6; q=60\n",
"\n",
"plot(0:0.01:1, betad, label=\"Beta($p,$q)\")"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"image/svg+xml": "\n\n"
},
"metadata": {},
"execution_count": 41
}
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": 42,
"source": [
"m=p/(q-1)\n",
"v=(p*(p + q - 1))/((q - 2)*(q - 1)^2)\n",
"\n",
"plot(0:.001:.5, ylim=1.0, betaprime, label=\"Beta Prime\")\n",
"\n",
"a = log(m^2/sqrt(v+m^2))\n",
"s = sqrt(log((v+m^2)/m^2))\n",
"lognormal(x) = (x*s*√(2*π))^-1 * exp(-(log(x) - a)^2/(2s^2))\n",
"\n",
"plot!(0:.001:.5, lognormal, label=\"LogNormal\")\n",
"\n",
"k=m^2/v; l=m/v\n",
"gammad(x) = (l^k)*x^(k-1)*exp(-l*x)/gamma(k)\n",
"\n",
"plot!(0:.001:.5, ylim=1.0, gammad, label=\"Gamma\")\n",
"\n",
"\n",
"vline!([m], label=\"Mean\")\n"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"image/svg+xml": "\n\n"
},
"metadata": {},
"execution_count": 42
}
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Gamma seems to be a better approximatio to the BetaPrime when the load and its variance is lower."
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": 43,
"source": [
"plot(0:.001:.5, x -> lognormal(x)-betaprime(x), label=\"Absolute error (LogNormal\")\n",
"plot!(0:.001:.5, x -> gammad(x)-betaprime(x), label=\"Absolute error (Gamma)\")"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"image/svg+xml": "\n\n"
},
"metadata": {},
"execution_count": 43
}
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"# Derivation of the transformation of beta distribution to beta prime distribution\n",
"\n",
"We want to apply the transformation $X \\rightarrow \\frac{1}{1-X}$ to the [beta distribution](https://en.wikipedia.org/wiki/Beta_distribution). Using the CDF [technique](https://math.stackexchange.com/a/13851) we have:\n",
"$$\n",
"\\begin{align*}\n",
"CDF_X(x) & = P(X \\leq x) \\rightarrow P(\\frac{1}{1-X} \\leq x) \\\\\n",
"& = P(1 \\leq x(1-X)) = P(X \\leq \\frac{x-1}{x}) \\\\\n",
"& = \\frac{1}{B(\\alpha,\\beta)} \\int_{0}^{\\frac{x-1}{x}} x^{\\alpha-1} (1-x)^{\\beta-1} dx\n",
"\\end{align*}\n",
"$$\n",
"\n",
"Differentiating both sides:\n",
"$$\n",
"\\begin{align*}\n",
"PDF_X(x) = \\frac{d}{dx} CDF_X(x) \n",
"& = \\frac{1}{B(\\alpha,\\beta)} (\\frac{x-1}{x})^{\\alpha-1} (1-\\frac{x-1}{x})^{\\beta-1} \\frac{d}{dx}(\\frac{x-1}{x}) \\\\\n",
"& = \\frac{1}{B(\\alpha,\\beta)} (x-1)^{\\alpha-1} x^{1-\\alpha} x^{1-\\beta} \\frac{1}{x^2} \\\\\n",
"& = \\frac{1}{B(\\alpha,\\beta)} (x-1)^{\\alpha-1} x^{-\\alpha-\\beta}\n",
"\\end{align*}\n",
"$$\n",
"\n",
"If we shift the distribution by one (representing the normalised minimum latency, i.e. a fixed service time) then $x \\rightarrow x-1$\n",
"$$\n",
"PDF_X(x) = \\frac{1}{B(\\alpha,\\beta)} x^{\\alpha-1} (x+1)^{-\\alpha-\\beta}\n",
"$$\n",
"Which is the beta prime distribution."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Inverse log approximation"
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": 44,
"source": [
"plot(0:0.01:0.99, x -> x/(1-x))\n",
"\n",
"plot!(0:0.01:0.99, x -> -1/log(x))"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"image/svg+xml": "\n\n"
},
"metadata": {},
"execution_count": 44
}
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": 45,
"source": [
"plot(0:0.01:0.99, x -> (x/(1-x) - -1/log(x)))"
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"image/svg+xml": "\n\n"
},
"metadata": {},
"execution_count": 45
}
],
"metadata": {}
}
],
"metadata": {
"orig_nbformat": 4,
"language_info": {
"name": "julia",
"version": "1.6.0",
"mimetype": "application/julia",
"file_extension": ".jl"
},
"kernelspec": {
"display_name": "Julia 1.6.0",
"language": "julia",
"name": "julia-1.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}