Le problème du vendeur de journaux (à une période de temps)
\n",
"$\\newcommand\\indi[1]{{\\mathbf 1}_{\\displaystyle #1}}$\n",
"$\\newcommand\\inde[1]{{\\mathbf 1}_{\\displaystyle\\left\\{ #1 \\right\\}}}$\n",
"$\\newcommand{\\ind}{\\inde}$\n",
"$\\newcommand\\E{{\\mathbf E}}$\n",
"$\\newcommand\\Cov{{\\mathrm Cov}}$\n",
"$\\newcommand\\Var{{\\mathrm Var}}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Chaque matin, le vendeur doit décider d’un nombre de journaux à commander $u \\in \\mathbb{U}= \\{0, 1,\\ldots\\}$ au prix unitaire $c > 0$. La demande du jour est incertaine $w \\in \\mathbb{W} = \\{0, 1,\\ldots\\}$.\n",
"\n",
"Si à la fin de la journée il lui reste des invendus (coût unitaire $c_S \\in \\mathbb{R}$):\n",
"$c_S(u − w )_+ = c_S \\max (u − w,0)$\n",
"\n",
"On suppose que $c > − c_S$.\n",
"\n",
"Si à la fin de la journée il n’a pas pu faire face à la demande on associe un coût unitaire $c_M$. Le coût lié à la non satisfaction de la demande est:\n",
"$c_M (w − u )_+ = c_M \\max (w − u,0)$.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 249,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np;\n",
"from scipy.special import comb;\n",
"from scipy.special import gamma;\n",
"import matplotlib.pyplot as plt\n",
"import io;"
]
},
{
"cell_type": "code",
"execution_count": 252,
"metadata": {},
"outputs": [],
"source": [
"binomiale=1\n",
"discrete=2\n",
"poisson=3\n",
"\n",
"def choisir_loi(law):\n",
"\n",
" if law == binomiale: ## loi binomiale\n",
" n=100\n",
" p=0.5\n",
"\n",
" buf = io.StringIO()\n",
" buf.write(\"loi binomiale(%d,%5.2f)\" % (n,p))\n",
" title = buf.getvalue()\n",
" \n",
" wi = np.linspace(0,n,num=n+1) ## les valeurs possibles {0,1,...,n}\n",
" ## pdf(\"bin\",x,n,p) = n!/x!(n-x)! p^x (1-p)^n-x \n",
" nn = n*np.ones(n+1)\n",
" pi=comb(nn,wi, exact=False) * pow(p,wi)*pow(1-p,n-wi)\n",
"\n",
" mu=n*p; ## moyenne de la binomiale \n",
" mu1= sum(pi*wi) ; ## vérification \n",
" if abs(mu-mu1) > 1.e-8:\n",
" print(\"something wrong in binomial law expectation\")\n",
" wait = input(\"PRESS ENTER TO CONTINUE.\")\n",
" \n",
" # un echantillong de taille N\n",
" N=1000\n",
" #W=grand(1,N,\"bin\",n,p);\n",
" W = np.random.binomial(n, p, N)\n",
" \n",
" if law == discrete: ## une loi discrète \n",
" n=3\n",
" wi=np.array([30,50,80])\n",
" # A vous\n",
" # pi = choisir des probabilités\n",
" buf = io.StringIO()\n",
" buf.write(\"loi discrète sur %d valeurs\" % wi.size)\n",
" title = buf.getvalue()\n",
" # A vous : mu= ## la moyenne\n",
"\n",
" N=1000\n",
" # un echantillong de taille N\n",
" # A vous : tirer un échantillon de taille N\n",
" \n",
" if law == poisson: ## loi de Poisson de paramètre mu \n",
" n=100;\n",
" p=0.5;\n",
" mu=n*p;\n",
" wi=np.linspace(0,n,num=n+1); ## les valeurs possibles \n",
" pi=( pow(mu,wi) *np.exp(-mu)) / gamma(wi+1);\n",
" buf = io.StringIO();\n",
" buf.write(\"Poisson %f\" % mu);\n",
" title = buf.getvalue();\n",
" # A vous : \n",
" # calculer la moyenne theorique, \n",
" # tirer un echantillon de taille N, \n",
" # comparer avec la moyenne empirique\n",
"\n",
" return pi, wi, W, title\n",
"\n",
"def hist_plot(samples,support=[]):\n",
" if np.size(support)==0:\n",
" support = np.sort(list(dict.fromkeys(samples)))\n",
" histo=np.zeros(support.size);\n",
" for i in range(support.size):\n",
" histo[i]=np.size(np.where(samples==support[i]))/(W.size)\n",
" # On trace cet histogramme\n",
" plt.bar(support, histo,width=1)\n"
]
},
{
"cell_type": "code",
"execution_count": 257,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAEZpJREFUeJzt3X+s3Xddx/Hny9Z1gmGDrhhpV1qygg4Mv+o2FQhhzHSgFOOWFYnsjyWVyAL+inYBlm3hD2aMA8OCLmw6qmGDgXgD1UUZiWigtnMTKGNyNwa7dLqOjuHQMQpv/zjfJofDvb3fe3tub+/5PB/JTb/fz/dz7nl/8715nU8/53s+J1WFJKkNP7bcBUiSThxDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+NIskDyR59Tx9Xp7k3hNVkzQOq5e7AGmlqqrPAM9b7jqkhXCkL0kNMfSlY0iyJsl7khzsft6TZE137JVJZpa7RmkhDH3p2N4OnAe8CHghcA7wjmWtSDoOhr50bG8Erqmqh6vqEHA18JvLXJO0aIa+dGzPAr42tP+1rk1akQx96dgOAs8e2t/YtUkrkqEvHduHgHckWZfkDOBK4K+XuSZp0bxPXzq2dwFPAz7f7X+ka5NWpPglKpLUDqd3JKkhhr4kNcTQl6SGGPqS1JCT7u6dM844ozZt2rTcZUjSinLnnXc+UlXr5ut30oX+pk2b2L9//3KXIUkrSpKvzd/L6R1JaoqhL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWrISfeJXOlksWnXJ2dtf+Ddrz3BlUjj40hfkhpi6EtSQwx9SWqIc/rSkLnm8efq4/y+VhpH+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcT79NW8Pvfm93ms9+xrJXCkL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQ3qFfpJtSe5NMp1k1yzH1yS5tTu+N8mmrn1Tkv9Lcnf38+fjLV+StBDzfjgrySrgeuACYAbYl2Sqqr401O0y4NGqOivJDuBa4JLu2H1V9aIx1y1JWoQ+I/1zgOmqur+qngRuAbaP9NkO3Nxt3wacnyTjK1OSNA59lmFYDzw4tD8DnDtXn6o6kuQxYG13bHOSu4BvA++oqs+MPkGSncBOgI0bNy7oBKSThUsyaCXoE/qzjdirZ5+HgI1V9c0kLwU+nuT5VfXtH+pYdQNwA8DWrVtHf7c0dsez3o60kvWZ3pkBzhza3wAcnKtPktXAacDhqvpuVX0ToKruBO4Dnnu8RUuSFqdP6O8DtiTZnOQUYAcwNdJnCri0274IuKOqKsm67o1gkjwH2ALcP57SJUkLNe/0TjdHfzlwO7AKuKmqDiS5BthfVVPAjcDuJNPAYQYvDACvAK5JcgT4PvDmqjq8FCciSZpfr/X0q2oPsGek7cqh7SeAi2d53EeBjx5njZKkMfETuZLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNaTXffqSFsbF13SycqQvSQ1xpK9muLKm5Ehfkppi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIa6yKS0x19bXycSRviQ1xNCXpIY4vaOJ5henSD/Mkb4kNaRX6CfZluTeJNNJds1yfE2SW7vje5NsGjm+McnjSf5gPGVLkhZj3tBPsgq4HrgQOBt4Q5KzR7pdBjxaVWcB1wHXjhy/Dvj74y9XknQ8+oz0zwGmq+r+qnoSuAXYPtJnO3Bzt30bcH6SACR5PXA/cGA8JUuSFqtP6K8HHhzan+naZu1TVUeAx4C1SZ4K/BFw9bGeIMnOJPuT7D906FDf2iVJC9Qn9DNLW/XsczVwXVU9fqwnqKobqmprVW1dt25dj5IkSYvR55bNGeDMof0NwME5+swkWQ2cBhwGzgUuSvLHwOnAD5I8UVXvO+7KJUkL1if09wFbkmwGvgHsAH5jpM8UcCnwWeAi4I6qKuDlRzskuQp43MCXpOUzb+hX1ZEklwO3A6uAm6rqQJJrgP1VNQXcCOxOMs1ghL9jKYuWJC1Or0/kVtUeYM9I25VD208AF8/zO65aRH2SpDHyE7mS1BBDX5Ia4oJr0gnk2vpabo70Jakhhr4kNcTQl6SGGPqS1BDfyNXE8duypLk50pekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhri2jvSMvELVbQcHOlLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0JekhnjLpiaCX5Eo9dNrpJ9kW5J7k0wn2TXL8TVJbu2O702yqWs/J8nd3c9/JPm18ZYvSVqIeUM/ySrgeuBC4GzgDUnOHul2GfBoVZ0FXAdc27V/EdhaVS8CtgF/kcT/XUjSMukz0j8HmK6q+6vqSeAWYPtIn+3Azd32bcD5SVJV/1tVR7r2U4EaR9GSpMXpE/rrgQeH9me6tln7dCH/GLAWIMm5SQ4AXwDePPQiIEk6wfqEfmZpGx2xz9mnqvZW1fOBnweuSHLqjzxBsjPJ/iT7Dx061KMkSdJi9An9GeDMof0NwMG5+nRz9qcBh4c7VNU9wHeAF4w+QVXdUFVbq2rrunXr+lcvSVqQPqG/D9iSZHOSU4AdwNRInyng0m77IuCOqqruMasBkjwbeB7wwFgqlyQt2Lx30lTVkSSXA7cDq4CbqupAkmuA/VU1BdwI7E4yzWCEv6N7+MuAXUm+B/wA+O2qemQpTkSSNL9et09W1R5gz0jblUPbTwAXz/K43cDu46xRkjQm3jMvnQT8QhWdKK69I0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQ7xPXyuW35YlLZwjfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGuJ9+tJJxrX1tZQc6UtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkP8RK5WFL8tSzo+jvQlqSG9Qj/JtiT3JplOsmuW42uS3Nod35tkU9d+QZI7k3yh+/dV4y1fkrQQ84Z+klXA9cCFwNnAG5KcPdLtMuDRqjoLuA64tmt/BPjVqvo54FJg97gKlyQtXJ+R/jnAdFXdX1VPArcA20f6bAdu7rZvA85Pkqq6q6oOdu0HgFOTrBlH4ZKkhevzRu564MGh/Rng3Ln6VNWRJI8BaxmM9I/6deCuqvru6BMk2QnsBNi4cWPv4qVJ5zLLGrc+I/3M0lYL6ZPk+QymfH5rtieoqhuqamtVbV23bl2PkiRJi9En9GeAM4f2NwAH5+qTZDVwGnC4298A/C3wpqq673gLliQtXp/Q3wdsSbI5ySnADmBqpM8UgzdqAS4C7qiqSnI68Engiqr613EVLUlanHlDv6qOAJcDtwP3AB+uqgNJrknyuq7bjcDaJNPA7wFHb+u8HDgLeGeSu7ufZ479LCRJvfT6RG5V7QH2jLRdObT9BHDxLI97F/Cu46xRkjQmfiJXkhpi6EtSQwx9SWqIoS9JDTH0Jakhrqevk55r6Evj40hfkhriSF9aIVx8TePgSF+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQ1yGQSclF1mTloahL61ArsOjxXJ6R5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktSQXqGfZFuSe5NMJ9k1y/E1SW7tju9NsqlrX5vk00keT/K+8ZYuSVqoeT+Rm2QVcD1wATAD7EsyVVVfGup2GfBoVZ2VZAdwLXAJ8ATwTuAF3Y+kMfPTuVqIPsswnANMV9X9AEluAbYDw6G/Hbiq274NeF+SVNV3gH9Jctb4Stakcr0daen1md5ZDzw4tD/Ttc3ap6qOAI8Ba/sWkWRnkv1J9h86dKjvwyRJC9Qn9DNLWy2iz5yq6oaq2lpVW9etW9f3YZKkBeoT+jPAmUP7G4CDc/VJsho4DTg8jgIlSePTJ/T3AVuSbE5yCrADmBrpMwVc2m1fBNxRVb1H+pKkE2PeN3Kr6kiSy4HbgVXATVV1IMk1wP6qmgJuBHYnmWYwwt9x9PFJHgCeBpyS5PXAL4/c+SNJOkF6fYlKVe0B9oy0XTm0/QRw8RyP3XQc9UmSxshP5EpSQ/y6RGmC+EEtzceRviQ1xJG+lpWfwpVOLEf6ktQQQ1+SGmLoS1JDnNOXJpR38mg2jvQlqSGGviQ1xNCXpIY4p68TznvzpeXjSF+SGuJIX2qAd/LoKEf6ktQQQ1+SGmLoS1JDnNPXCeEdO9LJwdCXGuObum1zekeSGmLoS1JDnN7RknEe/+TnVE97HOlLUkMMfUlqiNM7kgCnelph6GusnMeXTm6GvqQfMfri7ch/chj6Om6O7iefUz+To1foJ9kGvBdYBXygqt49cnwN8EHgpcA3gUuq6oHu2BXAZcD3gbdW1e1jq17LxqBvly8AK9u8oZ9kFXA9cAEwA+xLMlVVXxrqdhnwaFWdlWQHcC1wSZKzgR3A84FnAf+U5LlV9f1xn4iWnkGvUb4ArDx9RvrnANNVdT9AkluA7cBw6G8Hruq2bwPelyRd+y1V9V3gq0mmu9/32fGUr6VguGsxFvp344vE8ugT+uuBB4f2Z4Bz5+pTVUeSPAas7do/N/LY9aNPkGQnsLPbfTzJvb2qn9sZwCPH+TtWktbOF9o754k731x7zMMTd77zGMf5PrtPpz6hn1naqmefPo+lqm4AbuhRSy9J9lfV1nH9vpNda+cL7Z2z5zvZTuT59vlE7gxw5tD+BuDgXH2SrAZOAw73fKwk6QTpE/r7gC1JNic5hcEbs1MjfaaAS7vti4A7qqq69h1J1iTZDGwB/m08pUuSFmre6Z1ujv5y4HYGt2zeVFUHklwD7K+qKeBGYHf3Ru1hBi8MdP0+zOBN3yPAW07QnTtjmypaIVo7X2jvnD3fyXbCzjeDAbkkqQWusilJDTH0JakhExf6SbYluTfJdJJdy13PuCU5M8mnk9yT5ECSt3Xtz0jyj0m+0v379OWudZySrEpyV5JPdPubk+ztzvfW7iaDiZDk9CS3Jflyd51/oYHr+7vd3/MXk3woyamTdI2T3JTk4SRfHGqb9Zpm4M+6DPt8kpeMs5aJCv2hJSMuBM4G3tAtBTFJjgC/X1U/C5wHvKU7x13Ap6pqC/Cpbn+SvA24Z2j/WuC67nwfZbAUyKR4L/APVfUzwAsZnPfEXt8k64G3Alur6gUMbhg5upzLpFzjvwK2jbTNdU0vZHCn4xYGH1p9/zgLmajQZ2jJiKp6Eji6ZMTEqKqHqurfu+3/YRAI6xmc581dt5uB1y9PheOXZAPwWuAD3X6AVzFY8gMm6HyTPA14BYM74qiqJ6vqW0zw9e2sBn6i+5zPU4CHmKBrXFX/zODOxmFzXdPtwAdr4HPA6Ul+ely1TFroz7ZkxI8s+zApkmwCXgzsBX6qqh6CwQsD8Mzlq2zs3gP8IfCDbn8t8K2qOtLtT9J1fg5wCPjLbjrrA0meygRf36r6BvAnwNcZhP1jwJ1M7jU+aq5ruqQ5Nmmh32vZh0mQ5CeBjwK/U1XfXu56lkqSXwEerqo7h5tn6Top13k18BLg/VX1YuA7TNBUzmy6ueztwGYGq/E+lcEUx6hJucbzWdK/70kL/SaWfUjy4wwC/2+q6mNd838f/S9g9+/Dy1XfmP0S8LokDzCYrnsVg5H/6d1UAEzWdZ4BZqpqb7d/G4MXgUm9vgCvBr5aVYeq6nvAx4BfZHKv8VFzXdMlzbFJC/0+S0asaN189o3APVX1p0OHhpfCuBT4uxNd21KoqiuqakNVbWJwPe+oqjcCn2aw5AdM1vn+F/Bgkud1Tecz+ET7RF7fzteB85I8pfv7PnrOE3mNh8x1TaeAN3V38ZwHPHZ0GmgsqmqifoDXAP8J3Ae8fbnrWYLzexmD/+p9Hri7+3kNg3nuTwFf6f59xnLXugTn/krgE932cxis4zQNfARYs9z1jfE8XwTs767xx4GnT/r1Ba4Gvgx8EdgNrJmkawx8iMH7Fd9jMJK/bK5rymB65/ouw77A4K6msdXiMgyS1JBJm96RJB2DoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5Ia8v/UpLwF02hEdwAAAABJRU5ErkJggg==\n",
"text/plain": [
"