{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Beta分布的矩估计\n", "如果$x_i \\sim Beta(\\alpha, \\beta)$，由于$E(x_i)=\\frac{\\alpha}{\\alpha + \\beta}$，而$E(x_i^2)=\\frac{\\alpha ^2}{(\\alpha + \\beta)^2}+\\frac{\\alpha \\beta}{(\\alpha + \\beta)^2 (\\alpha+\\beta+1)}$，从而我们的矩估计即联立：$$\\frac{\\hat{\\alpha}}{\\hat{\\alpha} + \\hat{\\beta}}=\\bar{x}$$ $$\\frac{\\hat{\\alpha} ^2}{(\\hat{\\alpha} + \\hat{\\beta})^2}+\\frac{\\hat{\\alpha} \\hat{\\beta}}{(\\hat{\\alpha} + \\hat{\\beta})^2 (\\hat{\\alpha}+\\hat{\\beta}+1)}=\\overline{x^{2}}$$ 即可得到矩估计。在这里，我们将联立方程问题转化为一个最优化问题，即最小化： $$\\min_{\\hat{\\alpha},\\hat{\\beta}}\\left[\\frac{\\hat{\\alpha}}{\\hat{\\alpha} + \\hat{\\beta}}-\\bar{x} \\right]^2+\\left[\\frac{\\hat{\\alpha} ^2}{(\\hat{\\alpha} + \\hat{\\beta})^2}+\\frac{\\hat{\\alpha} \\hat{\\beta}}{(\\hat{\\alpha} + \\hat{\\beta})^2 (\\hat{\\alpha}+\\hat{\\beta}+1)}-\\overline{x^{2}}\\right]^2$$\n", "我们将会重复抽样、估计M=500次，并根据这500次的结果计算矩估计量的偏差（bias）、标准误(standard error)以及均方误差（mean sqrared error）。" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bias = [ 0.0392831 0.01443088]\n", "s.e. = [ 0.35045361 0.10847059]\n", "RMSE = [ 0.3526484 0.10942632]\n" ] } ], "source": [ "import numpy as np\n", "from numpy import random as nprd\n", "from scipy.optimize import minimize\n", "import scipy as sc\n", "\n", "def sampling(a,b,N):\n", " x=nprd.beta(a,b,N)\n", " return x\n", "\n", "def estimate(x):\n", " meanx=np.mean(x)\n", " x2=[xi**2 for xi in x]\n", " meanx2=np.mean(x2)\n", " def obj(theta):\n", " return (theta[0]/(theta[0]+theta[1])-meanx)**2 + ((theta[0]/(theta[0]+theta[1]))**2+(theta[0]*theta[1])/((theta[0]+theta[1])**2*(theta[0]+theta[1]+1))-meanx2)**2\n", " res=minimize(obj, np.array([1,1]), method='nelder-mead', options={'xtol': 1e-4, 'disp': False})\n", " return res\n", "\n", "M=500 ## simulation times\n", "N=200 ## sample size\n", "a=3\n", "b=1 ## true value\n", "RESULT=np.zeros((M,2), np.float64)\n", "for m in range(M):\n", " x=sampling(a,b,N)\n", " res=estimate(x)\n", " RESULT[m]=res.x\n", "\n", "MEAN_RESULT=np.average(RESULT, 0)\n", "BIAS=MEAN_RESULT-np.array([a,b])\n", "STD=np.std(RESULT, 0)\n", "MSE2=np.array([i**2 for i in STD])+np.array([i**2 for i in BIAS])\n", "MSE=np.array([np.sqrt(i) for i in MSE2])\n", "print(\"Bias = \", BIAS)\n", "print(\"s.e. = \", STD)\n", "print(\"RMSE = \", MSE)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Beta分布的极大似然估计\n", "由于Beta分布的对数似然函数为$$\\ln \\left( \\alpha, \\beta | x \\right)=\\sum_{i=1}^N \\left[ -\\ln (Beta(\\alpha,\\beta))+(\\alpha-1) \\ln (x_i) + (\\beta-1)\\ln (1-x_i) \\right]$$\n", "最大化似然函数，或者最小化负的似然函数，即可得到极大似然估计。" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bias = [ 0.0677773 0.0104206]\n", "s.e. = [ 0.31265354 0.0886861 ]\n", "RMSE = [ 0.31991561 0.08929621]\n" ] } ], "source": [ "import numpy as np\n", "from numpy import random as nprd\n", "from scipy.optimize import minimize\n", "import scipy as sc\n", "\n", "def sampling(a,b,N):\n", " x=nprd.beta(a,b,N)\n", " return x\n", " \n", "def estimate(x):\n", " def log_likelihood(theta):\n", " likeli=np.array([-1*np.log(sc.special.beta(theta[0],theta[1]))+(theta[0]-1)*np.log(xi)+(theta[1]-1)*np.log(1-xi) for xi in x])\n", " return -1*np.mean(likeli)\n", " res=minimize(log_likelihood, np.array([1,1]), method='nelder-mead', options={'xtol': 1e-4, 'disp': False})\n", " return res\n", "\n", "M=500 ## simulation times\n", "N=200 ## sample size\n", "a=3\n", "b=1 ## true value\n", "RESULT=np.zeros((M,2), np.float64)\n", "for m in range(M):\n", " x=sampling(a,b,N)\n", " res=estimate(x)\n", " RESULT[m]=res.x\n", "\n", "MEAN_RESULT=np.average(RESULT, 0)\n", "BIAS=MEAN_RESULT-np.array([a,b])\n", "STD=np.std(RESULT, 0)\n", "MSE2=np.array([i**2 for i in STD])+np.array([i**2 for i in BIAS])\n", "MSE=np.array([np.sqrt(i) for i in MSE2])\n", "print(\"Bias = \", BIAS)\n", "print(\"s.e. = \", STD)\n", "print(\"RMSE = \", MSE)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "通过比较，发现尽管Bias有大有效，但是极大似然估计的确比矩估计更有效。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 区间估计\n", "## 小样本正态总体\n", "如果$x_i\\sim N\\left( \\mu, \\sigma^2 \\right)$，那么$\\bar{x}\\sim N\\left(\\mu,\\frac{\\sigma^2}{N}\\right)$，或者，如果我们用样本方差替代总体方差，那么有：$$\\frac{\\bar{x}-\\mu}{\\sqrt{\\frac{s^2}{N}}}\\sim t(N-1)$$\n", "因而，对于总体均值$\\mu$的95%置信水平的区间估计为$$\\left[\\bar{x}-t_{1-\\alpha}(N-1)\\sqrt{\\frac{s^2}{N}},\\bar{x}+t_{1-\\alpha}(N-1)\\sqrt{\\frac{s^2}{N}}\\right]$$\n", "下面我们将重复1000次区间估计，理论上，在95%置信水平下，应该有950次区间估计包含真值，我们将看到这一结论是否成立。" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The prob. of included= 0.943\n" ] } ], "source": [ "import numpy as np\n", "from numpy import random as nprd\n", "from scipy import special as func\n", "\n", "def sampling(mu, sigma2, N):\n", " x=nprd.normal(mu,np.sqrt(sigma2),N)\n", " return x\n", "\n", "## true value\n", "mu=10\n", "sigma2=30\n", "N=10 # sample size\n", "## iteration times\n", "M=1000\n", "## confidence level\n", "alpha=0.05\n", "## results\n", "included=np.zeros(M)\n", "for i in range(M):\n", " x=sampling(mu,sigma2,10)\n", " xmean=np.mean(x)\n", " xstd=np.std(x)\n", " lower=xmean-func.stdtrit(N-1,1-alpha/2)*xstd/np.sqrt(N)\n", " upper=xmean+func.stdtrit(N-1,1-alpha/2)*xstd/np.sqrt(N)\n", " # 如果包含真值\n", " if mu>=lower and mu<=upper:\n", " included[i]=1\n", "print(\"The prob. of included=\",np.mean(included))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 区间估计\n", "## 大样本总体未知\n", "如果$x_i\\sim \\left( \\mu, \\sigma^2 \\right)$，那么当样本量足够大时，$\\bar{x}\\sim N\\left(\\mu,\\frac{\\sigma^2}{N}\\right)$，或者，如果我们用样本方差替代总体方差，那么当样本量足够大时，有：$$\\frac{\\bar{x}-\\mu}{\\sqrt{\\frac{s^2}{N}}}\\sim N\\left(0,1\\right)$$\n", "下面我们将重复1000次区间估计，理论上，在95%置信水平下，应该有950次区间估计包含真值，我们将看到这一结论是否成立。" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Sample size=10, the prob. of included= 0.909\n", "Sample size=50, the prob. of included= 0.951\n", "Sample size=100, the prob. of included= 0.95\n", "Sample size=500, the prob. of included= 0.942\n", "Sample size=1000, the prob. of included= 0.95\n" ] } ], "source": [ "import numpy as np\n", "from numpy import random as nprd\n", "from scipy import special as func\n", "\n", "def sampling(N):\n", " ## 产生一组样本，以0.5的概率为z+3，0.5的概率为z-3，其中z~N(0,1)，因而期望为0，方差为10\n", " d=nprd.rand(N)<0.5\n", " z=nprd.randn(N)\n", " x=np.array([z[i]+3 if d[i] else z[i]-3 for i in range(N)])\n", " return x\n", "\n", "## true value\n", "N=[10,50,100,500,1000] # sample size\n", "## iteration times\n", "M=1000\n", "## confidence level\n", "alpha=0.05\n", "for n in N:\n", " ## results\n", " included=np.zeros(M)\n", " for i in range(M):\n", " x=sampling(n)\n", " xmean=np.mean(x)\n", " xstd=np.std(x)\n", " lower=xmean-func.ndtri(1-alpha/2)*xstd/np.sqrt(n)\n", " upper=xmean+func.ndtri(1-alpha/2)*xstd/np.sqrt(n)\n", " # 如果包含真值\n", " #print(mu, lower, upper, mu>=lower, mu<=upper)\n", " if 0>=lower and 0<=upper:\n", " included[i]=1\n", " print(\"Sample size=%s, the prob. of included=\" % n, np.mean(included))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "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.7" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false } }, "nbformat": 4, "nbformat_minor": 2 }