Skip to content

Instantly share code, notes, and snippets.

@sens

sens/bayes-estimation.ipynb

Last active December 12, 2022 22:37
Show Gist options
  • Save sens/e7975c1aa7cb2a346c0780caa9000c3f to your computer and use it in GitHub Desktop.
Save sens/e7975c1aa7cb2a346c0780caa9000c3f to your computer and use it in GitHub Desktop.
Download ZIP
Calculations for Bayesian estimation of location and scale parameters
Display the source blob
Display the rendered blob
Raw
bayes-estimation.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "7305b972-8b23-498c-89c9-2a91b0310b78",
"metadata": {},
"source": [
"# Bayesian estimation of location and scale parameters from normally distributed samples\n",
"\n",
"In this note we look at the algebra underlying estimation in independently distributed normally distributed observations. Assume we have $n$ iid observations from a normal distribution with mean $\\mu$ and variance $\\phi$. We want to estimate those parameters. To decrease the algebra, we make use of the fact that the sample mean and the sample variance are sufficient statistics for this data. Thus the joint distribution of the data can be expressed as a function of those two statistics."
]
},
{
"cell_type": "markdown",
"id": "0a2542b1-3ff0-4130-a225-00d93d7106fc",
"metadata": {},
"source": [
"## Distribution definitions\n",
"\n",
"First, we define the pdf of the distributions we will use in this note. I have also added some checks to make sure that they are correctly defined."
]
},
{
"cell_type": "markdown",
"id": "518798bd-bccd-4948-ac00-6eb34da919f1",
"metadata": {},
"source": [
"### Normal pdf\n",
"\n",
"The pdf of a normal distribution with location parameter (mean) $\\mu$ and variance (scale) $\\phi$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "73b42e40-fed5-44cf-a564-bab8cadd9323",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{0}$}{\\it normalpdf}\\left(x , \\mu , \\varphi\\right):=\\frac{\\exp \\left(\\frac{-\\left(x-\\mu\\right)^2}{2\\,\\varphi}\\right)}{\\sqrt{2\\,\\pi\\,\\varphi}}\\]"
],
"text/plain": [
" 2\n",
" - (x - mu)\n",
" exp(-----------)\n",
" 2 phi\n",
"(%o0) normalpdf(x, mu, phi) := ----------------\n",
" sqrt(2 %pi phi)"
],
"text/x-maxima": [
"normalpdf(x,mu,phi):=exp((-(x-mu)^2)/(2*phi))/sqrt(2*%pi*phi)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"normalpdf(x,\\mu,\\phi) := exp(-(x-\\mu)^2/(2*\\phi)) / sqrt(2*%pi*\\phi);"
]
},
{
"cell_type": "markdown",
"id": "69805640-2d30-4bea-aeed-bb829f1b379e",
"metadata": {},
"source": [
"We verify that it integrates to 1."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "fd7efc52-2abf-4e9c-ae69-a4da0fff5c23",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{1}$}\\left[ \\varphi>0 \\right] \\]"
],
"text/plain": [
"(%o1) [phi > 0]"
],
"text/x-maxima": [
"[phi > 0]"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{2}$}1\\]"
],
"text/plain": [
"(%o2) 1"
],
"text/x-maxima": [
"1"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"assume(\\phi>0);\n",
"integrate(normalpdf(x,\\mu,\\phi),x,-inf,inf);"
]
},
{
"cell_type": "markdown",
"id": "0990239c-32c8-4e16-9897-6dd51b757f0e",
"metadata": {},
"source": [
"We can also get the Hermite polynomials by repeatedly differentiaing the pdf and dividing by the pdf. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "7c9798a8-4f5f-48f2-a16a-de0f9dad454b",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{3}$}-x\\]"
],
"text/plain": [
"(%o3) - x"
],
"text/x-maxima": [
"-x"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{4}$}x^2-1\\]"
],
"text/plain": [
" 2\n",
"(%o4) x - 1"
],
"text/x-maxima": [
"x^2-1"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{5}$}3\\,x-x^3\\]"
],
"text/plain": [
" 3\n",
"(%o5) 3 x - x"
],
"text/x-maxima": [
"3*x-x^3"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{6}$}x^4-6\\,x^2+3\\]"
],
"text/plain": [
" 4 2\n",
"(%o6) x - 6 x + 3"
],
"text/x-maxima": [
"x^4-6*x^2+3"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diff(normalpdf(x,0,1),x,1)/normalpdf(x,0,1);\n",
"diff(normalpdf(x,0,1),x,2)/normalpdf(x,0,1), ratsimp;\n",
"diff(normalpdf(x,0,1),x,3)/normalpdf(x,0,1), ratsimp;\n",
"diff(normalpdf(x,0,1),x,4)/normalpdf(x,0,1), ratsimp;"
]
},
{
"cell_type": "markdown",
"id": "37844c3b-bc89-41e4-a2bc-20bf3778422e",
"metadata": {},
"source": [
"### $\\chi^2$ pdf\n",
"\n",
"The pdf of a $\\chi^2$ distribution with scale parameter $\\phi$ and $d$ degrees of freedom."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "72020d38-aa86-4719-9061-d91c59522875",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{7}$}{\\it chisqpdf}\\left(x , d , \\varphi\\right):=\\frac{\\left(\\frac{x}{\\varphi}\\right)^{\\frac{d}{2}-1}\\,\\exp \\left(-\\frac{x}{2\\,\\varphi}\\right)}{\\varphi\\,2^{\\frac{d}{2}}\\,\\Gamma\\left(\\frac{d}{2}\\right)}\\]"
],
"text/plain": [
" x d/2 - 1 x\n",
" (---) exp(- -----)\n",
" phi 2 phi\n",
"(%o7) chisqpdf(x, d, phi) := -------------------------\n",
" d/2 d\n",
" phi 2 gamma(-)\n",
" 2"
],
"text/x-maxima": [
"chisqpdf(x,d,phi):=((x/phi)^(d/2-1)*exp(-x/(2*phi)))/(phi*2^(d/2)*gamma(d/2))"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"chisqpdf(x,d,\\phi) := (x/\\phi)^(d/2-1)*exp(-(x/(2*\\phi)))/(\\phi*2^(d/2)*gamma(d/2));"
]
},
{
"cell_type": "markdown",
"id": "4d639994-c4a7-45ce-a2f3-7af2ecbacdee",
"metadata": {},
"source": [
"We verify that the pdf integrates to 1."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "71f452fc-4513-4f46-a8bf-b5d10a587d9b",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdin",
"output_type": "stream",
"text": [
" d\n",
"Is - an integer?\n",
" 2\n",
" no;\n"
]
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{9}$}1\\]"
],
"text/plain": [
"(%o9) 1"
],
"text/x-maxima": [
"1"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"assume(d>0)$\n",
"integrate(chisqpdf(x,d,\\phi),x,0,inf);"
]
},
{
"cell_type": "markdown",
"id": "68bd2a78-9b99-4158-84ba-3698b2c7a651",
"metadata": {},
"source": [
"Expected value of the distribution."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "346e0317-79a4-4de2-af16-23c1e2f5392c",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdin",
"output_type": "stream",
"text": [
" d + 2\n",
"Is ----- an integer?\n",
" 2\n",
" no;\n"
]
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{11}$}d\\,\\varphi\\]"
],
"text/plain": [
"(%o11) d phi"
],
"text/x-maxima": [
"d*phi"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"SB-KERNEL:REDEFINITION-WITH-DEFUN: redefining MAXIMA::SIMP-UNIT-STEP in DEFUN\n",
"SB-KERNEL:REDEFINITION-WITH-DEFUN: redefining MAXIMA::SIMP-POCHHAMMER in DEFUN\n"
]
}
],
"source": [
"gamma_expand:true$\n",
"integrate(x*chisqpdf(x,d,\\phi),x,0,inf);"
]
},
{
"cell_type": "markdown",
"id": "ae78ff95-74aa-4abc-ad70-bfbe755d3250",
"metadata": {},
"source": [
"### Inverse $\\chi^2$ pdf\n",
"\n",
"This is a conjugate distribution for the $\\chi^2$ distribution and is useful for Bayesian inference."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "0d9f71b3-22a2-43de-9a78-9f3c8c64a3da",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{12}$}{\\it invchisqpdf}\\left(x , \\nu , \\tau\\right):=\\frac{\\left(\\frac{\\tau\\,\\nu}{2}\\right)^{\\frac{\\nu}{2}}\\,\\exp \\left(\\frac{-\\nu\\,\\tau}{2\\,x}\\right)}{\\Gamma\\left(\\frac{\\nu}{2}\\right)\\,x^{1+\\frac{\\nu}{2}}}\\]"
],
"text/plain": [
" tau nu nu/2 - nu tau\n",
" (------) exp(--------)\n",
" 2 2 x\n",
"(%o12) invchisqpdf(x, nu, tau) := --------------------------\n",
" nu 1 + nu/2\n",
" gamma(--) x\n",
" 2"
],
"text/x-maxima": [
"invchisqpdf(x,nu,tau):=(((tau*nu)/2)^(nu/2)*exp((-nu*tau)/(2*x)))\n",
" /(gamma(nu/2)*x^(1+nu/2))"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"invchisqpdf(x,\\nu,\\tau) := (\\tau*\\nu/2)^(\\nu/2) * exp(-\\nu*\\tau/(2*x))/(gamma(\\nu/2)*x^(1+\\nu/2));"
]
},
{
"cell_type": "markdown",
"id": "fec01ffc-b304-465a-8a59-1b38555ee9e6",
"metadata": {},
"source": [
"It integrates to 1."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "86e01aee-8426-431c-b53c-5062c593c374",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{14}$}1\\]"
],
"text/plain": [
"(%o14) 1"
],
"text/x-maxima": [
"1"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"assume(\\nu>0,\\tau>0)$\n",
"integrate(invchisqpdf(x,\\nu,\\tau),x,0,inf);"
]
},
{
"cell_type": "markdown",
"id": "edb77bda-11e6-455a-a104-73ba663f2e2f",
"metadata": {},
"source": [
"The expected value."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "b9c1db10-1fc7-4860-adfc-98c1f0df5986",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{17}$}\\frac{\\Gamma\\left(\\frac{\\nu}{2}-1\\right)\\,\\nu\\,\\tau}{2\\,\\Gamma\\left(\\frac{\\nu}{2}\\right)}\\]"
],
"text/plain": [
" nu\n",
" gamma(-- - 1) nu tau\n",
" 2\n",
"(%o17) --------------------\n",
" nu\n",
" 2 gamma(--)\n",
" 2"
],
"text/x-maxima": [
"(gamma(nu/2-1)*nu*tau)/(2*gamma(nu/2))"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gamma_expand:false$\n",
"assume(nu>2)$\n",
"postmean : integrate(x*invchisqpdf(x,\\nu,\\tau),x,0,inf);"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "94a661b9-0338-4455-89de-58ea162cca65",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{18}$}\\frac{\\nu\\,\\tau}{\\nu-2}\\]"
],
"text/plain": [
" nu tau\n",
"(%o18) ------\n",
" nu - 2"
],
"text/x-maxima": [
"(nu*tau)/(nu-2)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"subst(gamma(\\nu/2-1)*(\\nu/2-1),gamma(\\nu/2),postmean), ratsimp;"
]
},
{
"cell_type": "markdown",
"id": "5ebb636d-572d-4f28-b46b-6a94b612356c",
"metadata": {},
"source": [
"The mode."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "801c5a7d-2dcf-49fd-8ded-2c10df43a478",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{19}$}\\left[ x=\\frac{\\nu\\,\\tau}{\\nu+2} \\right] \\]"
],
"text/plain": [
" nu tau\n",
"(%o19) [x = ------]\n",
" nu + 2"
],
"text/x-maxima": [
"[x = (nu*tau)/(nu+2)]"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"solve(radcan(diff(log(invchisqpdf(x,\\nu,\\tau)),x))=0,x);"
]
},
{
"cell_type": "markdown",
"id": "f394aa03-a120-4051-8db6-2ff99abe902f",
"metadata": {},
"source": [
"## MLE for mean with known variance\n",
"\n",
"We now consider the simplest problem in the book. Since the sufficient statistic is the sample mean, we can write the likelihood of the data as a function of the likelihood of the mean which has mean $\\mu$ and variance $\\phi/n$. We denote by $Y$ the sample mean."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "8d7d0d54-0b35-4db0-8067-1d06f36f642e",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{20}$}-\\frac{\\log \\varphi}{2}-\\frac{\\mu^2\\,n}{2\\,\\varphi}+\\frac{Y\\,\\mu\\,n}{\\varphi}-\\frac{Y^2\\,n}{2\\,\\varphi}+\\frac{\\log n}{2}-\\frac{\\log \\pi}{2}-\\frac{\\log 2}{2}\\]"
],
"text/plain": [
" 2 2\n",
" log(phi) mu n Y mu n Y n log(n) log(%pi) log(2)\n",
"(%o20) (- --------) - ----- + ------ - ----- + ------ - -------- - ------\n",
" 2 2 phi phi 2 phi 2 2 2"
],
"text/x-maxima": [
"(-log(phi)/2)-(mu^2*n)/(2*phi)+(Y*mu*n)/phi-(Y^2*n)/(2*phi)+log(n)/2\n",
" -log(%pi)/2-log(2)/2"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"loglik1 : expand(radcan(log(normalpdf(Y,\\mu,\\phi/n))));"
]
},
{
"cell_type": "markdown",
"id": "e96160dc-ae8e-4ccf-a9f3-5cc8cccbdec7",
"metadata": {},
"source": [
"The score function."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "9c706b15-d9f2-47c5-9348-d7d240ae8f38",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{21}$}\\frac{Y\\,n}{\\varphi}-\\frac{\\mu\\,n}{\\varphi}\\]"
],
"text/plain": [
" Y n mu n\n",
"(%o21) --- - ----\n",
" phi phi"
],
"text/x-maxima": [
"(Y*n)/phi-(mu*n)/phi"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"score1 : diff(loglik1,\\mu);"
]
},
{
"cell_type": "markdown",
"id": "16388a66-903c-47e9-9712-a92383e22269",
"metadata": {},
"source": [
"The MLE, by soving the score equation."
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "58c94bc2-3cf0-48a5-808f-41c319679a56",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{22}$}\\left[ \\mu=Y \\right] \\]"
],
"text/plain": [
"(%o22) [mu = Y]"
],
"text/x-maxima": [
"[mu = Y]"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"solve(score1=0,\\mu);"
]
},
{
"cell_type": "markdown",
"id": "f292b5ab-6c42-4d4a-b4ec-a50ec8f0c93b",
"metadata": {},
"source": [
"The Fisher information."
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "14a465e1-c777-4535-8021-d686fb6a64ce",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{23}$}\\frac{n}{\\varphi}\\]"
],
"text/plain": [
" n\n",
"(%o23) ---\n",
" phi"
],
"text/x-maxima": [
"n/phi"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"-diff(loglik1,\\mu,2);"
]
},
{
"cell_type": "markdown",
"id": "7bc8719a-3fa0-4e01-8fad-b096815a6c66",
"metadata": {},
"source": [
"## MLE for mean and variance\n",
"\n",
"Let us denote by $Y$ the sample mean and by $S$ the sum of squares of deviations from the mean. The latter has a $\\chi^2$ distribution with $n-1$ degrees of freedom."
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "a616a927-1af8-471d-ac3e-39893ea5bd5a",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{24}$}-\\frac{n\\,\\varphi\\,\\log \\varphi+\\left(-\\log n+\\left(\\log 2-\\log S\\right)\\,n+2\\,\\log \\Gamma\\left(\\frac{n-1}{2}\\right)+3\\,\\log S+\\log \\pi\\right)\\,\\varphi+\\left(\\mu^2-2\\,Y\\,\\mu+Y^2\\right)\\,n+S}{2\\,\\varphi}\\]"
],
"text/plain": [
"(%o24) - (n phi log(phi) + ((- log(n)) + (log(2) - log(S)) n\n",
" n - 1 2 2\n",
" + 2 log(gamma(-----)) + 3 log(S) + log(%pi)) phi + (mu - 2 Y mu + Y ) n + S)\n",
" 2\n",
"/(2 phi)"
],
"text/x-maxima": [
"-(n*phi*log(phi)+((-log(n))+(log(2)-log(S))*n+2*log(gamma((n-1)/2))+3*log(S)\n",
" +log(%pi))\n",
" *phi+(mu^2-2*Y*mu+Y^2)*n+S)\n",
" /(2*phi)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"loglik2 : log( normalpdf(Y,\\mu,\\phi/n) ) + log( chisqpdf(S,n-1,\\phi) ), radcan;"
]
},
{
"cell_type": "markdown",
"id": "f8d41a33-379e-4330-ab53-688ce4e0e226",
"metadata": {},
"source": [
"The score functions."
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "1c4f86f1-f58e-4164-9c8f-a4155e11a7cc",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{25}$}-\\frac{\\left(\\mu-Y\\right)\\,n}{\\varphi}\\]"
],
"text/plain": [
" (mu - Y) n\n",
"(%o25) - ----------\n",
" phi"
],
"text/x-maxima": [
"-((mu-Y)*n)/phi"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"score2mu : diff(loglik2,\\mu), ratsimp;"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "bf74642c-a526-45c4-b14b-58deca9f8e93",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{26}$}-\\frac{n\\,\\varphi+\\left(-\\mu^2+2\\,Y\\,\\mu-Y^2\\right)\\,n-S}{2\\,\\varphi^2}\\]"
],
"text/plain": [
" 2 2\n",
" n phi + ((- mu ) + 2 Y mu - Y ) n - S\n",
"(%o26) - -------------------------------------\n",
" 2\n",
" 2 phi"
],
"text/x-maxima": [
"-(n*phi+((-mu^2)+2*Y*mu-Y^2)*n-S)/(2*phi^2)"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"score2phi : ratsimp(diff(loglik2,\\phi)), ratsimp;"
]
},
{
"cell_type": "markdown",
"id": "3abf436f-226f-430c-834c-4939adbf4401",
"metadata": {},
"source": [
"When we solve for the MLE we get the usual solutions. Note the denominator for the estimate of $\\phi$. It is one more than the degrees of freedom of $S$."
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "084ed133-bad3-4700-b524-1a69a6ef8a96",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{27}$}\\left[ \\left[ \\mu=Y , \\varphi=\\frac{S}{n} \\right] \\right] \\]"
],
"text/plain": [
" S\n",
"(%o27) [[mu = Y, phi = -]]\n",
" n"
],
"text/x-maxima": [
"[[mu = Y,phi = S/n]]"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"solve([score2mu = 0, score2phi = 0],[\\mu,\\phi]);"
]
},
{
"cell_type": "markdown",
"id": "68ac48c6-b045-4486-8c9b-e98a4450859a",
"metadata": {},
"source": [
"## Bayesian estimation of mean with known variance: non-informative prior\n",
"\n",
"We use the flat prior (Lebesgue measure) as prior for the mean. The posterior distribution is proportional to the likelihood. "
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "d78ee6f1-a383-4eff-ae05-e8af1f1474f8",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{28}$}\\frac{\\sqrt{n}\\,e^ {- \\frac{\\left(Y-\\mu\\right)^2\\,n}{2\\,\\varphi} }}{\\sqrt{2}\\,\\sqrt{\\pi}\\,\\sqrt{\\varphi}}\\]"
],
"text/plain": [
" 2\n",
" (Y - mu) n\n",
" - -----------\n",
" 2 phi\n",
" sqrt(n) %e\n",
"(%o28) ---------------------------\n",
" sqrt(2) sqrt(%pi) sqrt(phi)"
],
"text/x-maxima": [
"(sqrt(n)*%e^-(((Y-mu)^2*n)/(2*phi)))/(sqrt(2)*sqrt(%pi)*sqrt(phi))"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"post1 : normalpdf(Y,\\mu,\\phi/n);"
]
},
{
"cell_type": "markdown",
"id": "57326eb6-7737-4fdf-95ed-0a50693f2199",
"metadata": {},
"source": [
"The posterior mean is just the sample mean $Y$."
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "fbbd1818-9ac2-4b86-8858-612ae68f6a19",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{29}$}\\left[ n>0 \\right] \\]"
],
"text/plain": [
"(%o29) [n > 0]"
],
"text/x-maxima": [
"[n > 0]"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{30}$}0\\]"
],
"text/plain": [
"(%o30) 0"
],
"text/x-maxima": [
"0"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Principal Value\n"
]
}
],
"source": [
"assume(n>0);\n",
"integrate(\\mu*exp(logpost1),\\mu,-inf,inf);"
]
},
{
"cell_type": "markdown",
"id": "8dc28183-0d63-470f-a7bd-117895344e4c",
"metadata": {},
"source": [
"## Bayesian estimation of mean with known variance: informative prior\n",
"\n",
"Consider a proper (informative) prior with mean $\\mu_0$ and variance $\\phi_0$. The unscaled posterior is"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "e1ba3cb7-24d3-4fa8-90d4-35e667d1ef25",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{31}$}\\frac{\\sqrt{n}\\,e^ {- \\frac{\\left(\\mu^2-2\\,Y\\,\\mu+Y^2\\right)\\,n\\,\\varphi_{0}+\\left(\\mu_{0}^2-2\\,\\mu\\,\\mu_{0}+\\mu^2\\right)\\,\\varphi}{2\\,\\varphi\\,\\varphi_{0}} }}{2\\,\\pi\\,\\sqrt{\\varphi}\\,\\sqrt{\\varphi_{0}}}\\]"
],
"text/plain": [
" 2 2 2 2\n",
" (mu - 2 Y mu + Y ) n phi0 + (mu0 - 2 mu mu0 + mu ) phi\n",
" - --------------------------------------------------------\n",
" 2 phi phi0\n",
" sqrt(n) %e\n",
"(%o31) --------------------------------------------------------------------\n",
" 2 %pi sqrt(phi) sqrt(phi0)"
],
"text/x-maxima": [
"(sqrt(n)*%e^-(((mu^2-2*Y*mu+Y^2)*n*phi0+(mu0^2-2*mu*mu0+mu^2)*phi)\n",
" /(2*phi*phi0)))\n",
" /(2*%pi*sqrt(phi)*sqrt(phi0))"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"unscaledpost2 : normalpdf(Y,\\mu,\\phi/n) * normalpdf(\\mu,\\mu0,\\phi0), radcan;"
]
},
{
"cell_type": "markdown",
"id": "4ece9e43-52e4-4f7a-b201-8995f7f749f8",
"metadata": {},
"source": [
"We integrate out the parameter to get the marginal distribution of the data."
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "5de3dfd3-30ef-426e-860f-2c329a669bc8",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{34}$}\\frac{\\sqrt{n}\\,e^ {- \\frac{\\left(\\mu_{0}^2-2\\,Y\\,\\mu_{0}+Y^2\\right)\\,n}{2\\,n\\,\\varphi_{0}+2\\,\\varphi} }}{\\sqrt{2}\\,\\sqrt{\\pi}\\,\\sqrt{n\\,\\varphi_{0}+\\varphi}}\\]"
],
"text/plain": [
" 2 2\n",
" (mu0 - 2 Y mu0 + Y ) n\n",
" - -----------------------\n",
" 2 n phi0 + 2 phi\n",
" sqrt(n) %e\n",
"(%o34) ------------------------------------\n",
" sqrt(2) sqrt(%pi) sqrt(n phi0 + phi)"
],
"text/x-maxima": [
"(sqrt(n)*%e^-(((mu0^2-2*Y*mu0+Y^2)*n)/(2*n*phi0+2*phi)))\n",
" /(sqrt(2)*sqrt(%pi)*sqrt(n*phi0+phi))"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"assume(\\phi>0)$\n",
"assume(\\phi0>0)$\n",
"marg2 : integrate(unscaledpost2,\\mu,-inf,inf);"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "7c4bf78c-a818-44e7-9e2e-7673c5bccae7",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{35}$}\\frac{\\sqrt{n\\,\\varphi_{0}+\\varphi}\\,e^{\\frac{\\left(\\mu_{0}^2-2\\,Y\\,\\mu_{0}+Y^2\\right)\\,n}{2\\,n\\,\\varphi_{0}+2\\,\\varphi}-\\frac{\\left(\\mu^2-2\\,Y\\,\\mu+Y^2\\right)\\,n\\,\\varphi_{0}+\\left(\\mu_{0}^2-2\\,\\mu\\,\\mu_{0}+\\mu^2\\right)\\,\\varphi}{2\\,\\varphi\\,\\varphi_{0}}}}{\\sqrt{2}\\,\\sqrt{\\pi}\\,\\sqrt{\\varphi}\\,\\sqrt{\\varphi_{0}}}\\]"
],
"text/plain": [
" 2 2\n",
" (mu0 - 2 Y mu0 + Y ) n\n",
"(%o35) (sqrt(n phi0 + phi) expt(%e, -----------------------\n",
" 2 n phi0 + 2 phi\n",
" 2 2 2 2\n",
" (mu - 2 Y mu + Y ) n phi0 + (mu0 - 2 mu mu0 + mu ) phi\n",
" - --------------------------------------------------------))\n",
" 2 phi phi0\n",
"/(sqrt(2) sqrt(%pi) sqrt(phi) sqrt(phi0))"
],
"text/x-maxima": [
"(sqrt(n*phi0+phi)*%e^(((mu0^2-2*Y*mu0+Y^2)*n)/(2*n*phi0+2*phi)\n",
" -((mu^2-2*Y*mu+Y^2)*n*phi0+(mu0^2-2*mu*mu0+mu^2)*phi)\n",
" /(2*phi*phi0)))\n",
" /(sqrt(2)*sqrt(%pi)*sqrt(phi)*sqrt(phi0))"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"post2 : unscaledpost2/marg2;"
]
},
{
"cell_type": "markdown",
"id": "e7ad20a2-6d2e-4fa3-a14b-11aa1e1bc068",
"metadata": {},
"source": [
"This still looks ugly, but we can get the posterior mean by integration."
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "3b6bec0f-041e-4793-96db-9c25719ebd09",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{36}$}\\frac{Y\\,n\\,\\varphi_{0}+\\mu_{0}\\,\\varphi}{n\\,\\varphi_{0}+\\varphi}\\]"
],
"text/plain": [
" Y n phi0 + mu0 phi\n",
"(%o36) ------------------\n",
" n phi0 + phi"
],
"text/x-maxima": [
"(Y*n*phi0+mu0*phi)/(n*phi0+phi)"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"postmean2 : ratsimp(integrate(\\mu*post2,\\mu,-inf,inf));"
]
},
{
"cell_type": "markdown",
"id": "341c4bfe-f56a-45a9-8c68-0e22cd6ee9c2",
"metadata": {},
"source": [
"Define $$\\alpha = \\frac{1/\\phi_0}{( 1/\\phi_0 + n/\\phi )}.$$"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "b7070d34-8d4f-4d79-9359-179c7f8916bd",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{37}$}\\left[ \\varphi_{0}=-\\frac{\\left(\\alpha-1\\right)\\,\\varphi}{\\alpha\\,n} \\right] \\]"
],
"text/plain": [
" (alpha - 1) phi\n",
"(%o37) [phi0 = - ---------------]\n",
" alpha n"
],
"text/x-maxima": [
"[phi0 = -((alpha-1)*phi)/(alpha*n)]"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phiexpr : solve((1/\\phi0)/(1/\\phi0+n/\\phi) = \\alpha, \\phi0);"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "e2271d17-bd80-4901-ab76-d37b4fb2c16c",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{38}$}\\alpha\\,\\mu_{0}+Y\\,\\left(1-\\alpha\\right)\\]"
],
"text/plain": [
"(%o38) alpha mu0 + Y (1 - alpha)"
],
"text/x-maxima": [
"alpha*mu0+Y*(1-alpha)"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"collectterms(ratsimp(subst(rhs(phiexpr[1]),\\phi0,postmean2)),Y);"
]
},
{
"cell_type": "markdown",
"id": "206dc753-2796-4c87-b187-458c8dcd09fc",
"metadata": {},
"source": [
"We can see that the posterior mean is a linear combination of the prior mean and the sample mean."
]
},
{
"cell_type": "markdown",
"id": "a5094e6d-6501-471d-b832-d935dfd04fd8",
"metadata": {},
"source": [
"## Bayesian estimation of mean and variance: non-informative prior\n",
"\n",
"We now consider a \"non-informative\" prior for both mean and variance. The non-informative prior for variance is inversely proportional to $\\phi$. The posterior distribution is proportional to:"
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "306b67b4-8a51-4455-b978-811506a0c672",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{39}$}\\frac{S^{\\frac{n-1}{2}-1}\\,\\sqrt{n}\\,2^{-\\frac{n-1}{2}-\\frac{1}{2}}\\,\\varphi^{-\\frac{n-1}{2}-\\frac{3}{2}}\\,e^{-\\frac{\\left(Y-\\mu\\right)^2\\,n}{2\\,\\varphi}-\\frac{S}{2\\,\\varphi}}}{\\sqrt{\\pi}\\,\\Gamma\\left(\\frac{n-1}{2}\\right)}\\]"
],
"text/plain": [
" n - 1 n - 1 n - 1\n",
" ----- - 1 (- -----) - 1/2 (- -----) - 3/2\n",
" 2 2 2\n",
"(%o39) (S sqrt(n) 2 phi\n",
" 2\n",
" (Y - mu) n S\n",
" (- -----------) - -----\n",
" 2 phi 2 phi n - 1\n",
" %e )/(sqrt(%pi) gamma(-----))\n",
" 2"
],
"text/x-maxima": [
"(S^((n-1)/2-1)*sqrt(n)*2^((-(n-1)/2)-1/2)*phi^((-(n-1)/2)-3/2)\n",
" *%e^((-((Y-mu)^2*n)/(2*phi))-S/(2*phi)))\n",
" /(sqrt(%pi)*gamma((n-1)/2))"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"unscaledpost3 : normalpdf(Y,\\mu,\\phi/n) * chisqpdf(S,n-1,\\phi) * (1/\\phi);"
]
},
{
"cell_type": "markdown",
"id": "2c660cdc-6d1f-4fdc-bf72-7d5d0202ad0a",
"metadata": {},
"source": [
"We get the marginal distribution by integrating out. Then we get the posterior distribution."
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "8ff7e1b7-2fcd-4452-b16d-69832ef18f43",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{42}$}\\frac{S^{-\\frac{n}{2}+\\frac{n-1}{2}-\\frac{1}{2}}\\,\\Gamma\\left(\\frac{n}{2}-\\frac{1}{2}\\right)\\,2^{\\frac{n}{2}-\\frac{n-1}{2}-\\frac{1}{2}}}{\\Gamma\\left(\\frac{n-1}{2}\\right)}\\]"
],
"text/plain": [
" n - 1 n - 1\n",
" (- n/2) + ----- - 1/2 n/2 - ----- - 1/2\n",
" 2 n 1 2\n",
" S gamma(- - -) 2\n",
" 2 2\n",
"(%o42) ------------------------------------------------------\n",
" n - 1\n",
" gamma(-----)\n",
" 2"
],
"text/x-maxima": [
"(S^((-n/2)+(n-1)/2-1/2)*gamma(n/2-1/2)*2^(n/2-(n-1)/2-1/2))/gamma((n-1)/2)"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"assume(n>1)$\n",
"assume(S>0)$\n",
"marg3 : integrate(integrate(unscaledpost3,\\mu,-inf,inf),\\phi,0,inf);"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "c5fecf6d-84ac-496f-878a-98961aac2267",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{43}$}\\frac{S^{\\frac{n}{2}-\\frac{1}{2}}\\,\\sqrt{n}\\,\\varphi^{-\\frac{n-1}{2}-\\frac{3}{2}}\\,e^{-\\frac{\\left(Y-\\mu\\right)^2\\,n}{2\\,\\varphi}-\\frac{S}{2\\,\\varphi}}}{\\sqrt{\\pi}\\,\\Gamma\\left(\\frac{n}{2}-\\frac{1}{2}\\right)\\,2^{\\frac{n}{2}}}\\]"
],
"text/plain": [
" 2\n",
" n - 1 (Y - mu) n S\n",
" (- -----) - 3/2 (- -----------) - -----\n",
" n/2 - 1/2 2 2 phi 2 phi\n",
" S sqrt(n) phi %e\n",
"(%o43) ---------------------------------------------------------------\n",
" n 1 n/2\n",
" sqrt(%pi) gamma(- - -) 2\n",
" 2 2"
],
"text/x-maxima": [
"(S^(n/2-1/2)*sqrt(n)*phi^((-(n-1)/2)-3/2)\n",
" *%e^((-((Y-mu)^2*n)/(2*phi))-S/(2*phi)))\n",
" /(sqrt(%pi)*gamma(n/2-1/2)*2^(n/2))"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"post3 : unscaledpost3/marg3;"
]
},
{
"cell_type": "markdown",
"id": "7a145382-62d1-4972-82d2-47b38f18c4b8",
"metadata": {},
"source": [
"Let us get the mode of the posterior distribution (the analog of MLE)."
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "20d7f311-402f-43ed-beaf-4feed705993b",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{44}$}-\\frac{\\left(\\mu-Y\\right)\\,n}{\\varphi}\\]"
],
"text/plain": [
" (mu - Y) n\n",
"(%o44) - ----------\n",
" phi"
],
"text/x-maxima": [
"-((mu-Y)*n)/phi"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{45}$}-\\frac{\\left(n+2\\right)\\,\\varphi+\\left(-\\mu^2+2\\,Y\\,\\mu-Y^2\\right)\\,n-S}{2\\,\\varphi^2}\\]"
],
"text/plain": [
" 2 2\n",
" (n + 2) phi + ((- mu ) + 2 Y mu - Y ) n - S\n",
"(%o45) - -------------------------------------------\n",
" 2\n",
" 2 phi"
],
"text/x-maxima": [
"-((n+2)*phi+((-mu^2)+2*Y*mu-Y^2)*n-S)/(2*phi^2)"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"postscore31 : ratsimp(diff(log(post3),\\mu));\n",
"postscore32 : ratsimp(diff(log(post3),\\phi));"
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "6a46387b-eea4-42ed-85db-23a13142577a",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{46}$}\\left[ \\left[ \\mu=Y , \\varphi=\\frac{S}{n+2} \\right] \\right] \\]"
],
"text/plain": [
" S\n",
"(%o46) [[mu = Y, phi = -----]]\n",
" n + 2"
],
"text/x-maxima": [
"[[mu = Y,phi = S/(n+2)]]"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"solve([postscore31=0,postscore32=0],[\\mu,\\phi]);"
]
},
{
"cell_type": "markdown",
"id": "aab882ca-1416-45df-83fb-1499f8a4a6da",
"metadata": {},
"source": [
"Notice that the estimate of the location parameter is still the sample mean. The estimate of the scale parameter has divisor $n+2$ instead of $n$ for the MLE. Essentially, the prior is shrinking the estimate for the scale parameter a bit more towards zero.\n",
"\n",
"Now, instead of getting the MAP (maximum a posteriori) estimate we can go ahead and get the posterior distribution of the scale parameter by integrating out $\\mu$."
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "1c5546f0-51fd-458b-9a33-201c3844e210",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{47}$}\\frac{S^{\\frac{n-1}{2}}\\,e^ {- \\frac{S}{2\\,\\varphi} }}{\\Gamma\\left(\\frac{n-1}{2}\\right)\\,2^{\\frac{n-1}{2}}\\,\\varphi^{\\frac{n+1}{2}}}\\]"
],
"text/plain": [
" n - 1 S\n",
" ----- - -----\n",
" 2 2 phi\n",
" S %e\n",
"(%o47) ----------------------------\n",
" n - 1 n + 1\n",
" ----- -----\n",
" n - 1 2 2\n",
" gamma(-----) 2 phi\n",
" 2"
],
"text/x-maxima": [
"(S^((n-1)/2)*%e^-(S/(2*phi)))/(gamma((n-1)/2)*2^((n-1)/2)*phi^((n+1)/2))"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"postphi : ratsimp(integrate(post3,\\mu,-inf,inf)), radcan;"
]
},
{
"cell_type": "markdown",
"id": "98339fc2-bf3f-42ec-b60b-d788139a16c3",
"metadata": {},
"source": [
"This is the same as an inverse $\\chi^2$ distribution with $n-1$ degrees of freedom and scale $S/(n-1)$."
]
},
{
"cell_type": "code",
"execution_count": 34,
"id": "c6d736cb-c802-492f-8b8c-92d2a3c6ac74",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{48}$}1\\]"
],
"text/plain": [
"(%o48) 1"
],
"text/x-maxima": [
"1"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ratsimp(postphi/invchisqpdf(\\phi,n-1,S/(n-1)));"
]
},
{
"cell_type": "markdown",
"id": "3e739fce-82c4-410e-89a0-2db61b446821",
"metadata": {},
"source": [
"The mean of this distribution is:"
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "5b5c96ba-92f7-4caf-907d-cbe5e4ed8ace",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{49}$}\\frac{S}{n-3}\\]"
],
"text/plain": [
" S\n",
"(%o49) -----\n",
" n - 3"
],
"text/x-maxima": [
"S/(n-3)"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(S/(n-1))*(n-1)/(n-1-2);"
]
},
{
"cell_type": "markdown",
"id": "0c79f5bb-3e7c-4a82-a356-abca5f17feb0",
"metadata": {},
"source": [
"The mode of this distribution (the marginal posterior of $\\phi$) is"
]
},
{
"cell_type": "code",
"execution_count": 36,
"id": "ae254da6-53f9-41fc-8807-d96a263660a6",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{50}$}\\frac{S}{n+1}\\]"
],
"text/plain": [
" S\n",
"(%o50) -----\n",
" n + 1"
],
"text/x-maxima": [
"S/(n+1)"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(S/(n-1))*(n-1)/(n-1+2);"
]
},
{
"cell_type": "markdown",
"id": "5fa7e14e-ea42-4c07-9096-ac64c282847c",
"metadata": {},
"source": [
"Note that the mode of the joint posterior was attained at $$\\frac{S}{n+2}.$$"
]
},
{
"cell_type": "markdown",
"id": "55aa0dc1-9079-4cfd-b64d-60e78417bb8a",
"metadata": {},
"source": [
"## Bayesian estimation of mean and variance: informative prior\n",
"\n",
"Now we use a non-informative prior for the mean and an informative, inverse $\\chi^2$ prior for the variance with degrees of freedom $n_0$ and scale parameter$\\phi_0$."
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "cffc1025-db6e-41e4-8d96-cc86ab583335",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{51}$}\\frac{S^{\\frac{n-1}{2}-1}\\,\\sqrt{n}\\,n_{0}^{\\frac{n_{0}}{2}}\\,2^{-\\frac{n_{0}}{2}-\\frac{n-1}{2}-\\frac{1}{2}}\\,\\varphi^{-\\frac{n_{0}}{2}-\\frac{n-1}{2}-\\frac{3}{2}}\\,\\varphi_{0}^{\\frac{n_{0}}{2}}\\,e^{-\\frac{n_{0}\\,\\varphi_{0}}{2\\,\\varphi}-\\frac{\\left(Y-\\mu\\right)^2\\,n}{2\\,\\varphi}-\\frac{S}{2\\,\\varphi}}}{\\sqrt{\\pi}\\,\\Gamma\\left(\\frac{n-1}{2}\\right)\\,\\Gamma\\left(\\frac{n_{0}}{2}\\right)}\\]"
],
"text/plain": [
" n - 1 n - 1\n",
" ----- - 1 (- n0/2) - ----- - 1/2\n",
" 2 n0/2 2\n",
"(%o51) (S sqrt(n) n0 2\n",
" 2\n",
" n - 1 n0 phi0 (Y - mu) n S\n",
" (- n0/2) - ----- - 3/2 (- -------) - ----------- - -----\n",
" 2 n0/2 2 phi 2 phi 2 phi\n",
" phi phi0 %e )\n",
" n - 1 n0\n",
"/(sqrt(%pi) gamma(-----) gamma(--))\n",
" 2 2"
],
"text/x-maxima": [
"(S^((n-1)/2-1)*sqrt(n)*n0^(n0/2)*2^((-n0/2)-(n-1)/2-1/2)\n",
" *phi^((-n0/2)-(n-1)/2-3/2)*phi0^(n0/2)\n",
" *%e^((-(n0*phi0)/(2*phi))-((Y-mu)^2*n)/(2*phi)-S/(2*phi)))\n",
" /(sqrt(%pi)*gamma((n-1)/2)*gamma(n0/2))"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"post4 : normalpdf(Y,\\mu,\\phi/n) * chisqpdf(S,n-1,\\phi) * invchisqpdf(\\phi,n0,\\phi0);"
]
},
{
"cell_type": "code",
"execution_count": 38,
"id": "51b26d05-08cb-41a2-be29-fa3ee5f869bf",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{54}$}\\frac{S^{\\frac{n-3}{2}}\\,n_{0}^{\\frac{n_{0}}{2}}\\,\\varphi_{0}^{\\frac{n_{0}}{2}}\\,e^ {- \\frac{n_{0}\\,\\varphi_{0}+S}{2\\,\\varphi} }}{\\Gamma\\left(\\frac{n-1}{2}\\right)\\,\\Gamma\\left(\\frac{n_{0}}{2}\\right)\\,2^{\\frac{n_{0}+n-1}{2}}\\,\\varphi^{\\frac{n_{0}+n+1}{2}}}\\]"
],
"text/plain": [
" n - 3 n0 phi0 + S\n",
" ----- - -----------\n",
" 2 n0/2 n0/2 2 phi\n",
" S n0 phi0 %e\n",
"(%o54) ------------------------------------------------\n",
" n0 + n - 1 n0 + n + 1\n",
" ---------- ----------\n",
" n - 1 n0 2 2\n",
" gamma(-----) gamma(--) 2 phi\n",
" 2 2"
],
"text/x-maxima": [
"(S^((n-3)/2)*n0^(n0/2)*phi0^(n0/2)*%e^-((n0*phi0+S)/(2*phi)))\n",
" /(gamma((n-1)/2)*gamma(n0/2)*2^((n0+n-1)/2)*phi^((n0+n+1)/2))"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"assume(n0>0)$\n",
"a : ratsimp(integrate(post4,\\mu,-inf,inf))$\n",
"a : radcan(a);\n",
"a: subst(n1,n0+n+1,a)$\n",
"a: subst(n1-2,n0+n-1,a)$\n",
"a: subst(S1,n0*\\phi0+S,a)$"
]
},
{
"cell_type": "code",
"execution_count": 39,
"id": "54af1820-73ca-4d4e-b169-bd129006b8a6",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{60}$}\\frac{S^{\\frac{n-3}{2}}\\,n_{0}^{\\frac{n_{0}}{2}}\\,\\Gamma\\left(\\frac{n_{1}-2}{2}\\right)\\,\\varphi_{0}^{\\frac{n_{0}}{2}}}{S_{1}^{\\frac{n_{1}-2}{2}}\\,\\Gamma\\left(\\frac{n-1}{2}\\right)\\,\\Gamma\\left(\\frac{n_{0}}{2}\\right)}\\]"
],
"text/plain": [
" n - 3\n",
" -----\n",
" 2 n0/2 n1 - 2 n0/2\n",
" S n0 gamma(------) phi0\n",
" 2\n",
"(%o60) ------------------------------------\n",
" n1 - 2\n",
" ------\n",
" 2 n - 1 n0\n",
" S1 gamma(-----) gamma(--)\n",
" 2 2"
],
"text/x-maxima": [
"(S^((n-3)/2)*n0^(n0/2)*gamma((n1-2)/2)*phi0^(n0/2))\n",
" /(S1^((n1-2)/2)*gamma((n-1)/2)*gamma(n0/2))"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"assume(n1-2>0)$\n",
"assume(S1>0)$\n",
"integrate(a,\\phi,0,inf), radcan;"
]
},
{
"cell_type": "markdown",
"id": "5367d6f8-de41-424e-af0b-a71024644598",
"metadata": {},
"source": [
"Still looks ugly."
]
},
{
"cell_type": "markdown",
"id": "72fe1acf-4860-4c17-9903-0569be68c1d4",
"metadata": {},
"source": [
"Wel solve for the posterior mode."
]
},
{
"cell_type": "code",
"execution_count": 40,
"id": "b6a5a58b-abe1-4a5c-815d-29383c3ae15e",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{61}$}-\\frac{\\left(\\mu-Y\\right)\\,n}{\\varphi}\\]"
],
"text/plain": [
" (mu - Y) n\n",
"(%o61) - ----------\n",
" phi"
],
"text/x-maxima": [
"-((mu-Y)*n)/phi"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{62}$}\\frac{n_{0}\\,\\varphi_{0}+\\left(-n_{0}-n-2\\right)\\,\\varphi+\\left(\\mu^2-2\\,Y\\,\\mu+Y^2\\right)\\,n+S}{2\\,\\varphi^2}\\]"
],
"text/plain": [
" 2 2\n",
" n0 phi0 + ((- n0) - n - 2) phi + (mu - 2 Y mu + Y ) n + S\n",
"(%o62) ----------------------------------------------------------\n",
" 2\n",
" 2 phi"
],
"text/x-maxima": [
"(n0*phi0+((-n0)-n-2)*phi+(mu^2-2*Y*mu+Y^2)*n+S)/(2*phi^2)"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"postscore41 : ratsimp(diff(log(post4),\\mu));\n",
"postscore42 : ratsimp(diff(log(post4),\\phi));"
]
},
{
"cell_type": "code",
"execution_count": 41,
"id": "c6a1ef03-d879-4c75-a562-94b82472afc7",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{63}$}\\left[ \\left[ \\mu=Y , \\varphi=\\frac{n_{0}\\,\\varphi_{0}+S}{n_{0}+n+2} \\right] \\right] \\]"
],
"text/plain": [
" n0 phi0 + S\n",
"(%o63) [[mu = Y, phi = -----------]]\n",
" n0 + n + 2"
],
"text/x-maxima": [
"[[mu = Y,phi = (n0*phi0+S)/(n0+n+2)]]"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"solve([postscore41=0,postscore42=0],[\\mu,\\phi]);"
]
},
{
"cell_type": "markdown",
"id": "f672bc91-16e0-45d7-a140-8b7623557ddb",
"metadata": {},
"source": [
"Note that this has the same form as that of the mode for the non-informative prior case. The denominator is prior plus data degrees of freedom plus 3. The numerator is the prior plus data sum of squares."
]
},
{
"cell_type": "markdown",
"id": "0933aadd-69be-43c3-99cf-98152a5759cb",
"metadata": {},
"source": [
"## Linear regression\n",
"\n",
"Let us consider the linear regression case where the observations are normally distributed with mean:\n",
"\n",
"$$E(y_i) = \\beta_0 + \\beta_1 x_i,$$\n",
"and $$V(y_i) = \\phi.$$\n",
"\n",
"We can then write the log likelihood as follows:"
]
},
{
"cell_type": "code",
"execution_count": 42,
"id": "0e9f92e7-8d33-47ad-8a03-c19bae7c13e3",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{64}$}-\\frac{n\\,\\log \\left(\\pi\\,\\varphi\\right)}{2}-\\frac{\\sum_{i=1}^{n}{\\left(y_{i}-\\beta_{1}\\,x_{i}-\\beta_{0}\\right)^2}}{2\\,\\varphi}\\]"
],
"text/plain": [
" n\n",
" ====\n",
" \\ 2\n",
" > (y - beta1 x - beta0)\n",
" / i i\n",
" ====\n",
" n log(%pi phi) i = 1\n",
"(%o64) (- --------------) - ------------------------------\n",
" 2 2 phi"
],
"text/x-maxima": [
"(-(n*log(%pi*phi))/2)-('sum((y[i]-beta1*x[i]-beta0)^2,i,1,n))/(2*phi)"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"loglik : sum( -(y[i]-\\beta0 -beta1*x[i])^2/(2*\\phi),i,1,n ) - (n/2)*log(%pi*\\phi);"
]
},
{
"cell_type": "markdown",
"id": "a031156c-de7e-4c0a-b6ac-778f2c4dd91f",
"metadata": {},
"source": [
"We declare the sum function to be linear to enable simplifications."
]
},
{
"cell_type": "code",
"execution_count": 43,
"id": "a52feaf6-102a-4434-9e62-c1b17ae42f70",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{65}$}\\mathbf{done}\\]"
],
"text/plain": [
"(%o65) done"
],
"text/x-maxima": [
"done"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"declare (sum, linear);"
]
},
{
"cell_type": "markdown",
"id": "ed7e93c3-730a-45d6-b1cd-eec36d6890d5",
"metadata": {},
"source": [
"We take derivatives with respect to the parameters of interest."
]
},
{
"cell_type": "code",
"execution_count": 44,
"id": "360cef42-b5d8-4724-afe1-94e3f5e04eb5",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{66}$}-\\frac{\\beta_{0}\\,n}{\\varphi}+\\frac{\\sum_{i=1}^{n}{y_{i}}}{\\varphi}-\\frac{\\beta_{1}\\,\\sum_{i=1}^{n}{x_{i}}}{\\varphi}\\]"
],
"text/plain": [
" n n\n",
" ==== ====\n",
" \\ \\\n",
" > y beta1 > x\n",
" / i / i\n",
" ==== ====\n",
" beta0 n i = 1 i = 1\n",
"(%o66) (- -------) + -------- - --------------\n",
" phi phi phi"
],
"text/x-maxima": [
"(-(beta0*n)/phi)+('sum(y[i],i,1,n))/phi-(beta1*'sum(x[i],i,1,n))/phi"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{67}$}\\frac{\\sum_{i=1}^{n}{x_{i}\\,y_{i}}}{\\varphi}-\\frac{\\beta_{1}\\,\\sum_{i=1}^{n}{x_{i}^2}}{\\varphi}-\\frac{\\beta_{0}\\,\\sum_{i=1}^{n}{x_{i}}}{\\varphi}\\]"
],
"text/plain": [
" n n n\n",
" ==== ==== ====\n",
" \\ \\ 2 \\\n",
" > x y beta1 > x beta0 > x\n",
" / i i / i / i\n",
" ==== ==== ====\n",
" i = 1 i = 1 i = 1\n",
"(%o67) ----------- - -------------- - --------------\n",
" phi phi phi"
],
"text/x-maxima": [
"('sum(x[i]*y[i],i,1,n))/phi-(beta1*'sum(x[i]^2,i,1,n))/phi\n",
" -(beta0*'sum(x[i],i,1,n))/phi"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{68}$}\\frac{\\sum_{i=1}^{n}{\\left(y_{i}-\\beta_{1}\\,x_{i}-\\beta_{0}\\right)^2}}{2\\,\\varphi^2}-\\frac{n}{2\\,\\varphi}\\]"
],
"text/plain": [
" n\n",
" ====\n",
" \\ 2\n",
" > (y - beta1 x - beta0)\n",
" / i i\n",
" ====\n",
" i = 1 n\n",
"(%o68) ------------------------------ - -----\n",
" 2 2 phi\n",
" 2 phi"
],
"text/x-maxima": [
"('sum((y[i]-beta1*x[i]-beta0)^2,i,1,n))/(2*phi^2)-n/(2*phi)"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"scoreb0 : expand(diff(loglik,\\beta0));\n",
"scoreb1 : expand(diff(loglik,\\beta1));\n",
"scorephi : (diff(loglik,\\phi));"
]
},
{
"cell_type": "markdown",
"id": "5f1376df-b021-4382-bb15-2da42b664bdb",
"metadata": {},
"source": [
"To simplify the formulas, we assume, without loss of generality that the $x$'s have been standardized so that the have sum zero and sum of squares n. Then we have some further simplification of these expressions."
]
},
{
"cell_type": "code",
"execution_count": 45,
"id": "5a71187e-ecd7-4f77-9408-3f9539234b07",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{69}$}\\frac{\\sum_{i=1}^{n}{y_{i}}}{\\varphi}-\\frac{\\beta_{0}\\,n}{\\varphi}\\]"
],
"text/plain": [
" n\n",
" ====\n",
" \\\n",
" > y\n",
" / i\n",
" ====\n",
" i = 1 beta0 n\n",
"(%o69) -------- - -------\n",
" phi phi"
],
"text/x-maxima": [
"('sum(y[i],i,1,n))/phi-(beta0*n)/phi"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{71}$}\\frac{\\sum_{i=1}^{n}{x_{i}\\,y_{i}}}{\\varphi}-\\frac{\\beta_{1}\\,n}{\\varphi}\\]"
],
"text/plain": [
" n\n",
" ====\n",
" \\\n",
" > x y\n",
" / i i\n",
" ====\n",
" i = 1 beta1 n\n",
"(%o71) ----------- - -------\n",
" phi phi"
],
"text/x-maxima": [
"('sum(x[i]*y[i],i,1,n))/phi-(beta1*n)/phi"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"scoreb0 : subst(0,sum(x[i],i,1,n),scoreb0);\n",
"scoreb1 : subst(0,sum(x[i],i,1,n),expand(scoreb1))$\n",
"scoreb1 : subst(n,sum(x[i]^2,i,1,n),expand(scoreb1));"
]
},
{
"cell_type": "code",
"execution_count": 46,
"id": "ae9f6bde-5b85-4681-9fcb-c9cd605ae698",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{72}$}\\left[ \\beta_{0}=\\frac{\\sum_{i=1}^{n}{y_{i}}}{n} , \\beta_{1}=\\frac{\\sum_{i=1}^{n}{x_{i}\\,y_{i}}}{n} , \\varphi=\\frac{\\sum_{i=1}^{n}{\\left(y_{i}-\\beta_{1}\\,x_{i}-\\beta_{0}\\right)^2}}{n} \\right] \\]"
],
"text/plain": [
" n n\n",
" ==== ====\n",
" \\ \\\n",
" > y > x y\n",
" / i / i i\n",
" ==== ====\n",
" i = 1 i = 1\n",
"(%o72) [beta0 = --------, beta1 = -----------, \n",
" n n\n",
" n\n",
" ====\n",
" \\ 2\n",
" > (y - beta1 x - beta0)\n",
" / i i\n",
" ====\n",
" i = 1\n",
" phi = ------------------------------]\n",
" n"
],
"text/x-maxima": [
"[beta0 = ('sum(y[i],i,1,n))/n,beta1 = ('sum(x[i]*y[i],i,1,n))/n,\n",
" phi = ('sum((y[i]-beta1*x[i]-beta0)^2,i,1,n))/n]"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"linsolve([expand(scoreb0)=0,expand(scoreb1)=0,scorephi=0],[\\beta0,\\beta1,phi]);"
]
},
{
"cell_type": "markdown",
"id": "9bf6ed9a-b3a4-4c1d-a238-639d6d56dc19",
"metadata": {},
"source": [
"It can also be shown that these three estimates are sufficient statistics and so we define three statistics for the linear regression case.\n",
"\n",
"$$B0=\\frac{1}{n}\\sum y_i$$\n",
"$$B1=\\frac{1}{n}\\sum y_i x_i$$\n",
"$$S=\\sum (y_i - B0 - B1 x_i)^2$$\n",
"\n",
"It is easily shown that $B0$, $B1$ and $S$ are independent. The first two are normally distributed with mean $\\beta_0$ and $\\beta_1$ respectively and variance $\\phi/n$. The third term has a $\\chi^2$ distribution with $n-2$ degrees of freedom, multiplied by $\\phi$. Thus, we can write the likelihood for the linear regression case as follows."
]
},
{
"cell_type": "code",
"execution_count": 47,
"id": "af21136d-2b88-4bc1-9873-978604a46a65",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{73}$}\\frac{S^{\\frac{n-2}{2}-1}\\,n\\,2^{-\\frac{n-2}{2}-1}\\,\\varphi^{-\\frac{n-2}{2}-1}\\,e^{-\\frac{\\left(B_{1}-\\beta_{1}\\right)^2\\,n}{2\\,\\varphi}-\\frac{\\left(B_{0}-\\beta_{0}\\right)^2\\,n}{2\\,\\varphi}-\\frac{S}{2\\,\\varphi}}}{\\pi\\,\\Gamma\\left(\\frac{n-2}{2}\\right)}\\]"
],
"text/plain": [
" n - 2 n - 2 n - 2\n",
" ----- - 1 (- -----) - 1 (- -----) - 1\n",
" 2 2 2\n",
"(%o73) (S n 2 phi\n",
" 2 2\n",
" (B1 - beta1) n (B0 - beta0) n S\n",
" (- ---------------) - --------------- - -----\n",
" 2 phi 2 phi 2 phi n - 2\n",
" %e )/(%pi gamma(-----))\n",
" 2"
],
"text/x-maxima": [
"(S^((n-2)/2-1)*n*2^((-(n-2)/2)-1)*phi^((-(n-2)/2)-1)\n",
" *%e^((-((B1-beta1)^2*n)/(2*phi))-((B0-beta0)^2*n)/(2*phi)\n",
" -S/(2*phi)))\n",
" /(%pi*gamma((n-2)/2))"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"linreglik : normalpdf(B0,\\beta0,\\phi/n) * normalpdf(B1,\\beta1,\\phi/n) * chisqpdf(S,n-2,\\phi);"
]
},
{
"cell_type": "code",
"execution_count": 48,
"id": "70c3227e-5664-4a4c-b7b4-2cacb1e1a781",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{74}$}\\frac{S^{\\frac{n-4}{2}}\\,n\\,e^ {- \\frac{\\left(\\beta_{1}^2-2\\,B_{1}\\,\\beta_{1}+\\beta_{0}^2-2\\,B_{0}\\,\\beta_{0}+B_{1}^2+B_{0}^2\\right)\\,n+S}{2\\,\\varphi} }}{\\pi\\,\\Gamma\\left(\\frac{n-2}{2}\\right)\\,2^{\\frac{n}{2}}\\,\\varphi^{\\frac{n}{2}}}\\]"
],
"text/plain": [
"(%o74) \n",
" 2 2 2 2\n",
" n - 4 (beta1 - 2 B1 beta1 + beta0 - 2 B0 beta0 + B1 + B0 ) n + S\n",
" ----- - -------------------------------------------------------------\n",
" 2 2 phi\n",
" S n %e\n",
" --------------------------------------------------------------------------\n",
" n - 2 n/2 n/2\n",
" %pi gamma(-----) 2 phi\n",
" 2"
],
"text/x-maxima": [
"(S^((n-4)/2)*n\n",
" *%e^-(((beta1^2-2*B1*beta1+beta0^2-2*B0*beta0+B1^2+B0^2)*n+S)\n",
" /(2*phi)))\n",
" /(%pi*gamma((n-2)/2)*2^(n/2)*phi^(n/2))"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"radcan(linreglik);"
]
},
{
"cell_type": "markdown",
"id": "063d9dab-960d-4cda-a5b9-34085a02956e",
"metadata": {},
"source": [
"With an informative inverse $\\chi^2$ prior as before, we get that the posterior distribution is proportional to:"
]
},
{
"cell_type": "code",
"execution_count": 49,
"id": "c2d595b8-307e-42cd-8b89-a82ac2e2ca3b",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{75}$}\\frac{S^{\\frac{n-4}{2}}\\,n\\,n_{0}^{\\frac{n_{0}}{2}}\\,\\varphi_{0}^{\\frac{n_{0}}{2}}\\,e^ {- \\frac{n_{0}\\,\\varphi_{0}+\\left(\\beta_{1}^2-2\\,B_{1}\\,\\beta_{1}+\\beta_{0}^2-2\\,B_{0}\\,\\beta_{0}+B_{1}^2+B_{0}^2\\right)\\,n+S}{2\\,\\varphi} }}{\\pi\\,\\Gamma\\left(\\frac{n-2}{2}\\right)\\,\\Gamma\\left(\\frac{n_{0}}{2}\\right)\\,2^{\\frac{n_{0}+n}{2}}\\,\\varphi^{\\frac{n_{0}+n+2}{2}}}\\]"
],
"text/plain": [
" n - 4\n",
" -----\n",
" 2 n0/2 n0/2\n",
"(%o75) (S n n0 phi0\n",
" 2 2 2 2\n",
" n0 phi0 + (beta1 - 2 B1 beta1 + beta0 - 2 B0 beta0 + B1 + B0 ) n + S\n",
" - -----------------------------------------------------------------------\n",
" 2 phi\n",
" %e )\n",
" n0 + n n0 + n + 2\n",
" ------ ----------\n",
" n - 2 n0 2 2\n",
"/(%pi gamma(-----) gamma(--) 2 phi )\n",
" 2 2"
],
"text/x-maxima": [
"(S^((n-4)/2)*n*n0^(n0/2)*phi0^(n0/2)\n",
" *%e^-((n0*phi0+(beta1^2-2*B1*beta1+beta0^2-2*B0*beta0+B1^2+B0^2)*n\n",
" +S)\n",
" /(2*phi)))\n",
" /(%pi*gamma((n-2)/2)*gamma(n0/2)*2^((n0+n)/2)*phi^((n0+n+2)/2))"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"post5 : linreglik * invchisqpdf(\\phi,n0,\\phi0), radcan;"
]
},
{
"cell_type": "markdown",
"id": "d0a020da-53e6-443e-b812-47de7898a862",
"metadata": {},
"source": [
"Differentiating we get:"
]
},
{
"cell_type": "code",
"execution_count": 50,
"id": "5b73ffed-e7e0-4d0b-90a7-0d234ec0d61a",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{76}$}-\\frac{\\left(\\beta_{0}-B_{0}\\right)\\,n}{\\varphi}\\]"
],
"text/plain": [
" (beta0 - B0) n\n",
"(%o76) - --------------\n",
" phi"
],
"text/x-maxima": [
"-((beta0-B0)*n)/phi"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{77}$}-\\frac{\\left(\\beta_{1}-B_{1}\\right)\\,n}{\\varphi}\\]"
],
"text/plain": [
" (beta1 - B1) n\n",
"(%o77) - --------------\n",
" phi"
],
"text/x-maxima": [
"-((beta1-B1)*n)/phi"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{78}$}\\frac{n_{0}\\,\\varphi_{0}+\\left(-n_{0}-n-2\\right)\\,\\varphi+\\left(\\beta_{1}^2-2\\,B_{1}\\,\\beta_{1}+\\beta_{0}^2-2\\,B_{0}\\,\\beta_{0}+B_{1}^2+B_{0}^2\\right)\\,n+S}{2\\,\\varphi^2}\\]"
],
"text/plain": [
" 2 2\n",
"(%o78) (n0 phi0 + ((- n0) - n - 2) phi + (beta1 - 2 B1 beta1 + beta0\n",
" 2 2 2\n",
" - 2 B0 beta0 + B1 + B0 ) n + S)/(2 phi )"
],
"text/x-maxima": [
"(n0*phi0+((-n0)-n-2)*phi+(beta1^2-2*B1*beta1+beta0^2-2*B0*beta0+B1^2+B0^2)*n\n",
" +S)\n",
" /(2*phi^2)"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"postscore50 : ratsimp(diff(log(post5),\\beta0));\n",
"postscore51 : ratsimp(diff(log(post5),\\beta1));\n",
"postscore52 : ratsimp(diff(log(post5),\\phi));"
]
},
{
"cell_type": "markdown",
"id": "9f7a7e61-a552-4f4a-9d61-296a4f6df6ed",
"metadata": {},
"source": [
"The solutions are:"
]
},
{
"cell_type": "code",
"execution_count": 51,
"id": "7ac44472-59f0-40e8-b630-566f31fc1e38",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{79}$}\\left[ \\left[ \\beta_{0}=B_{0} , \\beta_{1}=B_{1} , \\varphi=\\frac{n_{0}\\,\\varphi_{0}+S}{n_{0}+n+2} \\right] \\right] \\]"
],
"text/plain": [
" n0 phi0 + S\n",
"(%o79) [[beta0 = B0, beta1 = B1, phi = -----------]]\n",
" n0 + n + 2"
],
"text/x-maxima": [
"[[beta0 = B0,beta1 = B1,phi = (n0*phi0+S)/(n0+n+2)]]"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"solve([postscore50=0,postscore51=0,postscore52=0],[\\beta0,\\beta1,\\phi]);"
]
},
{
"cell_type": "markdown",
"id": "717fce3b-24e9-4c80-8528-3def763cdc69",
"metadata": {},
"source": [
"Note that the formula for the posterior mode for the variance parameter has a formula identical to the normal mean case.\n",
"\n",
"If we want to do maximum likelihood, we can do the following."
]
},
{
"cell_type": "code",
"execution_count": 52,
"id": "bf6cc077-7329-4d0b-a5db-1c06ce664f69",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{80}$}\\frac{S^{\\frac{n-4}{2}}\\,n\\,e^ {- \\frac{\\left(\\beta_{1}^2-2\\,B_{1}\\,\\beta_{1}+\\beta_{0}^2-2\\,B_{0}\\,\\beta_{0}+B_{1}^2+B_{0}^2\\right)\\,n+S}{2\\,\\varphi} }}{\\pi\\,\\Gamma\\left(\\frac{n-2}{2}\\right)\\,2^{\\frac{n}{2}}\\,\\varphi^{\\frac{n}{2}}}\\]"
],
"text/plain": [
"(%o80) \n",
" 2 2 2 2\n",
" n - 4 (beta1 - 2 B1 beta1 + beta0 - 2 B0 beta0 + B1 + B0 ) n + S\n",
" ----- - -------------------------------------------------------------\n",
" 2 2 phi\n",
" S n %e\n",
" --------------------------------------------------------------------------\n",
" n - 2 n/2 n/2\n",
" %pi gamma(-----) 2 phi\n",
" 2"
],
"text/x-maxima": [
"(S^((n-4)/2)*n\n",
" *%e^-(((beta1^2-2*B1*beta1+beta0^2-2*B0*beta0+B1^2+B0^2)*n+S)\n",
" /(2*phi)))\n",
" /(%pi*gamma((n-2)/2)*2^(n/2)*phi^(n/2))"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"post6 : linreglik, radcan;"
]
},
{
"cell_type": "code",
"execution_count": 53,
"id": "0edd202d-3451-4d42-846d-c9771636ed8f",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{81}$}-\\frac{\\left(\\beta_{0}-B_{0}\\right)\\,n}{\\varphi}\\]"
],
"text/plain": [
" (beta0 - B0) n\n",
"(%o81) - --------------\n",
" phi"
],
"text/x-maxima": [
"-((beta0-B0)*n)/phi"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{82}$}-\\frac{\\left(\\beta_{1}-B_{1}\\right)\\,n}{\\varphi}\\]"
],
"text/plain": [
" (beta1 - B1) n\n",
"(%o82) - --------------\n",
" phi"
],
"text/x-maxima": [
"-((beta1-B1)*n)/phi"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{83}$}-\\frac{n\\,\\varphi+\\left(-\\beta_{1}^2+2\\,B_{1}\\,\\beta_{1}-\\beta_{0}^2+2\\,B_{0}\\,\\beta_{0}-B_{1}^2-B_{0}^2\\right)\\,n-S}{2\\,\\varphi^2}\\]"
],
"text/plain": [
"(%o83) \n",
" 2 2 2 2\n",
" n phi + ((- beta1 ) + 2 B1 beta1 - beta0 + 2 B0 beta0 - B1 - B0 ) n - S\n",
" - -------------------------------------------------------------------------\n",
" 2\n",
" 2 phi"
],
"text/x-maxima": [
"-(n*phi+((-beta1^2)+2*B1*beta1-beta0^2+2*B0*beta0-B1^2-B0^2)*n-S)/(2*phi^2)"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"postscore60 : ratsimp(diff(log(post6),\\beta0));\n",
"postscore61 : ratsimp(diff(log(post6),\\beta1));\n",
"postscore62 : ratsimp(diff(log(post6),\\phi));"
]
},
{
"cell_type": "code",
"execution_count": 54,
"id": "8f810af9-d141-486a-9fde-f309dee57eed",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{84}$}\\left[ \\left[ \\beta_{0}=B_{0} , \\beta_{1}=B_{1} , \\varphi=\\frac{S}{n} \\right] \\right] \\]"
],
"text/plain": [
" S\n",
"(%o84) [[beta0 = B0, beta1 = B1, phi = -]]\n",
" n"
],
"text/x-maxima": [
"[[beta0 = B0,beta1 = B1,phi = S/n]]"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"solve([postscore60=0,postscore61=0,postscore62=0],[\\beta0,\\beta1,\\phi]);"
]
},
{
"cell_type": "markdown",
"id": "2dfda3ea-9a46-456f-9fac-d24d0237b8d4",
"metadata": {},
"source": [
"For a non-informative prior on the variance, we get:"
]
},
{
"cell_type": "code",
"execution_count": 55,
"id": "e7e9d374-78a0-48d9-8dac-96d7fc1efabb",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{85}$}\\frac{S^{\\frac{n-4}{2}}\\,n\\,e^ {- \\frac{\\left(\\beta_{1}^2-2\\,B_{1}\\,\\beta_{1}+\\beta_{0}^2-2\\,B_{0}\\,\\beta_{0}+B_{1}^2+B_{0}^2\\right)\\,n+S}{2\\,\\varphi} }}{\\pi\\,\\Gamma\\left(\\frac{n-2}{2}\\right)\\,2^{\\frac{n}{2}}\\,\\varphi^{\\frac{n+2}{2}}}\\]"
],
"text/plain": [
"(%o85) \n",
" 2 2 2 2\n",
" n - 4 (beta1 - 2 B1 beta1 + beta0 - 2 B0 beta0 + B1 + B0 ) n + S\n",
" ----- - -------------------------------------------------------------\n",
" 2 2 phi\n",
" S n %e\n",
" --------------------------------------------------------------------------\n",
" n + 2\n",
" -----\n",
" n - 2 n/2 2\n",
" %pi gamma(-----) 2 phi\n",
" 2"
],
"text/x-maxima": [
"(S^((n-4)/2)*n\n",
" *%e^-(((beta1^2-2*B1*beta1+beta0^2-2*B0*beta0+B1^2+B0^2)*n+S)\n",
" /(2*phi)))\n",
" /(%pi*gamma((n-2)/2)*2^(n/2)*phi^((n+2)/2))"
]
},
"execution_count": 55,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"post7 : linreglik/phi, radcan;"
]
},
{
"cell_type": "code",
"execution_count": 56,
"id": "c55df0a5-296b-496e-8f68-561719c48b37",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\[\\tag{${\\it \\%o}_{89}$}\\left[ \\left[ \\beta_{0}=B_{0} , \\beta_{1}=B_{1} , \\varphi=\\frac{S}{n+2} \\right] \\right] \\]"
],
"text/plain": [
" S\n",
"(%o89) [[beta0 = B0, beta1 = B1, phi = -----]]\n",
" n + 2"
],
"text/x-maxima": [
"[[beta0 = B0,beta1 = B1,phi = S/(n+2)]]"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"postscore70 : ratsimp(diff(log(post7),\\beta0))$\n",
"postscore71 : ratsimp(diff(log(post7),\\beta1))$\n",
"postscore72 : ratsimp(diff(log(post7),\\phi))$\n",
"solve([postscore70=0,postscore71=0,postscore72=0],[\\beta0,\\beta1,\\phi]);"
]
},
{
"cell_type": "markdown",
"id": "0f36542a-248a-4626-9fd1-01ba709d9643",
"metadata": {},
"source": [
"Once again, the formula for the variance is identical to that for iid observations from a normal distribution."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Maxima",
"language": "maxima",
"name": "maxima"
},
"language_info": {
"codemirror_mode": "maxima",
"file_extension": ".mac",
"mimetype": "text/x-maxima",
"name": "maxima",
"pygments_lexer": "maxima",
"version": "5.46.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment