Bayesian Statistics for Data Science - Greco
Frank Jung
23 December 2025
Two Coins
Suppose we have 2 different coins:
- Coin 1 is a fair coin with the P(Heads) = 0.5
- Coin 2 is a biased coin with the P(Heads) = 0.9
This is our prior distribution. We randomly select one of the coins (with
equal probability) and flip it once. It lands on heads. It is likely we have
chosen the biased on, but how sure are we? It’s more than 50%, but less than
100%.
We can use Bayes’ Rule to determine the posterior probability that we
selected Coin 2 given that we observed heads.
\[\begin{align*}
p(Coin 1 | Heads)
&= \frac {p(Heads \cap Coin 1)}
{p(Heads)} \\
\\
&= \frac {p(Heads \mid Coin 1) \; p(Coin 1)}
{p(Heads \cap Coin 1) + p(Heads \cap Coin 2)} \\
\\
&= \frac {p(Heads \mid Coin 1) \; p(Coin 1)}
{p(Heads \mid Coin 1) \; p(Coin 1) + p(Heads \mid Coin 2) \; p(Coin 2)} \\
\\
&= (0.5 * 0.5) / (0.5 * 0.5 + 0.9 * 0.5) \\
&\approx 0.3571429
\end{align*}\]
Where:
p_coin_1_given_heads <- (0.5 * 0.5) / (0.5 * 0.5 + 0.9 * 0.5)
p_coin_1_given_heads_percent <- round(p_coin_1_given_heads * 100, 3)
So the posterior probability that we selected Coin 1 given that we observed
heads is approximately 35.714%.
Therefore, the posterior probability that we selected Coin 2 given that we
observed heads is approximately 64.286%.
We can also calculate this directly using Bayes’ Rule:
p_coin_2_given_heads <- 0.9 * 0.5 / (0.9 * 0.5 + 0.5 * 0.5)
p_coin_2_given_heads_percent <- round(p_coin_2_given_heads * 100, 2)
The posterior probability that we selected Coin 2 given that we observed
heads is approximately 64.29%.
Multiple Observations
We can flip the coin multiple times and update our beliefs accordingly.
Suppose on our second flip we see tails, so our belief should shift back
towards the fair coin:
\[\begin{align*}
p(Coin 1 | HT)
&= \frac {p(HT \mid Coin 1) \; p(Coin 1)}
{p(HT \mid Coin 1) \; p(Coin 1) + p(HT \mid Coin 2) \; p(Coin 2)} \\
\\
&= (0.5 * 0.5 * 0.5) / (0.5 * 0.5 * 0.5 + 0.9 * 0.1 * 0.5) \\
&\approx 0.7352941
\end{align*}\]
p_coin_1_given_ht <- (0.5 * 0.5 * 0.5) / (0.5 * 0.5 * 0.5 + 0.9 * 0.1 * 0.5)
p_coin_1_given_ht_percent <- round(p_coin_1_given_ht * 100, 2)
The posterior probability of Coin 1 giving us Heads then Tails is approximately
73.53%.
We can also use the previous posterior as our new prior to calculate this
posterior. In this case it will be:
p_coin_1_given_ht <- 0.357 * 0.5 / (0.357 * 0.5 + 0.1 * 0.643)
p_coin_1_given_ht_percent <- round(p_coin_1_given_ht * 100, 2)
The posterior probability of Coin 1 giving us Heads then Tails is approximately
73.52%.
Coin Flipping
You have two coins: Coin 1 is a fair coin with equal probability of heads and
tails, and Coin 2 has a 75% probability of heads and a 25% probability of tails.
You choose a coin at random and flip it. If it lands heads, what is the
probability that you chose Coin 1?
We want to find \(P(C_1 | H)\):
\[\begin{align*}
P(C_1 | H)
&= \frac {P(H | C_1) P(C_1)} {P(H | C_1) P(C_1) + P(H | C_2) P(C_2)} \\
&= \frac {0.5 \cdot 0.5} {0.5 \cdot 0.5 + 0.75 \cdot 0.5} \\
&= \frac {2} {5} \\
&= 0.4
\end{align*}\]
p_c_1_given_h <- (0.5 * 0.5) / (0.5 * 0.5 + 0.75 * 0.5)
p_c_1_given_h_percent <- round(p_c_1_given_h * 100, 2)
The probability that you chose Coin 1 given that it landed heads is about
40%.
Medical Testing
A certain disease has a prevalence rate of 1% in the population. A diagnostic
test for this disease has a 98% true positive rate (sensitivity) and a 97% true
negative rate (specificity). If a person tests positive for the disease, what
is the probability that they actually have the disease?
We want to find \(P(D \mid +)\):
p_d_given_pos <- (0.98 * 0.01) / (0.98 * 0.01 + (1 - 0.97) * (1 - 0.01))
p_d_given_pos_percent <- round(p_d_given_pos * 100, 2)
So, the probability that a random person has the disease, given they have
tested positive, is approximately 24.81%.
Email Spam Filtering
An email spam filter has a 99% success rate in detecting spam emails and a 95%
success rate in correctly identifying legitimate emails. If 20% of all emails
are spam, what is the probability that an email flagged as spam is actually
spam?
We want to find \(P(S | F)\):
p_s_given_f <- (0.99 * 0.2) / (0.99 * 0.2 + (1 - 0.95) * (1 - 0.2))
p_s_given_f_percent <- round(p_s_given_f * 100, 2)
The probability that an email flagged as spam is actually spam is about
83.19%.
Manufacturing Defects
A factory has two production lines (Line 1 and Line 2) manufacturing the same
product. Line 1 produces 70% of the products, while Line 2 produces 30% of the
products. Line 1 has a defect rate of 4%, and Line 2 has a defect rate of 8%.
If a product is found to be defective, what is the probability it was produced
by Line 1?
We want to find \(P(L_1 | D)\):
p_l1_given_d <- (0.04 * 0.7) / (0.04 * 0.7 + 0.08 * 0.3)
p_l1_given_d_percent <- round(p_l1_given_d * 100, 2)
The probability that a defective product was produced by Line 1 is about
53.85%.
Internet Connectivity
Two internet service providers (ISP 1 and ISP 2) offer internet connectivity in
a neighbourhood. ISP 1 serves 80% of the households, while ISP 2 serves the
remaining 20%. The connection failure rates for ISP 1 and ISP 2 are 1% and 3%,
respectively. If a household experiences a connection failure, what is the
probability that they are using ISP 1?
We want to find \(P(ISP_1 \mid F)\):
p_i1_given_f <- (0.01 * 0.8) / (0.01 * 0.8 + 0.03 * 0.2)
p_i1_given_f_percent <- round(p_i1_given_f * 100, 2)
The probability that a household experiencing a connection failure is using ISP
1 is about 57.14%.
Beta Distribution
Plot the Beta distribution for \(\alpha = 10\) and \(\beta = 10\):
## Loading required package: ggplot2
Probability Calculation
What is the probability that \(\theta\) is between 0.4 and 0,6?
The prior probability that \(\theta\) is between 0.4 and 0.6 is
about 62.78%.
95% Credible Interval
To calculate the 95% credible interval for the Beta(10, 10) distribution:
lower_bound <- qbeta(0.025, 10, 10)
upper_bound <- qbeta(0.975, 10, 10)
Which gives us the interval
[0.2886, 0.7114]
Probability density function
The prior is a probability density function (PDF) of the Beta distribution:
\[\begin{align*}
\pi(x; \alpha, \beta)
&= \frac {\Gamma (\alpha + \beta)} {\Gamma(\alpha) \Gamma(\beta)} \; x^{\alpha - 1} (1 - x)^{\beta - 1} \\
&= \frac {1} {B(\alpha, \beta)} \; x^{\alpha - 1} (1 - x)^{\beta - 1}
\end{align*}\]
The beta function, \(B\) is a normalization constant to ensure that the total
probability is 1. In the above equations \(x\) is a realization—an observed value
that actually occurred—of a random variable \(X\).
The likelihood function could be modeled as a binomial distribution:
\[L(x; n, \theta) = \binom {n} {x} \; \theta^x (1 - \theta)^{n - x}\]
Recall that the posterior distribution is proportional to the product of the
prior and the likelihood. And if remove anything that does not depend on
\(\theta\) from the likelihood if we are interested in only the proportionality:
\[\begin{align*}
\text{Posterior}
&\propto \pi(\theta; \alpha, \beta) \times L(x; n, \theta) \\
&\propto \theta^{\alpha - 1} (1 - \theta)^{\beta - 1} \times \theta^x (1 - \theta)^{n - x} \\
&\propto \theta^{x + \alpha - 1} (1 - \theta)^{n - x + \beta - 1} \\
&\propto \text{Beta}(\alpha + x, \beta + n - x) \\
\end{align*}\]
Where the constant of proportionality is the beta function.
So we have a prior of Beta(\(\alpha\), \(\beta\)) and a posterior of
Beta(\(\alpha + x\), \(\beta + n - x\)).
\[\begin{align*}
\text{Prior} &\approx \text{Beta}(\alpha_{prior}, \beta_{prior}) \\
\text{Posterior} &\approx \text{Beta}(\alpha_{prior} + x, \beta_{prior} + n - x)
\end{align*}\]
Where:
- \(\alpha_{posterior} = \alpha_{prior} + x\) which adds the number of successes to \(\alpha\)
- \(\beta_{posterior} = \beta_{prior} + n - x\) which adds the number of failures to \(\beta\)
An interpretation here is that \(\alpha\) represents the number of successes,
while \(\beta\) represents the number of failures.
The posterior distribution describes our updated beliefs after observing the
data.
Beta-Binomial Survey Example
This is a Beta-Binomial example: Given a marketing campaign we want to learn a
the response rate to our emails. We will model the response rate \(\theta\) with a
Binomial random variable. A response is \(Y=1\) and no response is \(Y=0\).
We will use a beta prior on \(\theta\) to represent our beliefs about the response
rate before seeing any data.
We initially estimate that we will have a 50% response rate, so we set our prior
with:
samples <- 10000 # number of samples to draw
alpha_prior <- 5 # our prior belief of the response rate
beta_prior <- 5 # our prior belief of the non-response rate
We can visualise the prior distribution:
Let’s find the 95% credible interval for our parameter \(\theta\):
ci <- qbeta(c(0.025, 0.975), alpha_prior, beta_prior)
The 95% prior credible interval for our prior is between
0.212 and 0.788.
We have sent out 50 emails and 20 of them have responded.
n <- 50 # total emails sent
x <- 20 # total responses
We can update our prior to get the posterior distribution:
alpha_post <- alpha_prior + x
beta_post <- beta_prior + (n - x)
This gives us a posterior distribution of Beta(25, 35).
We can visualise the posterior distribution:
Let’s find the 95% posterior credible interval for our parameter \(\theta\):
ci_post <- qbeta(c(0.025, 0.975), alpha_post, beta_post)
The 95% posterior credible interval for our posterior is from
0.296 to 0.543.
If we were to send out another 50 emails, what is the predictive distribution
for the number of responses we would get? The predictive distribution accounts
for both the inherent variability in the data and the uncertainty in our
estimate of \(\theta\).
predictive_responses <- rbeta(samples, alpha_post, beta_post)
predictive_samples <- rbinom(samples, size = n, prob = predictive_responses)
We can visualise the predictive distribution:
This is a Beta-Binomial model, where we have a Beta prior on the response rate
\(\theta\) and a Binomial likelihood for the number of responses given \(\theta\).
In summary we have simulated a survey response scenario using a hierarchical
process (or compound distribution). We started by simulating the probability of
getting a response to our emails using a beta distribution, and then used that
probability of success to simulate from a binomial distribution.
The VGAM R package
provides functions to work with the Beta-Binomial distribution directly.
## Loading required package: VGAM
## Loading required package: stats4
## Loading required package: splines
Note that this is not the same as just picking from a binomial distribution
(see below). As the compound beta-binomial distribution above has more variance
than the binomial distribution alone, because it accounts for the uncertainty
in our estimate of \(\theta\).
Gamma Distribution
Average number of cars passing through an intersection
Suppose we want to learn about the average number of cars that pass an
intersection in a given day. We thinks its around 100 cars per day, but are not
certain. We put a prior on \(\theta\) using a Gamma distribution with shape=10
and rate=0.1:
\[ \lambda \sim \text{Gamma}(shape=10, rate=0.1) \]
Note that the mean = shape/rate = 10/0.1 = 100 cars per day.
Next we observe the intersection for n = 10 days. The number of cars observed
each day are:
observed_cars <- c(49, 53, 45, 58, 51, 39, 44, 42, 64, 40)
total_cars <- sum(observed_cars)
n_days <- length(observed_cars)
So our posterior distribution is:
shape_post <- shape + sum(observed_cars)
rate_post <- rate + n_days
This gives us a posterior distribution of Gamma(495, 10.1).
The posterior mean is shape/rate = 495/10.1 =
49.01 cars per day.
Credible Interval
shape <- 30 # alpha
rate <- 5 # beta
Given a Gamma(30, 5) distribution, what is the 95% credible
interval?
# Calculate 95% credible interval
ci_lower <- qgamma(0.025, shape = shape, rate = rate)
ci_upper <- qgamma(0.975, shape = shape, rate = rate)
The 95% credible interval is [4.048, 8.33].
Mean
For a Gamma distribution we can get the mean using the formula:
\[ \text{Mean} = \frac {shape} {rate} \]
mean <- shape / rate
The mean of the Gamma(30, 5) distribution is approximately
6.
Mode
The mode of a Gamma distribution can be calculated using the formula:
\[ \text{Mode} = \frac {shape - 1} {rate} \quad \text{for } shape \ge 1, \text{ and } 0 \text{ for } shape < 1 \]
mode <- (shape - 1) / rate
Here, the mode of the Gamma(30, 5) distribution is approximately
5.8.
Customers arriving at a store
We have a store and are studying the number of customers arriving per day. We
observe the number of customers arriving for 10 days and want to infer the
average number of customers arriving per day.
We can model the daily arrivals using a Poisson distribution where rate
parameter \(\lambda\) represents the average number of arrivals per day. Since
\(\lambda\) is unknown, we place a Gamma prior on it. The Gamma distribution is a
suitable prior for \(\lambda\) because it is always positive and can represent a
wide range of shapes depending on its parameters.
We think around 10 customers arrive per day, but we are not certain. This means
the expected value of our prior should be 10:
\[\frac {\alpha} {\beta} = 10\]
As we are uncertain, then standard deviation should be relatively high (maybe
also around 10):
\[\sqrt {\frac {\alpha} {\beta^2}} = 10\]
Solving these two equations gives us \(\alpha = 1\) and \(\beta = 0.1\).
alpha <- 1
beta <- 0.1
We can visualise the prior distribution with a histogram:
We think this prior distribution is reasonable as we are confident that we have
fewer then 40 customers per day.
ci <- qgamma(c(0.025, 0.975), shape = alpha, rate = beta)
The 95% credible interval for our prior is between
0.25 and 36.89 customers per day.
We simulate customer arrivals for 10 days.
set.seed(123) # for reproducibility
n <- 10 # number of days
true_lambda <- 5 # true average number of arrivals per day
arrivals <- rpois(n, lambda = true_lambda)
arrivals
## [1] 4 7 4 8 9 2 5 8 5 5
We now can update the parameters and plot a histogram of the posterior
distribution:
alpha_post <- alpha + sum(arrivals) # add total arrivals to alpha
beta_post <- beta + n # add number of days to beta
This gives us a posterior distribution of Gamma(58, 10.1).
Let’s find the 95% credible interval for our posterior distribution:
ci_post <- qgamma(c(0.025, 0.975), shape = alpha_post, rate = beta_post)
The 95% credible interval for our posterior distribution is between
4.36 and 7.31 customers per day.
This is a predictive distribution for future observations, accounting for the
inherent daily uncertainty and our uncertainty about the true average arrival
rate \(\lambda\).
predictive_lambda <- rgamma(10000, shape = alpha_post, rate = beta_post)
future_arrivals <- rpois(10000, lambda = predictive_lambda)
hist(future_arrivals,
main = "Predictive Distribution of Daily Customers",
xlab = "Number of Customers",
ylab = "Frequency",
col = "steelblue",
border = "steelblue"
)
We can also use the negative binomial distribution to model the predictive
distribution directly:
require(VGAM)
param1 <- alpha_post
param2 <- beta_post / (1 + beta_post)
future_arrivals_nb <- rnbinom(10000, size = alpha_post, prob = param2)
hist(future_arrivals_nb,
main = "Predictive Distribution of Daily Customers (Negative Binomial)",
xlab = "Number of Customers",
ylab = "Frequency",
col = "steelblue",
border = "steelblue"
)
The point estimates for the average daily customers are:
mean_lambda <- alpha_post / beta_post
median_lambda <- qgamma(0.5, shape = alpha_post, rate = beta_post)
mode_lambda <- (alpha_post - 1) / beta_post
The posterior customers per day:
- mean is approximately 5.74
- median is approximately 5.71
- mode is approximately 5.64
Normal Distribution
Let’s assume we’re conducting quality control in a factory. We have a machine
that produces candy, and the weight of these candies follows a normal
distribution.
We know the standard deviation of the weight is 2 grams, but the machine has a
control to adjust the average weight.
The setting on the machine has been lost, so we don’t know the average weight
it is set for.
We weigh 10 pieces of candy, and we want to infer the average weight from these
measurements.
I think the weight of the average candy is around 100 grams (that’s what
the dial was set at before!), but I’m very uncertain.
Prior
The prior hyper-parameters are:
mu_prior <- 100
sigma_prior <- 50
The prior distribution is Normal(100, 50), which looks
like:
ci <- qnorm(c(0.025, 0.975), mu_prior, sigma_prior)
The 95% prior credible interval for the average weight is between
2 and 198 grams.
Posterior
Next we will simulate the data of the weights of 10 observed candies. In real
life, this data would be given to us, we are simulating it so that we don’t
have to read in files of data, and to demonstrate that the true mean which is
unknown to us is 50.
We will simulate 10 candies using the Normal distribution with a mean of 50
grams and a standard deviation of 2 grams (Normal(50, 2)).
set.seed(123) # for reproducibility
n <- 10 # number of candy pieces
true_mu <- 50 # true average weight of candy (in grams)
sigma <- 2 # known standard deviation (in grams)
x <- rnorm(n, true_mu, sigma) # observed weights
The observed weights are: 48.88, 49.54, 53.12, 50.14, 50.26, 53.43, 50.92, 47.47, 48.63, 49.11 grams.
Given this data, we will now update our parameters, and then visualize the
posterior distribution.
mu_post <- (sigma^2 * mu_prior + n * mean(x) * sigma_prior^2) / (sigma^2 + n * sigma_prior^2)
sigma_post <- sqrt((sigma^2 * sigma_prior^2) / (sigma^2 + n * sigma_prior^2))
The posterior distribution is
Normal(50.16, 0.63).
ci_post <- qnorm(c(0.025, 0.975), mu_post, sigma_post)
This has a 95% credible interval between 48.92 and
51.4 grams.
Predictive Distribution
If we were to draw future observations, what would they look like?
The prior predictive distribution is given by:
\[ \hat{y} \sim \text{Normal}(\mu_{prior}, \sqrt{\sigma^2 + \sigma_{prior}^2}) \]
Where as the posterior predictive distribution is given by:
\[ \hat{y} \sim \text{Normal}(\mu_{post}, \sqrt{\sigma^2 + \sigma_{post}^2}) \]
We can simulate this by first drawing values of \(\mu\) from the posterior, and
then drawing future observations from a normal distribution with those \(\mu\)
values and the known standard deviation \(\sigma\). (This is a hierarchical or
compound process.)
mu_draws <- rnorm(100000, mu_post, sigma_post)
future_obs <- rnorm(1000000, mean = mu_draws, sd = sigma)
We can also do this by direct simulation:
future_obs_direct <- rnorm(1000000, mean = mu_post, sd = sqrt(sigma_post^2 + sigma^2))
future_ci_direct <- quantile(future_obs_direct, c(0.025, 0.975))
Which gives us a prediction interval for future observations between
46.04 and
54.27 grams.
Unknown Initial Standard Deviation
If we don’t know the standard deviation, we can place a prior on it as well.
The prior will be related to the Gamma distribution.
The predictive distribution will then be a Student’s t-distribution, which has
fatter tails than the normal distribution, reflecting our increased uncertainty
about the parameters.
Simple Linear Regression
In this section we will explore a more complicated scenario’s. We will be
looking at simple linear regression:
- What is the model (likelihood) that generates the data?
- What are the parameters we want to estimate?
- How do we express our knowledge about the parameters with our prior?
- How do we update our beliefs about the parameters with the data (posterior)?
We have discussed three types of scenarios so far:
- Counts of successes or failures (Bernoulli/Binomial)
- Count data (Poisson)
- Continuous or Quantitative data (Normal)
In each of these scenarios we tried to estimate a parameter (\(\theta\),
\(\lambda\), \(\mu\)) that described the data, \(X_1, X_2, \ldots, X_n\).
We will now look at a slightly more complicated scenario and see how we can use
Bayesian methods.
The previous examples involved univariate data (a single variable). We will now
look at bivariate data (two variables) that are correlated. The two variables
are paired. For example height and weight of individuals.
We can model the relationship between the two variables using simple linear
regression:
\[ y_i = \beta_0 + \beta_1 x_i + \epsilon_i \]
Where $ (x_i, y_i) $ are the paired observations, $ _i $ is the error
term, and \(\beta_0\) is the intercept and \(\beta_1\) is the slope of the line.
If \(x\) was height and \(y\) was weight, then \(\beta_1\) would represent the change
in weight for a unit change in height.
In simple linear regression we assume
\(\epsilon_i \sim \text{Normal}(0, \sigma^2)\).
For now we will assume that \(\sigma^2\) is known, but commonly it is unknown and
we will use a inverse gamma prior on it.
We will use a Normal prior on the regression coefficients \(\beta_0\) and
\(\beta_1\). The posterior distribution will also be Normal, as the Normal
distribution is conjugate to itself. We could also put a prior on the variance.
Then we will have a Normal-Inverse-Gamma prior on the parameters.
If \(\sigma^2\) is unknown, then the posterior distribution for the coefficients
will be a Student’s t-distribution.
You can also use an uninformative prior on the coefficients. For example you
could use a uniform distribution on \(\theta\) in a Beta distribution.
Web Site Visits and Purchases
For this example we will assume there is a linear relationship between the
amount a persons spends and the time they spend on a website.
We will simulate the data for this example with:
- the true intercept \(\beta_0 = 10\)
- the true slope \(\beta_1 = 20\)
- the true variance \(\sigma = 5\)
set.seed(123) # for reproducibility
n <- 100 # number of observations
beta_0_true <- 10 # true intercept
beta_1_true <- 20 # true slope
sigma_true <- 5 # true standard deviation
x <- runif(n, min = 0, max = 5) # predictor variable (time on site in hours)
y <- rnorm(n, mean = beta_0_true + beta_1_true * x, sd = sigma_true) # response variable (amount spent)
Plot the simulated data:
We want to do is learn about the parameters \(\beta_0\) and \(\beta_1\) given the
data.
## Loading required package: brms
## Loading required package: Rcpp
## Loading 'brms' package (version 2.23.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
## Attaching package: 'brms'
## The following objects are masked from 'package:VGAM':
##
## acat, cratio, cumulative, dfrechet, dirichlet, exponential,
## frechet, geometric, lognormal, multinomial, negbinomial, pfrechet,
## qfrechet, rfrechet, s, sratio
## The following object is masked from 'package:stats':
##
## ar
## Compiling Stan program...
## Start sampling
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 7e-06 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.008 seconds (Warm-up)
## Chain 1: 0.006 seconds (Sampling)
## Chain 1: 0.014 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 2e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.009 seconds (Warm-up)
## Chain 2: 0.007 seconds (Sampling)
## Chain 2: 0.016 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 2e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.008 seconds (Warm-up)
## Chain 3: 0.006 seconds (Sampling)
## Chain 3: 0.014 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 2e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.008 seconds (Warm-up)
## Chain 4: 0.006 seconds (Sampling)
## Chain 4: 0.014 seconds (Total)
## Chain 4:
summary(fit)
## Family: gaussian
## Links: mu = identity
## Formula: y ~ x
## Data: data_df (Number of observations: 100)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 9.87 0.97 7.94 11.75 1.00 4079 3244
## x 19.89 0.33 19.25 20.54 1.00 4515 3090
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 4.79 0.34 4.18 5.53 1.00 3826 2670
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit)
The posterior mean for the intercept \(\beta_0\) is approximately
9.87,
with a 95% credible interval of
(7.94,
11.75).
The posterior mean for the slope \(\beta_1\) is approximately
19.89,
with a 95% credible interval of
(19.25,
20.54).
The posterior mean for the variance \(\sigma\) is approximately
4.79,
with a 95% credible interval of
(4.18,
5.53).
Certificate