Assgn3 Submitted 200204_Aryan

Assignment_2/.vscode/settings.json

@@ -0,0 +1,3 @@
+{
+    "jupyter.jupyterServerType": "local"
+}

Assignment_3/.ipynb_checkpoints/200204_Aryan-checkpoint.ipynb

@@ -0,0 +1,293 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "430b2fea",
+   "metadata": {},
+   "source": [
+    "# Question 1\n",
+    "For Standard Cauchy Distributions: $\\pi(x)=\\frac{1}{\\pi(x^{2}+1)}$\\\n",
+    "For Standard Normal Distributions: $g(x)=e^{\\frac{-x^{2}}{2}}\\frac{1}{\\sqrt{2\\pi}}$\n",
+    "$$\n",
+    "\\frac{\\pi(x)}{g(x)}=\\frac{\\sqrt{2\\pi}e^{\\frac{-x^{2}}{2}}}{\\pi(x^{2}+1)}\\\n",
+    "\\frac{\\pi(x)}{g(x)}=\\frac{\\sqrt{2}e^{\\frac{-x^{2}}{2}}}{\\sqrt{\\pi}(x^{2}+1)}\\\n",
+    "$$\n",
+    "$$sup\\frac{\\pi(x)}{g(x)}=\\sqrt{\\frac{2}{\\pi}}$$\n",
+    "Also, since we are estimating the mean, then $h(x)=1$ implying that $Var(h(x))=0$.\\\n",
+    "Since this is finite, using Theorem 2 of W4L11, we obtain that the variance of the estimator will be finite."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "90c64e88",
+   "metadata": {},
+   "source": [
+    "# Question 2"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "46e56fce",
+   "metadata": {},
+   "source": [
+    "## Question 2.1\n",
+    "Yes, the weighted importance estimator should also have a finite variance.The reasoning is as follows:\\\n",
+    "Since $sup\\frac{f(x)}{g(x)}<\\infty$, we can say that $sup\\frac{f(x)}{\\tilde{g(x)}}<\\infty$, as $\\tilde{g(x)} \\propto g(x)$ and the proportionality constant is known to be finite."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d873c79d",
+   "metadata": {},
+   "source": [
+    "## Question 2.2\n",
+    "I feel that the benefit of Importance Sampling means that in each iteration, we get an acceptable sample, however for Accept-Reject, this is not true. We generally expect to get an acceptable sample in more than one iteration.\\\n",
+    "This implies that Importance Sampling is faster to perform."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "74c3dd26",
+   "metadata": {},
+   "source": [
+    "# Question 3\n",
+    "\n",
+    "Let us define the following:\n",
+    "$$ \n",
+    "\\tilde{\\pi}(x)=e^{\\frac{-x^{2}}{2}}\\prod_{i=1}^{n}(1+\\frac{(y_{i}-x)^{2}}{v})^{-(v+1)/2}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "bbd08524",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "normalcdf (generic function with 1 method)"
+      ]
+     },
+     "execution_count": 1,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "using Distributions\n",
+    "using Random\n",
+    "Random.seed!(1)\n",
+    "function sampleY(v,n)\n",
+    "    global i=0\n",
+    "    global Y=[]\n",
+    "    global v=5\n",
+    "    d=TDist(v)\n",
+    "    while(i<n)\n",
+    "        global i=i+1\n",
+    "        y=rand(d,1)\n",
+    "        append!(Y,y)\n",
+    "    end\n",
+    "return Y\n",
+    "end\n",
+    "\n",
+    "function PI(Y,x,v,n)\n",
+    "    global m=exp(-x*x/2)\n",
+    "    global i1=0\n",
+    "    while(i1<n)\n",
+    "        global i1=i1+1\n",
+    "        global m=m*(1+(Y[i1]-x)^2/v)^(-(v+1)/2)\n",
+    "    end\n",
+    "    return m\n",
+    "end\n",
+    "\n",
+    "function normalcdf(x)\n",
+    "    a=exp((-x*x/2))/((2*pi)^(0.5))\n",
+    "    return a\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "0ee3e60a",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "secondMoment (generic function with 1 method)"
+      ]
+     },
+     "execution_count": 2,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "function firstMoment(n,v,num)\n",
+    "    Y=sampleY(v,n)\n",
+    "    global i3=0\n",
+    "    global top=0\n",
+    "    global bot=0\n",
+    "    global flag=0\n",
+    "    while (flag==0)\n",
+    "#         println(i3)\n",
+    "        x=randn()\n",
+    "        g=normalcdf(x)\n",
+    "        p=PI(Y,x,v,n)\n",
+    "#         println(\"p:\",p,\"g:\",g)\n",
+    "        top=top+(x*p)/g\n",
+    "        bot=bot+p/g\n",
+    "        i3=i3+1\n",
+    "        if(i3>=num)\n",
+    "            global flag=1\n",
+    "            break\n",
+    "        end\n",
+    "    end\n",
+    "    return top/bot\n",
+    "end\n",
+    "\n",
+    "function secondMoment(n,v,num)\n",
+    "    Y=sampleY(v,n)\n",
+    "    global i2=0\n",
+    "    global top2=0\n",
+    "    global bot2=0\n",
+    "    while (i2<num)\n",
+    "        i2=i2+1\n",
+    "        x=randn()\n",
+    "        g=normalcdf(x)\n",
+    "        p=PI(Y,x,v,n)\n",
+    "        top2=top2+x*x*p/g\n",
+    "        bot2=bot2+p/g\n",
+    "    end\n",
+    "    return top2/bot2\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "7f122d36",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "For v=5\n",
+      "Expectation:0.3195498376470879\n",
+      "Variance:0.048214074599329894\n"
+     ]
+    }
+   ],
+   "source": [
+    "println(\"For v=5\")\n",
+    "expectation=firstMoment(50,5,10000)\n",
+    "println(\"Expectation:\",expectation)\n",
+    "println(\"Variance:\",expectation^2-secondMoment(50,5,10000))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "14d4471f",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "For v=1\n",
+      "Expectation:0.04382631588330899\n",
+      "Variance:-0.18224589091390006\n"
+     ]
+    }
+   ],
+   "source": [
+    "println(\"For v=1\")\n",
+    "expectation=firstMoment(50,1,10000)\n",
+    "println(\"Expectation:\",expectation)\n",
+    "println(\"Variance:\",expectation^2-secondMoment(50,1,10000))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "c08f649e",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "For v=2\n",
+      "Expectation:0.20644879379356426\n",
+      "Variance:0.00898394357852978\n"
+     ]
+    }
+   ],
+   "source": [
+    "println(\"For v=2\")\n",
+    "expectation=firstMoment(50,2,10000)\n",
+    "println(\"Expectation:\",expectation)\n",
+    "println(\"Variance:\",(expectation)^2-secondMoment(50,2,10000))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "142ee15c",
+   "metadata": {},
+   "source": [
+    "We can observe that variance is finite"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fa9535c2",
+   "metadata": {},
+   "source": [
+    "# Question 4"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "520ce363",
+   "metadata": {},
+   "source": [
+    "Knowing the following:\n",
+    "$$\n",
+    "\\pi (\\lambda | D) \\propto \\pi(D|\\lambda)\\pi(\\lambda)\\\\\n",
+    "\\pi(\\lambda)=\\frac{\\lambda ^{\\alpha -1}e^{-\\beta\\lambda}\\beta^{\\alpha}}{\\Gamma{\\alpha}}\\\\\n",
+    "\\pi(D|\\lambda)=\\prod_{i=1}^{n}\\frac{\\lambda^{Y_{i}}e^{-\\lambda}}{Y_{i}!}\\\\\n",
+    "\\pi (\\lambda | D) \\propto \\lambda ^{\\alpha -1}e^{-\\beta\\lambda}*\\prod_{i=1}^{n}\\frac{\\lambda^{Y_{i}}e^{-\\lambda}}{Y_{i}!}\\\\\n",
+    "\\pi (\\lambda | D) \\propto \\lambda ^{\\alpha +\\sum_{i=1}^{n}Y_{i}-1}e^{-\\beta\\lambda}*\\prod_{i=1}^{n}\\frac{e^{-\\lambda}}{Y_{i}!}\\\\\n",
+    "\\pi (\\lambda | D) \\propto \\lambda ^{\\alpha +\\sum_{i=1}^{n}Y_{i}-1}e^{-\\beta\\lambda-n\\lambda}*\\prod_{i=1}^{n}\\frac{1}{Y_{i}!}\\\\\n",
+    "\\pi (\\lambda | D) \\propto \\lambda ^{\\alpha +\\sum_{i=1}^{n}Y_{i}-1}e^{-\\beta\\lambda-n\\lambda}\\\\\n",
+    "$$\n",
+    "$\\pi (\\lambda | D)$ represents a Gamma Distribution: $Gamma(\\alpha +\\sum_{i=1}^{n}Y_{i},\\beta+n)$\n",
+    "This becomes our posterior for $\\lambda$.\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "@webio": {
+   "lastCommId": null,
+   "lastKernelId": null
+  },
+  "kernelspec": {
+   "display_name": "Julia 1.7.2",
+   "language": "julia",
+   "name": "julia-1.7"
+  },
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.7.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}