{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# How to generate random real numbers with a power-law distribution?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notes from Appendix D of:
\n",
"Clauset, A., Shalizi, C. R., & Newman, M. E. J. (2007). Power-law distributions in empirical data. SIAM Review, 51(4), 43. https://doi.org/10.1109/ICPC.2008.18"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, generate set of random real numbers $r$ that are uniformaly distributed in the interval $[0, 1)$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from math import floor, ceil"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"r = np.random.uniform(low=0.0, high=1.0, size=1000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's plot to see if they're uniformly distributed."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"NumBins = 10\n",
"MyBins = np.linspace(0, 1, NumBins+1)\n",
"\n",
"r_hist, r_bins = np.histogram(r, bins=MyBins, density=True)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFWCAYAAACre65zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAATVUlEQVR4nO3df6xkd3nf8c+DXYfwy1Tspq38g3VS07JyIzlaOUSkxSmktUGyVRVSu6WNU5eNkppUMkW4SetSRyoEGqVNawKmRE6ixj9CVboNSx0lMQLRGHnBYMV23W6ME28d4g0xphUxxuLpHzPA9eWud+7du3Pmu/f1kq40Z+a7c589vuu3zrkzZ6q7AwCstudMPQAAcHyCDQADEGwAGIBgA8AABBsABiDYADCA06ce4GTYtWtX79mzZ+oxAGBTPvWpT/1xd+/e6LFTMth79uzJoUOHph4DADalqn7/WI85JQ4AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGIBgA8AABBsABnBKfvgHwCrbc92Hpx7hGR5+5+umHoEFOMIGgAEINgAMQLABYAB+hz04vwsD2BkcYQPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQAD8D7sBazSe529z3nnWKWfu8TPHkzNETYADMARNrAlzgDAcjnCBoABCDYADECwAWAAgg0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADmDTYVfWLVfVYVf3uMR6vqvr5qjpcVfdW1fcse0YAWAVTH2HfnOSSZ3n80iTnz7/2J/mFJcwEACtn0mB398eS/MmzLLk8yS/3zF1JXlxVf2E50wHA6pj6CPt4zkryyJrtI/P7AGBHWfVg1wb39YYLq/ZX1aGqOnT06NGTPBYALNeqB/tIknPWbJ+d5NGNFnb3Td29r7v37d69eynDAcCyrHqwDyT5B/NXi78iyRPd/YdTDwUAy3b6lN+8qm5JcnGSXVV1JMm/TPJnkqS735vkYJLXJjmc5MtJfmSaSQFgWpMGu7uvPM7jneQfL2kcAFhZq35KHACIYAPAECY9Jc6pb891H556hGd4+J2vm3oEgC1xhA0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAfh4TXa0Vfr4Tx/9ub38t926Vdp3yXj772RxhA0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAQg2AAxAsAFgAIINAAMQbAAYwOlTDwAAm7Hnug9PPcI3PPzO1y3teznCBoABCDYADECwAWAAkwe7qi6pqger6nBVXbfB4+dW1Z1VdU9V3VtVr51iTgCY0qTBrqrTktyY5NIke5NcWVV71y3750lu7+4Lk1yR5D3LnRIApjf1EfZFSQ5390Pd/VSSW5Ncvm5NJ3nR/PaZSR5d4nwAsBKmflvXWUkeWbN9JMn3rlvz9iS/UVVvTvL8JK9ZzmgAsDqmPsKuDe7rddtXJrm5u89O8tokv1JV3zJ3Ve2vqkNVdejo0aMnYVQAmM7UwT6S5Jw122fnW095X53k9iTp7t9J8twku9Y/UXff1N37unvf7t27T9K4ADCNqYN9d5Lzq+q8qjojsxeVHVi35g+SvDpJqurlmQXbITQAO8qkwe7up5Nck+SOJA9k9mrw+6rqhqq6bL7sLUneVFWfTXJLkqu6e/1pcwA4pU39orN098EkB9fdd/2a2/cneeWy5wKAVTL1KXEAYAGCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAQg2AAxAsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMADBBoABCDYADECwAWAAgg0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAQg2AAzg9M0srqpK8pokP5jkryU5N8muJH+a5LEkn0ny20kOdPf/2d5RAWDnWijYVfW8JD+R5Eczi3TNH3oys1B/e5LvTPJdSf52kn9XVf8tyc929//Y7qEBYKc57inxqvqRJP87yb/O7Ej6X2V2hP3i7n5ed5/d3S/JLP57k/zDJP85yaVJPl5Vt1XVuSfrLwAAO8EiR9gfSPKhJO/o7ruPtai7O8n/nH/dXFUvSvLDSa5LclWSG054WgDYoRYJ9r7u/vRmn7i7v5Tk31fV+5Ps2eyfBwC+6bjB3kqs1/35JzM76gYAtmjTb+uqqoeq6saTMQwAsLGtvA97d5IntnsQAODYthLs+zJ7+9a2qKpLqurBqjpcVdcdY80PVdX9VXVfVf3qdn1vABjFpi6cMvfzSf5jVX13d997It+8qk5LcmNmbxM7kuTuqjrQ3fevWXN+kn+W5JXd/XhVfceJfE8AGNFWgn0kyW8m+URVvS/J3Uk+n6TXL+zujx3nuS5Kcri7H0qSqro1yeVJ7l+z5k1Jbuzux+fP+dgWZgaAoW0l2B/NLM6V5NpsEOo1TjvOc52V5JE120eSfO+6NS9Lkqr6xPz53t7d/339E1XV/iT7k+Tcc12nBYBTy1aCfUOePdKbURvct/65T09yfpKLk5yd2dXTLujuLz7jD3XflOSmJNm3b992zQcAK2HTwe7ut2/j9z+S5Jw122cneXSDNXd191eTfK6qHsws4Me86hoAnGqm/njNu5OcX1XnVdUZSa5IcmDdmg8l+YEkqapdmZ0if2ipUwLAxCYNdnc/neSaJHckeSDJ7d19X1XdUFWXzZfdkeQLVXV/kjuTvLW7vzDNxAAwjeOeEp9/TOa/6O7PbPbJq+rbkvx4kie7+xc2WtPdB5McXHff9Wtud2Yvbrt2s98fAE4Vixxh/+Ukn6qqj1TV36mq5x7vD1TVy6vqHUk+l+RnkvzfE5wTAHa0RV50tjfJP0nyk0n+RpKnqurTSQ4l+cMkjyd5bpKXZBb3V2T2dq1K8htJ/ml3/+72jw4AO8cin9b11ST/pqrek+TvJbk6swuefN/Xl+SZb886muT9Sd5zoldCAwBmFn5bV3d/ObMQv7+qXpRZsM/N7Mj6T5M8luTe7r7vZAwKADvZVi6cku7+Umav3gYAlmBLwU6SqnpBkr+V5MIkZ2b2kZv3JPkv3f3/tmc8ACDZYrCr6g1J3pvkxXnm7687yb+tqh/t7g9uw3wAQLYQ7Kr6wSS3JPlakl/O7MNAPp/kz2d2RbK/m+SWqvpid//m9o0KADvXVo6wr0/ylSR/tbs/ve6xX6qq/5DkY/N1gg0A22Arlya9MMltG8Q6SdLdh5LcnuR7TmQwAOCbthLsr2R2wZRn8+h8HQCwDbYS7I8n+f7jrHllZqfFAYBtsJVgvy3JX6mqd1bV89c+UFXPr6p3JbkgyXXbMSAAsLUXnb0tyb1J3ppk//y64n+U5M9l9nvrMzM7un5b1dp3fKW7++oTGxcAdqatBPuqNbdfnOSvb7DmVfOvtTqz65ADAJu0lWCft+1TAADPatPB7u7fPxmDAADHtpUXnQEASybYADAAwQaAAQg2AAxAsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMADBBoABCDYADECwAWAAgg0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADmDzYVXVJVT1YVYer6rpnWff6quqq2rfM+QBgFUwa7Ko6LcmNSS5NsjfJlVW1d4N1L0zyE0k+udwJAWA1TH2EfVGSw939UHc/leTWJJdvsO6nk7wryZPLHA4AVsXUwT4rySNrto/M7/uGqrowyTnd/evLHAwAVsnUwa4N7utvPFj1nCQ/l+Qtx32iqv1VdaiqDh09enQbRwSA6U0d7CNJzlmzfXaSR9dsvzDJBUk+WlUPJ3lFkgMbvfCsu2/q7n3dvW/37t0ncWQAWL6pg313kvOr6ryqOiPJFUkOfP3B7n6iu3d1957u3pPkriSXdfehacYFgGlMGuzufjrJNUnuSPJAktu7+76quqGqLptyNgBYJadPPUB3H0xycN191x9j7cXLmAkAVs3Up8QBgAUINgAMQLABYACCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAQg2AAxAsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMADBBoABCDYADECwAWAAgg0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAQg2AAxAsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMIDJg11Vl1TVg1V1uKqu2+Dxa6vq/qq6t6p+q6peOsWcADClSYNdVacluTHJpUn2JrmyqvauW3ZPkn3d/d1JPpjkXcudEgCmN/UR9kVJDnf3Q939VJJbk1y+dkF339ndX55v3pXk7CXPCACTmzrYZyV5ZM32kfl9x3J1ko9s9EBV7a+qQ1V16OjRo9s4IgBMb+pg1wb39YYLq96YZF+Sd2/0eHff1N37unvf7t27t3FEAJje6RN//yNJzlmzfXaSR9cvqqrXJPmpJK/q7q8saTYAWBlTH2HfneT8qjqvqs5IckWSA2sXVNWFSd6X5LLufmyCGQFgcpMGu7ufTnJNkjuSPJDk9u6+r6puqKrL5sveneQFSX6tqj5TVQeO8XQAcMqa+pR4uvtgkoPr7rt+ze3XLH0oAFgxU58SBwAWINgAMADBBoABCDYADECwAWAAgg0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGIBgA8AABBsABiDYADAAwQaAAQg2AAxAsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMADBBoABCDYADECwAWAAgg0AAxBsABiAYAPAAAQbAAYg2AAwAMEGgAEINgAMQLABYACCDQADEGwAGMDkwa6qS6rqwao6XFXXbfD4t1XVbfPHP1lVe5Y/JQBMa9JgV9VpSW5McmmSvUmurKq965ZdneTx7v6LSX4uyc8sd0oAmN7UR9gXJTnc3Q9191NJbk1y+bo1lyf5pfntDyZ5dVXVEmcEgMlNHeyzkjyyZvvI/L4N13T300meSPKSpUwHACuiunu6b171hiR/s7v/0Xz77ye5qLvfvGbNffM1R+bbvzdf84V1z7U/yf755l9K8uAS/gqbtSvJH089xKDsu62z706M/bd19t3mvbS7d2/0wOnLnmSdI0nOWbN9dpJHj7HmSFWdnuTMJH+y/om6+6YkN52kObdFVR3q7n1TzzEi+27r7LsTY/9tnX23vaY+JX53kvOr6ryqOiPJFUkOrFtzIMkPz2+/Pslv95SnBQBgApMeYXf301V1TZI7kpyW5Be7+76quiHJoe4+kOQDSX6lqg5ndmR9xXQTA8A0pj4lnu4+mOTguvuuX3P7ySRvWPZcJ8lKn7Jfcfbd1tl3J8b+2zr7bhtN+qIzAGAxU/8OGwBYgGBvM5daPTEL7L9rq+r+qrq3qn6rql46xZyr6Hj7bs2611dVV5VX784tsu+q6ofmP3v3VdWvLnvGVbbAv9tzq+rOqrpn/m/3tVPMObzu9rVNX5m9cO73knxnkjOSfDbJ3nVrfjzJe+e3r0hy29Rzr8rXgvvvB5I8b377x+y/xffdfN0Lk3wsyV1J9k099yp8Lfhzd36Se5L82fn2d0w996p8Lbj/bkryY/Pbe5M8PPXcI345wt5eLrV6Yo67/7r7zu7+8nzzrszeu89iP3tJ8tNJ3pXkyWUOt+IW2XdvSnJjdz+eJN392JJnXGWL7L9O8qL57TPzrdfbYAGCvb1cavXELLL/1ro6yUdO6kTjOO6+q6oLk5zT3b++zMEGsMjP3cuSvKyqPlFVd1XVJUubbvUtsv/enuSNVXUks3cFvTls2uRv6zrFbHSkvP5l+Ius2akW3jdV9cYk+5K86qRONI5n3XdV9ZzMPu3uqmUNNJBFfu5Oz+y0+MWZndX5eFVd0N1fPMmzjWCR/Xdlkpu7+2er6vsyu7bGBd39tZM/3qnDEfb22sylVvNsl1rdoRbZf6mq1yT5qSSXdfdXljTbqjvevnthkguSfLSqHk7yiiQHvPAsyeL/bv9rd3+1uz+X2WcVnL+k+VbdIvvv6iS3J0l3/06S52Z2nXE2QbC3l0utnpjj7r/5ad33ZRZrv0f8pmfdd939RHfv6u493b0ns9//X9bdh6YZd6Us8u/2Q5m94DFVtSuzU+QPLXXK1bXI/vuDJK9Okqp6eWbBPrrUKU8Bgr2N5r+T/vqlVh9IcnvPL7VaVZfNl30gyUvml1q9Nskx336z0yy4/96d5AVJfq2qPlNV6//HsCMtuO/YwIL77o4kX6iq+5PcmeStve4TA3eqBfffW5K8qao+m+SWJFc5UNk8VzoDgAE4wgaAAQg2AAxAsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMADBBoABCDawkKraU1VdVTdX1cuq6raqeqyqvlZVF089H5zqfLwmsFnfleSTSf5Xkv+U5NuTfGnSiWAHEGxgs74/yTu6+yenHgR2Eh/+ASykqvYk+VySP0ryUp9FDsvld9jAZn1WrGH5BBvYrM9PPQDsRIINbJbfo8EEBBsABiDYADAAwQaAAQg2AAzA+7ABYACOsAFgAIINAAMQbAAYgGADwAAEGwAGINgAMADBBoABCDYADECwAWAAgg0AA/j/qGDMkq4xKH4AAAAASUVORK5CYII=\n",
"text/plain": [
"