{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Когнитивные технологии\n", "\n", "*Алла Тамбовцева*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Генерирование выборок из различных распределений. Правдоподобие выборки." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Генерирование выборок из распространённых распределений" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Получить выборку заданного объёма из определённого распределения с некоторыми параметрами можно с помощью того же модуля `stats` из библиотеки `scipy`." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import scipy.stats as st" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Для этого нужно создать случайную величину с заданным распределением, сохранить её в переменную и применить к ней метод для создания выборки. Для примера создадим выборку из биномиального распределения. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "X = st.binom(n=10, p=0.7) # случайная величина\n", "X.rvs(8) # выборка объёма 8" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Понятно, что таким образом можно создать выборку из любого дискретного или непрерывного распределения." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0, 0, 4, 2, 0, 3, 1, 4, 2, 3])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# из распределения Пуассона\n", "Y = st.poisson(2)\n", "Y.rvs(10)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 3.28938486, 7.19821919, -1.48464074, -0.89904451, -0.24306394,\n", " 8.51202274, -1.65676397, 5.70952848, 7.7097788 , 0.33977637,\n", " -0.10745563, 4.93665913])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# из нормального распределения с a=2, sigma=3\n", "st.norm(2, 3).rvs(12)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "На всякий случай, вспомним, что для обеспечения воспроизводимости, нужно зафиксировать стартовую точку псевдослучайного алгоритма, который генерирует выборку (*seed*):" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 3.41430549, -1.57292708, 6.29812091, 1.06204431, -0.1617662 ,\n", " 4.66148882, 4.57876524, 0.09042949, 2.04708912, -4.72805486,\n", " 5.45010717, 4.97583807])" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from scipy import random as rd\n", "\n", "rd.seed(1234)\n", "st.norm(2, 3).rvs(12)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 3.41430549, -1.57292708, 6.29812091, 1.06204431, -0.1617662 ,\n", " 4.66148882, 4.57876524, 0.09042949, 2.04708912, -4.72805486,\n", " 5.45010717, 4.97583807])" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rd.seed(1234)\n", "st.norm(2, 3).rvs(12) # выборка не меняется" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Генерирование выборок из произвольных распределений" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Теперь рассмотрим более интересную задачу – сгенерируем произвольную дискретную случайную величину на основе массива значений и вероятностей, а затем возьмём из неё выборку. " ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "x = np.array([-1, 0, 5, 8])\n", "p = np.array([0.1, 0.4, 0.25, 0.25])\n", "\n", "W = st.rv_discrete(name='dist', values=(x, p)) # назовём распределение dist" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 8, 8, 0, -1, 0, 5, 8, 0, -1, 5, 8, 0, 0, 8, 5])" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "W.rvs(size=15) # выборка из 15 наблюдений" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Кстати, ранее это мы не обсуждали, но таким образом можно создать любую случайную величину и определить её характеристики:" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.15" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "W.expect()" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array(3.15), array(12.4275))" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "W.stats()" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.25" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "W.pmf(5) # P(W=5)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.75" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "W.cdf(5) # P(W<=5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "А теперь сгенерируем непрерывную случайную величину аналогичным образом, только вместо массива вероятностей будет функция плотности вероятности. Пусть \n", "\n", "$$\n", "f(x) = \n", "\\begin{cases}\n", "0.5 x, x \\in [0, 2] \\\\\n", "0, x \\notin [0, 2]\n", "\\end{cases}.\n", "$$\n", "\n", "Создадим класс `my_cont_dist`, который будет описывать случайную величину с непрерывным распределением и задавать её функцию распределения:" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [], "source": [ "class my_cont_dist(st.rv_continuous):\n", " def _pdf(self, x):\n", " return 0.5 * x\n", "\n", "CV = my_cont_dist(a=0, b=2, name='cont_dist')" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.3333333333333335" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "CV.expect() # математическое ожидание " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Построим график функции плотности на отрезке $x \\in [0, 2]$:" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [], "source": [ "x = np.linspace(0, 2)\n", "fx = CV.pdf(x)" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4VGX+/vH3h9CrNAWBUAQVsIBGbAiyoBQpP9uK2EkIXVFRURQRLKCrosIidV1dUVdsEwi9SC9BFBBFERQiIEEpAgIheX5/ZPA7ZgMZYCYnmblf15XLmXOembl9ONw5nJk5x5xziIhIZCnkdQAREQk9lbuISARSuYuIRCCVu4hIBFK5i4hEIJW7iEgEUrmLiEQglbuISARSuYuIRKDCXr1wpUqVXK1atbx6eRGRAmnVqlW7nHOVcxvnWbnXqlWLlJQUr15eRKRAMrOfghmnwzIiIhFI5S4iEoFU7iIiEUjlLiISgVTuIiIRKNdyN7OJZrbTzNYdZ72Z2etmttHM1pjZJaGPKSIiJyOYPfe3gDYnWN8WqOf/SQRGn34sERE5HbmWu3NuAfDbCYZ0At52WZYBZ5hZ1VAFFBGJFAfTD/LYrMf4aU9QH1U/LaE45l4N2BpwP9W/7H+YWaKZpZhZSlpaWgheWkSkYJi7eS4Xjr6QF5e8SPL3yWF/vVCUu+WwLMerbjvnxjrn4pxzcZUr5/rtWRGRAm/PoT1083Wj5dstKWSFmH/PfHpe1jPsrxuK0w+kAjUC7lcHtoXgeUVECjTfBh89p/Zkx/4dPHrVowy+djAlipTIk9cOxZ67D7jb/6mZK4C9zrntIXheEZECaeeBnXSe3JlO73eiYomKLE9YzvDrhudZsUMQe+5m9h5wLVDJzFKBp4EiAM65N4FkoB2wETgI3BeusCIi+ZlzjnfXvssD0x9g/5H9DG0xlEevfpSiMUXzPEuu5e6cuz2X9Q7oHbJEIiIF0Na9W+kxtQfJ3ydzRfUrmNBxAg0qN/Asj2en/BURiQSZLpMxKWN4bPZjZLgMRrQeQZ8mfYgpFONpLpW7iMgp+v7X70lISmDBTwtoVacVY9uPpXb52l7HAlTuIiIn7WjmUV5Z+gpPz3+aYjHFmNBxAvc1ug+znD4Z7g2Vu4jISfhqx1fE++JZtX0VN55/I6PajaJqmfz3pXyVu4hIEA4fPcyzC55l2OJhVChRgQ9v/ZCb69+cr/bWA6ncRURysXTrUuJ98Xyz6xvuvvhuXrn+FSqWrOh1rBNSuYuIHMf+I/t5cu6TvL78daqXrc60O6bRpu6JTpKbf6jcRURyMOuHWSROSeTHPT/S+7LevNDyBcoUK+N1rKCp3EVEAuz+Yzf9Z/Zn4pcTObfiuSy4dwHX1LzG61gnTeUuIuL3yTef0Cu5F2kH0ni86eMMaj6I4oWLex3rlKjcRSTq7di/g77T+jJ5/WQaVWnE1C5TuaRqwb5iqMpdRKKWc4531rxDv+n9OJh+kOf/9jz9r+pPkZgiXkc7bSp3EYlKP+35ie5TujPjhxlcVeMqJnScwPmVzvc6Vsio3EUkqmS6TEavHM2AOQNwzvFG2zfodVkvClkoLm+Rf6jcRSRqbNi1gYSkBBZtWUTrc1ozpv0Yap5R0+tYYaFyF5GIl56RzstLX2bw/MGULFKStzq9xd0X351vTx0QCip3EYloq7evJt4Xz+odq7mlwS280fYNqpSu4nWssFO5i0hEOnT0EEM/H8rwxcOpVLISH/39I26qf5PXsfKMyl1EIs7iLYuJ98Wz4dcN3NfoPl6+/mXKlyjvdaw8pXIXkYix/8h+npjzBCNXjCS2XCwz7pzB9edc73UsT6jcRSQizPxhJolJiWzZu4W+TfryXMvnKF20tNexPKNyF5EC7bc/fuPhmQ/z1pdvcX6l81l430Kujr3a61ieU7mLSIH10fqP6J3cm10HdzHwmoE82ezJAnuir1BTuYtIgbP99+30mdaHj7/5mMZVGjP9zuk0qtLI61j5ispdRAoM5xxvffkWD818iD/S/2BYy2E8fNXDFC6kKstOMyIiBcLm3ZtJnJLI7E2zuSb2GsZ3HM+5Fc/1Ola+pXIXkXwtIzODUStH8ficxylkhRjVbhQ94npE3Im+Qk3lLiL51jdp3xDvi2dp6lLa1m3Lm+3fJLZcrNexCgSVu4jkO+kZ6by4+EWGLBhC6aKleefGd7jjwjsi+kRfoaZyF5F8ZdW2VXT1dWXNL2v4e8O/80bbNziz1JlexypwgjpoZWZtzGyDmW00swE5rI81s3lmttrM1phZu9BHFZFI9kf6HwyYPYDLx19O2oE0PrntEz645QMV+ynKdc/dzGKAUcB1QCqw0sx8zrn1AcOeBP7rnBttZg2AZKBWGPKKSARa8NMCEnwJfP/b9yQ0TuCl61/ijOJneB2rQAtmz70JsNE5t8k5dwR4H+iUbYwDyvpvlwO2hS6iiESqfYf30Xtqb5q/1ZyjmUeZfddsxnUcp2IPgWCOuVcDtgbcTwUuzzZmMDDTzPoCpYBWIUknIhFr2vfT6D6lO6n7UnnwigcZ2mIopYqW8jpWxAhmzz2nt6ddtvu3A28556oD7YB3zP73Q6hmlmhmKWaWkpaWdvJpRaTA23VwF3d9chftJrWjTLEyLIlfwiutX1Gxh1gwe+6pQI2A+9X538Mu8UAbAOfcUjMrDlQCdgYOcs6NBcYCxMXFZf8FISIRzDnHh+s/pE9yH3Yf2s1TzZ5i4DUDKVa4mNfRIlIw5b4SqGdmtYGfgc5Al2xjtgAtgbfMrD5QHNCuuYgAsO33bfSa2ovPNnxG3NlxzO44m4vOusjrWBEt13J3zh01sz7ADCAGmOic+9rMhgApzjkf8DAwzsweJOuQzb3OOe2Zi0Q55xwTVk+g/8z+HM44zEvXvUS/K/rpRF95IKgZds4lk/XxxsBlgwJurwd0dnwR+dOm3ZvoltSNuZvn0rxmc8Z3HE/dCnW9jhU19OtTREIqIzOD15e/zsC5AylcqDCjbxhN4qWJOtFXHlO5i0jIfL3za+J98Sz/eTk31LuBN9u/SfWy1b2OFZVU7iJy2o5kHGHYomE8u+BZyhUvx6SbJtH5gs460ZeHVO4iclpW/rySeF88a3eupcuFXRjRegSVS1X2OlbUU7mLyCk5mH6Qp+c9zSvLXqFq6ar4OvvocF4Hr2OJn8pdRE7a/B/n0y2pGxt/20j3S7szvNVwyhUv53UsCaByF5Gg7T20l8dmP8aYVWM4p/w5zL17Li1qt/A6luRA5S4iQZn63VS6T+nO9v3b6X9lf55p8Qwli5T0OpYch8pdRE4o7UAa/Wb0Y9LaSVxw5gV8fNvHNKnWxOtYkguVu4jkyDnH++ve5/7p97P30F6eufYZBjQdQNGYol5HkyCo3EXkf6TuS6Xn1J5M+W4KTao1YULHCVxw5gVex5KToHIXkT9lukzGfzGeR2Y9QnpGOq9c/wr3X34/MYVivI4mJ0nlLiIAbPxtI92SujH/x/m0qNWCcR3GcU6Fc7yOJadI5S4S5Y5mHmXEshE8Ne8pisYUZVyHccQ3jtepAwo4lbtIFFv7y1riffGs3LaSjud15J/t/km1stW8jiUhoHIXiUKHjx7m+YXP8/yi5ylfvDwf3PIBtza4VXvrEUTlLhJllqUuI94Xz/q09dxx4R2MaDOCSiUreR1LQkzlLhIlDhw5wFPznmLEshFUK1uNqV2m0q5eO69jSZio3EWiwJxNc+iW1I3NezbTM64nw1oNo2yxsl7HkjBSuYtEsD2H9vDIzEcYv3o89SrU4/N7P6dZzWZex5I8oHIXiVCfffsZPaf2ZOeBnTx61aMMvnYwJYqU8DqW5BGVu0iE2XlgJ/dPu58Pvv6Ai8+6mKTbk7j07Eu9jiV5TOUuEiGcc7y79l0emP4A+4/s59kWz/Lo1Y9SJKaI19HEAyp3kQiwde9WekztQfL3yVxZ/UomdJxA/cr1vY4lHlK5ixRgmS6TMSljeHT2o2S6TEa0HkGfJn10oi9RuYsUVN/9+h0JvgQWbllIqzqtGNt+LLXL1/Y6luQTKneRAuZo5lFeWfoKT89/muKFizOx40TubXSvTh0gf6FyFylAvtrxFfG+eFZtX8WN59/IqHajqFqmqtexJB9SuYsUAIePHubZBc8ybPEwKpSowH9v+S+3NLhFe+tyXCp3kXxuydYlJPgS+GbXN9xz8T28fP3LVCxZ0etYks+p3EXyqf1H9jNwzkDeWPEGNcrVYPod02ldt7XXsaSAKBTMIDNrY2YbzGyjmQ04zpi/m9l6M/vazCaFNqZIdJn1wywuHH0hb6x4g96X9WZdz3Uqdjkpue65m1kMMAq4DkgFVpqZzzm3PmBMPeBx4Grn3G4zOzNcgUUi2e4/dvPwzIf515f/4ryK57HgvgU0jW3qdSwpgII5LNME2Oic2wRgZu8DnYD1AWO6AaOcc7sBnHM7Qx1UJNJ98s0n9EruRdqBNB5v+jiDmg+ieOHiXseSAiqYcq8GbA24nwpcnm3MuQBmthiIAQY756ZnfyIzSwQSAWJjY08lr0jE2bF/B32n9WXy+sk0qtKI5C7JNK7a2OtYUsAFU+45fdbK5fA89YBrgerAQjO7wDm35y8Pcm4sMBYgLi4u+3OIRBXnHG9/9TYPzniQg+kHef5vz9P/qv460ZeERDDlngrUCLhfHdiWw5hlzrl0YLOZbSCr7FeGJKVIhPlpz090n9KdGT/M4OoaVzO+43jOr3S+17EkggTzaZmVQD0zq21mRYHOgC/bmE+BFgBmVomswzSbQhlUJBJkukxGrhhJw382ZNGWRbzR9g0W3LdAxS4hl+ueu3PuqJn1AWaQdTx9onPuazMbAqQ453z+ddeb2XogA3jEOfdrOIOLFDQbdm0g3hfP4q2LaX1Oa8a0H0PNM2p6HUsilDnnzaHvuLg4l5KS4slri+Sl9Ix0Xl76MoPnD6ZkkZK82vpV7r74bp06QE6Jma1yzsXlNk7fUBUJo9XbVxPvi2f1jtXcXP9mRrYbSZXSVbyOJVFA5S4SBoeOHmLI50N4cfGLVCpZiY/+/hE31b/J61gSRVTuIiG2aMsiEnwJbPh1A/c1uo+Xr3+Z8iXKex1LoozKXSREfj/8O0/MeYJRK0cRWy6WGXfO4Ppzrvc6lkQplbtICEzfOJ3uU7qzde9W+jTpw/Mtn6d00dJex5IopnIXOQ2/HvyVh2Y+xNtfvc35lc5nUddFXFXjKq9jiajcRU6Fc46PvvmI3sm9+e2P33jymicZ2GygTvQl+YbKXeQkbf99O72Te/PJt59wadVLmXnnTC6ucrHXsUT+QuUuEiTnHP/68l88NOMhDmccZljLYTx81cMULqS/RpL/aKsUCcLm3ZtJnJLI7E2zaVazGeM6jOPciud6HUvkuFTuIieQkZnByBUjeWLuE8RYDKNvGE3ipYkUsqCuUCniGZW7yHGsT1tPgi+BpalLaVu3LWPaj6FGuRq5P1AkH1C5i2STnpHO8MXDGbpgKKWLluadG9/hjgvv0Im+pEBRuYsEWLVtFV19XVnzyxpua3gbr7d9nTNL6XrvUvCo3EWAP9L/YPD8wfxj6T84q9RZfHrbp3Q6v5PXsUROmcpdot6CnxaQ4Evg+9++J6FxAi9d/xJnFD/D61gip0XlLlFr3+F9DJg9gNEpo6l9Rm1m3zWblnVaeh1LJCRU7hKVkr9PpseUHqTuS+XBKx5kaIuhlCpayutYIiGjcpeosuvgLh6c8SD/WfMfGlRuwJL4JVxR/QqvY4mEnMpdooJzjg/Xf0if5D7sPrSbQc0G8cQ1T1CscDGvo4mEhcpdIt6237fRa2ovPtvwGXFnxzG742wuOusir2OJhJXKXSKWc44JqyfQf2Z/Dmcc5qXrXqLfFf10oi+JCtrKJSJt2r2JbkndmLt5Ls1rNmd8x/HUrVDX61gieUblLhElIzOD15e/zsC5AylcqDBj2o8h4ZIEnehLoo7KXSLGup3riPfFs+LnFdxQ7wbebP8m1ctW9zqWiCdU7lLgHck4wgsLX+C5hc9Rrng5Jt00ic4XdNaJviSqqdylQFv580q6+rqybuc6br/gdl5r8xqVS1X2OpaI51TuUiAdTD/IoHmDeHXZq1QtXRVfZx8dzuvgdSyRfEPlLgXO/B/nk+BL4IfdP9D90u4MbzWccsXLeR1LJF9RuUuBsffQXh6d9ShjvxjLOeXPYe7dc2lRu4XXsUTypaA+H2Zmbcxsg5ltNLMBJxh3i5k5M4sLXUQRSNqQRMN/NmT86vH0v7I/a3quUbGLnECue+5mFgOMAq4DUoGVZuZzzq3PNq4McD+wPBxBJTqlHUjjgekP8N6697jgzAv4+LaPaVKtidexRPK9YPbcmwAbnXObnHNHgPeBnC5RMxR4ETgUwnwSpZxzTFo7ifqj6jN5/WSeufYZViWuUrGLBCmYcq8GbA24n+pf9iczawzUcM5NCWE2iVKp+1Lp+H5H7vj4DupWqMvq7qsZ1HwQRWOKeh1NpMAI5g3VnL4J4v5caVYIeBW4N9cnMksEEgFiY2ODSyhRI9NlMm7VOB6Z9QgZLoNXW79K3yZ9iSkU43U0kQInmHJPBWoE3K8ObAu4Xwa4AJjv/0ZgFcBnZh2dcymBT+ScGwuMBYiLi3OI+G38bSMJvgQ+/+lzWtZuydgOY6lTvo7XsUQKrGDKfSVQz8xqAz8DnYEux1Y65/YClY7dN7P5QP/sxS6Sk6OZRxmxbARPzXuKYjHFGNdhHPGN43XqAJHTlGu5O+eOmlkfYAYQA0x0zn1tZkOAFOecL9whJTKt+WUN8b54Ural0Om8Tvzzhn9ydpmzvY4lEhGC+hKTcy4ZSM62bNBxxl57+rEkkh0+epjnFj7HC4teoHzx8nxwywfc2uBW7a2LhJC+oSp5alnqMuJ98axPW8+dF93JiNYjqFiyotexRCKOyl3yxIEjB3hy7pO8tvw1qpWtxtQuU2lXr53XsUQilspdwm7Opjl0S+rG5j2b6RnXk2GthlG2WFmvY4lENJW7hM3uP3bTf2Z/Jn45kXoV6vH5vZ/TrGYzr2OJRAWVu4TFZ99+Rs+pPdl5YCePXf0YTzd/mhJFSngdSyRqqNwlpH7Z/wv3T7+f/379Xy466yKSbk/i0rMv9TqWSNRRuUtIOOd4d+27PDD9AfYf2c+zLZ7l0asfpUhMEa+jiUQllbucti17t9BjSg+mbZzGldWvZELHCdSvXN/rWCJRTeUupyzTZfJmyps8NvsxMl0mr7V5jd6X9daJvkTyAZW7nJLvfv2OBF8CC7cspFWdVoxtP5ba5Wt7HUtE/FTuclKOZh7l5SUv8/T8rE+/TOw4kXsb3atTB4jkMyp3CdpXO76iq68rX2z/ghvPv5FR7UZRtUxVr2OJSA5U7pKrQ0cP8eyCZxm+eDgVS1Rk8q2TubnBzV7HEpETULnLCS3ZuoR4Xzzf7vqWey6+h1dav0KFEhW8jiUiuVC5S472H9nPwDkDeWPFG9QoV4Npd0yjTd02XscSkSCp3OV/zPxhJolJify09yf6XNaH51s+T5liZbyOJSInQeUuf9r9x24emvkQb335FudVPI+F9y2kaWxTr2OJyClQuQsAH3/zMb2Te5N2II3Hmz7OoOaDKF64uNexROQUqdyj3I79O+iT3IePvvmIRlUakdwlmcZVG3sdS0ROk8o9Sjnn+PdX/+ahGQ9xMP0gL7R8gYevfFgn+hKJECr3KPTjnh/pPqU7M3+YSdPYpozvMJ7zKp3ndSwRCSGVexTJdJmMWjGKx+c8jpkxsu1Iel7Wk0JWyOtoIhJiKvco8e2ub0nwJbB462Jan9OaMe3HUPOMml7HEpEwUblHuPSMdF5a8hLPfP4MpYqU4t//79/cddFdOtGXSIRTuUewL7Z/Qbwvni93fMktDW5hZNuRnFX6LK9jiUgeULlHoD/S/2DI50N4aclLVC5VmY///jE31r/R61gikodU7hFm0ZZFxPvi+e7X7+jaqCv/uP4flC9R3utYIpLHVO4R4vfDv/P4nMcZtXIUtc6oxay7ZtGqTiuvY4mIR1TuEWD6xukkJiWSui+Vfpf3Y+jfhlK6aGmvY4mIh1TuBdivB3/lwRkP8s6ad6hfqT6Luy7myhpXeh1LRPIBlXsB5Jzjw/Uf0ie5D7sP7eapZk8x8JqBFCtczOtoIpJPBPXVRDNrY2YbzGyjmQ3IYf1DZrbezNaY2Rwz07djwmT779u56b83cdvk24gtF0tKtxSGtBiiYheRv8i13M0sBhgFtAUaALebWYNsw1YDcc65i4DJwIuhDhrtnHNMXD2R+qPqM33jdF5s9SLLEpZxcZWLvY4mIvlQMIdlmgAbnXObAMzsfaATsP7YAOfcvIDxy4A7Qxky2m3evZnEKYnM3jSbZjWbMa7DOM6teK7XsUQkHwum3KsBWwPupwKXn2B8PDAtpxVmlggkAsTGxgYZMXplZGYwcsVInpj7BDEWw+gbRpN4aaJO9CUiuQqm3HM6CYnLcaDZnUAc0Dyn9c65scBYgLi4uByfQ7KsT1tPgi+BpalLaVu3LWPaj6FGuRpexxKRAiKYck8FAlulOrAt+yAzawUMBJo75w6HJl70Sc9IZ/ji4QxdMJQyRcvwnxv/Q5cLu+hEXyJyUoIp95VAPTOrDfwMdAa6BA4ws8bAGKCNc25nyFNGiZRtKcT74lnzyxpua3gbr7d9nTNLnel1LBEpgHItd+fcUTPrA8wAYoCJzrmvzWwIkOKc8wEvAaWBD/17mFuccx3DmDuiHEw/yOD5g3l56ctUKV2FT2/7lE7nd/I6logUYEF9ick5lwwkZ1s2KOC2TmJyij7/8XMSkhLY+NtGul3SjReve5Ezip/hdSwRKeD0DVWP7Du8j8dmPcabq96kTvk6zLl7Dn+r/TevY4lIhFC5e2Dqd1PpMbUH237fxkNXPMSQFkMoVbSU17FEJIKo3PNQ2oE0+s3ox6S1k2hYuSGTb53M5dVP9JUBEZFTo3LPA845Pvj6A/pO68veQ3t5uvnTPHHNExSNKep1NBGJUCr3MPt538/0nNqTpO+SuOzsy5jQcQIXnnWh17FEJMKp3MPEOcf4L8bTf1Z/0jPS+cd1/6DfFf2IKRTjdTQRiQIq9zD44bcf6JbUjXk/zuPaWtcyrsM46lao63UsEYkiKvcQysjM4LXlr/Hk3CcpElOEMe3HkHBJgk70JSJ5TuUeIut2riPeF8+Kn1fQ/tz2jL5hNNXLVvc6lohEKZX7aTqScYQXFr7Acwufo1zxcky6aRKdL+isE32JiKdU7qdhxc8riPfFs27nOrpc2IURrUdQuVRlr2OJiKjcT8XB9IM8NfcpRiwfQdXSVUm6PYn257b3OpaIyJ9U7idp3uZ5JCQlsGn3Jrpf2p3hrYZTrng5r2OJiPyFyj1Iew/t5ZFZjzDui3GcU/4c5t49lxa1W3gdS0QkRyr3ICRtSKLH1B7s2L+D/lf255kWz1CySEmvY4mIHJfK/QTSDqRx//T7eX/d+1x45oV8etunXFbtMq9jiYjkSuWeA+cc7617j/un3c++w/t45tpnGNB0gE70JSIFhso9m9R9qfSc2pMp303h8mqXM6HjBBqe2dDrWCIiJ0Xl7pfpMhm3ahyPzHqEDJfBq61fpW+TvjrRl4gUSCp3yLp+aVI35v84n5a1WzK2w1jqlK/jdSwRkVMW1eV+NPMoI5aN4Kl5T1EsphjjO4yna+OuOnWAiBR4UVvua35ZQ7wvnpRtKXQ8ryOjbxjN2WXO9jqWiEhIRF25Hz56mOcWPscLi16gfPHyfHDLB9za4FbtrYtIRImqcl+Wuox4Xzzr09Zz10V38WrrV6lYsqLXsUREQi4qyv3AkQM8OfdJXlv+GtXLVie5SzJt67X1OpaISNhEfLnP3jSbxKRENu/ZTM+4ngxrNYyyxcp6HUtEJKwittz3HNrDwzMeZuKXE6lXoR6f3/s5zWo28zqWiEieiMhy//TbT+k1tRc7D+xkwNUDGNR8ECWKlPA6lohInomocv9l/y/0ndaXD9d/SKMqjZjSZQqXVL3E61giInkuIsrdOcd/1vyHfjP6sf/Ifp7723M8ctUjFIkp4nU0ERFPFApmkJm1MbMNZrbRzAbksL6YmX3gX7/czGqFOujxbNm7hXaT2nH3p3dzXsXz+LL7lzxxzRMqdhGJarmWu5nFAKOAtkAD4HYza5BtWDyw2zlXF3gVGB7qoNllukxGrRhFw382ZOFPC3m9zessvG8h9SvXD/dLi4jke8EclmkCbHTObQIws/eBTsD6gDGdgMH+25OBkWZmzjkXwqx/2rBrAwlJCSzasojr6lzH2A5jqXVGrXC8lIhIgRRMuVcDtgbcTwUuP94Y59xRM9sLVAR2hSJkoImrJ9Jrai9KFCnBvzr9i3suvkenDhARySaYcs+pObPvkQczBjNLBBIBYmNjg3jp/3VuxXNpf257RrYbSZXSVU7pOUREIl0w5Z4K1Ai4Xx3YdpwxqWZWGCgH/Jb9iZxzY4GxAHFxcad0yKZpbFOaxjY9lYeKiESNYD4tsxKoZ2a1zawo0BnwZRvjA+7x374FmBuu4+0iIpK7XPfc/cfQ+wAzgBhgonPuazMbAqQ453zABOAdM9tI1h5753CGFhGREwvqS0zOuWQgOduyQQG3DwG3hjaaiIicqqC+xCQiIgWLyl1EJAKp3EVEIpDKXUQkAqncRUQikHn1cXQzSwN+OsWHVyIMpzYIAeU6Ocp18vJrNuU6OaeTq6ZzrnJugzwr99NhZinOuTivc2SnXCdHuU5efs2mXCcnL3LpsIyISARSuYuIRKCCWu5jvQ5wHMp1cpTr5OXXbMp1csKeq0AecxcRkRMrqHvuIiJyAvmu3E/nYtxm9rh/+QYza53HuR4ys/VmtsbM5phZzYB1GWb2pf8n++mSw53rXjNLC3j9hIB195jZ9/6fe7I/Nsy5Xg3I9J2Z7QlYF875mmhmO81s3XHWm5m97s+9xswuCVgXlvkKItMd/ixrzGyJmV0csO5HM1vrn6uUUGU6iWzXmtnegD+vQQHrTrgNhDnXIwGZ1vlv/m+1AAAEi0lEQVS3qQr+dWGZMzOrYWbzzOwbM/vazB7IYUzebV/OuXzzQ9YphX8A6gBFga+ABtnG9ALe9N/uDHzgv93AP74YUNv/PDF5mKsFUNJ/u+exXP77+z2cr3uBkTk8tgKwyf/f8v7b5fMqV7bxfck6lXRY58v/3M2AS4B1x1nfDphG1tXFrgCW58F85ZbpqmOvRdaF6pcHrPsRqOThfF0LTDndbSDUubKN7UDWNSbCOmdAVeAS/+0ywHc5/H3Ms+0rv+25/3kxbufcEeDYxbgDdQL+7b89GWhpZuZf/r5z7rBzbjOw0f98eZLLOTfPOXfQf3cZWVesCrdg5ut4WgOznHO/Oed2A7OANh7luh14L0SvfULOuQXkcJWwAJ2At12WZcAZZlaVMM5Xbpmcc0v8rwl5t20de+3c5ut4TmfbDHWuPNm+nHPbnXNf+G//DnxD1vWlA+XZ9pXfyj2ni3Fnn5y/XIwbOHYx7mAeG85cgeLJ+u18THEzSzGzZWb2/0KU6WRy3ez/J+BkMzt2ycR8MV/+w1e1gbkBi8M1X8E4XvZwztfJyL5tOWCmma2yrGsUe+FKM/vKzKaZWUP/snwxX2ZWkqyS/ChgcdjnzLIOFzcGlmdblWfbV1AX68hDp3Mx7qAu0n2Kgn5uM7sTiAOaByyOdc5tM7M6wFwzW+uc+yGPciUB7znnDptZD7L+1fO3IB8bzlzHdAYmO+cyApaFa76C4cX2FRQza0FWuQdeRPhq/1ydCcwys2/9e7V55Quyvg6/38zaAZ8C9cgH8+XXAVjsnAvcyw/rnJlZabJ+mfRzzu3LvjqHh4Rl+8pve+4nczFu7K8X4w7mseHMhZm1AgYCHZ1zh48td85t8/93EzCfrN/oeZLLOfdrQJZxwKXBPjacuQJ0Jts/mcM4X8E4XvZwzleuzOwiYDzQyTn367HlAXO1E/iE0B2KDIpzbp9zbr//djJQxMwq4fF8BTjR9hXyOTOzImQV+7vOuY9zGJJ321eo31Q4zTckCpP1RkJt/u9NmIbZxvTmr2+o/td/uyF/fUN1E6F7QzWYXI3JegOpXrbl5YFi/tuVgO8J0RtLQeaqGnD7RmCZ+783cDb785X3366QV7n8484j680ty4v5CniNWhz/DcIb+OsbXivCPV9BZIol6z2kq7ItLwWUCbi9BGgTyrkKIluVY39+ZJXkFv/cBbUNhCuXf/2xHb9SeTFn/v/vt4ERJxiTZ9tXSDeCEE1QO7LeZf4BGOhfNoSsvWGA4sCH/o19BVAn4LED/Y/bALTN41yzgV+AL/0/Pv/yq4C1/o17LRCfx7leAL72v/484PyAx3b1z+NG4L68zOW/PxgYlu1x4Z6v94DtQDpZe0vxQA+gh3+9AaP8udcCceGeryAyjQd2B2xbKf7ldfzz9JX/z3hgKOcqyGx9AravZQT8AsppG8irXP4x95L1IYvAx4Vtzsg6XOaANQF/Vu282r70DVURkQiU3465i4hICKjcRUQikMpdRCQCqdxFRCKQyl1EJAKp3EVEIpDKXUQkAqncRUQi0P8HOtcZFZieW34AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from matplotlib import pyplot as plt\n", "% matplotlib inline\n", "\n", "plt.plot(x, fx, 'g')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Правдоподобие выборок" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Правдоподобие выборки $x_1, x_2, \\ldots, x_n$, взятой из дискретного распределения, определяется так:\n", "\n", "$$\n", "L(x_1, x_2, \\ldots, x_n) = \\prod\\limits_{i=1}^{n}p_i\n", "$$\n", "\n", "Правдоподобие выборки $x_1, x_2, \\ldots, x_n$, взятой из непрерывного распределения, определяется так:\n", "\n", "$$\n", "L(x_1, x_2, \\ldots, x_n) = \\prod\\limits_{i=1}^{n}f(x_i), \\text{где $f$ – функция распределения.}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Когда Python генерирует псевдослучайные выборки из распределений, в результате получаются выборки с высоким правдоподобием. Убедимся в этом. Для начала напишем функцию, которая будет вычислять правдоподобие выборки из какого-нибудь распределения. Для примера возьмём биномиальное распределение." ] }, { "cell_type": "code", "execution_count": 74, "metadata": {}, "outputs": [], "source": [ "def binom_likelihood(n, p, sample):\n", " X = st.binom(n, p)\n", " L = 1\n", " for x in sample:\n", " P = X.pmf(x)\n", " L *= P\n", " return L" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Функция выше перемножает все вероятности значений в поданной на вход выборке и возвращает значение правдоподобия. Из-за того, что перемножаться всё время будут значения меньше $1$, результаты будут очень малы и сравнивать значения правдоподобия будет неудобно. Поэтому давайте поправим функцию так, чтобы она заодно возвращала натуральный логарифм правдоподобия." ] }, { "cell_type": "code", "execution_count": 83, "metadata": {}, "outputs": [], "source": [ "def binom_likelihood(n, p, sample):\n", " X = st.binom(n, p)\n", " L = 1\n", " for x in sample:\n", " P = X.pmf(x)\n", " L *= P\n", " return L, np.log(L)" ] }, { "cell_type": "code", "execution_count": 84, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.00024873671122804943, -8.299115605551385)" ] }, "execution_count": 84, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sample1 = st.binom(10, 0.5).rvs(5)\n", "binom_likelihood(10, 0.5, sample1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Кажется, что значение правдоподобия очень мало. Но это относительно. Предложим выборку со значением, которое не может принадлежать случайной величине с биномиальным распределением с $n=10$, например, $12$:" ] }, { "cell_type": "code", "execution_count": 85, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:7: RuntimeWarning: divide by zero encountered in log\n", " import sys\n" ] }, { "data": { "text/plain": [ "(0.0, -inf)" ] }, "execution_count": 85, "metadata": {}, "output_type": "execute_result" } ], "source": [ "binom_likelihood(10, 0.5, [3, 6, 12, 4, 0]) # L=0 - выборка неправдоподобна совсем" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Или выборки с нетипичными (маловероятными) значениями:" ] }, { "cell_type": "code", "execution_count": 86, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(2.1582735598713085e-09, -19.953957212680532)" ] }, "execution_count": 86, "metadata": {}, "output_type": "execute_result" } ], "source": [ "binom_likelihood(10, 0.5, [3, 2, 10, 9, 8])" ] }, { "cell_type": "code", "execution_count": 88, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(3.996802888650572e-12, -26.24552635223885)" ] }, "execution_count": 88, "metadata": {}, "output_type": "execute_result" } ], "source": [ "binom_likelihood(10, 0.5, [0, 1, 10, 1, 2])" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }