Posterior counterparts to paired, exact, and contingency-table tests
Problem Set 1 begins with a null hypothesis and calculates the probability of a test statistic under a null reference distribution. This assignment begins with a likelihood and a prior distribution. Bayes’ rule updates the prior to a posterior distribution. Posterior probabilities are probabilities over parameter values conditional on the observed data. They are not p values.
The exercises follow the same sequence as Problem Set 1. A normal location model corresponds to the paired t analysis. A robust Student t location model serves the same inferential role as the Wilcoxon analysis. A beta–Bernoulli model corresponds to the exact sign analysis. Separate beta–binomial models correspond to the two-row association analyzed by Fisher’s exact test. A posterior predictive independence model corresponds to the chi-squared analysis. The reading-time data become relevant in Exercises 1–4. The syntactic data become relevant in Exercises 5–6.
All six exercises are required. The points sum to 100. Every parameter is scalar. Bold symbols denote collections of independent and identically distributed observations, not vector-valued parameters.
Download the complete Problem Set 2 package. Extract the package and open ps2.ipynb from the extracted ps2-assignment directory. Select the R kernel. The package contains the R notebook, the assignment source, the README, the preparation scripts, and the redistributable teaching tables.
The files are also available separately:
Use <- for assignment, use snake_case object names, and record a seed before every random calculation. Every numerical answer must be produced by visible code.
Each derivation is divided into short steps. Complete only the displayed line or lines in each step. The surrounding equations supply the notation and the algebraic operation. You do not need to reproduce them.
In an R block, replace every ___ with the requested expression or value. Code without a blank is supplied as scaffolding and should be retained.
# Run the notebook from the extracted ps2-assignment directory.
stopifnot(
file.exists("README.md"),
file.exists("data/provo-content-word-positions.csv"),
file.exists("data/ud_wh_dependencies.csv")
)
library(dplyr)
library(ggplot2)
library(tidyr)The data for Exercises 1–4 come from the Provo Corpus. One row represents one content-word position. The response mean_log_ffd is the mean log first-fixation duration over observed readers who fixated that position. Thus the analysis is conditional on fixation and targets an observed-reader mean, not individual reading behavior. Text_ID identifies the passage and Word_POS identifies the word class.
The comparison is paired by passage. Within each passage, average mean_log_ffd over noun positions and over verb positions. Then define
\delta_p = \overline{\log(\operatorname{FFD})}_{p,\mathrm{noun}} - \overline{\log(\operatorname{FFD})}_{p,\mathrm{verb}}.
Thus \exp(\delta_p) is the ratio of the noun-position geometric mean FFD to the verb-position geometric mean FFD in passage p.
Load data/provo-content-word-positions.csv. Retain noun and verb positions. Verify that each passage contains at least four positions from each class. Summarize within passage and class, then join by Text_ID to obtain one row per passage. Do not align the summaries by row position.
positions <- read.csv("data/provo-content-word-positions.csv")
str(positions)
sort(unique(positions$Word_POS))
noun_verb_positions <- positions |>
filter(Word_POS %in% c("Noun", "Verb"))
position_counts <- noun_verb_positions |>
count(Text_ID, Word_POS)
stopifnot(all(position_counts$n >= 4))
passage_class_means <- noun_verb_positions |>
group_by(Text_ID, Word_POS) |>
summarize(mean_log_ffd = mean(mean_log_ffd), .groups = "drop")
passage_pairs <- passage_class_means |>
pivot_wider(
names_from = Word_POS,
values_from = mean_log_ffd
) |>
rename(
mean_log_ffd_noun = Noun,
mean_log_ffd_verb = Verb
) |>
mutate(delta_log_ffd = mean_log_ffd_noun - mean_log_ffd_verb)
stopifnot(nrow(passage_pairs) == 55, !anyNA(passage_pairs$delta_log_ffd))Calculate the number of passages, the mean, the standard deviation, the median, and the number of negative, zero, and positive differences. Plot the differences with zero marked.
passage_summary <- passage_pairs |>
summarize(
passage_count = n(),
mean_delta = ___,
sd_delta = ___,
median_delta = ___,
negative_count = ___,
zero_count = ___,
positive_count = ___
)
passage_summary
ggplot(passage_pairs, aes(x = delta_log_ffd)) +
geom_dotplot(binwidth = .01) +
geom_vline(xintercept = 0, linetype = "dashed") +
labs(x = "Noun minus verb mean log FFD", y = "Count")Report the requested values in the subtraction order noun minus verb.
Replace this sentence with the sample size and requested summaries.
Let Y_p:\Omega\to\mathbb R be the random noun minus verb difference for passage p, and let \delta_p be its observed value. Let M:\Omega\to\mathbb R be the population mean difference and \Sigma:\Omega\to\mathbb R_{>0} the among-passage standard deviation. Use
Y_p\mid M=\mu,\Sigma=\sigma \sim \mathcal N(\mu,\sigma^2)
with scalar priors
M\sim\mathcal N(0,.05^2) \qquad\text{and}\qquad \Sigma\sim\operatorname{Exponential}(20).
In the course notation, the second argument of \mathcal N(\mu,\sigma^2) is the variance. In brms, normal(0, .05) uses the standard deviation .05.
The construction of a posterior from a likelihood and prior is developed in posterior distributions. Conjugate priors develops the beta–Bernoulli case. The normal–normal derivation below is not supplied there. You must carry out that derivation before fitting the model.
This derivation establishes a conjugate update that is not derived in the notes. Let P be the number of passages and define
\overline Y = \frac{1}{P}\sum_{p=1}^{P}Y_p, \qquad \overline\delta = \frac{1}{P}\sum_{p=1}^{P}\delta_p.
Conditional on a fixed value \Sigma=\sigma,
\overline Y\mid M=\mu,\Sigma=\sigma \sim \mathcal N\left(\mu,\frac{\sigma^2}{P}\right).
Write m_0=0 and s_0=.05. The quantity to derive is the scalar conditional density
f_{M\mid\Sigma,\overline Y}(\mu\mid\sigma,\overline\delta).
The normal kernel is \exp\{-(\text{observation}-\text{mean})^2/(2\,\text{variance})\}. Substitute observation \overline\delta, mean \mu, and variance \sigma^2/P. Complete one line.
Target: f_{\overline Y\mid M,\Sigma}(\overline\delta\mid\mu,\sigma) as a kernel in \mu.
Your LaTeX
\begin{aligned} f_{\overline Y\mid M,\Sigma}(\overline\delta\mid\mu,\sigma) &\propto\text{enter the likelihood kernel in }\mu. \end{aligned}
Use the same normal kernel with mean m_0 and variance s_0^2. Complete one line.
Target: f_M(\mu) as a kernel in \mu.
Your LaTeX
\begin{aligned} f_M(\mu)&\propto\text{enter the prior kernel in }\mu. \end{aligned}
Multiply the two kernels. Their exponents add. Expand the two squares and discard \overline\delta^2 and m_0^2 because they do not depend on \mu. The remaining terms have the form
\exp\left\{-\frac12(A\mu^2-2B\mu)\right\}.
Collect the coefficient of \mu^2 into A and the coefficient paired with 2\mu into B. Complete two lines.
Your LaTeX
\begin{aligned} A&=\text{enter the coefficient of }\mu^2,\\ B&=\text{enter the coefficient paired with }2\mu. \end{aligned}
Use A\mu^2-2B\mu=A(\mu-B/A)^2-B^2/A. The normal kernel therefore has mean B/A and variance 1/A. Substitute the expressions for A and B from Step 3. Complete two lines.
Your LaTeX
\begin{aligned} m_P&=\text{enter the posterior mean in the original quantities},\\ s_P^2&=\text{enter the posterior variance in the original quantities}. \end{aligned}
Insert the two quantities from Step 4 into the normal distribution. Complete one line.
M\mid\Sigma=\sigma,\overline Y=\overline\delta \sim\mathcal N\left(\text{enter the mean},\text{ enter the variance}\right).
Define the event that the population mean difference is positive. Use full set-builder notation. Then standardize the posterior normal distribution to express its probability with the standard normal cumulative distribution function \Phi.
Define the event in the first line. Standardizing zero under \mathcal N(m_P,s_P^2) gives the second line. Complete two lines.
Your LaTeX
\begin{aligned} B_M^+&=\{\omega\in\Omega\mid\text{enter the condition on }M(\omega)\},\\ \mathbb P(B_M^+\mid\Sigma=\sigma,\overline Y=\overline\delta) &=\text{enter the expression in }m_P,s_P,\text{ and }\Phi. \end{aligned}
A prior predictive distribution is the distribution of replicated data obtained by drawing parameters from their priors and data from the likelihood. Draw at least 10,000 values of \mu, \sigma, and one replicated difference. Plot the replicated differences and summarize \exp(Y_{\mathrm{rep}}).
set.seed(___)
prior_draw_count <- 10000
mu_prior <- rnorm(prior_draw_count, mean = 0, sd = .05)
sigma_prior <- rexp(prior_draw_count, rate = 20)
y_rep_prior <- rnorm(
prior_draw_count,
mean = mu_prior,
sd = sigma_prior
)
ggplot(tibble(y_rep_prior), aes(x = y_rep_prior)) +
geom_histogram(bins = 60) +
geom_vline(xintercept = 0, linetype = "dashed")
quantile(exp(y_rep_prior), probs = c(.025, .25, .5, .75, .975))Fit the normal model with brms. Use four chains, at least 2,000 iterations per chain, the cmdstanr backend, and a recorded seed.
library(brms)
library(posterior)
detected_cores <- parallel::detectCores()
core_count <- if (is.na(detected_cores)) 1L else min(4L, detected_cores)
normal_fit <- brm(
delta_log_ffd ~ 1,
data = passage_pairs,
family = gaussian(),
prior = c(
prior(normal(0, .05), class = Intercept),
prior(exponential(20), class = sigma)
),
chains = 4,
cores = core_count,
iter = 2000,
seed = ___,
backend = "cmdstanr",
refresh = 0
)Report \widehat R, bulk and tail effective sample size, and divergent transitions for the intercept and sigma. The Markov-chain checks put these diagnostics in the order used here.
normal_diagnostics <- summarise_draws(
as_draws(normal_fit),
"rhat",
"ess_bulk",
"ess_tail"
)
normal_diagnostics |>
filter(variable %in% c("b_Intercept", "sigma"))
normal_nuts <- nuts_params(normal_fit)
divergence_count <- sum(
normal_nuts$Parameter == "divergent__" &
normal_nuts$Value == 1
)
divergence_countFrom the posterior draws, calculate the median and 95% equal-tailed credible interval for M and \exp(M). Estimate \mathbb P(B_M^+\mid\boldsymbol\delta) as the proportion of intercept draws greater than zero, where \boldsymbol\delta=(\delta_1,\ldots,\delta_P).
normal_draws <- as_draws_df(normal_fit)
mu_draw <- normal_draws$b_Intercept
mu_summary <- c(
median = median(mu_draw),
lower_95 = quantile(mu_draw, .025),
upper_95 = quantile(mu_draw, .975)
)
ratio_draw <- exp(mu_draw)
ratio_summary <- c(
median = median(ratio_draw),
lower_95 = quantile(ratio_draw, .025),
upper_95 = quantile(ratio_draw, .975)
)
probability_mu_positive <- mean(mu_draw > 0)
mu_summary
ratio_summary
probability_mu_positiveReport the computational checks, posterior summaries, and posterior event probability.
Replace this sentence with the requested normal-model result.
The Wilcoxon test reduces the influence of numerical spacing by using ranks. A Bayesian Student t model takes a different route: it retains the measured differences but assigns more probability to extreme observations than a normal likelihood does. It is not a Bayesian version of the signed-rank test. It is the corresponding robust location analysis in this assignment.
Replace the normal likelihood in Exercise 2 with
Y_p\mid M=\mu,\Sigma=\sigma,\Nu=\nu \sim t_\nu(\mu,\sigma),
where M and \Sigma retain the priors from Exercise 2 and the scalar degrees-of-freedom parameter has prior \Nu\sim\operatorname{Gamma}(2,.1). Here the gamma distribution uses shape 2 and rate .1, so its mean is 2/.1=20. The Student t distribution approaches a normal distribution as \nu increases and has heavier tails for smaller \nu.
Fit this model with family = student() and prior(gamma(2, .1), class = nu). Use the same chain settings and seed policy as Exercise 2.
student_fit <- brm(
delta_log_ffd ~ 1,
data = passage_pairs,
family = student(),
prior = c(
prior(normal(0, .05), class = Intercept),
prior(exponential(20), class = sigma),
prior(gamma(2, .1), class = nu)
),
chains = 4,
cores = core_count,
iter = 2000,
seed = ___,
backend = "cmdstanr",
refresh = 0
)
student_diagnostics <- summarise_draws(
as_draws(student_fit),
"rhat",
"ess_bulk",
"ess_tail"
)
student_diagnostics |>
filter(variable %in% c("b_Intercept", "sigma", "nu"))
student_nuts <- nuts_params(student_fit)
sum(
student_nuts$Parameter == "divergent__" &
student_nuts$Value == 1
)Calculate the posterior median and 95% interval for M, \exp(M), and \Nu. Calculate \mathbb P(B_M^+\mid\boldsymbol\delta) from the posterior draws.
student_draws <- as_draws_df(student_fit)
student_mu_draw <- student_draws$b_Intercept
nu_draw <- student_draws$nu
student_summary <- tibble(
model = "Student t",
quantity = c("M", "exp(M)", "Nu"),
median = c(
median(student_mu_draw),
median(exp(student_mu_draw)),
median(nu_draw)
),
lower_95 = c(___, ___, ___),
upper_95 = c(___, ___, ___)
)
student_probability_positive <- mean(student_mu_draw > 0)
student_summary
student_probability_positiveReport the requested values beside the corresponding normal-model values.
Replace this sentence with the normal and robust location results.
Define Z_p:\Omega\to\{0,1\} by
Z_p(\omega)= \begin{cases} 1 & \text{if }Y_p(\omega)>0,\\ 0 & \text{if }Y_p(\omega)\le0. \end{cases}
Let \Pi:\Omega\to[0,1] be the population probability of a positive passage difference. Use
Z_p\mid\Pi=\pi\sim\operatorname{Bernoulli}(\pi) \qquad\text{and}\qquad \Pi\sim\operatorname{Beta}(2,2).
The conjugate-prior notes derive the beta–Bernoulli update. Apply that update here.
Let k_{\mathrm{obs}}=\sum_{p=1}^{P}z_p.
The first beta parameter receives the number of ones, and the second receives the number of zeros. Add those counts to the prior parameters (2,2). Complete one line.
Your LaTeX
\Pi\mid\mathbf Z=\mathbf z \sim\operatorname{Beta}(\text{enter the updated parameters in terms of }P\text{ and }k_{\mathrm{obs}}).
Use the code from Exercise 1 to obtain P and k_{\mathrm{obs}}, then substitute them. Complete one line.
\Pi\mid\mathbf Z=\mathbf z \sim\operatorname{Beta}(\text{enter the two numerical parameters}).
Calculate the posterior mean, median, and 95% equal-tailed interval. Define B_\Pi^+=\{\omega\in\Omega\mid\Pi(\omega)>.5\} and calculate \mathbb P(B_\Pi^+\mid\mathbf z) with pbeta(..., lower.tail = FALSE). Verify it with at least 50,000 independent posterior draws and report the Monte Carlo standard error.
passage_count <- nrow(passage_pairs)
positive_count <- sum(passage_pairs$delta_log_ffd > 0)
nonpositive_count <- passage_count - positive_count
alpha_posterior <- 2 + ___
beta_posterior <- 2 + ___
pi_summary <- c(
mean = alpha_posterior / (alpha_posterior + beta_posterior),
median = qbeta(.5, alpha_posterior, beta_posterior),
lower_95 = qbeta(.025, alpha_posterior, beta_posterior),
upper_95 = qbeta(.975, alpha_posterior, beta_posterior)
)
exact_probability <- pbeta(
.5,
alpha_posterior,
beta_posterior,
lower.tail = FALSE
)
set.seed(___)
simulation_count <- 50000
pi_draw <- rbeta(simulation_count, alpha_posterior, beta_posterior)
mc_probability <- mean(pi_draw > .5)
mcse <- sqrt(mc_probability * (1 - mc_probability) / simulation_count)
pi_summary
c(exact_probability = exact_probability, simulation = mc_probability, mcse = mcse)Report k_{\mathrm{obs}}, the posterior summaries, the event probability, and its Monte Carlo standard error.
Replace this sentence with the requested sign-model result.
The data for Exercises 5–6 come from the training split of Universal Dependencies English EWT. The preparation scripts retain sentences containing exactly one selected what, which, or who token annotated as nsubj or obj. Each row of data/ud_wh_dependencies.csv represents one retained sentence. The relevant columns are sentence_id, wh_word, and deprel.
The categorical analysis in PS1 used a conditional null distribution over 2×2 tables. The Fisher’s exact-test notes develop that construction. Here we instead place a beta prior on each scalar row probability and compare the resulting posterior distributions. The calculation does not condition on both observed margins, and its posterior event probability is not Fisher’s p value.
For this exercise, retain what and who. Let N_w be the number of retained sentences with word w, and let K_w be the number annotated obj. Use the two scalar models
K_w\mid\Pi_w=\pi_w\sim\operatorname{Binomial}(N_w,\pi_w) \qquad\text{and}\qquad \Pi_w\sim\operatorname{Beta}(1,1),
where w\in\{\mathrm{what},\mathrm{who}\}.
Load the data, verify that sentence_id is unique, and construct the two-by-two table with rows what, who and columns nsubj, obj.
wh_dependencies <- read.csv("data/ud_wh_dependencies.csv")
str(wh_dependencies)
stopifnot(!anyDuplicated(wh_dependencies$sentence_id))
what_who_data <- wh_dependencies |>
filter(wh_word %in% c("what", "who"))
what_who_table <- with(
what_who_data,
table(
factor(wh_word, levels = c("what", "who")),
factor(deprel, levels = c("nsubj", "obj"))
)
)
row_totals <- rowSums(what_who_table)
object_counts <- what_who_table[, "obj"]
what_who_table
row_totals
object_countsApply the beta–binomial update separately to the two scalar parameters. For each row, the first beta parameter is 1+K_w and the second is 1+N_w-K_w. Substitute the counts from the table. Complete two lines.
Your LaTeX
\begin{aligned} \Pi_{\mathrm{what}}\mid\text{data} &\sim\operatorname{Beta}(\text{enter the parameters}),\\ \Pi_{\mathrm{who}}\mid\text{data} &\sim\operatorname{Beta}(\text{enter the parameters}). \end{aligned}
Draw at least 100,000 independent values from each posterior. Pair draws by simulation index and calculate
\Delta_\Pi=\Pi_{\mathrm{what}}-\Pi_{\mathrm{who}}
and
\operatorname{OR} = \frac{\Pi_{\mathrm{what}}/(1-\Pi_{\mathrm{what}})} {\Pi_{\mathrm{who}}/(1-\Pi_{\mathrm{who}})}.
Calculate the posterior median and 95% interval for both quantities. Define B_\Delta^+=\{\omega\in\Omega\mid\Delta_\Pi(\omega)>0\} and estimate its posterior probability and Monte Carlo standard error.
alpha_what <- 1 + object_counts["what"]
beta_what <- 1 + row_totals["what"] - object_counts["what"]
alpha_who <- 1 + object_counts["who"]
beta_who <- 1 + row_totals["who"] - object_counts["who"]
set.seed(___)
association_draw_count <- 100000
pi_what_draw <- rbeta(association_draw_count, alpha_what, beta_what)
pi_who_draw <- rbeta(association_draw_count, alpha_who, beta_who)
delta_pi_draw <- ___
odds_ratio_draw <- ___
delta_summary <- quantile(delta_pi_draw, c(.025, .5, .975))
odds_ratio_summary <- quantile(odds_ratio_draw, c(.025, .5, .975))
probability_delta_positive <- mean(delta_pi_draw > 0)
delta_mcse <- sqrt(
probability_delta_positive * (1 - probability_delta_positive) /
association_draw_count
)
delta_summary
odds_ratio_summary
c(probability = probability_delta_positive, mcse = delta_mcse)Report the two observed row counts and all requested posterior summaries.
Replace this sentence with the requested two-row association result.
Repeat the calculation with which and who. Define the difference as \Pi_{\mathrm{which}}-\Pi_{\mathrm{who}} and orient the odds ratio with the which odds in the numerator. Use a new pair of scalar beta posteriors rather than a vector-valued parameter.
which_who_data <- wh_dependencies |>
filter(wh_word %in% c("which", "who"))
which_who_table <- with(
which_who_data,
table(
factor(wh_word, levels = c("which", "who")),
factor(deprel, levels = c("nsubj", "obj"))
)
)
which_who_totals <- rowSums(which_who_table)
which_who_objects <- which_who_table[, "obj"]
alpha_which <- 1 + which_who_objects["which"]
beta_which <-
1 + which_who_totals["which"] - which_who_objects["which"]
alpha_who_second <- 1 + which_who_objects["who"]
beta_who_second <-
1 + which_who_totals["who"] - which_who_objects["who"]
set.seed(___)
comparison_draw_count <- 100000
pi_which_draw <- rbeta(
comparison_draw_count,
alpha_which,
beta_which
)
pi_who_second_draw <- rbeta(
comparison_draw_count,
alpha_who_second,
beta_who_second
)
delta_which_who_draw <- ___
odds_ratio_which_who_draw <- ___
quantile(delta_which_who_draw, c(.025, .5, .975))
quantile(odds_ratio_which_who_draw, c(.025, .5, .975))
mean(delta_which_who_draw > 0)Put the what–who and which–who posterior summaries in one table.
Replace this sentence with the requested comparison of the two associations.
Treat the retained what and who sentences as a finite population. For each n\in\{20,40,80,160\}, draw 500 simple random samples without replacement. For each sample, update the two beta posteriors and record the posterior mean of \Delta_\Pi and its 95% interval width. Use one recorded seed.
Use 5,000 paired posterior draws inside each replicate to estimate the interval width. The posterior mean of \Delta_\Pi may be calculated directly as the difference between the two beta means. A sample may omit one word class. The beta prior still yields a proper posterior in that case: use object count zero and row total zero for the omitted class.
what_who_population <- wh_dependencies |>
filter(wh_word %in% c("what", "who"))
analyze_sentence_sample <- function(
sample_size,
replicate_id,
posterior_draw_count = 5000
) {
sampled_sentences <- what_who_population |>
slice_sample(n = sample_size, replace = FALSE)
sampled_table <- with(
sampled_sentences,
table(
factor(wh_word, levels = c("what", "who")),
factor(deprel, levels = c("nsubj", "obj"))
)
)
sampled_totals <- rowSums(sampled_table)
sampled_objects <- sampled_table[, "obj"]
alpha_what <- 1 + sampled_objects["what"]
beta_what <- 1 + sampled_totals["what"] - sampled_objects["what"]
alpha_who <- 1 + sampled_objects["who"]
beta_who <- 1 + sampled_totals["who"] - sampled_objects["who"]
posterior_mean_delta <-
alpha_what / (alpha_what + beta_what) -
alpha_who / (alpha_who + beta_who)
pi_what <- rbeta(posterior_draw_count, alpha_what, beta_what)
pi_who <- rbeta(posterior_draw_count, alpha_who, beta_who)
delta_draw <- ___
delta_interval <- quantile(delta_draw, c(.025, .975))
tibble(
sample_size = sample_size,
replicate_id = replicate_id,
posterior_mean_delta = posterior_mean_delta,
interval_width = ___
)
}
set.seed(___)
sampling_plan <- expand_grid(
sample_size = c(20, 40, 80, 160),
replicate_id = seq_len(500)
)
sampling_results <- bind_rows(Map(
analyze_sentence_sample,
sampling_plan$sample_size,
sampling_plan$replicate_id
))For each n, calculate the mean and standard deviation of the posterior-mean estimates and the median interval width. Plot the empirical sampling distributions of the posterior means by n.
sampling_summary <- sampling_results |>
group_by(sample_size) |>
summarize(
mean_posterior_mean = mean(posterior_mean_delta),
sd_posterior_mean = sd(posterior_mean_delta),
median_interval_width = median(interval_width),
.groups = "drop"
)
sampling_summary
ggplot(
sampling_results,
aes(x = factor(sample_size), y = posterior_mean_delta)
) +
geom_boxplot() +
labs(
x = "Sentences sampled",
y = "Posterior mean of the object-probability difference"
)Report the sample-size summary table.
Replace this sentence with the sample-size summary table.
Use all three word forms. Let g\in\{\mathrm{what},\mathrm{which},\mathrm{who}\} index rows, N_g the row total, and K_g the object count. Under independence, all rows share the scalar object probability \Pi_0:
K_g\mid\Pi_0=\pi_0 \sim \operatorname{Binomial}(N_g,\pi_0), \qquad \Pi_0\sim\operatorname{Beta}(1,1).
A posterior predictive distribution is obtained by drawing \Pi_0 from its posterior and then drawing a replicated table from the likelihood.
Construct the observed three-by-two table in row order what, which, who and column order nsubj, obj. Calculate its expected counts under independence and its Pearson statistic using the formula from the chi-squared-test notes.
three_word_table <- with(
wh_dependencies,
table(
factor(wh_word, levels = c("what", "which", "who")),
factor(deprel, levels = c("nsubj", "obj"))
)
)
row_totals <- rowSums(three_word_table)
column_totals <- colSums(three_word_table)
grand_total <- sum(three_word_table)
expected_counts <- outer(row_totals, column_totals) / grand_total
pearson_contributions <- ___
x_squared_observed <- sum(pearson_contributions)
three_word_table
expected_counts
x_squared_observedGenerate at least 100,000 posterior predictive tables. For each replicate, draw one value of \Pi_0, draw one object count for each fixed row total, construct the replicated table, and calculate its Pearson statistic. The function below returns NA_real_ if a replicated column total is zero, because the corresponding expected counts would also be zero. Such a replicate is extraordinarily unlikely here, but the rule makes the calculation defined.
object_counts <- three_word_table[, "obj"]
alpha_zero <- 1 + sum(object_counts)
beta_zero <- 1 + sum(row_totals - object_counts)
set.seed(___)
predictive_draw_count <- 100000
pi_zero_draw <- rbeta(predictive_draw_count, alpha_zero, beta_zero)
replicated_statistic <- function(pi_zero) {
replicated_objects <- rbinom(
length(row_totals),
size = row_totals,
prob = pi_zero
)
replicated_table <- cbind(
nsubj = row_totals - replicated_objects,
obj = replicated_objects
)
replicated_columns <- colSums(replicated_table)
if (any(replicated_columns == 0)) {
return(NA_real_)
}
replicated_expected <- outer(
rowSums(replicated_table),
replicated_columns
) / sum(replicated_table)
sum(___)
}
replicated_x_squared <- vapply(
pi_zero_draw,
replicated_statistic,
numeric(1)
)
replicated_x_squared <- replicated_x_squared[!is.na(replicated_x_squared)]Define the posterior predictive tail event in full set-builder notation.
The event contains replicated statistics at least as large as the observed statistic. Complete one line.
Your LaTeX
B_{X^2}^{\mathrm{rep}} =\{\omega\in\Omega\mid\text{enter the discrepancy condition}\}.
Let M be the number of retained replications. Average an indicator of the event over the M statistics. Complete one line.
\mathbb P(B_{X^2}^{\mathrm{rep}}\mid\mathbf K=\mathbf k) \approx\text{enter the Monte Carlo average}.
Plot the replicated statistics with the observed statistic marked. Calculate the posterior predictive tail probability and its Monte Carlo standard error. Calculate the observed Pearson residuals and identify the two cells with the largest absolute residuals.
ggplot(tibble(replicated_x_squared), aes(x = replicated_x_squared)) +
geom_histogram(bins = 50) +
geom_vline(xintercept = x_squared_observed, color = "red")
tail_indicator <- replicated_x_squared >= x_squared_observed
tail_probability <- mean(tail_indicator)
tail_mcse <- sqrt(
tail_probability * (1 - tail_probability) /
length(tail_indicator)
)
pearson_residuals <- ___
residual_table <- as.data.frame(as.table(pearson_residuals)) |>
mutate(absolute_residual = abs(Freq)) |>
arrange(desc(absolute_residual))
c(probability = tail_probability, mcse = tail_mcse)
slice_head(residual_table, n = 2)Report the observed statistic, posterior predictive tail probability, Monte Carlo standard error, and two largest absolute residuals with their signs.
Replace this sentence with the requested posterior predictive result.
Text_ID rather than aligned by row position.