--- title: "STK4021/9021 Compulsory Exercise 2026 - Worked solution" author: "Odd Kolbjørnsen" date: "`r format(Sys.Date(), '%d %B %Y')`" output: html_document: toc: true toc_depth: 3 number_sections: false df_print: paged fontsize: 18px; pdf_document: toc: true number_sections: false fontsize: 18px header-includes: - \usepackage{amsmath} - \usepackage{booktabs} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE, fig.width=7, fig.height=4.5, fig.align="center") # Packages used only for simulation and stable multivariate-normal calculations. if (!requireNamespace("MASS", quietly=TRUE)) stop("Install package MASS") # Text files must be in the same directory as this Rmd file. set.seed(4021) A <- as.matrix(read.table("A.txt", header=FALSE)) d00 <- scan("seis00.txt", quiet=TRUE) d01 <- scan("seis01.txt", quiet=TRUE) d02 <- scan("seis02.txt", quiet=TRUE) n <- nrow(A) idx <- seq_len(n) mu <- 0 sigma_m2 <- 1 ell <- 20 sigma_e2 <- 0.04 Delta_mu <- 5 R <- outer(idx, idx, function(i,j) exp(-3*abs(i-j)/ell)) Sigma_m <- sigma_m2 * R Sigma_e <- sigma_e2 * diag(n) mu0 <- rep(mu,n) # Numerically stable helpers: avoid forming matrix inverses explicitly. logdet_spd <- function(S) { 2*sum(log(diag(chol(S)))) } log_mvn <- function(x, mean, Sigma) { r <- as.vector(x - mean) quadratic_form <- as.numeric(crossprod(r, solve(Sigma, r))) -0.5 *( length(x)*log(2*pi) + logdet_spd(Sigma) + quadratic_form ) } normalize_logweights <- function(logw) { logw <- as.numeric(logw) z <- exp(logw - max(logw)) z / sum(z) } gaussian_update <- function(d, A, mu_pri_m, Sigma_pri_m, Sigma_lik_e) { S_d <- A %*% Sigma_pri_m %*% t(A) + Sigma_lik_e K <- Sigma_pri_m %*% t(A) %*% solve(S_d) mu_post_m <- as.vector( mu_pri_m + K %*% (d - A %*% mu_pri_m)) Sigma_post_m <- Sigma_pri_m - K %*% A %*% Sigma_pri_m # Ensure numerical symmetry Sigma_post_m <- (Sigma_post_m + t(Sigma_post_m)) / 2 list( mu_post_m = mu_post_m, Sigma_post_m = Sigma_post_m ) } zone_mean <- function(T,B, base_mu=mu, drop=Delta_mu) { z <- rep(base_mu,n) z[T:B] <- base_mu-drop z } # Grid of all admissible non-zero zones. zone_grid <- expand.grid(T=30:60, B=40:70) zone_grid <- zone_grid[zone_grid$T <= zone_grid$B,] zone_grid$V <- zone_grid$B-zone_grid$T+1 rownames(zone_grid) <- NULL zone_posterior_fixed <- function(d, A,Sigma_m, Sigma_e,prior=NULL) { Sd <- A %*% Sigma_m%*% t(A)+Sigma_e if (is.null(prior)) prior <- rep(1/nrow(zone_grid),nrow(zone_grid)) loglik <- vapply(seq_len(nrow(zone_grid)), function(k) { mm <- zone_mean(zone_grid$T[k],zone_grid$B[k]) log_mvn(d,as.vector(A%*%mm),Sd) }, numeric(1)) p <- normalize_logweights(log(prior)+loglik) transform(zone_grid, loglik=loglik, posterior=p) } plot_zone_heatmap <- function(tab, main) { Z <- matrix(NA_real_,31,31,dimnames=list(T=30:60,B=40:70)) Z[cbind(tab$T-29,tab$B-39)] <- tab$posterior maxLim = max(log(tab$posterior))-0.5 minLim = maxLim-4.5 image(30:60,40:70,Z,xlab="Top T",ylab="Base B",main=main, col=hcl.colors(40,"YlOrRd",rev=TRUE)) contour(30:60,40:70,Z,add=TRUE,drawlabels=FALSE,levels=exp(seq(minLim, maxLim, by = 2))) } summarize_zone <- function(tab) { pT <- rowsum(tab$posterior,tab$T) pB <- rowsum(tab$posterior,tab$B) pV <- rowsum(tab$posterior,tab$V) list(T=pT,B=pB,V=pV, means=c(T=sum(tab$T*tab$posterior), B=sum(tab$B*tab$posterior), V=sum(tab$V*tab$posterior)), MAP=tab[which.max(tab$posterior),c("T","B","V","posterior")]) } ``` # Part 1: Gaussian Bayesian inversion and simulation ## Part 1(a): Joint and posterior distributions The model is \[ \mathbf d=\mathbf A\mathbf m+\boldsymbol\varepsilon,\qquad \mathbf m\sim N(\boldsymbol\mu_m,\boldsymbol\Sigma_m),\qquad \boldsymbol\varepsilon\sim N(\mathbf0,\boldsymbol\Sigma_e), \] where $\mathbf m$ and $\boldsymbol\varepsilon$ are independent. Since both vectors are multivariate Gaussian and independent, they form the joint multivariate Gaussian vector \[ \begin{pmatrix}\mathbf m\\\boldsymbol\varepsilon\end{pmatrix} \sim N\left( \begin{pmatrix}\boldsymbol\mu_m\\\mathbf 0\end{pmatrix}, \begin{pmatrix}\boldsymbol\Sigma_m & \mathbf 0\\ \mathbf 0 & \boldsymbol\Sigma_e\end{pmatrix} \right). \] The required vector is the linear transformation \[ \begin{pmatrix}\mathbf m\\\mathbf d\end{pmatrix} = \underbrace{\begin{pmatrix}\mathbf I & \mathbf 0\\ \mathbf A & \mathbf I\end{pmatrix}}_{\mathbf H} \begin{pmatrix}\mathbf m\\\boldsymbol\varepsilon\end{pmatrix}. \] A linear transformation of a multivariate Gaussian vector is again multivariate Gaussian. Its mean and covariance are therefore obtained as $\mathbf H\boldsymbol\mu$ and $\mathbf H\boldsymbol\Sigma\mathbf H^T$. Hence \[ \boldsymbol\mu_d=E(\mathbf d)=\mathbf A\boldsymbol\mu_m, \quad \boldsymbol\Sigma_d=\operatorname{Var}(\mathbf d)=\mathbf A\boldsymbol\Sigma_m\mathbf A^T+ \boldsymbol\Sigma_e, \] and \[ \operatorname{Cov}(\mathbf d,\mathbf m)=\mathbf A\boldsymbol\Sigma_m, \qquad \operatorname{Cov}(\mathbf m,\mathbf d)=\boldsymbol\Sigma_m\mathbf A^T. \] Therefore \[ \begin{pmatrix}\mathbf m\\\mathbf d\end{pmatrix} \sim N\left[ \begin{pmatrix}\boldsymbol\mu_m\\\mathbf A\boldsymbol\mu_m\end{pmatrix}, \begin{pmatrix} \boldsymbol\Sigma_m & \boldsymbol\Sigma_m\mathbf A^T\\ \mathbf A\boldsymbol\Sigma_m & \boldsymbol\Sigma_d \end{pmatrix}\right]. \] Conditioning a joint Gaussian gives \[ \mathbf m\mid\mathbf d\sim N(\boldsymbol\mu_{m\mid d}, \boldsymbol\Sigma_{m\mid d}), \] with \[ \boldsymbol\mu_{m\mid d}=\boldsymbol\mu_m+ \boldsymbol\Sigma_m\mathbf A^T\boldsymbol\Sigma_d^{-1} (\mathbf d-\mathbf A\boldsymbol\mu_m), \] \[ \boldsymbol\Sigma_{m\mid d}=\boldsymbol\Sigma_m- \boldsymbol\Sigma_m\mathbf A^T\boldsymbol\Sigma_d^{-1} \mathbf A\boldsymbol\Sigma_m. \] ## Part 1(b): Prior summaries and simulations Here $\mu_{m,i}=0$ and \[ (\Sigma_m)_{ij}=\exp\{-3|i-j|/20\}. \] The pointwise 95% prior interval is $0\pm1.96$, but the components are strongly correlated at short vertical separation. ```{r prior, fig.width=8, fig.height=5} prior_sd <- sqrt(diag(Sigma_m)) prior_lower <- mu0 - 1.96 * prior_sd prior_upper <- mu0 + 1.96 * prior_sd mprior <- MASS::mvrnorm(5, mu0, Sigma_m) plot(idx, mu0, type="n", ylim=range(prior_lower, prior_upper, mprior), xlab="Profile component i", ylab=expression(m[i]), main="Prior mean, pointwise 95% interval, and five simulations") polygon(c(idx, rev(idx)), c(prior_lower, rev(prior_upper)), col=adjustcolor("steelblue", .25), border=NA) matlines(idx, t(mprior), lty=1, lwd=1, col=hcl.colors(5, "Dark 3")) lines(idx, mu0, lwd=2, col="black") legend("topright", c("Prior mean", "Pointwise 95% interval", "Simulations"), col=c("black", adjustcolor("steelblue", .5), hcl.colors(5, "Dark 3")[1]), lwd=c(2, 8, 1), bty="n", cex=.85) ``` ## Part 1(c): Posterior for `seis00.txt` ```{r posterior00, fig.width=8, fig.height=5} post00 <- gaussian_update(d00, A, mu0, Sigma_m, Sigma_e) psd <- sqrt(pmax(diag(post00$Sigma_post_m), 0)) post_lower <- post00$mu_post_m - 1.96 * psd post_upper <- post00$mu_post_m + 1.96 * psd mpost <- MASS::mvrnorm(5, post00$mu_post_m, post00$Sigma_post_m) plot(idx, post00$mu_post_m, type="n", ylim=range(post_lower, post_upper, mpost), xlab="Profile component i", ylab=expression(m[i]), main="Posterior mean, pointwise 95% interval, and five simulations") polygon(c(idx, rev(idx)), c(post_lower, rev(post_upper)), col=adjustcolor("tomato", .25), border=NA) matlines(idx, t(mpost), lty=1, lwd=1, col=hcl.colors(5, "Dark 3")) lines(idx, post00$mu_post_m, lwd=2, col="black") legend("topright", c("Posterior mean", "Pointwise 95% interval", "Simulations"), col=c("black", adjustcolor("tomato", .5), hcl.colors(5, "Dark 3")[1]), lwd=c(2, 8, 1), bty="n", cex=.85) ``` ```{r covariance-images, fig.width=8, fig.height=7} plot_covariance <- function(S, main, zlim, cols) { layout(matrix(c(1, 2), nrow=1), widths=c(5, 1)) par(mar=c(4, 4, 3, 1)) image(idx, idx, S, xlab="Component i", ylab="Component j", main=main, col=cols, zlim=zlim, asp=1, useRaster=TRUE) par(mar=c(4, 1, 3, 4)) z <- seq(zlim[1], zlim[2], length.out=length(cols)+1) image(1, z, matrix(z, nrow=1), col=cols, axes=FALSE, xlab="", ylab="", zlim=zlim, useRaster=TRUE) axis(4, las=1, cex.axis=.8) mtext("Covariance", side=4, line=2.5) layout(1) } cov_zlim <- range(Sigma_m, post00$Sigma_post_m) cov_cols <- hcl.colors(60, "Blue-Red 3") plot_covariance(Sigma_m, "Prior covariance", cov_zlim, cov_cols) plot_covariance(post00$Sigma_post_m, "Posterior covariance", cov_zlim, cov_cols) ``` The posterior covariance is smaller than the prior covariance in data-informed directions. It need not be uniformly small because the forward operator smooths the profile and is weakly sensitive to some directions. ## Part 1(d): Posterior distribution of exceedance volume For $c=-1$, \[ V=h(\mathbf m)=\sum_i |m_i-c|I(m_i=abs(r_obs)))/(length(r_perm)+1) hist(r_perm,breaks=30,xlab="Lag-one correlation",main="Permutation null distribution",xlim=c(-1,1)) abline(v=r_obs,col="red",lwd=3) c(observed=r_obs,p_value=p_perm) ``` ```{r permutation more samples} maxCor=cor(sort(d00)[-length(d00)],sort(d00)[-1]) r_perm <- replicate(10000,lag1(sample(d00,replace = FALSE))) # We include the observed permutation as one case. # This makes us avoid estimated p-value of 0 p_perm <- (1+sum(abs(r_perm)>=abs(r_obs)))/(length(r_perm)+1) hist(r_perm,breaks=30,xlab="Lag-one correlation",main="Permutation null distribution",xlim=c(-1,1)) abline(v=r_obs,col="red",lwd=3) c(observed=r_obs,p_value=p_perm,max_r=maxCor) ``` A small permutation p-value rejects exchangeability: the position of an observation in the trace contains information through serial dependence. # Part 2: Model selection and decision making For a zone $(T,B)$, define $\boldsymbol\mu_m^{T,B}$ by reducing components $T,\ldots,B$ by $\Delta\mu=5$. The covariance is unchanged. Thus \[ \mathbf d\mid T,B\sim N(\mathbf A\boldsymbol\mu_m^{T,B},\boldsymbol\Sigma_d), \quad \boldsymbol\Sigma_d=\mathbf A\boldsymbol\Sigma_m\mathbf A^T+\boldsymbol\Sigma_e. \] For a discrete prior $p(T,B)$, \[ p(T,B\mid\mathbf d)=\frac{p(\mathbf d\mid T,B)p(T,B)} {\sum_{(t,b)\in\mathcal G}p(\mathbf d\mid t,b)p(t,b)}. \] ## Part 2(a): Posterior top, base and thickness We first inspect all three observed traces on common axes. ```{r observed-traces, fig.width=8, fig.height=5} trace_cols <- c("black", "steelblue", "tomato") matplot(idx, cbind(d00,d01,d02), type="l", lty=1, lwd=1.5, col=trace_cols, xlab="Trace component i", ylab=expression(d[i]), main="Observed seismic traces") legend("topright", c("seis00","seis01","seis02"), col=trace_cols, lty=1, lwd=1.5, bty="n") ``` For every admissible $(T,B)$, the stored `loglik` is \[ \log p(\mathbf d\mid T,B)=\log\phi_n\!\left( \mathbf d;\mathbf A\boldsymbol\mu_m^{T,B}, \mathbf A\boldsymbol\Sigma_m\mathbf A^T+\boldsymbol\Sigma_e\right). \] With a uniform prior on the admissible pairs, normalizing these marginal likelihoods gives the posterior over $(T,B)$. ```{r zones, fig.width=10, fig.height=5} z01 <- zone_posterior_fixed(d01,A,Sigma_m,Sigma_e) z02 <- zone_posterior_fixed(d02,A,Sigma_m,Sigma_e) plot_zone_loglik <- function(tab, main) { Z <- matrix(NA_real_,31,31,dimnames=list(T=30:60,B=40:70)) Z[cbind(tab$T-29,tab$B-39)] <- tab$loglik filled.contour(30:60,40:70,Z, color.palette=function(n) hcl.colors(n,"YlOrRd",rev=TRUE), xlab="Top T",ylab="Base B",main=main, plot.axes={axis(1);axis(2);contour(30:60,40:70,Z,add=TRUE, drawlabels=TRUE,col="grey25",lwd=.7)}, key.title=title(main="log p(d|T,B)",cex.main=.8)) } plot_zone_loglik(z01,"seis01: log marginal likelihood") plot_zone_loglik(z02,"seis02: log marginal likelihood") par(mfrow=c(1,2)) plot_zone_heatmap(z01,"seis01: p(T,B | d)") plot_zone_heatmap(z02,"seis02: p(T,B | d)") par(mfrow=c(1,1)) s01 <- summarize_zone(z01); s02 <- summarize_zone(z02) s01$means; s01$MAP s02$means; s02$MAP plot_marginals <- function(s,main) { par(mfrow=c(1,3)) plot(as.numeric(rownames(s$T)),s$T[,1],type="h",xlab="T",ylab="Probability",main=paste(main,"top")) plot(as.numeric(rownames(s$B)),s$B[,1],type="h",xlab="B",ylab="Probability",main=paste(main,"base")) plot(as.numeric(rownames(s$V)),s$V[,1],type="h",xlab="V",ylab="Probability",main=paste(main,"thickness")) par(mfrow=c(1,1)) } plot_marginals(s01,"seis01") plot_marginals(s02,"seis02") ``` The marginal prior is not uniform for $T$, $B$, or $V$: different values occur in different numbers of admissible pairs. ## Part 2(b): Including the zero-thickness model It would be inappropriate to assign the no-zone model the same prior mass as each individual $(T,B)$ pair: the non-zero hypothesis is represented by many models and would then receive a large prior probability merely because it was discretized finely. This is a point of discussion with subject matter experts. Perhaps equal chance of zone vs no zone. Which gives the prior odds between no zone and some zone equal to one. \[ P(M_0)=0.5,\qquad P(M_1)=0.5, \qquad P(T,B\mid M_1)=1/|\mathcal G|. \] ```{r zero_model} zero_posterior <- function(d,A,Sigma_m,Sigma_e,zone_grid) { ## Sd <- A%*%Sigma_m%*%t(A)+Sigma_e ll0 <- log_mvn(d,as.vector(A%*%mu0),Sd) ll1 <- vapply(seq_len(nrow(zone_grid)),function(k) log_mvn(d,as.vector(A%*%zone_mean(zone_grid$T[k],zone_grid$B[k])),Sd),numeric(1)) logprior <- c(log(.5),rep(log(.5/nrow(zone_grid)),nrow(zone_grid))) pp <- normalize_logweights(logprior+c(ll0,ll1)) list(p_zero=pp[1],zone=transform(zone_grid,posterior=pp[-1]),all=pp) } zm01 <- zero_posterior(d01,A,Sigma_m,Sigma_e,zone_grid) zm02 <- zero_posterior(d02,A,Sigma_m,Sigma_e,zone_grid) c(seis01=zm01$p_zero,seis02=zm02$p_zero) ``` ## Part 2(c): Investment decision Let action $a=1$ mean invest and $a=0$ mean do not invest. Profit is \[ U(1,V)=V-10,\qquad U(0,V)=0, \] and loss may be defined as $L(a,V)=-U(a,V)$. The Bayes action invests if \[ E[V-10\mid\mathbf d]>0\quad\Longleftrightarrow\quad E[V\mid\mathbf d]>10. \] The zero model has $V=0$ and must be included in this expectation. ```{r decision} decision_summary <- function(zm) { EV <- sum(zm$zone$V*zm$zone$posterior) c(P_zero=zm$p_zero,E_V=EV,E_profit_invest=EV-10, invest=as.numeric(EV>10)) } rbind(seis01=decision_summary(zm01),seis02=decision_summary(zm02)) ``` For communication to decision makers, we can report the posterior probability of no zone, expected thickness (if zone), probability that thickness exceeds the break-even value 10, and expected profit from investing, in addition to the binary recommendation. ```{r decision_uncertainty} decision_risk <- function(zm) { p_gt10 <- sum(zm$zone$posterior[zm$zone$V>10]) c(P_zero=zm$p_zero,P_V_gt_10=p_gt10, expected_thickness=sum(zm$zone$V*zm$zone$posterior)) } rbind(seis01=decision_risk(zm01),seis02=decision_risk(zm02)) ``` # Part 3: Uncertain parameters ## Part 3(a): Uncertain common mean Write \[ \mathbf m=\mu\mathbf1+\boldsymbol\eta,\qquad \mu\sim N(\mu_0,\tau_\mu^2),\qquad \boldsymbol\eta\sim N(\mathbf0,\boldsymbol\Sigma_m), \] independently. Conditional on $\mu$, $E(\mathbf m\mid\mu)=\mu\mathbf1$ and $\operatorname{Var}(\mathbf m\mid\mu)=\boldsymbol\Sigma_m$. Hence \[ E(\mathbf m)=\mu_0\mathbf1, \] \[ \begin{aligned} \operatorname{Var}(\mathbf m) &=E\{\operatorname{Var}(\mathbf m\mid\mu)\} +\operatorname{Var}\{E(\mathbf m\mid\mu)\}\\ &=\boldsymbol\Sigma_m+\tau_\mu^2\mathbf1\mathbf1^T. \end{aligned} \] For the data, \[ E(\mathbf d)=\mathbf A\mu_0\mathbf1, \qquad \operatorname{Var}(\mathbf d)= \mathbf A\operatorname{Var}(\mathbf m)\mathbf A^T+\boldsymbol\Sigma_e. \] Thus the fixed-mean and uncertain-mean covariances are \[ \begin{array}{ll} \text{fixed mean:}& \operatorname{Var}(\mathbf m)=\boldsymbol\Sigma_m,\quad \operatorname{Var}(\mathbf d)=\mathbf A\boldsymbol\Sigma_m\mathbf A^T+\boldsymbol\Sigma_e,\\[1mm] \text{uncertain mean:}& \operatorname{Var}(\mathbf m)=\boldsymbol\Sigma_m+\tau_\mu^2\mathbf1\mathbf1^T,\\ &\operatorname{Var}(\mathbf d)=\mathbf A(\boldsymbol\Sigma_m+ \tau_\mu^2\mathbf1\mathbf1^T)\mathbf A^T+\boldsymbol\Sigma_e. \end{array} \] ```{r uncertain-mean-covariances, fig.width=10, fig.height=9} tau_mu2 <- .25 one_n <- rep(1,n) Sigma_m_fixedmean <- Sigma_m Sigma_m_uncmean <- Sigma_m + tau_mu2*tcrossprod(one_n) Sigma_d_fixedmean <- A%*%Sigma_m_fixedmean%*%t(A)+Sigma_e Sigma_d_uncmean <- A%*%Sigma_m_uncmean%*%t(A)+Sigma_e plot_cov_panel <- function(S,main,zlim) image(idx,idx,S,xlab="Component i", ylab="Component j",main=main,col=hcl.colors(60,"Blue-Red 3"), zlim=zlim,asp=1,useRaster=TRUE) par(mfrow=c(2,2),mar=c(4,4,3,1)) plot_cov_panel(Sigma_m_fixedmean,"Cov(m): fixed mean",range(Sigma_m_fixedmean,Sigma_m_uncmean)) plot_cov_panel(Sigma_m_uncmean,"Cov(m): uncertain mean",range(Sigma_m_fixedmean,Sigma_m_uncmean)) plot_cov_panel(Sigma_d_fixedmean,"Cov(d): fixed mean",range(Sigma_d_fixedmean,Sigma_d_uncmean)) plot_cov_panel(Sigma_d_uncmean,"Cov(d): uncertain mean",range(Sigma_d_fixedmean,Sigma_d_uncmean)) par(mfrow=c(1,1)) rbind(m=c(fixed=mean(diag(Sigma_m_fixedmean)),uncertain=mean(diag(Sigma_m_uncmean))), d=c(fixed=mean(diag(Sigma_d_fixedmean)),uncertain=mean(diag(Sigma_d_uncmean)))) plot(one_n, type = "l", lwd = 2,col = "blue", ylim = c(-0.2, 1.2), xlab = "Index",ylab = "Value", main = expression(paste("Comparison of ", 1, " and A%*%1"))) lines(A %*% one_n,col = "red",lwd = 2,lty = 2) legend( "topright", legend = c("one_n", "A %*% one_n"), col =c("blue", "red"), lwd = 2, lty = c(1, 2), bty = "n") ``` Uncertainty in the mean adds the rank-one term $\tau_\mu^2\mathbf1\mathbf1^T$ to the prior covariance of $\mathbf m$ and the rank-one term $\tau_\mu^2(\mathbf A\mathbf1)(\mathbf A\mathbf1)^T$ to the marginal covariance of $\mathbf d$. Note taht the variance of the data $\Sigma_d$ do not change, even though the variance of the parameter $\mathbf m$ increases. This is because $\mathbf A\mathbf1 = \mathbf0$ ```{r uncertain_mean} z01_uncmean <- zone_posterior_fixed(d01,A,Sigma_m_uncmean,Sigma_e) par(mfrow=c(1,2)); plot_zone_heatmap(z01,"Fixed common mean") plot_zone_heatmap(z01_uncmean,"Uncertain common mean"); par(mfrow=c(1,1)) rbind(fixed=summarize_zone(z01)$MAP, uncertain_mean=summarize_zone(z01_uncmean)$MAP) ``` The additional coherent background variation do not change the data distribution, so there will not be a change in the marginal likelihood $p(\mathbf d |T,B)$. ## Part 3(b): Uncertain variance scale The exercise uses \[ \boldsymbol\Sigma_m=\sigma_m^2\mathbf R, \qquad \sigma_e^2=\nu\sigma_m^2, \qquad \nu=0.04. \] Thus $\nu=\sigma_e^2/\sigma_m^2$ is known. No new $\kappa$ is needed. Let $s=\sigma_m^2$. Conditional on $s$ and mean model $M$, \[ \mathbf d\mid s,M\sim N_n(\mathbf A\boldsymbol\mu_m^M,s\mathbf Q), \qquad \mathbf Q=\mathbf A\mathbf R\mathbf A^T+\nu\mathbf I. \] Choose the shape-scale inverse-gamma prior \[ s\sim\operatorname{IG}(a_0,b_0),\qquad p(s)=\frac{b_0^{a_0}}{\Gamma(a_0)}s^{-(a_0+1)}e^{-b_0/s}. \] Its mean is $b_0/(a_0-1)$. We use $a_0=2.5$, $b_0=1.5$, giving prior mean 1. With $\mathbf r_M=\mathbf d-\mathbf A\boldsymbol\mu_m^M$, \[ p(\mathbf d\mid s,M)\propto |\mathbf Q|^{-1/2}s^{-n/2} \exp\left\{-\frac{\frac12\mathbf r_M^T\mathbf Q^{-1}\mathbf r_M}{s}\right\}. \] Multiplication by the prior gives \[ p(s\mid\mathbf d,M)\propto s^{-(a_0+n/2+1)} \exp\left\{-\frac{b_0+\frac12\mathbf r_M^T\mathbf Q^{-1}\mathbf r_M}{s}\right\}. \] Therefore \[ \boxed{s\mid\mathbf d,M\sim\operatorname{IG}(a_M,b_M)},\qquad a_M=a_0+n/2, \quad b_M=b_0+\tfrac12\mathbf r_M^T\mathbf Q^{-1}\mathbf r_M. \] ```{r scale_posterior, fig.width=8, fig.height=5} # Default # a0 <- 2.5; b0 <- 1.5; nu <- .04 # to see effect of prior and integration a0 <- 11; b0 <- 30; nu <- .04 Q <- A%*%R%*%t(A)+nu*diag(n) ig_density <- function(x,a,b) b^a/gamma(a)*x^(-a-1)*exp(-b/x) scale_posterior <- function(d,A,a0,b0,Q,mean_m=mu0) { r <- as.vector(d-A%*%mean_m) an <- a0+n/2 bn <- b0+as.numeric(crossprod(r,solve(Q,r)))/2 c(shape=an,scale=bn,mean=bn/(an-1),mode=bn/(an+1), q025=1/qgamma(.975,shape=an,rate=bn), q975=1/qgamma(.025,shape=an,rate=bn)) } # example constant mean (M_0) and data d01 from file seis01.txt # the Posterior will change with Model and data. scale01 <- scale_posterior(d01,A,a0,b0,Q,mu0) scale01 xmax <- max(1/qgamma(.001,a0,rate=b0), 1/qgamma(.001,scale01["shape"],rate=scale01["scale"])) x <- seq(1e-4,xmax,length.out=1500) y0 <- ig_density(x,a0,b0) y1 <- ig_density(x,scale01["shape"],scale01["scale"]) plot(x,y0,type="l",lwd=2,col="steelblue",xlab=expression(sigma[m]^2), ylab="Density",ylim=c(0,max(y0,y1)),main=expression("Prior and posterior for "*sigma[m]^2)) lines(x,y1,lwd=2,col="tomato") abline(v=c(b0/(a0-1),scale01["mean"]),col=c("steelblue","tomato"),lty=2) legend("topright",c("Prior","Posterior","Prior mean","Posterior mean"), col=c("steelblue","tomato","steelblue","tomato"),lty=c(1,1,2,2),lwd=c(2,2,1,1),bty="n") ``` The prior scale $b_0$ is common to every candidate. However, the posterior scale \[ b_M=b_0+\tfrac12\mathbf r_M^T\mathbf Q^{-1}\mathbf r_M \] changes with $(T,B)$ because the fitted mean, and therefore the residual, changes. ```{r integrated_scale, fig.width=10, fig.height=5} zone_posterior_scale_integrated <- function(d,a0,b0,A,Q,zone_grid) { # computes b_M for each model, bm-> b_model # also compute the marginal distribution for 3c, i.e. lw-> loglik # also compute the normalized posterior. given presence. an <- a0+n/2 bm <- vapply(seq_len(nrow(zone_grid)),function(k) { mm <- zone_mean(zone_grid$T[k],zone_grid$B[k]) r <- as.vector(d-A%*%mm) b0+as.numeric(crossprod(r,solve(Q,r)))/2 },numeric(1)) lw <- -an*log(bm) transform(zone_grid,b_model=bm,loglik=lw, posterior=normalize_logweights(lw)) } z01_scale <- zone_posterior_scale_integrated(d01,a0,b0,A,Q,zone_grid) BZ <- matrix(NA_real_,31,31,dimnames=list(T=30:60,B=40:70)) BZ[cbind(z01_scale$T-29,z01_scale$B-39)] <- z01_scale$b_model filled.contour(30:60,40:70,BZ, color.palette=function(n) hcl.colors(n,"Viridis"), xlab="Top T",ylab="Base B",main=expression("Posterior scale "*b[M]*" across models, data set 1"), plot.axes={axis(1);axis(2);contour(30:60,40:70,BZ,add=TRUE, drawlabels=TRUE,col="white",lwd=.7)}, key.title=title(main=expression(b[M]),cex.main=.9)) z02_scale <- zone_posterior_scale_integrated(d02,a0,b0,A,Q,zone_grid) BZ2 <- matrix(NA_real_,31,31,dimnames=list(T=30:60,B=40:70)) BZ2[cbind(z02_scale$T-29,z02_scale$B-39)] <- z02_scale$b_model filled.contour(30:60,40:70,BZ2, color.palette=function(n) hcl.colors(n,"Viridis"), xlab="Top T",ylab="Base B",main=expression("Posterior scale "*b[M]*" across models, data set 2"), plot.axes={axis(1);axis(2);contour(30:60,40:70,BZ2,add=TRUE, drawlabels=TRUE,col="white",lwd=.7)}, key.title=title(main=expression(b[M]),cex.main=.9)) ``` The $b_M$ contour illustrates directly that $b_M$ is model-dependent. Smaller $b_M$ means a smaller generalized residual sum of squares and a larger marginal likelihood. Integrating the scale replaces the fixed-scale Gaussian likelihood by a Student likelihood and generally reduces the penalty for moderately large residuals. ## Part 3(c): Integrating scale uncertainty in the zone posterior For $M=(T,B)$, integrate $s$ from the normal--inverse-gamma model: \[ \begin{aligned} p(\mathbf d\mid M) &=\int p(\mathbf d\mid s,M)p(s)\,ds\\ &=\frac{\Gamma(a_0+n/2)}{\Gamma(a_0)(2\pi b_0)^{n/2}|\mathbf Q|^{1/2}} \left(1+\frac{\mathbf r_M^T\mathbf Q^{-1}\mathbf r_M}{2b_0}\right)^{-(a_0+n/2)}. \end{aligned} \] Hence \[ \boxed{\mathbf d\mid M\sim t_{2a_0}\left( \mathbf A\boldsymbol\mu_m^M,\frac{b_0}{a_0}\mathbf Q\right)}. \] The degrees of freedom are $2a_0$. The displayed matrix is the Student scale matrix, not its covariance. For $a_0>1$, \[ \operatorname{Var}(\mathbf d\mid M)=\frac{b_0}{a_0-1}\mathbf Q. \] The Student distribution is heavier-tailed than a Gaussian with the same centre and covariance because it averages over uncertainty in $\sigma_m^2$. Equivalently, consider the relation : \[p(d|M)=\frac{p(d|s,M)p(s|M)}{p(s|d,M)}\] We then get: \[ p(\mathbf d\mid M)\propto b_M^{-(a_0+n/2)}, \] since $\mathbf Q$ is common to all candidates. Proof: For a given mean model \(M\), define \[ c_M=\frac12\mathbf r_M^T\mathbf Q^{-1}\mathbf r_M. \] Then \[ b_M=b_0+c_M, \qquad a_M=a_0+\frac n2. \] The likelihood, including its normalizing constant, is \[ p(\mathbf d\mid s,M) = (2\pi)^{-n/2} |\mathbf Q|^{-1/2} s^{-n/2} \exp\left(-\frac{c_M}{s}\right). \] The prior for \(s\) is \[ p(s\mid M)=p(s) = \frac{b_0^{a_0}}{\Gamma(a_0)} s^{-(a_0+1)} \exp\left(-\frac{b_0}{s}\right), \] The posterior distribution is \[ s\mid\mathbf d,M\sim\operatorname{IG}(a_M,b_M), \] and hence its density is \[ p(s\mid\mathbf d,M) = \frac{b_M^{a_M}}{\Gamma(a_M)} s^{-(a_M+1)} \exp\left(-\frac{b_M}{s}\right). \] Substituting the likelihood, the prior, and posterior densities into the relation above gives \begin{align*} p(\mathbf d\mid M) &= \frac{ (2\pi)^{-n/2} |\mathbf Q|^{-1/2} s^{-n/2} \exp(-c_M/s) \, \dfrac{b_0^{a_0}}{\Gamma(a_0)} s^{-(a_0+1)} \exp(-b_0/s) }{ \dfrac{b_M^{a_M}}{\Gamma(a_M)} s^{-(a_M+1)} \exp(-b_M/s) } \\[4pt] &= (2\pi)^{-n/2} |\mathbf Q|^{-1/2} \frac{b_0^{a_0}\Gamma(a_M)} {\Gamma(a_0)b_M^{a_M}} \, s^{-n/2-(a_0+1)+(a_M+1)} \times \exp\left( -\frac{c_M+b_0-b_M}{s} \right). \end{align*} Since \[ a_M=a_0+\frac n2 \text{ and } b_M=b_0+c_M, \] the powers of \(s\) cancel, and the exponential term also cancels: \[ -\frac n2-(a_0+1)+(a_M+1)=0,\quad \text{and} \quad c_M+b_0-b_M=0. \] Therefore, \[ p(\mathbf d\mid M) = (2\pi)^{-n/2} |\mathbf Q|^{-1/2} \frac{b_0^{a_0}\Gamma(a_0+n/2)} {\Gamma(a_0)} b_M^{-(a_0+n/2)}. \] Because \(a_0\), \(b_0\), \(n\), and \(\mathbf Q\) are common to all mean models under consideration, all factors apart from \(b_M\) are constant with respect to \(M\). Consequently, \[ \boxed{ p(\mathbf d\mid M) \propto b_M^{-(a_0+n/2)} =b_M^{-(a_M)}} \] ```{r Model code uncertainty in scale} par(mfrow=c(1,2)); plot_zone_heatmap(z01,"Fixed scale") plot_zone_heatmap(z01_scale,"Scale integrated out"); par(mfrow=c(1,1)) rbind(fixed=summarize_zone(z01)$MAP, integrated=summarize_zone(z01_scale)$MAP) z01_scale[which.min(z01_scale$b_model),c("T","B","V","b_model")] z01_scale[which.max(z01_scale$b_model),c("T","B","V","b_model")] range(z01_scale$b_model) ``` # Part 4: SVD reduction (STK9021 only) ## Part 4(a): Singular-value decomposition Let \[ A=U\Lambda V^T \] be the singular value decomposition. For square $A$, $U$ and $V$ are orthogonal, $U^TU=V^TV=I$, and $\Lambda$ is diagonal with nonnegative singular values $\lambda_1\ge\cdots\ge\lambda_n\ge0$. The columns $u_i$ and $v_i$ are orthonormal bases for the data and parameter spaces, respectively, and satisfy \[ Av_i=\lambda_i u_i, \qquad A^Tu_i=\lambda_i v_i. \] The corresponding expansions are \[ \mathbf m=V\widetilde{\mathbf m}=\sum_i\widetilde m_i v_i, \quad \widetilde m_i=v_i^T\mathbf m, \] \[ \mathbf d=U\widetilde{\mathbf d}=\sum_i\widetilde d_i u_i, \quad \widetilde d_i=u_i^T\mathbf d. \] ```{r svd} SA <- svd(A,nu=n,nv=n) U <- SA$u; lambda <- SA$d; V <- SA$v plot(lambda,type="b",log="y",xlab="Index",ylab="Singular value (log scale), and noise level", main="Singular values of A, and noise level (standard deviation)") lines(c(1,100),sqrt(c(0.04,0.04)),lty=2) par(mfrow=c(5,2),mar=c(2,3,2,1)) for (i in 1:5) { plot(V[,i],type="l",main=paste0("u",i),xlab="",ylab="", ylim=c(-0.3,0.3)) plot(U[,i],type="l",main=paste0("v",i),xlab="",ylab="", ylim=c(-0.3,0.3)) } #Also vectors coresponding to smallest eigen values par(mfrow=c(5,2),mar=c(2,3,2,1)) for (i in 97:101) { plot(V[,i],type="l",main=paste0("u",i),xlab="",ylab="", ylim=c(-0.3,0.3)) plot(U[,i],type="l",main=paste0("v",i),xlab="",ylab="", ylim=c(-0.3,0.3)) } par(mfrow=c(1,1)) # Numerical checks of the SVD identities. c(reconstruction=max(abs(A-U%*%diag(lambda)%*%t(V))), U_orthogonality=max(abs(crossprod(U)-diag(n))), V_orthogonality=max(abs(crossprod(V)-diag(n)))) ``` ## Part 4(b): Transformation to independent equations Define \[ \widetilde{\mathbf d}=U^T\mathbf d, \quad \widetilde{\boldsymbol\varepsilon}=U^T\boldsymbol\varepsilon, \quad \widetilde{\mathbf m}=V^T\mathbf m. \] Premultiplication of $\mathbf d=A\mathbf m+\boldsymbol\varepsilon$ by $U^T$ gives \[ \widetilde{\mathbf d}=U^T U\Lambda V^T\mathbf m+U^T\boldsymbol\varepsilon =\Lambda\widetilde{\mathbf m}+\widetilde{\boldsymbol\varepsilon}. \] Thus, componentwise, \[ \widetilde d_i=\lambda_i\widetilde m_i+\widetilde\varepsilon_i. \] If $\boldsymbol\varepsilon\sim N(0,\sigma_e^2I)$, then \[ \widetilde{\boldsymbol\varepsilon}=U^T\boldsymbol\varepsilon \sim N(0,U^T\sigma_e^2IU)=N(0,\sigma_e^2I). \] The transformed errors are therefore independent Gaussian variables with the same variance. ## Part 4(c): Interpretation The original inverse problem is equivalent to $n$ scalar observation equations. Large singular values identify well-observed parameter-space directions; small singular values identify directions that are weakly informed by the data and hence remain prior-sensitive. ```{r svd_data_check} dtilde <- crossprod(U,d00) mtilde_post <- crossprod(V,post00$mu_post_m) head(data.frame(i=1:n,d_tilde=dtilde, lambda=lambda,m_tilde=mtilde_post),10) ``` # Reproducibility information ```{r session-info} sessionInfo() ```