1 Introduction: The Purpose and Scope of This Element

The Bayesian philosophy about statistical inference is much older than the better-established traditional Frequentist/Likelihoodist paradigm, starting with Reference BayesBayes (1763) and considered by Laplace, Gauss, and others. Scholars such as Jeffreys, Zellner, Savage, de Finetti, and Lindley reactivated interest in Bayesian methods in the middle of the last century in response to observed deficiencies in classical techniques. Unfortunately many of the specifications developed by these early Bayesians, while superior in theoretical foundation, led to mathematically intractable forms. This problem has been famously solved in recent years by a revolution in statistical computing techniques. Yet regrettably these remarkable developments have not permeated all data-analytic academic fields including the social and behavioral sciences in particular. This may be because of a lack of effective introductory material.

The Bayesian paradigm is ideally suited to the type of data analysis performed by social scientists because it recognizes the mobility of population parameters, incorporates prior knowledge that researchers possess, and updates estimates as new data are observed. Because most empirical work in the social sciences is observational rather than experimental, subjects are rarely cooperative, and systematic effects are often more elusive than in other fields, the Bayesian approach to modeling uncertainty in parameter estimates provides a more robust and realistic picture of the data generating process.

Most people would be surprised to hear that there are different philosophical views within statistics. In fact, this led to a major split in the discipline through most of the twentieth century. A big reason for this schism was that most of the giants of the adolescent age of statistics in the late nineteenth and early twentieth century were openly hostile towards the Bayesian paradigm. Another reason was (notice the past tense) that computing resources were not available in this period for producing useful Bayesian results for the types of models that statisticians and others wanted to specify. Neither of these two issues are important now. It is still useful to understand how the typology of statistics still permeates written and spoken discussions today.

The first generation of modern statisticians were Frequentists from the idea of performing some experiment or test multiple times using a long stream of independent observations from the same data generating source. If a frequentist wanted to determine the probability of heads when presented with a new coin, the procedure would be to flip it into a sand pit hundreds or thousands of times and record the long-term results. This comes from the classic Neyman/Pearson/Wald setup in the late nineteenth century where the orthodox view is sampling is infinite (or as long as desired) and decision rules can be sharp. They also mostly studied physical phenomenon and therefore viewed many target population parameters as fixed by nature where variation occurred from less than perfectly accurate instruments and observers. This idea of replicated tests is a key part of frequentist statistics and led to classic frequentist procedures such as the $1-\alpha$ confidence interval: an interval that over $100(1-\alpha )$ % of replications contains the true value of the parameter on average. To this day the dichotomy between pure frequentism and the regular practice of social science statistics makes this definition confusing to both students and practitioners since it is literally built on the idea of repeating the exact same experiment multiple times, which few of us are fortunate enough to be able to do. Since variability to a frequentist rests in the data and the parameters are assumed fixed points by nature, the key probabilistic quantity of interest is the probability of some function of the data given some hypothesis: $p\left(f\right(\text{data}\left)|H_A\right)$ . The objective is a point estimate and the associated variance of this estimate that informs a formal procedure, setting some $\alpha$ level set in advance in a distributional test. In this setup $H_A$ is accompanied by a complementary hypothesis $H_B$ , meaning that the state of the world is explained by either. This sets up clear definition and calculation of Type I and Type II errors. The frequentist accepts $H_A$ if $p\left(f\right(\text{data}\left)|H_A\right)<\alpha$ and accepts $H_B$ if $p\left(f\right(\text{data}\left)|H_B\right)\ge \alpha$ . That sentence should feel uncomfortable to every practicing empirical social scientist based on the use of the word “accept.” This is the critical difference between actual frequentists and others who do not have a perceived infinite flow of independent identically distributed data and sharp complementary hypotheses. Nearly all social scientists have limited, often one-off, datasets that are contextual in time and space and cannot be replicated 19 more times to achieve a true $0.95$ confidence interval. This why nearly everyone we know in the social sciences who calls themselves a frequentist is not a frequentist.

Why were the leading figures of early statistics openly hostile to the Bayesian approach? Fisher in particular objected to the uniform distribution as a choice of prior distribution, and most of them disliked the idea of inverting the order of conditional probability with Bayes’ Law, which will be discussed at length in the next section. But there were other fundamental philosophical differences that ran counter to frequentist thinking. From the Bayes/Laplace/de Finetti tradition, all unknown quantities are treated probabilistically by assigning them a distribution or probability value, and the state of the world can always be updated by conditioning on new information as it arrives. In the absence of an unending stream of data, the data that do arrive are treated as observed and fixed by the limited sampling process. The most caustic arguments, though, centered around the interpretation of probability. To a Bayesian probability is not from a frequent, long-run, series of experiments but instead based on the “degree of belief” given necessarily limited information before or after a specific dataset arrives. This has been called a “subjective” interpretation of probability but this is a poor description since what the Bayesian interpretation really is based on is what does the currently existing evidence imply about an unknown parameter, say $\theta$ , using probability as the mechanism of description. So the key quantity of interest is $p\left(\theta |f\right(\text{data}\left)\right)$ , which is a distribution. Here $f\left(\text{data}\right)$ is any applicable function of the data from a model or other treatment. Therefore instead of resorting to mechanical hypothesis testing a Bayesian can ask questions such as what is the probability that $\theta$ is greater than zero? Or what is the probability that a specific treatment changes the status of the treated group over the control group? Bayesians construct these inferential statements by taking a prior probability statement about the quantity of interest $p\left(\theta \right)$ and conditioning it on a new set of data to produce an updated version, $p\left(\theta |f\right(\text{data}\left)\right)$ , providing new knowledge if the data are informative.

The overall history of Bayesian statistics as described earlier is an almost soap opera like tale of scientific sociology. The Reverend Thomas Bayes died in 1761 without publishing his article “An Essay towards Solving a Problem in the Doctrine of Chances,” which was then submitted on his behalf by his friend Richard Price (they are both buried in the small Bunhill Cemetery in central London, which you can visit). Some have conjectured that Bayes did not believe in or did not have full confidence in this work, therefore making Bayes not “Bayesian.” This idea was later revisited in the late nineteenth through the middle twentieth centuries in particular by scholars such as Jeffreys, de Finetti, Good, Savage, Lindley, and Zellner. Their work was difficult because many of the giants of the time in statistics such as Pearson (Karl), Neyman and Pearson (Egon), Wald, and, of course, Fisher were openly hostile to the idea of Bayesian inference, particularly obtaining $p\left(\theta |f\right(\text{data}\left)\right)$ using uniformly distributed prior distributions. Fascinating accounts can be found in Reference StiglerStigler (1982), Reference Stigler(1983), and Reference DaleDale (2012). The non-Bayesian charges included: priors are always assigned subjectively, inverting likelihood functions is illogical, and it is easy to specify Bayesian models such that unknown parameter estimation is difficult or impossible. The latter insult was the only one that had an element of truth to it, but this was solved in 1990 through the introduction of Markov chain Monte Carlo to use computational power rather than brute analytical derivation. We now exist in an era where there are no intellectual or analytical impediments to specifying and estimating Bayesian models.

This discussion here so far has left out the bulk of as-practiced social science statistics, which uses limited data collection resources to find evidence for an effect through statistics calculated once in space and time. This is generally from the central challenge in studying humans socially, politically, or biomedically: we are not fixed parameters of nature, and groups of humans interacting with each other are even less so. So because of the instability of human behavior, datasets tend to be limited in generality, and contextual to specific phenomenon being studied at that moment. It would be terrific to go back and get substantially more survey data relevant to the 2016 US presidential election to understand why predictions were so wrong, but we cannot do that even with unlimited funding. Armed with a single point estimate and uncertainty around it for some phenomenon of interest (maximum likelihood estimates as described in Section 3), social scientists are clearly closer to the Bayesian paradigm. In fact, all such models are special cases of Bayesian models with a uniform prior distribution (Section 4), and as the data size gets very large they are identical for any nonpathological choice of the prior distribution (also Section 4). So in fact the overwhelming majority of empirical social scientists are actually just Bayesians who do not yet know it. Hopefully this discussion is motivation to read on.

It is important to provide some pedagogical notes at this beginning point in the Bayesian journey provided here. First, one cannot do meaningful Bayesian inference in the social sciences or elsewhere without some calculus operations. This is perhaps one of the traditional impediments for wide acceptance and use of Bayesian statistics. So there will be some basic calculus in the statistical theory presented in a very gentle way here. It will be as direct and simple as possible with an emphasis on a general understanding of the procedures rather than detailed mathematical expositions. In the twenty-first century practicing empirical social scientists, data scientists, and others performing daily data analysis tasks do not need to do routine calculus calculations because almost invariably the computer does it for them. So the key task is to understand what these operations do rather than master the specific mathematical processes.

There are two primary calculus tools: differentiation and integration. Differentiation takes a function and determines what is the instantaneous rate of change for that function at a specific point. If you are driving a car you have some velocity but in the course of traveling you have changes in that velocity. For instance, after the light turns green you hit the gas and are not yet going very fast but your instantaneous rate of change is very high and positive because you are speeding up. Conversely, if you are cruising at a relatively steady state on the highway your instantaneous rate of change is very low even at high speed. Then as you see a traffic light or obstacle ahead you start to slow down using the break pedal and your instantaneous rate of change becomes large and negative. We often use derivatives to make a statement like this on functions to find points of interest like where this instantaneous rate of change is zero, which is a local maxima or minima of a curvilinear form. Integrals are instead about space. They can measure how much area is there under a curve to the x-axis between two points on that x-axis, like a slice of the curve above zero. This is very useful when dealing with probability functions and general Bayesian inference as we will shortly find out. The key point is to not worry so much about the notation, which some people find odd or intimidating looking, but to think broadly about what is being accomplished with these operations.

Second, we strongly emphasize using computer simulation to produce desired results rather than analytical calculations, which are of course traditional in twentieth-century Bayesian statistics. It turns out that nearly all of the quantities of interest can be produced faster and with less effort by requiring a machine to do repetitive work. This theme is also central to progress in Bayesian statistics that occurred right at the end of the twentieth century and led to enormous leaps in our ability to specify and estimate Bayesian models of interest. We will repeatedly make use of the simulation idea that quantities of interest, like estimated parameters, can be effectively described by generating many samples from the relevant distribution and then treating these values like data and doing simple summaries. This idea goes back to the dawn of computing and is vitally important in all areas of statistics in the twenty-first century.

All code and data used throughout this element are stored in the GitHub repository: https://github.com/jgill22/Bayesian.Social.Science.Statistics, and can also be conveniently executed online through the Code Ocean capsule located at: https://codeocean.com/capsule/8772484.

2 Basic Probability Principles and Bayes Law

Knowing basic probability theory is a requirement for understanding statistics, and even more so Bayesian statistics. Therefore in this section we introduce the required fundamentals and interested readers can consult Reference GillGill (2014) for more details. Probability is nothing more than a standard way to describe uncertainty. When unknown or future events are assigned a small positive number near zero we think of them as being very unlikely to occur. Conversely when these events are assigned a number less than but near one we think of them as being very likely to occur. When we perform the experiment of flipping a fair coin then each event/outcome, heads versus tails, is equally likely to happen so we assign them both the probability of occurrence of 0.5. Implied by this brief discussion is a set of rules that map the occurrence of events to their probability of happening. Start with a set of nonoverlapping discrete events, $A_1,A_2,\dots ,A_k$ that make up a sample space, $S$ , which is the collection of all things that can happen for this experiment. The probabilities associated with these events are denoted $p\left(A_1\right),p\left(A_2\right),\dots ,p\left(A_k\right)$ , which are just functions in the spirit of $f\left(x\right)=x^2$ except we use $p\left(\right)$ instead of $f\left(\right)$ as reminder that these are probability functions. Returning to the example of flipping a fair coin we have $A_1=heads$ , and $A_2=tails$ , and $p\left(A_1\right)=p\left(A_2\right)=0.5$ . First we stipulate that:

• the probability of any realizable event is between zero and one: $p\left(A_i\right)\in [0:1]$ for all of the $A_i$ in $i=1,\dots k$ in the sample space $S$ .

Next we require that one of the events must occur:

• some event happens with probability one: $p\left(S\right)=1$ .

The next rule is a little more tricky and states that if we are considering whether any of a sub-collection of events can occur, say $A_1$ or $A_3$ or $A_7$ , then the associated probability of this set is the sum of their individual probabilities. More formally:

• The probability of unions (collections) of $n$ “pairwise disjoint” (that is, non-overlapping) events is the sum of their individual probabilities: $p\left(\bigcup _{i=1}^n{A}_i\right)=\sum _{i=1}^{n}p\left(A_i\right)$

(where $\bigcup$ means “union”). So $p(A_1\bigcup A_3\bigcup A_7)=p\left(A_1\right)+p\left(A_3\right)+p\left(A_7\right)$ . To further illustrate this last rule return to flipping a fair coin from above: “the probability that we get heads or tails must be one, $p(A_1\bigcup A_2)=p\left(heads\right)+p\left(tails\right)=1$ . In total these rules are called the Kolmogorov probability axioms (Reference KolmogorovKolmogorov, 1933), and they codify the modern definition of a probability function. Before this definition was institutionalized by mathematicians, authors would present their version of these statements at the beginning of the paper and prove or derive some objective based on their version, which meant that readers had to exert extra effort to understand the work than was actually necessary.

The next important principle to understand is conditional probability, which is simply the idea that prior information may change the probability of interest. For example if we are interested the probability of rolling a with a far die, and we know in advance that the outcome is an even number, then sample space reduces from all six outcomes to the set so now instead of $\frac{1}{6}$ . In more formal notation for event $A$ , the unconditional (marginal) probability is $p\left(A\right)$ and the conditional probability given a prior event $B$ is $p\left(A|B\right)$ , “the probability of event $A$ given event $B$ has occurred.” If $B$ is meaningful in the calculation of $p\left(A\right)$ then $p\left(A\right)\ne p\left(A|B\right)$ and we would not want to ignore this information. This is calculated mathematically by:

$p\left(A|B\right)=\frac{p(A\cap B)}{p\left(B\right)},$(2.1)

given that $p\left(B\right)\ne 0$ . Here $\cap$ in the numerator indicates “joint” or “intersection” meaning that the two events both occur, “ $A$ and $B$ .” Henceforth $\cap$ will be replaced with a comma, as in $p(A,B)$ , which is more common in statistical notation for the joint probability of $A$ and $B$ . For the die rolling experiment we know that because the only way both could occur is if we roll a . Since we know that the probability of rolling an even number is $\frac{1}{2}$ , then the calculation of the conditional probability previously is calculated by:

$p\left(A|B\right)=\frac{p(A,B)}{p\left(B\right)}=\frac{1/6}{1/2}=\frac{1}{3}.$

Obviously this principle applies to more elaborate settings as we will see with Bayesian inference. The important idea to remember is that if we have prior information relevant to a probability calculation we would like to perform, then it will affect the resulting probability number when conditioned upon.

Conditional probability is order dependent meaning that the probability of event $A$ given that event $B$ has occurred is not the same as the probability of event $B$ given that event $A$ has occurred:

$p\left(A|B\right)\ne p\left(B|A\right)\text{and therefore from above}\frac{p(A,B)}{p\left(B\right)}\ne \frac{p(B,A)}{p\left(A\right)},$

although

$p(A,B)=p(B,A)$

because order does not matter for joint probability statements just like “and” does not imply an order in English language statements. Given this fact we can write the two conditional probability statements accordingly and rearrange:

$\begin{array}{cccc}p\left(A|B\right)& =\frac{p(A,B)}{p\left(B\right)}& p\left(B|A\right)& =\frac{p(A,B)}{p\left(A\right)}\\ p\left(A|B\right)p\left(B\right)& =p(A,B)& p\left(B|A\right)p\left(A\right)& =p(A,B).\end{array}$

Therefore:

$p\left(A|B\right)p\left(B\right)=p\left(B|A\right)p\left(A\right)⟶p\left(A|B\right)=\frac{p\left(A\right)}{p\left(B\right)}p\left(B|A\right),$(2.2)

meaning that we can switch the order of conditioning by multiplying with the ratio of marginal distributions in the order of the left-hand-side. This is so important that it has a special name, Bayes’ Law from Reference BayesBayes (1763), although we know that Laplace (1811) independently discovered this idea around the same time. In fact, this statement is the core of Bayesian inference in statistics as will be covered in detail in Section 4. The important idea to remember is that order is critical in conditional probability statements, and there is even a name for misunderstanding this principle called the “prosecutor’s fallacy” since courtroom arguments sometimes make this exact mistake.

Consider the following example from the early months of the Covid pandemic as reported in The New York Times (August 6, 2020) for the city of New York that previous spring at the height of their pandemic emergency. The reported probabilities in the city population for having had Covid ( $Cov$ ) and therefore having antibodies, along with test accuracy, are for randomly chosen resident:

• the estimated probability of Covid infection: $p\left(Cov\right)=0.10$ (prevalence)

• the probability of correct positive ( $Pos$ ) classification: $p\left(Pos|Cov\right)=0.875$ (sensitivity)

• the probability of correct negative $\left(Po{s}^{\complement }\right)$ classification: $p\left(Po{s}^{\complement }|Co{v}^{\complement }\right)=0.975$ (specificity), where the “ $\complement$ ” in the exponent denotes the complement of the associated probability, meaning $p\left(even{t}^{\complement }\right)=1-p\left(event\right)$ .

• Now suppose we want $p\left(Cov|Pos\right)$ , from:

  $p\left(Cov|Pos\right)=\frac{p\left(Cov\right)}{p\left(Pos\right)}p\left(Pos|Cov\right),$
  which is the probability of actual Covid infection given a positive test result.

To make this Bayes’ Law calculation we need the marginal $p\left(Pos\right)$ , which we do not immediately have. It can be obtained from the Law of Total Probability, which states that here there are only two ways to test positive (since this example has only two outcomes), testing positive with having Covid and testing positive not having had Covid. Using this fact, the definition of conditional probability, and complementation, we can calculate this marginal:

$\begin{array}{cc}p\left(Pos\right)& =p(Pos,Cov)+p(Pos,Co{v}^{\complement })\\ & \text{[from the Law of Total Probability]}\\ & =p\left(Pos|Cov\right)p\left(Cov\right)+p\left(Pos|Co{v}^{\complement }\right)p\left(Co{v}^{\complement }\right)\\ & \text{[from turning joints into conditionals times marginals]}\\ & =p\left(Pos|Cov\right)p\left(Cov\right)+[1-p(Po{s}^{\complement }|Co{v}^{\complement }\left)\right]p\left(Co{v}^{\complement }\right)\\ & \text{[using complementation in the second term]}\\ & =\left(0.875\right)\left(0.10\right)+(1-0.975)(1-0.10)=\mathrm{0.11.}\end{array}$

Now we have all of the ingredients for the Bayes’ Law calculation:

$p\left(Cov|Pos\right)=\frac{p\left(Cov\right)}{p\left(Pos\right)}p\left(Pos|Cov\right)=\frac{0.10}{0.11}\left(0.875\right)=\mathrm{0.7954545.}$

This means that the probability that this individual has had Covid given a positive test classification is approximately 0.80. From the principle of complementation it also means that about 20% of those with a positive test classification never had the disease. These are not numbers that would make epidemiologists or policy makers particularly happy, but this was an extraordinary time and location for this disease.

Consider using computer simulation to ask probability questions. In traditional basic statistics courses it is common to see exercises that ask for probabilities over sub-regions of the support of a distribution. Most commonly this is tail values of a normal distribution (z-scores), but we can make it as general as we want. The key method introduced here, Monte Carlo simulation dating back to the earliest era of digital computing, is based on generating data according to a known distribution and manipulating the values arithmetically to compute desired quantities. Since computers can easily generate large numbers of specified random variables, the Law of Large Numbers and the Central Limit Theorem are applicable and we can replace analytical human calculations with easier computational work. Suppose we want to calculate the density of a normal distribution with mean 3 and standard deviation 2 that is above zero? This is $p_{normal}(y>0|\mu =3,\sigma =2)$ . A recipe is to generate one million (1M) values from this distribution, count the number of values above zero, and divide that by the number generated. This is three lines of R or Python code (a single line if one wants to be clever), and returns the value $0.933$ . This could be done more directly with specific functions in these languages, and this is an over-simplified use of Monte Carlo simulation. We provide this example here in the adjacent code boxes as a way to introduce simulation based calculation of probabilities with samples, which is a critical and necessary tool in more complex settings for Bayesian inference.

R Code for Simulated Probability Calculations
n.sims <- 1000000
y <- rnorm(n.sims,mean=3,sd=2)
length(y[y>0])/n.sims

Python Code for Simulated Probability Calculations
import numpy as np
n_sims = 1000000
y = np.random.normal(3,2,n_sims)
sum(y>0)/n_sims

The purpose of this section has been to introduce (or review) the necessary basic probability theory for continuing our studies in Bayesian inference. Probability is a deep area of study in mathematics and is pervasive in statistics, data science, and machine learning. Many readers will eventually want to continue studies in this area.

3 What Is a Likelihood Function and Why Care

Many statistics are produced by closed-form mathematical expressions that are well understood in both their calculation and their subsequent properties. The best known form of these is the calculation of the estimate of the vector of regression parameters for a multivariate linear model, which is simply $\hat{\beta }=(X^{\prime }X{)}^{-1}{X}^{\prime }y$ where $X$ is the matrix of explanatory variables down columns with a leading column of 1s and $y$ is the vector of outcomes (bold indicates a vector or matrix structure). This is literally the most studied expression in statistics and we know everything about its properties. Another class of statistics are produced by assuming a distributional property, which leads to a functional form to be numerically maximized in order to produce estimates with known and optimal properties. The latter parametric approach is the topic of this section and a vast number of models in modern statistics need to use this approach to produce useful results for the researcher.

Random variables are defined by their probability mass functions for the discrete case (PMF), or their probability density functions for the continuous case (PDF). These distributions are probability functions, in the sense of the last section, that describe assumed variation in some random variable, say $Y$ over a specific range of values (support), and are typically conditional on parameters, say $\theta$ . For the discrete case we specify $p(Y=y)$ for the probability that the random variable $Y$ takes on some specific realization $y$ . So in the case of flipping a fair coin we say $p(Y=1)=0.5$ and $p(Y=0)=0.5$ , where $1$ indicates heads and $0$ indicates tails (this numerical assignment is standard but arbitrary). In the continuous case we consider $f\left(y\right)$ and often focus intervals on the real line ( $R$ ) since exact numbers do not make sense as specific outcomes since they have probability zero of occurring on the real line. When the random variable $Y$ is conditioned on parameters, that can be known or unknown, then we want to use conditional probability as discussed in the last section. Suppose that the random variable under study, Y, is distributed Exponential: $p\left(\text{Y}|\theta \right)=\theta e^{-\theta \text{Y}}$ (the “rate” version) with support $(0:\infty )$ . This is a conditional probability statement since the distribution of $Y$ is overtly dependent on the value of $\theta$ as in the definitional form in (2.1). Note that this is a PDF since values of $Y$ are defined over the real line from zero to positive infinity, and this makes the exponential distribution useful for modeling time to an event. Several versions with different values of $\theta$ are shown in Figure 1.

Figure 1 The exponential distribution for different $\theta$ values

For the moment consider a “generic” conditional distribution for the random variable $Y$ , $p\left(Y|\theta \right)$ , but we observe an independent sample of $n$ of them from the same data generating process conditioned on the same $\theta$ : $y_1,y_2,\dots ,y_n$ . This is called an “independent, identically distributed” set of data, usually abbreviated IID or iid since they come from the same distribution but are produced without being conditional on each other in any way. Classic data that do not have this property come from a time series where typically serially observed values are conditional on previous observations, and conditions (conditional parameters) are likely to change over time like stock prices or agricultural output.

As a quick illustration, return to the coin flipping example but suppose we flipped the fair coin three times. The observed set of heads and tails is clearly an IID dataset, denoted here as $y=(y_1,y_2,y_3)$ with partial PMF $p(Y=y),y=0,1$ . For just one of these events the probability of an observation is fully described by the PMF. Now ask the question: what is the probability of observing the data $(H,T,H)$ ? It is intuitively given by:

$p(H,T,H)=p(Y=1)\times p(Y=0)\times p(Y=1)=\left(0.5\right)\left(0.5\right)\left(0.5\right)=0.125$

(with standard encoding of heads and tails). We know that it is a fair coin so let us specify a conditional parameter to remind us of this fact: $p(Y=1|\theta =0.5)=0.5$ . Since only two things can happen we also know for a fact that $p(Y=0|\theta =0.5)=0.5$ , but this is redundant information and we only need the first statement to have full information. Returning to the generic description of the observed data, $y=(y_1,y_2,y_3)$ , we can more generally state that:

$p(y|\theta =0.5)=p(y_1|\theta =0.5)p(y_2|\theta =0.5)p(y_3|\theta =0.5)=\prod _{i=1}^{3}p(y_i|\theta =0.5),$

where the product notation ( $\prod$ ) is a convenience to make such statements with arbitrary and possibly large sample sizes manageable notationally. Because it is a fair coin, $\theta =0.5$ , with only two sides, the probability of observing any specific $n$ sized sample will always be equal to $(0.5{)}^n$ , which is oversimplified for our intentions. A slightly more interesting variant would be a “biased” coin with, say, $p(Y=1|\theta =0.75)$ . Now the probability of observing the sequence from above changes to:

$p(H,T,H)=p(Y=1)\times p(Y=0)\times p(Y=1)=\left(0.75\right)\left(0.25\right)\left(0.75\right)=\mathrm{0.140625.}$

Note that this numeric value is higher than before. This is because it is more likely to draw a head and we did draw more heads than tails.

Now consider that we do not know the value of $\theta$ before sampling the data but we are willing to impose a parametric assumption with some form of $p(Y=y|\theta )$ . For the exponentially distributed random variable earlier we specified that $p\left(\text{Y}|\theta \right)=\theta e^{-\theta \text{Y}}$ . In most statistical modeling we make such a parametric assumption at the beginning of the analysis based on previous information, standard practice, or exploratory criteria (although there exist “nonparametric” modeling approaches, also called “semiparametric,” in practice). Of course not knowing the true value of $\theta$ is not only more realistic, it adds substantially to the challenge of analyzing the data and calculating probabilities as we did previously. So now for an arbitrary set of $n$ data values:

$p\left(y|\theta \right)=p\left(y_1|\theta \right)p\left(y_2|\theta \right)\cdots p\left(y_n|\theta \right)=\prod _{i=1}^{n}p\left(y_i|\theta \right).$(3.1)

This would be just as easy to mathematically manipulate as the earlier examples, even for large $n$ , except that we do not know $\theta$ . A typical statistical problem is to observe the sample of $y$ values and use them to estimate some unknown population parameters such as $\theta$ . In fact this is the exact definition of “inference” in statistics. For example, most basic statistics texts start with using a sample of $y$ values to calculate $\bar{y}$ to estimate the true population mean, ${\mu }_y$ . Importantly, once the vector of $y$ values is observed then it is fixed for our purposes. We can go get coffee for ten minutes, or go on vacation for two weeks, and when we return and look on our computer’s hard drive it will still be the same (barring something catastrophic of course). So what is unknown (and possibly random) is $\theta$ . In two of the most important papers in the history of statistics Fisher (Reference FisherFisher, 1922, Reference Fisher1925), considered this situation and decided to first change the notation of (3.1) to reflect that it is the $y$ that is known and the $\theta$ that is target of our inquiry:

$L\left(\theta |y\right)=\prod _{i=1}^{n}p\left(y_i|\theta \right).$(3.2)

This is a notational sleight of hand since we have discussed conditional probability as going in the opposite direction and in the last section showed that order critically matters in conditional probability, but it is expressed this way to emphasize what is known and what is unknown. So now attention is focused on estimating the best value of the unknown $\theta$ given the observed and fixed vector $y$ . So what we want now is the value of $\theta$ , denoted $\hat{\theta }$ (“theta hat”), that is most likely to have generated the observed data given an assumed parametric form. For this reason the “L” in (3.2) stands for “Likelihood” and (3.2) is called a likelihood function. But it is important to remember that technically a likelihood function is just the product of the $n$ PMFs or PDFs of the observed IID data.

Now our job is to find the best $\hat{\theta }$ for (3.2), that is, the value of $y$ that is “most likely to have generated the observed data.” This is called the Maximum Likelihood Estimate (MLE). For of the most common forms of $p\left(y|\theta \right)$ used in practice the form of $L\left(\theta |y\right)$ is concave to the x-axis as shown in Figure 2. This means that the function has a single unique maximum value of $L\left(\theta |y\right)$ for a given observed dataset. Here is where the big intellectual leap occurs. This unique value, $\hat{\theta }$ , is the single value that is most likely to have produced the observed data, indicated by the dotted line in Figure 2. Therefore if the sample data is truly representative of the population data and the parametric assumption is correct, then this is the best estimate of the unknown population value of $\theta$ . In addition, this estimator has been proven to have optimal properties under common circumstances (Reference BirnbaumBirnbaum, 1962). This is a very subtle and deep idea and most readers (at all levels!) do not fully grasp the theoretical importance right away.

Figure 2 A generic likelihood function

Perhaps the best way to add intuition is to work through a detailed and realistic example. Consider a standard model for evaluating data that are counts: nonnegative integers without an upper bound. Another way of thinking of such count data is in terms of durations: the time waiting for some event of interest. If the probability of an event is proportional to the length of the wait, then the number of events in a given time period can be modeled with the Poisson distribution:

$p\left(y|\theta \right)=\frac{e^{-\theta }{\theta }^y}{y!},y\in I^{+},\theta \in R^{+}.$(3.3)

which is read as the “the probability that exactly $y$ events occur given parameter $\theta$ .” Notationally $y\in I^{+}$ indicates that the support of the outcome variable $y$ is over the nonnegative integers, and $\theta \in R^{+}$ indicates that the support of $\theta$ is over the positive real line, $R^{+}$ . The assumption of proportionality is usually quite reasonable because over longer periods of time the event has more “opportunities” to occur. Here $\theta$ is called the intensity parameter and gives the mean (expected) number of events and the (expected) variance. The concept of expectation is described in detail in Section 6. To use the Poisson probability model for counts, we need to assume:

1. 1. Non-Simultaneity: Two events cannot occur at exactly the same time.

2. 2. IID: Events in different time segments are independent and identically distributed.

3. 3. Proportionality: For small time periods, the probability of an event is proportional to the length of time passed in the period so far, and is not dependent on the number of previous events in this period.

This is a basic model for counts and there are many variants when the aforementioned assumptions are violated and for other complexities.

Now assume that we have a vector of counts, $y$ , and we want the MLE of the unknown $\theta$ parameter. The likelihood function is created from the joint distribution of the observed data as shown earlier:

$L\left(\theta |y\right)=\prod _{i=1}^n\frac{e^{-\theta }{\theta }^{y_i}}{y_i!}=\frac{e^{-\theta }{\theta }^{y_1}}{y_1!}\frac{e^{-\theta }{\theta }^{y_2}}{y_2!}\cdots \frac{e^{-\theta }{\theta }^{y_n}}{y_n!}=e^{-n\theta }{\theta }^{\sum y_i}{\left(\prod _{i=1}^n{y}_i!\right)}^{-1}.$

Suppose now that we have the count data: $y=(5,1,1,1,0,0,3,2,3,4)$ , then the likelihood function from plugging in these values previously is:

$L\left(\theta |y\right)=\frac{e^{-10\theta }{\theta }^{20}}{207360},$

which is the probability of observing this exact sample. With likelihood function calculations is often easier to deal the logarithm:

$\begin{array}{c}\log L\left(\theta \right|y)=\ell (\theta \left|y\right)=log\left(e^{-n\theta }{\theta }^{\sum yi}{\left(\prod _{i-1}^{n}yi!\right)}^{-1}\right)\\ =-n\theta +\sum _{i-1}^{n}yi\log \left(\theta \right)-\log \left(\prod _{i=1}^{n}yi!\right)\end{array}$(3.4)

where the “ $logL$ ” and “ $\ell$ ” notation are both commonly used. For our small example this is numerically:

$\ell \left(\theta |y\right)=-10\theta +20log\left(\theta \right)-\underset{12.242}{\underbrace{log(207360{)}_{}}}.$(3.5)

Importantly, for the family of functions that we will use the likelihood function and the log-likelihood function have the same mode (maximum of the function) for $\theta$ , and they are both guaranteed to be concave to the x-axis. This is illustrated for this example in Figure 3. The property of having only one mode and no minima is important because our calculus tool (taking the derivative, setting equal to zero) finds points with a minima or maxima without regard for which, so in this case we know what we are finding.

Figure 3 Poisson likelihood function and log-likelihood function

So now let’s use freshman calculus to find out the maximum of the log likelihood function. This is at the $\theta$ point when first derivative of $\ell \left(\theta |y\right)$ equals zero, where the dotted line in Figure 2 touches the curve. This means that we take the first derivative with regard to the parameter of interest denoted $\frac{d}{d\theta }$ , set the resulting equation equal to zero, and solve for $\theta$ . Here the statement $\frac{d}{d\theta }\ell \left(\theta |y\right)\equiv 0$ is unsurprisingly called the Likelihood Equation. For the numerical example we start with $\ell \left(\theta |y\right)$ and take the derivative (which uses the exponent rule in calculus: $\frac{d}{dx}\left(\text{some constant}\right)x^k=\left(\text{some constant}\right)k{x}^{k-1}$ ), and setting equal to zero gives:

$\frac{d}{d\theta }\ell \left(\theta |y\right)=\frac{d}{d\theta }(-10\theta +20log(\theta )-12.242)=-10+20{\theta }^{-1}\equiv 0,$(3.6)

where we use the additional rules that $\frac{d}{dx}\left(\text{some constant alone}\right)=0$ , and $\frac{d}{dx}log\left(x\right)=x^{-1}$ . So that $20{\theta }^{-1}=10$ , and therefore the MLE is $\hat{\theta }=2$ (note the hat indicating that this is now an estimate). This is the most likely $\theta$ value with the Poisson PMF to have generated the data $y=(5,1,1,1,0,0,3,2,3,4)$ . Note that the mathematical expression $\equiv 0$ means “set equal to zero,” as opposed to $=0$ which means “is equal to zero.”

In the more general sense with arbitrary data we can perform the same process to get a generic form of the Poisson MLE:

$\begin{array}{c}\ell \left(\theta \right|y)=-n\theta +\sum _{i=1}^{n}yi\log \left(\theta \right)-\log \left(\prod _{i=1}^{n}yi!\right)\\ \left[takethefirstderivativewithrespectto\theta \right]\\ \frac{d}{d\theta }\ell \left(\theta \right|y)=-n+\frac{1}{\theta }\sum _{i=1}^{n}yi-0\equiv 0\\ \left[\text{algebraically rearrange}\right]\\ \hat{\theta }=\frac{1}{n}\sum _{i=1}^{n}yi=\overline{y}\end{array}$

Note that it is not true that the MLE is always the data mean; this is just a special property of the Poisson PMF. So to summarize the general process, perform the following steps:

1. 1. identify the PMF or PDF applicable to the outcome

2. 2. create the likelihood function from the joint distribution of the observed data

3. 3. change to the log-likelihood for convenience

4. 4. take the first derivative with respect to the parameter of interest

5. 5. set this equation equal to zero

6. 6. solve algebraically for the MLE.

Of course with more complicated forms this process is done in software such as with generalized linear models regression (logit, probit, gamma, etc.) through a process called Iteratively Weighted Least Squares; see Reference Gill and TorresGill and Torres (2019).

R Code to for the Poisson MLE Example
# A POISSON LIKELIHOOD AND LOG-LIKELIHOOD FUNCTION
llhfunc<-function(X,p,do.log=TRUE) {
    d <- rep(X,length(p))
    q.vec <- rep(length(X),length(p))
    p.vec <- rep(p,q.vec)
    d.mat <- matrix(dpois(d,p.vec,log=do.log),
        ncol=length(p))
    if (do.log==TRUE) apply(d.mat,2,sum)
    else apply(d.mat,2,prod)
}

y.vals<-c(1,3,1,5,2,6,8,11,0,0)

# EXAMPLE RUN FOR TWO POSSIBLE VALUES OF THETA: 4 AND 30
llhfunc(y.vals,c(4,30))

# USE THE R CORE FUNCTION FOR OPTIMIZING,
# par=STARTING VALUES,
# control=list(fnscale=-1) INDICATES A MAXIMIZATION,
# bfgs=QUASI-NEWTON ALGORITHM
mle <- optim(par=1,fn=llhfunc,X=y.vals,
    control=list(fnscale=-1),method="BFGS")

# MAKE A PRETTY GRAPH OF THE LOG AND NON-LOG VERSIONS
ruler <- seq(from=.01, to=20, by= .01)
poisson.ll <- llhfunc(y.vals,ruler)
poisson.l <- llhfunc(y.vals,ruler,do.log=FALSE)

par(oma=c(3,3,1,1),mar=c(0,0,0,0),mfrow=c(2,1))
plot(ruler,poison.l,type="l",xaxt="n",lwd=3)
text(mean(ruler),mean(poison.l),
    "Poisson Likelihood Function")
plot(ruler,poison.ll,,type="l",lwd=3)
text(mean(ruler)+5,mean(poison.ll)/2,
    "Poisson Log-Likelihood Function")

Python Code to for the Poisson MLE Example
import numpy as np
from scipy.stats import poisson
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# POISSON LIKELIHOOD AND LOG_LIKELIHOOD FUNCTION
def llhfunc(X, p, do_log=True):
    lX, lp = len(X), len(p)
    d = np.tile(X, lp)
    u = np.repeat(p, lX)
    p_pmf = [poisson.pmf(d[i:i+lX], u[i:i+lX])
             for i in range(0, lX*lp, lX)]
    d_mat = np.log(p_pmf).T if do_log else \
            np.array(p_pmf).T
    return np.sum(d_mat, axis=0) if do_log else \
           np.prod(d_mat, axis=0)

y_vals = np.array([1, 3, 1, 5, 2, 6, 8, 11, 0, 0])
# EXAMPLE RUN FOR TWO POSSIBLE VALUES OF THETA: 4 AND 30
llhfunc(y_vals,np.array([4, 30]))

# USE minimize() FROM scipy.optimize
mle = minimize(fun=lambda p: -llhfunc(y_vals, p),
               x0=1,
               method="BFGS")
# fun=lambda p: -llhfunc(y_vals, p): DEFINE FUNCTION
# x0=1: STARTING VALUES
# method="BFGS": QUASI-NEWTON ALGORITHM

ruler = np.arange(.01, 20.01, .01)
poison_ll = llhfunc(y_vals, ruler)
poison_l = llhfunc(y_vals, ruler, do_log=False)

fig, axs = plt.subplots(2, 1, figsize=(8, 12))
axs[0].plot(ruler, poison_l, linewidth=3)
axs[0].annotate('Poisson Likelihood Function',
                xy=(np.mean(ruler),
                    np.mean(poison_l)),
                xytext=(np.mean(ruler)-2,
                        np.mean(poison_l)))
axs[0].tick_params(which='both',
                   bottom=False,
                   labelbottom=False)
axs[1].plot(ruler, poison_ll, linewidth=3)
axs[1].set_xlabel('Support of $\Theta$')
axs[1].annotate('Poisson Log-Likelihood Function',
                xy=(np.mean(ruler),
                    np.mean(poison_ll)),
                xytext=(np.mean(ruler)+2,
                        np.mean(poison_ll)+20))
plt.subplots_adjust(hspace=0)
plt.show()

We can also measure the uncertainty of the MLE around this point. Recall that all reported statistics should be accompanied by an associated measure of uncertainty, usually a standard error. To obtain this measurement we will return to basic derivative calculus. When working with functions the first derivative measures slope of the tangent line at given points (zero was of interest in our case before) and the second derivative measures “curvature” of the function at a given point. The second derivative is obtained by taking a function’s first derivative and performing the same operation again on that: taking $\frac{d}{d\theta }\ell \left(\theta |y\right)$ and repeating the differentiation process, $\frac{d}{d\theta }\left(\frac{d}{d\theta }\ell \right(\theta |y\left)\right)$ . Getting this curvature is important because the more peaked the log-likelihood function is at the MLE, the more “certain” the data are about this estimator. Conversely the more diffuse the log-likelihood function around the MLE, the less the data are saying about the this estimate.

Going from curvature to a standard error of the MLE is relatively simple. The square root of the negative inverse of the expected value of the second derivative is the standard error of the MLE. We will save the details of expected value for Section 6, but the operation is very simple here. Calculate the first derivative:

$\frac{d}{d\theta }\ell \left(\theta |y\right)=-n+\frac{1}{\theta }\sum _{i=1}^n{x}_i$

and then take another (now second) derivative of this form:

$\frac{d^2}{d{\theta }^2}\ell \left(\theta |y\right)=\frac{d}{d\theta }\left(\frac{d}{d\theta }\ell \left(\theta |y\right)\right)=-{\theta }^{-2}\sum _{i=1}^n{x}_i$

again using the exponential rule in calculus. The expected value (estimate) of $\theta$ is the MLE (data mean for the Poisson case) allowing us to replace $\hat{\theta }$ with $\bar{y}$ , so:

$\text{Var}\left(\hat{\theta }\right)={\left({\hat{\theta }}^{-2}\sum _{i=1}^n{x}_i\right)}^{-1}=\frac{{\hat{\theta }}^2}{\sum _{i=1}^n{x}_i}=\frac{{\bar{y}}^2}{n\bar{y}}=\frac{\bar{y}}{n}.$

Now we have the MLE for the general Poisson case, $\bar{y}$ , and its associated standard error, $(y/n{)}^{-1/2}$ . In our numerical example the standard error is then simply $\sqrt{2/10}=0.4472$ .

4 The Core of Bayesian Inference: Prior Times Likelihood

So far everything in this exposition has been leading up to this section, which is easily the most important one. How are Bayesian models constructed, and how is Bayesian inference performed? Here we will look at the way prior information and data are combined to provide inferences that are a balance between the two.

We start with broad high-level language giving the philosophical principles of Bayesian inference before proceeding to the mechanics and details of the process. The first, and perhaps most important, principle is that all unknowns are treated probabilistically, meaning they are assigned a conditional distributional statement $p(\cdot |\cdot )$ , or a probability value $[0:1]$ . This includes unknown parameters, missing data values, and even model choices. Thus uncertainty is always a probabilistic phenomenon. Specifically data enters the inference using probability models in the form of likelihood functions, which, recall are just the joint probability of the observed data using an assumed parametric form. Then all parameters (variables of interest or requirement) are assigned distributions with all relevant information prior to considering the observed data at hand: PMFs or PDFs depending on the level of measurement.

The Bayesian inferential process proceeds by using prior information combined data information to produce an updated distribution for each parameter of interest with the principle of inverse probability through conditional distributions presented in Section 2 as Bayes’ Law. This provides a full probabilistic description that can be evaluated in many ways. The overall process can be then described in three steps:

1. 1. Specify a probability model for unknown parameter values that includes some prior knowledge about the parameters if available.

2. 2. Update knowledge about the unknown parameters by conditioning this probability model on observed data.

3. 3. Evaluate the fit of the model to the data and the sensitivity of the conclusions to the assumptions.

We will, of course, go over these steps in detail in the process of this monograph. For now readers should have an appreciation for the considerable difference between Bayesian inference and the maximum likelihood process for estimation described in Section 3, even though both use a likelihood function.

4.1 The Core Bayesian Process

Bayesian inference is about updating. The core operation is taking a distribution that represents the current state of knowledge about some phenomenon of interest and updating it with recently acquired information to produce a new distribution that is more informed. Our starting point is (2.2) in Section 2, which stated $p\left(A|B\right)=\frac{p\left(A\right)}{p\left(B\right)}p\left(B|A\right)$ , and can be thought of as updating $A$ with new information $B$ from a version of $A$ with no conditioning, $P\left(A\right)$ . Now replace this $A$ with $\theta$ as some model parameter that we would like to estimate, and replace $B$ with a more direct expression of the data, $y$ . Therefore $p\left(B|A\right)$ is a joint distribution of the data conditional on the parameter, the likelihood function: $L\left(\theta |y\right)$ . Now Bayes’ Law can be expressed in the real quantities of interest:

$\pi \left(\theta |y\right)=\frac{p\left(\theta \right)L\left(\theta |y\right)}{p\left(y\right)},$(4.1)

where the $\pi \left(\right)$ notation is just a reminder that the left-hand-side expression is different than the $p\left(\theta \right)$ on the right-hand-side. So far this is very intuitive, but how do we interpret $p\left(y\right)$ ? This is the unconditional probability of observing the data that we have in hand right now. Remember to a Bayesian once the data are observed they are fixed forever and no longer have a random (stochastic) characteristic, as opposed to a frequentist who has an unending stream of random data. So in a Bayesian sense then $p\left(y\right)=1$ . So why is it even there if it is just an unimportant constant? One answer is that it is a constant that ensures that $\pi \left(\theta |y\right)$ sums or integrates to one so that it is a proper probability statement in the sense of the Kolmogorov probability axioms discussed in Section 2. So specifically, what would happen if we left it out is that $\pi \left(\theta |y\right)$ would not sum (PMF) or integrate (PDF) to $1$ as a proper probability function should. If this denominator is actually not the number $1$ , then what is it? It is whatever number that in divisor makes the right-hand-side sum or integrate to one. It is easier for now to consider only the continuous case where all of the functions in (4.1) above are PDFs. If we want $p\left(\theta \right)L\left(\theta |y\right)$ to integrate to $1$ then we need to divide it by the area it occupies under the curve over the support of $\theta$ on the x-axis labeled $\Theta$ as shown hypothetically in Figure 4. Remember $y$ is fixed here so only $\theta$ has a distributional property. This area is calculated by the integral:

Figure 4 Integrated likelihood

${\int }_{\Theta }p\left(\theta \right)L\left(\theta |y\right)d\theta =f\left(y\right),$(4.2)

where: $\int$ is a curvy version of $\sum$ for an interval measured variable, $\Theta$ is the support of the variable $\theta$ on the x-axis, and $d\theta$ is a reminder that this calculation is done with respect to the variable $\theta$ (“the variable of integration”). Some of this notation feels like redundant information but it is a hint that integrals can be much more complex and nuanced with many variables and different types of support (Section 6). So if $p\left(\theta \right)L\left(\theta |y\right)$ integrated to $1.4$ , say, we would just divide it by $1.4$ to make $\pi \left(\theta |y\right)$ a proper probability statement. So actually performing the integral (typically done computationally, if at all) would produce a number that we would use accordingly.

All of this discussion of the denominator leads to two forms of the key Bayesian inferential expression. First there is:

$\pi \left(\theta |y\right)=\frac{p\left(\theta \right)L\left(\theta |y\right)}{{\int }_{\Theta }p\left(\theta \right)L\left(\theta |y\right)d\theta },$(4.3)

which is the primary definitional form that says that the posterior distribution is obtained with Bayes’ Law from the prior distribution times the likelihood function with a scaling factor in the denominator. So we have a balance or compromise of prior information and data information that become posterior information. Now that we know that the denominator’s main purpose is simply to make the posterior distribution a proper one and that this purely numeric information can be recovered at any time, it is convenient to drop it and use a proportional statement ( $\propto$ ) for the relationship:

$\pi \left(\theta |y\right)\propto p\left(\theta \right)L\left(\theta |y\right),$(4.4)

which is read as:

$\text{Posterior Probability}\propto \text{Prior Probability}\times \text{Likelihood Function}.$

This is the form that you see most often because it contains the core of the relationship and simply omits an inconvenience that can be recovered later after key calculations. The omitted denominator, called the “integrated likelihood,” the “marginal likelihood,” the “marginal probability of the data,” or the “predictive probability of the data,” seems unimportant here but it is actually very useful in other contexts, including the model comparison tool Bayes Factor.

From all of this we see the key fundamental principle of Bayesian inference: start with a prior distribution, $p\left(\theta \right)$ , that is unconditional on the data at hand and reflects prior belief in the distribution stipulated, multiply it with the likelihood function, $L\left(\theta |y\right)$ , where the data enters the equation, to produce the posterior distribution, $\pi \left(\theta |y\right)$ , which is a distribution that reflects our latest knowledge about the phenomenon of interest. Notice that the result is a distribution not a point estimate, meaning that uncertainty is built-in to the results as shown by how peaked or flattened this posterior distribution appears. One powerful feature of this process is that in the future if additional data ( $y_2$ ) are observed then we can treat the current posterior based on the earlier data ( $y_1$ ) as a prior distribution in (4.4), and use the new data in a new likelihood function, and create an updated posterior to reflect the latest knowledge about $\theta$ :

$\pi \left(\theta |y_2\right)\propto \pi \left(\theta |y_1\right)L\left(\theta |y_2\right).$(4.5)

Plus, the final posterior created in this process is exactly the same as if we had both datasets ( $y_1,y_2$ ) at the same time. This updating process can also be repeated as many times as data arrival permits with the same principle applying.

4.2 Mathematical Example of Posterior Calculation

To show how this process works with a specific model and actual social science data, return to the Poisson model from Section 2 where we developed the likelihood function for count data using a Poisson specification to provide inference for an unknown intensity parameter $\theta$ . Recall that the derived likelihood function for this setup with an $n$ -length sample $y$ was:

$L\left(\theta |y\right)=exp[-n\theta ]{\theta }^{\sum y_i}{\left(\prod _{i=1}^n{y}_i!\right)}^{-1}.$(4.6)

A convenient and flexible prior distribution for $\theta$ with the right support is the gamma distribution:

$p(\theta |\alpha ,\beta )={\beta }^{\alpha }\Gamma (\alpha {)}^{-1}{\theta }^{\alpha -1}exp[-\beta \theta ],\alpha ,\beta ,\theta >0$(4.7)

given a shape parameter $\alpha$ , and rate parameter $\beta$ (the latter alternately expressed as a scale parameter $1/\beta$ ). Also, the Gamma function ( $\Gamma \left(\right)$ earlier is the generalization of the factorial function ( $k!=k\times (k-1)\times (k-2)\cdots 2\times 1$ ) that can be applied to noninteger values( $\Gamma \left(k\right)={\int }_0^{\infty }{t}^{k-1}{e}^{-t}dt,k>0$ ). The exponential distribution shown in Figure 1 is actually a special case of the gamma distribution where the $\alpha$ parameter is set to one. As discussed in Section 4, it is the researcher’s responsibility to set the two “hyperparameter” values ( $\alpha$ and $\beta$ ) according to theory, convention, or desired vagueness. The posterior distribution for $\theta$ is given by the calculation starting from (4.4):

$\begin{array}{cc}\pi \left(\theta |y\right)& \propto p\left(\theta \right)L\left(\theta |y\right)\\ & \text{[plug (4.7) and (4.6) in that order]}\\ & ={\beta }^{\alpha }\Gamma (\alpha {)}^{-1}{\theta }^{\alpha -1}exp[-\beta \theta ]\times exp[-n\theta ]{\theta }^{\sum y_i}{\left(\prod _{i=1}^n{y}_i!\right)}^{-1}\\ & \text{[strip off all constants using proportionality]}\\ & \propto {\theta }^{\alpha -1+\sum _{i=1}^n{y}_i}exp[-(\beta +n\left)\theta \right].\end{array}$(4.8)

The second step is justified because we are already using proportionality by ignoring the denominator of Bayes’ Law, and therefore these constants add no important information. They will be recovered in a later step to make posterior statements “proper” in the sense that $\pi \left(\theta |y\right)$ integrates to 1 as probability function should. From the second step we can also see that ignoring the constants for the moment makes the form of the posterior more clear and less cluttered. Notationally these Bayesian expressions are often given without all of the fixed parameters on the side of the conditional if the dependency is obvious as it is with (4.8). This expression can be more completely described as $p(\theta |y,\alpha ,\beta )$ on the left side of $\propto$ , but the conditionality is clear without this. In more complex specifications listing all of these hyperparameter and fixed parameter values often muddies the notation.

To further improve our intuition we can define:

${\alpha }^{†}=\alpha +\sum _{i=1}^n{y}_i{\beta }^{†}=\beta +n,$

which means that we can express the posterior distribution in the last line of (4.8) as:

$\pi \left(\theta |y\right)\propto {\theta }^{{\alpha }^{†}}exp[-({\beta }^{†}\left)\theta \right],$

which is the kernel (the part of a distribution absent constants) of a gamma distribution different than the gamma prior given by (4.7), which means that if use we proportionality in the other direction (adding back the constants) we get:

$\pi \left(\theta |y\right)=\left({\beta }^{†}{)}^{{\alpha }^{†}}\Gamma \right({\alpha }^{†}{)}^{-1}{\theta }^{{\alpha }^{†}-1}exp[-{\beta }^{†}\theta ].$(4.9)

This means that the posterior distribution of $\theta$ is a gamma form with parameters ${\alpha }^{†}$ and ${\beta }^{†}$ previously, and that we now know everything about it. So we can report its mean, median, variance, quantile variances, and more. For instance the mean and variance of gamma distributed random variable $X$ are:

• $E\left[X\right]=\frac{\alpha }{\beta }$ , rate version.

• $\text{Var}\left[X\right]=\frac{\alpha }{{\beta }^2}$ , rate version.

The first expression earlier for the expected value will be derived in Section 6 as an example of an expected value calculation. It is also now a trivial exercise to plot it in R or Python. This is a particularly elegant result partly because it is a gamma to gamma conjugate relationship through the Poisson likelihood function. Conjugacy is a property of a set of specific Bayesian models where pairing the distributional form of the prior distribution with a specific likelihood function means that the posterior distribution has the same form as the prior, with different parameterization of course.

We can also think about the expected value (mean) of this posterior gamma distribution for the parameter $\theta$ in an informed way with some modest algebra using the aforementioned definition:

$E\left[\theta |y\right]=\overline{\theta |y}=\frac{{\alpha }^{†}}{{\beta }^{†}}=\frac{\alpha +\sum _{i=1}^n{y}_i}{\beta +n}=\left(\frac{n}{\beta +n}\right)\bar{y}+\left(\frac{\beta }{\beta +n}\right)\left(\frac{\alpha }{\beta }\right).$(4.10)

The concept of expected value will be discussed at length in Section 6, but for now think of it as simply the mean of the distribution of the random variable $\theta$ given the data $y$ . This rewrite of the posterior mean into a weighed sum of the data mean and the prior mean shows some important principles here and is present in all Bayesian inference. Consider earlier what happens as $n$ gets very large moving towards infinity. The second term in the sum goes away because $n$ is only in the denominator of the bracketed component (the second weight), and in the first term the bracketed component (the first weight) goes to one for any reasonable choice of the number $\beta$ . So asymptotically (in the limit as the data size gets bigger) the posterior mean for $\theta$ converges to the data mean, $\bar{y}$ , contained in the first term, and the prior mean, $\alpha /\beta$ , in the second term becomes irrelevant. So the data wins over the prior in the limit, and this is true in every Bayesian specification no matter how complex. It also makes intuitive sense: if we have massive amounts of (presumably high quality) data then prior information should matter less or not at all. Conversely, when the data size is limited to a modest number we want to be able to rely on high quality prior information. Furthermore, this is exactly how science is supposed to work: a large quantity of reliable new information should change our knowledge about some phenomenon of interest relative to prior information, and modest amounts of new information may not do so.

4.3 Empirical Example of Posterior Calculation

Consider describing active shooter incidents by year in the United States from 2000 to 2021. The report cautiously notes that the definition of “mass shooter incident” is not fully agreed upon and therefore uses instead the FBI definition of active shooter incidents: “one or more individuals actively engaged in the killing or attempting to kill people in a populated area” to produce the count data reproduced in Table 1 with $n=22$ , and $\sum y_i=429$ . It turns out that now there is very little calculation to be done since we have a full recipe for the posterior of interest from (4.9) with the parameter definitions in (4.2).

Table 1 Active shooter incidents in the United States

Year	2000	2001	2002	2003	2004	2005	2006	2007
Count	3	10	7	12	5	11	12	14
Year	2008	2009	2010	2011	2012	2013	2014	2015
Count	9	19	27	13	20	18	19	19
Year	2016	2017	2018	2019	2020	2021
Count	19	31	30	30	40	61

Note: Counts by Year in the United States.

Source: Pew Research Center, April 6, 2023, from FBI statistics.

Our first job is to determine a reasonable prior form for $\theta$ , which means selecting the two parameters in (4.7). A reasonably diffuse and conservative choice is $\alpha =14$ and $\beta =2$ for a prior mean and variance of $7$ and $3.5$ , respectively, reflecting a much lower level of active shooter incidents before the year 2000. So stipulating this prior is an implicit Bayesian hypothesis test that these incidents have increased in the last two decades over previous eras. This prior decision leads to the simple conjugate update described in (4.9) on 25 to produce the gamma distribution posterior for $\theta$ :

$G_{\text{prior}}(14,2)⟶G_{\text{posterior}}(14+429,2+22)=G(443,24),$

meaning that the updated posterior mean and variance of $\theta$ are $18.46$ and $0.77$ respectively. This progression from prior through data to posterior is shown in Figure 5. The effect of the data on the prior is not only to move the center of the distribution but also to decrease the variability indicating that the data are informative in addition to what is contained in the prior. We can also compare quantiles using the R function qgamma or the Python function gdtr, as shown in Table 2.

Figure 5 Gamma Poisson model prior to posterior update

Table 2 Quantiles from the prior and posterior, active shooter model

Quantile	0.05	0.25	0.50	0.75	0.95
Prior Distribution	4.232	5.664	6.834	8.155	10.334
Posterior Distribution	17.040	17.859	18.444	19.042	19.924

R Code for the Gamma Poisson Update
y <- c(3,10,7,12,5,11,12,14,9,19,27,13,20,18,19,19,19,
       31, 30,30,40,61)
par(mar=c(6,6,2,2), cex.lab=1.5)
ruler <- seq(from=0,to=25,length=300)
alpha <-14; beta <- 2
plot(ruler,dgamma(ruler,alpha,beta),type="l",
    ylim=c(0,0.5),lwd=3, col="grey70",
    ylab="Posterior Density",xlab="Support")
lines(ruler,dgamma(ruler,alpha+sum(y),beta+length(y)),
    lwd=4)
text(7,0.25,"Prior",col="grey70",cex=1.5,adj=0.5)
text(18.46,0.485,"Posterior",col="grey10",cex=1.5,
    adj=0.5)

Python Code for the Gamma Poisson Update
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import gamma

y = np.array([3,10,7,12,5,11,12,14,9,19,27,13,
              20,18,19,19,19,31,30,30,40,61])
ruler = np.linspace(0, 25, num=300)
alpha = 14; beta = 2

# GAMMA DISTRIBUTION IN scipy.stats USES SCALE FOR
# THE SECOND PARAMETER: 1/rate
d_val_prior = gamma.pdf(ruler, alpha, scale=1/beta)
d_val_post = gamma.pdf(ruler, alpha+sum(y),
                       scale=1/(beta+len(y)))
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(ruler, d_val_prior, color='grey',
        linewidth=3, label="Prior")
ax.plot(ruler, d_val_post, color='black',
        linewidth=4, label="Posterior")
ax.set_xlabel('Support')
ax.set_ylabel('Density')
ax.text(6.5, 0.235, 'Prior', color='grey',
        fontsize=12, ha='center')
ax.text(18.46, 0.475, 'Posterior', color='black',
        fontsize=12, ha='center')
plt.ylim([0, 0.51])
plt.show()

5 Prior Probabilities and the Progression of Human Knowledge

Prior to posterior inference is how scientific knowledge accumulates, and prior information abounds in the social sciences. We have hundreds of years of self-aware study of the human condition and human interactions, both qualitative and quantitatively. In fact it seems daft to ignore deep literature, common wisdom, recorded history, and researcher intuition when creating statistical models to understand social science phenomena. Prior specifications can come from previous studies, conflicting theories in a given literature, researcher experience, nonparametrics and other data-oriented sources, diagnostic objectives, expert judgments, and even previous posterior distributions.

As discussed, the Bayesian process of inference starts with assigning prior distributions for unknown parameters. Prior distributions are necessary in Bayesian models. These unknown parameter distributions are operationalized with observed explanatory variables in a simple model. Such prior distributions range from very informative descriptions based on previous research in the field to purposefully vague and uncertain forms that reflect high levels uncertainty or possibly previous ignorance. It is important to notice that the prior distribution is not an inconvenience imposed on the researcher by the treatment of unknown quantities. It is instead an opportunity to include existing knowledge systematically in the model. Such prior information can include: qualitative, narrative, statistical, and intuitive information.

The use of prior distributions has been controversial at various times in the history of statistics. This is partly because guidance on the selection of priors is less firm than with other parts of the model. However, there are often strong arguments for particular forms of the prior: little or vague knowledge often justifies a diffuse or even a uniform prior, certain probability models logically lead to particular forms of the prior (conjugacy), and the prior allows researchers to include additional information collected outside the current study. This controversy is also because early leading statisticians were opposed to the use of “subjective” information in specifications, thus this word became a derogatory term for Bayesian practitioners. However, it is a completely misguided charge since all decisions made in the process of model development, estimation, and reporting are subjective. It is subjective which data are used; it is subjective which variables are selected; it is subjective which model is selected; it is subjective which likelihood function is specified; it is subjective which significance level is selected (note that the standard levels are not theoretically derived, they are merely “conveniences”); it is subjective which software is used; it is subjective which estimation procedure is specified; it is subjective how the results are presented; and it is subjective where and how findings are distributed. So every statistical model ever developed is “subjective.”

There are also technical reasons to allay concerns about the use of prior distributions in statistical models. The maximum likelihood estimate (Section 3) is equal to the Bayesian posterior mode with the appropriately bounded uniform prior. This means that a uniform prior can be found to equate the point estimates, although Bayesians generally prefer to give more information than that through distributional descriptions. A likelihood based model with no prior and a corresponding Bayesian model are asymptotically equivalent for any nonpathological choice of prior: as the sample size increases, the likelihood progressively dominates the prior, as demonstrated for one particular case in Section 4.2. Because of the Central Limit Theorem, both estimates are also identically normally distributed for large data. These properties are explored in technical detail in Reference Diaconis and FreedmanDiaconis and Freedman (1986).

5.1 Conjugate Model Specifications

The phenomenon of conjugacy happens when the distributional form of the prior distribution flows down to the posterior distribution given a specific likelihood function:

${\pi }_f\left(\theta |y\right)=p_f\left(\theta \right)L_{\mapsto f}\left(\theta |y\right),$(5.1)

where “ $\mapsto$ ” means “maps to” and the “ $f$ ” designation means that $\pi \left(\right)$ and $p\left(\right)$ are the same PDF or PMF but with different parameterizations due to the influence of $L\left(\right)$ . We have already seen one example of a conjugate model setup in Section 4 where a gamma distribution prior and a Poisson likelihood produced a gamma distribution posterior distribution. Conjugacy is not a requirement of Bayesian specifications, but it is especially important in the history of Bayesian statistics because before the advent of advanced computational tools in the 1990s it was a mathematically guaranteed way to be able to produce tractable posterior forms. Up until the 1990s it was easy to specify a desired Bayesian (usually regression) model where the joint posterior distribution could not be integrated to produce marginal distributions for a results summary. The most common conjugate relationships are given in Table 3. See Reference GillGill (2014), Appendix B for mathematical details on these forms. An excellent book-long exposition that derives these conjugate models and related forms is Reference ZellnerZellner (1996).

Table 3 Conjugate prior-likelihood pairings for Bayesian models

Likelihood Function Form	Prior/Posterior Distribution
Bernoulli	Beta
Binomial	Beta
Multinomial	Dirichlet
Negative Binomial	Beta
Poisson	Gamma
Gamma (including ${\chi }^2$ ,
and Exponential)	Gamma
Normal for $\mu$	Normal
Normal for ${\sigma }^2$	Inverse Gamma
Pareto for $\alpha$	Gamma
Pareto for $\beta$	Pareto
Uniform	Pareto

An easy, obvious, and convenient conjugate setup is the normal-normal for the mean parameter and the normal-inverse gamma for the variance parameter. This is also important because of the popularity of normal specifications in standard practice and the effect of the Central Limit Theorem. Suppose that we have $n$ -length $y$ , a vector of IID normally distributed data with population mean $-\infty <\mu <\infty$ and population variance ${\sigma }^2>0$ . Then the two-variable normal likelihood function for these data is:

$L(\mu ,{\sigma }^2|y)=(2\pi {\sigma }^2{)}^{-\frac{n}{2}}\prod _{i=1}^n\text{exp}\left(-\frac{1}{2{\sigma }^2}(y_i-\mu {)}^2\right),$(5.2)

constructed as shown in Section 3. As per Table 3 we first specify an inverse gamma distribution for the unknown ${\sigma }^2$ :

$\begin{array}{cc}p({\sigma }^2|a,b)& =\frac{b^a}{\Gamma \left(a\right)}({\sigma }^2{)}^{-(a+1)}exp[-b/{\sigma }^2]\\ & \propto ({\sigma }^2{)}^{-(a+1)}exp\left(-b/{\sigma }^2\right),\end{array}$(5.3)

where the form in the second line is just the kernel of the PDF for ${\sigma }^2$ stripping off the constants with proportionality. If a variable $X$ is distributed gamma, then $1/X$ is distributed inverse gamma, and this distribution is the conjugate form in the normal-normal setup. Here $a>0$ and $b>0$ are “hyperparameter values” of the prior that are set by the researcher. An additional utility of this example is in further demonstrating the usefulness of proportionality in simplifying Bayesian calculations. The normal kernel prior distribution for $\mu$ is given by:

$p(\mu |{\sigma }^2,m)\propto ({\sigma }^2{)}^{-\frac{1}{2}}exp\left(-\frac{1}{2{\sigma }^2}(\mu -m{)}^2\right),$(5.4)

where $m$ is the prior mean parameter for $\mu$ set by the researcher. The new wrinkle here is that this prior is actually called a “semi-conjugate” form since it needs to be conditioned on the ${\sigma }^2$ parameter in order for the conjugacy property to hold. The joint posterior distribution is produced in the same way that we did in Section 4, except that there are now two variables instead of one in the model, meaning that the posterior distribution is a “joint posterior” distribution for unknown $\mu$ and ${\sigma }^2$ and is thus slightly more complicated than the posterior form in (4.3). This joint posterior distribution is produced from multiplying the two prior distribution kernels times the likelihood function:

$\pi (\mu ,{\sigma }^2|y)\propto p(\mu |{\sigma }^2,m)\times p({\sigma }^2|a,b)\times L(\mu ,\sigma |y).$(5.5)

Now substitute the expressions (5.2), (5.3), and (5.4) into the definition of the posterior distribution earlier and collect like terms into a $({\sigma }^2{)}^{\text{something}}$ component and an $exp\left[\text{something else}\right]$ component to produce:

$\pi (\mu ,{\sigma }^2|y)\propto ({\sigma }^2{)}^{-\alpha -\frac{n}{2}-\frac{3}{2}}exp\left(-\frac{1}{{\sigma }^2}b-\frac{1}{2{\sigma }^2}\sum _{i=1}^n(y_i-\mu {)}^2-\frac{1}{2{\sigma }^2}(\mu -m{)}^2\right).$(5.6)

This joint form is large and awkward but can easily be marginalized (see details in Reference GillGill (2014)) with integration as explained further in Section 6. It is important not to be intimidated by the length of this posterior expression earlier: it is just a series of additive and multiplicative states in two exponentials multiplied together. Performing the marginalization first for ${\sigma }^2$ means integrating out the $\mu$ parameter:

$\begin{array}{cc}\pi \left({\sigma }^2|y\right)& ={\int }_{-\infty }^{\infty }\pi (\mu ,{\sigma }^2|y)d\mu \\ & \propto ({\sigma }^2{)}^{-a-\frac{n}{2}-\frac{3}{2}}exp\left(-\frac{1}{{\sigma }^2}\left(b+\frac{1}{2}\sum _{i=1}^n{y}_i^2-\frac{1}{2}n{\bar{y}}^2\right)\right).\end{array}$(5.7)

This feels complicated! The integral averages over uncertainty from the $\mu$ parameter leaving just a distributional expression for ${\sigma }^2$ , and we will discuss integration extensively in Section 6. Actually this expression earlier is neater than it seems. Looking back at the form of (5.3) we see that this marginal posterior is the kernel of an inverse gamma for ${\sigma }^2$ described by:

$p\left({\sigma }^2|y\right)\propto ({\sigma }^2{)}^{-\text{(something}+1)}\times exp[-\text{(something else)}/{\sigma }^2]$

where

$\begin{array}{cc}\text{something}& =a+\frac{n}{2}+\frac{1}{2}\end{array}$

$\begin{array}{cc}\text{something else}& =b+\frac{1}{2}\sum _{i=1}^n{y}_i^2-\frac{1}{2}n{\bar{y}}^2,\end{array}$

which defines the parameters of another inverted gamma distribution! So now we know that the posterior distribution of ${\sigma }^2$ is:

${\sigma }^2|y\sim IG\left(a+\frac{n}{2}+\frac{1}{2},b+\frac{1}{2}\sum _{i=1}^n{y}_i^2-\frac{1}{2}n{\bar{y}}^2\right).$(5.8)

Remember that no matter how awkward the parameterization is (although this one is not too bad), once we have the posterior distribution fully described we know everything about it and can describe it to readers any way we want. For instance the mean and variance of the generic inverse gamma PDF for $X$ are given by:

• $E\left[X\right]=\frac{b}{a-1},a>1$ .

• $\text{Var}\left[X\right]=\frac{b^2}{(a-1{)}^2(a-2)},a>2$ .

The marginal (really conditional because of the semi-conjugacy) posterior distribution for $\mu$ can be obtained by using the definition of conditional probability ( $p\left(A|B\right)=p(A,B)/p\left(B\right)$ ) and then some basic algebra as follows:

$\pi (\mu |\sigma ,y)=\frac{\pi (\mu ,{\sigma }^2|y)}{\pi \left({\sigma }^2|y\right)}={\sigma }^{-2}exp\left(-\frac{n}{2{\sigma }^2}\left({\mu }^2-2\frac{n\bar{y}+m}{n}\mu +\frac{n{\bar{y}}^2+m^2}{n}\right)\right)$(5.9)

(the $y$ just “comes along for the ride” here in the conditional probability statement). We can use the same “something” approach as before in looking at this form and see that this is the kernel of another normal distribution, according to:

$\mu |{\sigma }^2,y\sim N\left(\frac{n\bar{y}+m}{n},\frac{{\sigma }^2}{n}\right).$(5.10)

Here that the prior dependence on ${\sigma }^2$ flows through to the posterior for $\mu$ , which is why this is called semi-conjugate. Note also that as $n$ gets very large the posterior mean of $\mu$ converges to the data mean, which is intuitive but can be shown more technically:

$\underset{n\to \infty }{\text{lim}}\left(\frac{n\bar{y}+m}{n}\right)=\underset{n\to \infty }{\text{lim}}\left(\frac{n}{n}\bar{y}\right)+\underset{n\to \infty }{\text{lim}}\left(\frac{m}{n}\right)=\bar{y}.$(5.11)

After splitting up the first fraction into two components we see that the ratio of $n$ values becomes $1$ in the first term (even though they are infinity, they are the same “flavor” of infinity), and that the $n$ in the denominator of the second term eliminates that one. This is another clear illustration of the Bayesian principle that the data always “win” in the limit.

5.2 Analysis of Salary Data for Data Scientists

As an example of the normal-normal conjugate Bayesian model, consider the annual US salary numbers in thousands for Data Scientist positions listed for hire at Glassdoor in April 2023 with the sample in Table 4. This is the first 20 listed out of 8,880 since that is the most one can see without a lengthy registration process. Suppose we want to use this modest sample to make inferences about the national picture for Data Science pay. Salary figures in a specific occupation are often approximately normally distributed, even though a small sample of $n=20$ may not appear to be so. The sample mean is 172.1619 and the sample standard deviation is 26.7349. The hyperparameter values are set at $m=170,a=2,b=550$ to somewhat resemble sample values and provide a diffuse and conservative picture of variability. Calculation of the marginal posterior distributions directly follows from the recipes in (5.8) and (5.10). The results are shown in Figure 6. For the $\mu$ parameter we see that the prior is influential in the production of the posterior, which makes sense when the data size is only $n=20$ . This is an important illustration of the sample size principle in Bayesian inference: when the data are modest we will rely heavily on the specification of the prior and should be very introspective about its choice. For the ${\sigma }^2$ parameter we see that the data were more influential in determining the form of the posterior but not dramatically so.

Table 4 A year 2023 sample of salaries for data scientists in $1,000$ s

214.907	125.985	156.206	176.417	210.296	181.129	187.393
172.322	187.450	135.631	162.032	198.629	143.636	188.978
191.123	185.914	153.376	134.740	198.416	138.659

Figure 6 Normal-normal conjugate model for data science salaries

R Code for Normal-Normal Model
library(MCMCpack) # FOR dinvgamma(x, shape, scale = 1)
# DATA
salary <- scan("glassdoor.dat")
salary <- salary/1000; n <- length(salary)
# HYPERPRIOR VALUES
m <- 170; a <- 2; b <- 550

# POSTERIOR PARAMETERS
post.a <- a + n/2 + 1/2
post.b <- b + 0.5*sum(salary^2) - 0.5*n*mean(salary)^2
post.mu <- (n*mean(salary) + m)/n
post.var <- post.b/(post.a - 1)

# GRAPH
par(oma=c(1,1,1,1), mar=c(3,5,1,1),cex.lab=1.5,
    mfrow=c(2,1))
ruler <- seq(60,300,length=500)
prior.dens <- dnorm(ruler,m,sqrt(b/(a-1)))
plot(ruler,prior.dens,type="l",ylim=c(0,0.018),lwd=3,
    col="grey70", ylab=expression(paste(mu," Density")),
    xlab="")
post.dens <- dnorm(x=ruler,mean=post.mu,
    sd=sqrt(post.var))
lines(ruler,post.dens,lwd=3, col="grey30")
text(125,0.010,"Prior",col="grey70",cex=1.10,adj=0.5)
text(250,0.0050,"Posterior",col="grey10",cex=1.10,
    adj=0.5)
ruler <- seq(0,2000,length=500)
prior.dens <- dinvgamma(ruler,a,b)
plot(ruler,prior.dens,type="l",ylim=c(0,0.0025),
    xlim=c(-200,2000), lwd=3, col="grey70", xlab="",
    ylab=expression(paste(sigma^2," Density")))
post.dens <- dinvgamma(ruler,post.a,post.b)
lines(ruler,post.dens,lwd=3,col="grey30")
text(-50,0.0020,"Prior",col="grey70",cex=1.10,adj=0.5)
text(1060,0.0010,"Posterior",col="grey10",cex=1.10,
    adj=0.5)

Python Code for Normal-Normal Model
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm, invgamma

# DATA
salary = np.loadtxt("glassdoor.dat")
salary /= 1000; n = len(salary)

# HYPERPRIOR VALUES
m = 170; a = 2; b = 550

# POSTERIOR PARAMETERS
post_a = a + n/2 + 1/2
post_b = b + 0.5*sum(salary**2) - \
         0.5*n*np.mean(salary)**2
post_mu = (n*np.mean(salary) + m)/n
post_var = post_b/(post_a - 1)

# GRAPH
fig, axs = plt.subplots(2, 1, figsize=(8, 12))
ruler = np.linspace(60, 300, num=500)
prior_mu = norm.pdf(ruler, m, scale=np.sqrt(b/(a-1)))
post_mu = norm.pdf(ruler, post_mu, np.sqrt(post_var))
axs[0].plot(ruler, prior_mu, color='grey',
            linewidth=3, label='Prior')
axs[0].plot(ruler, post_mu, color='black',
            linewidth=3, label='Posterior')
axs[0].set_ylabel(r'$\mu$ Density')
axs[0].text(125, 0.010, 'Prior', color='grey',
            fontsize=12, ha='center')
axs[0].text(248, 0.0050, 'Posterior', color='black',
            fontsize=12, ha='center')
ruler = np.linspace(0, 2000, num=500)
prior_s2 = invgamma.pdf(ruler, a, scale=b)
post_s2 = invgamma.pdf(ruler, post_a, scale=post_b)
axs[1].plot(ruler, prior_s2, color='grey',
            linewidth=3, label='Prior')
axs[1].plot(ruler, post_s2, color='black',
            linewidth=3, label='Posterior')
axs[1].set_ylabel(r'$\sigma^2$ Density')
axs[1].text(0, 0.0020, 'Prior', color='grey',
            fontsize=12, ha='center')
axs[1].text(1050, 0.0010, 'Posterior', color='black',
            fontsize=12, ha='center')
plt.subplots_adjust(hspace=0.15)
plt.show()