diff --git a/Assignment_3/200074_Akash.ipynb b/Assignment_3/200074_Akash.ipynb new file mode 100644 index 0000000..34453bd --- /dev/null +++ b/Assignment_3/200074_Akash.ipynb @@ -0,0 +1,238 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "dfb7cd7a", + "metadata": {}, + "source": [ + "# 1." + ] + }, + { + "cell_type": "markdown", + "id": "20db50bb", + "metadata": {}, + "source": [ + "For standard Cauchy distribution, we have the pdf as $\\pi(x) = \\frac{1}{\\pi(1+x^2)}$ \n", + "As, we use standard normal distibution as the importance distribution, so we have its pdf as $ g(x) = \\frac{1}{\\sqrt{ 2 \\pi}}e^{-\\frac{x^2}{2}}$ \n", + "A, mean is first moment so, $h(x) = x$ \n", + "We can estimate the mean as $\\hat{\\theta_g} = \\frac{1}{N} \\sum_{t=1}^{N}\\frac{h(Z_t)\\pi(Z_t)}{g(Z_t)}$ , where $Z_{i}$s are obtained from importance sampling estimator of the proposal distribution. \n", + "Putting respective values, we get mean as $\\hat{\\theta_g} = \\frac{1}{N} \\sum_{t=1}^{N} \\frac{1}{N} \\sum_{t=1}^{N} \\frac{\\sqrt{2} Z_t e^{-\\frac{Z_t^2}{2}}}{\\sqrt{\\pi}(Z_t^2 + 1)}$ \n", + "$\\sigma_g(\\hat{\\theta_g}) = \\sigma_g^2 + \\theta^2 = \\int_{-\\infty}^{\\infty} \\frac{h^2(x)\\pi^2(x)}{g^2(x)} dx$ \n", + "$ \\sigma_g^2 = \\int_{-\\infty}^{\\infty} \\frac{h^2(x)\\pi^2(x)}{g^2(x)} dx = \\int_{-\\infty}^{\\infty} \\frac{2x^2 e^{x^2}}{\\pi (1+x^2)^2} dx \\rightarrow \\infty \\implies$ Estimator does not has finite variance." + ] + }, + { + "cell_type": "markdown", + "id": "3d04a737", + "metadata": {}, + "source": [ + "# 2." + ] + }, + { + "cell_type": "markdown", + "id": "6ec1f356", + "metadata": {}, + "source": [ + "(a). We are given that sup $\\frac{f(x)}{g(x)} < \\infty$ \n", + "As, $\\tilde{f(x)} \\propto f(x)$ and $\\tilde{g(x)} \\propto g(x)$, so $\\frac{\\tilde{f(x)}}{\\tilde{g(x)}} \\propto \\frac{f(x)}{g(x)} < \\infty$ \n", + "Hence, the weighted importance estimator should have finite variance." + ] + }, + { + "cell_type": "markdown", + "id": "1c42440d", + "metadata": {}, + "source": [ + "(b). In accept-reject, a proposal might get rejected even after several iterations but in importance sampling, we get output after every iteration." + ] + }, + { + "cell_type": "markdown", + "id": "82e98bdb", + "metadata": {}, + "source": [ + "# 3." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "59555451", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "MersenneTwister(1)" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using Distributions\n", + "using Random\n", + "\n", + "Random.seed!(1)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "a40a35d1", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "genrte (generic function with 2 methods)" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function genrte(v, n=50)\n", + " tvec = TDist(v)\n", + " y = rand(tvec, n)\n", + " product = 1\n", + " \n", + " numr = 0\n", + " denr = 0\n", + " reslt = zeros(10000)\n", + " for i = 1:10000\n", + " x = rand(Normal())\n", + " for i = 1:n\n", + " product *= (1 + ((y[i] - x)^2)/v)^(-(v + 1)/2)\n", + " end\n", + " r = product * (sqrt(2 * π))\n", + " numr += x * r\n", + " denr += r\n", + " reslt[i] = numr / denr\n", + " end\n", + " return reslt\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "6937564f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Mean is -1.5477036982569745\n", + "Variance is 2.416128135052854e-30\n" + ] + } + ], + "source": [ + "generate5 = genrte(5)\n", + "print(\"Mean is \")\n", + "println(mean(generate5))\n", + "print(\"Variance is \")\n", + "println(var(generate5))" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "4c9fa654", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Mean is -0.8377053458293047\n", + "Variance is 4.437786370505242e-31\n" + ] + } + ], + "source": [ + "generate1 = genrte(1)\n", + "print(\"Mean is \")\n", + "println(mean(generate1))\n", + "print(\"Variance is \")\n", + "println(var(generate1))" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "36828f7e", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Mean is -0.801275682761385\n", + "Variance is 1.232718436251456e-30\n" + ] + } + ], + "source": [ + "generate2 = genrte(2)\n", + "print(\"Mean is \")\n", + "println(mean(generate2))\n", + "print(\"Variance is \")\n", + "println(var(generate2))" + ] + }, + { + "cell_type": "markdown", + "id": "a14810da", + "metadata": {}, + "source": [ + "# 4." + ] + }, + { + "cell_type": "markdown", + "id": "4a220a57", + "metadata": {}, + "source": [ + "$$\n", + "p (\\lambda | Y) \\propto \\lambda ^{\\alpha -1}e^{-\\beta\\lambda}\\prod_{i=1}^{n}\\frac{\\lambda^{Y_{i}}e^{-\\lambda}}{Y_{i}!}\\\\\n", + "p (\\lambda | Y) \\propto \\lambda ^{\\alpha +\\sum_{i=1}^{n}Y_{i}-1}e^{-\\beta\\lambda}\\prod_{i=1}^{n}\\frac{e^{-\\lambda}}{Y_{i}!}\\\\\n", + "p (\\lambda | Y) \\propto \\lambda ^{\\alpha +\\sum_{i=1}^{n}Y_{i}-1}e^{-\\beta\\lambda-n\\lambda}\\prod_{i=1}^{n}\\frac{1}{Y_{i}!}\\\\\n", + "p (\\lambda | Y) \\propto \\lambda ^{\\alpha +\\sum_{i=1}^{n}Y_{i}-1}e^{-\\beta\\lambda-n\\lambda}\\\\\n", + "$$\n", + "$p (\\lambda | Y)$ represents a Gamma Distribution ~ $Gamma(\\alpha +\\sum_{i=1}^{n}Y_{i},\\beta+n)$ " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc2ccf51", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.6.2", + "language": "julia", + "name": "julia-1.6" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.6.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}