### ### AUTOR: Arnošt Komárek ### ### ÚČEL: Vážení lehkých objektů, část 1 ### - regresní model, příklad dle knihy ### Box, G. E., Tiao, G. C. (1973). Bayesian Inference in Statistical Analysis. ### John Wiley and Sons. ### ### LOG: 20101106 vytvořeno ### 20101119 přidána inference založená na přímé simulaci z aposteriorního rozdělení pomocí podminování ### 20101212 přidáno MCMC (Gibbsův algoritmus) ### 20141022 mírná revize ### 20181122 přidána možnost vlastního gamma apriorního rozdělení pro tau ### (pouze v části před MCMC) ### ### ============================================================================================================ ROOT <- "/home/komarek/teach/mff_2022/nmst431_Bayes/" library("colorspace") library("mvtnorm") ### Parametry gamma rozdělení pro tau g <- 0; h <- 0 #g <- 1; h <- 0.005 ### Data ### ------------------------------------------------------------------------------ y <- c(109,85, 114,121,140,122,125,129,98,134,133, 217,203,243,229,233,221,221) ## y[1:2] hmotnost objektu A ## y[3:11] hmotnost objektu B ## y[12:18] hmotnost obou objektů A i B xA <- c(rep(1, 2), rep(0, 9), rep(1, 7)) xB <- c(rep(0, 2), rep(1, 9), rep(1, 7)) Data <- data.frame(y = y, xA = xA, xB = xB) print(Data) ### Frekventistický normální lineární regresní model odhadnutý MNČ ### ------------------------------------------------------------------ fitLS <- lm(y ~ -1 + xA + xB, data = Data, x = TRUE) (sfitLS <- summary(fitLS)) b <- coef(fitLS) SSe <- deviance(fitLS) s2 <- sfitLS$sigma^2 ## = s^2 s <- sqrt(s2) varb <- vcov(fitLS) ## = s^2(X'X)^{-1} XtXinv <- (1/s2) * varb dfres <- fitLS$df.residual UKAZ <- data.frame(matrix(b, nrow = 1), SSe, dfres, s2, s) colnames(UKAZ) <- c(names(b), "SSe", "dfres", "s2", "s") print(UKAZ) print(varb) ### frekventistické intervaly spolehlivosti pro regresní koeficienty (spolehlivost 95 % není sdružená) (IS.MNC <- confint(fitLS, conf.level = 0.95)) ### rozumné rozsahy pro kreslení aposteriorních hustot #grid <- confint(fitLS, level=0.9999) ### Princip neurčitosti, případně tau ~ gamma(g, h) ### ------------------------------------------------------------------- ### parametry marginálního aposteriorního rozdělení regresních koeficientů loc1 <- b ## poloha MVT = aposteriorní střední hodnota cmat1 <- (SSe + 2*h)/(dfres + 2*g) * XtXinv ## scale matice MVT nu1 <- dfres + 2*g ## stupně volnosti MVT ### parametry marginálního rozdělení inverzního rozptylu (přesnosti) tau (c1 <- (dfres + 2*g)/2) (d1 <- (SSe + 2*h)/2) ### odhady pro tau (Etau <- c1/d1) ## aposteriorní střední hodnota (medtau <- qgamma(0.5, shape = c1, rate = d1)) ## aposteriorní medián (ETtau <- qgamma(c(0.025, 0.975), shape = c1, rate = d1)) ## 95% ET věrohodnostní interval ### numerické nalezení HPD intervalu pro tau ### (řešení soustavy dvou rovnic o dvou neznámých) fun.to.min.tau <- function(x) { tomin1 <- abs(dgamma(x[1], shape = c1, rate = d1) - dgamma(x[2], shape = c1, rate = d1))/100 tomin2 <- abs(pgamma(x[2], shape = c1, rate = d1) - pgamma(x[1], shape = c1, rate = d1) - 0.95) return(tomin1 + tomin2) } opt.tau <- optim(par = ETtau, fn = fun.to.min.tau) (HPDtau <- opt.tau$par) ### kontrola dgamma(HPDtau, shape = c1, rate = d1) pgamma(HPDtau[2], shape = c1, rate = d1) - pgamma(HPDtau[1], shape = c1, rate = d1) ### odhady pro sigma (medsigma <- 1/sqrt(medtau)) (ETsigma <- 1/sqrt(ETtau[2:1])) ### numerické nalezení HPD intervalu pro sigma fun.to.min.sigma <- function(x) { tomin1 <- abs((2/x[1]^3)*dgamma(1/x[1]^2, shape = c1, rate = d1) - (2/x[2]^3)*dgamma(1/x[2]^2, shape = c1, rate = d1)) tomin2 <- abs(pgamma(1/x[1]^2, shape = c1, rate = d1) - pgamma(1/x[2]^2, shape = c1, rate = d1) - 0.95) return(tomin1 + tomin2) } opt.sigma <- optim(par = ETsigma, fn = fun.to.min.sigma) (HPDsigma <- opt.sigma$par) ### kontrola (2/HPDsigma[1]^3) * dgamma(1/HPDsigma[1]^2, shape = c1, rate = d1) (2/HPDsigma[2]^3) * dgamma(1/HPDsigma[2]^2, shape = c1, rate = d1) pgamma(1/HPDsigma[1]^2, shape = c1, rate = d1) - pgamma(1/HPDsigma[2]^2, shape = c1, rate = d1) ### aposteriorní hustoty (marginální) pro bety lgrid <- 500 x <- fx <- fx2 <- list() for (i in 1:length(b)){ #x[[i]] <- seq(grid[i, 1], grid[i, 2], length=lgrid) x[[i]] <- seq(75, 145, length = lgrid) fx[[i]] <- dmvt(matrix(x[[i]], ncol = 1), delta = loc1[i], sigma = matrix(cmat1[i, i], nrow = 1), df = nu1, log = FALSE) x2 <- (x[[i]] - loc1[i])/sqrt(cmat1[i, i]) fx2[[i]] <- dt(x2, df=nu1) / sqrt(cmat1[i, i]) } ### sdružená aposteriorní hustota pro bety (dvourozměrné t-rozdělení) xbeta <- list(beta1 = seq(80, 115, length = 100), beta2 = seq(110, 135, length = 100)) xx <- cbind(rep(xbeta[[1]], length(xbeta[[2]])), rep(xbeta[[2]], each = length(xbeta[[1]]))) fbety <- matrix(dmvt(xx, delta = loc1, sigma = cmat1, df = nu1, log = FALSE), ncol = length(xbeta[[2]])) rownames(fbety) <- xbeta[[1]] colnames(fbety) <- xbeta[[2]] ### aposteriorní hustota (marginální) pro tau a sigma^2 xtau <- seq(0.001, 0.020, length=lgrid) xsig2 <- 1/xtau[lgrid:1] xsig <- sqrt(xsig2) fxtau <- dgamma(xtau, shape=c1, rate=d1) fxsig2 <- (1/xsig2^2) * fxtau[lgrid:1] ## viz věta o transformaci fxsig <- (2/xsig^3) * fxtau[lgrid:1] ## viz věta o transformaci ### obrázek: jednorozměrné marginální aposteriorní hustoty jednotlivých parametrů layout(matrix(c(1,1,2, 1,1,3), nrow=2, byrow=TRUE)) par(bty="n") plot(x[[1]], fx[[1]], ylim=c(0, max(c(fx[[1]], fx[[2]]))), type="l", col="green", xlab=expression(beta), ylab="Aposteriorní hustota", main="Hmotnost A a B") lines(x[[1]], fx2[[1]], col="darkgreen") ## zcela identicka cara lines(x[[2]], fx[[2]], col="orange") lines(x[[2]], fx2[[2]], col="red") ## zcela identicka cara plot(xtau, fxtau, type="l", col="blue", xlab=expression(tau), ylab="Aposteriorní hustota", main="Přesnost") plot(xsig, fxsig, type="l", col="darkblue", xlab=expression(sigma), ylab="Aposteriorní hustota", main="Směrodatná odchylka") ### obrázek: marginální aposteriorní hustoty tau a sigma #PDF("Vaha_lehkych_13", width = 10, height = 7) par(mfrow=c(1, 2), bty="n", mar=c(4, 4, 4, 0)+0.1) plot(xtau, fxtau, type="l", col="blue", xlab=expression(tau), ylab="Aposteriorní hustota", main="Reziduální přesnost", lwd=2) plot(xsig, fxsig, type="l", col="darkblue", xlab=expression(sigma), ylab="Aposteriorní hustota", main="Reziduální směrodatná odchylka", lwd=2) #dev.off() ### obrázek: dvourozměrná sdružená aposteriorní hustota (beta1, beta2) #PDF("Vaha_lehkych_14", width = (35/25)*7, height = 7) COLP <- rev(heat_hcl(25, c=c(80, 30), l=c(30, 90), power=c(1/5, 2))) par(mfrow=c(1, 1), mar=c(4, 4, 1, 1)+0.1, bty="n") image(xbeta[[1]], xbeta[[2]], fbety, xlab=expression(beta[1]), ylab=expression(beta[2]), col=COLP) contour(xbeta[[1]], xbeta[[2]], fbety, col="brown", add=TRUE) #dev.off() ### obrázek: marginální aposteriorní hustoty beta1 a beta2 #PDF("Vaha_lehkych_15", width = 10, height = 7) par(mfrow=c(1, 1), bty="n", mar=c(4, 4, 1, 1)+0.1) plot(x[[1]], fx[[1]], ylim=c(0, max(c(fx[[1]], fx[[2]]))), type="l", col="darkgreen", xlab=expression(beta), ylab="Aposteriorní hustota", main="", lwd=3) lines(x[[2]], fx[[2]], col="red3", lwd=3) legend(77, 0.1, legend=c("A", "B"), col=c("darkgreen", "red3"), lty=1, lwd=3) #dev.off() ### HPD i ET věrohodnostní intervaly pro regresní koeficienty ### aposteriorní pravděpodobnost je 95 % marginálně, pro každý regresní koeficient zvlášt HPD.presny <- matrix(NA, ncol=2, nrow=length(b)) rownames(HPD.presny) <- c("beta1", "beta2") colnames(HPD.presny) <- c("dolni", "horni") qq <- qt(c(0.025, 0.975), df=nu1) for (i in 1:length(b)){ HPD.presny[i,] <- qq * sqrt(cmat1[i, i]) + loc1[i] } print(HPD.presny) print(IS.MNC) ### sdružená věrohodnostní množina (elipsa) ivarb <- chol2inv(chol(varb)) U <- chol(ivarb) ## U'*U = s^{-2}*X'X Qp <- function(beta1, beta2){ u <- as.numeric(U %*% (c(beta1, beta2) - b)) Qp <- as.numeric(crossprod(u))/length(b) return(Qp) } lgrids <- 100 beta1 <- seq(80, 115, length=lgrids) beta2 <- seq(110, 137, length=lgrids) Qp12 <- matrix(NA, nrow=lgrids, ncol=lgrids) for (i in 1:lgrids) for (j in 1:lgrids) Qp12[i, j] <- Qp(beta1[i], beta2[j]) #PDF("Vaha_lehkych_16", width = (35/25)*7, height = 7) par(mfrow=c(1, 1), bty="n", mar=c(4, 4, 1, 1)+0.1) clevel <- c(0.5, 0.75, 0.9, 0.95) contour(beta1, beta2, Qp12, levels=qf(clevel, length(b), dfres), col="blue", labels=paste(clevel*100, " %", sep=""), xlab="Hmotnost A", ylab="Hmotnost B") contour(beta1, beta2, Qp12, levels=qf(0.95, length(b), dfres), col="red", labels="95 %", add=TRUE, lwd=2) abline(v=HPD.presny[1,], col="brown", lty=4, lwd=2) abline(h=HPD.presny[2,], col="brown", lty=4, lwd=2) points(b[1], b[2], col="orange", cex=2, pch=16) points(b[1], b[2], col="brown", cex=2, pch=13) #dev.off() ### Generování z aposteriorního rozdělení při apriorním rozdělení dle principu neurčitosti ### (přímá simulace založená na podminování) ### ----------------------------------------------------------------------------------------- ### QR rozklad matice X a odmocninová matice pro X'X (bude se hodit pro generování z vícerozměrného normálního rozdělení) ### X = QR, přičemž QQ' = H = X(X'X)^(-1)X' (Q má ve sloupcích bázi M(X), H je projekční matice do M(X)) ### ==> pro L = X'Q platí LL' = X'X, tj. L je odmocninová matice pro symetrickou pozitivně definitní matici X'X ### a dále, označíme-li jako iL = L^(-1), máme iL'iL = (X'X)^(-1), tj. iL' je odmocninová matice pro (X'X)^(-1) X <- fitLS$x Q <- qr.Q(fitLS$qr) L <- t(X) %*% Q L[abs(L) < 1e-15] <- 0 ### ochrana proti zaokrouhlovacím chybám iL <- solve(L) ### samotné generování ### Poznámka: v R není potřeba for cyklus! M <- 1000 #M <- 1e6 set.seed(221913282) tau <- rgamma(M, shape = (dfres + 2*g)/2, rate = (SSe + 2*h)/2) ## náhodný výběr z p(tau | y) sigma <- 1 / sqrt(tau) ## náhodný výběr z p(sigma | y) z <- matrix(rnorm(2*M, mean=0, sd=1), ncol=2) ## řádky obsahují náh. výběr z N(0, I_2) z <- z %*% iL ## nyní obsahují řádky náh. výběr z N(0, (X'X)^(-1)) ## kontrola apply(z, 2, mean) ## pro velké M musí vyjít blízké (0, 0) var(z) ## pro velké M musí vyjít blízké matici (X'X)^(-1) solve(t(X) %*% X) z <- z * cbind(sigma, sigma) ## nyní obsahují řádky náh. výběr z N(0, sigma^2*(X'X)^(-1)), přičemž sigma se liší řádek od řádku beta <- z + matrix(rep(b, M), nrow=M, byrow=TRUE) ## nyní obsahují řádky náh. výběr z N(b, sigma^2*(X'X)^(-1)), přičemž sigma se liší řádek od řádku colnames(beta) <- c("beta1", "beta2") ### scatterploty vygenerovaných hodnot ### (nekresli to pro velké M) PCH <- 21 BGC <- "skyblue" COL <- "blue4" if (M < 1e6){ #PDF(paste("Vaha_lehkych_02-", M, sep = ""), width = 10, height = 7) layout(matrix(c(1,1, 2,3), ncol = 2)) par(bty="n", mar = c(4, 4, 4, 0)+0.1) image(xbeta[[1]], xbeta[[2]], fbety, xlab = expression(beta[1]), ylab = expression(beta[2]), col = COLP) points(beta[,1], beta[,2], col = "grey30", pch=16, xlab = expression(beta[1]), ylab = expression(beta[2])) contour(xbeta[[1]], xbeta[[2]], fbety, col="brown", add=TRUE) plot(beta[,1], tau, pch = PCH, col = COL, bg = BGC, xlab = expression(beta[1]), ylab = expression(tau)) plot(beta[,2], tau, pch = PCH, col = COL, bg = BGC, xlab = expression(beta[2]), ylab = expression(tau)) #dev.off() #PDF(paste("Vaha_lehkych_03-", M, sep = ""), width = 10, height = 7) layout(matrix(c(1,1, 2,3), ncol=2)) par(bty = "n", mar = c(4, 4, 4, 0)+0.1) image(xbeta[[1]], xbeta[[2]], fbety, xlab = expression(beta[1]), ylab = expression(beta[2]), col=COLP) points(beta[,1], beta[,2], col="grey30", pch=16, xlab=expression(beta[1]), ylab=expression(beta[2])) contour(xbeta[[1]], xbeta[[2]], fbety, col="brown", add=TRUE) plot(beta[,1], sigma, pch = PCH, col = COL, bg = BGC, xlab = expression(beta[1]), ylab = expression(sigma)) plot(beta[,2], sigma, pch = PCH, col = COL, bg = BGC, xlab = expression(beta[2]), ylab = expression(sigma)) #dev.off() } ### srovnání histogramů vygenerovaných hodnot se skutečnými aposteriorními hustotami ylim.beta <- c(0, 0.11) #PDF(paste("Vaha_lehkych_01-", M, sep = ""), width = 10, height = 7) par(mfrow=c(2, 2), bty="n", mar=c(4, 4, 4, 0)+0.1) hist(beta[, 1], prob=TRUE, xlim=c(75, 125), ylim=ylim.beta, xlab=expression(beta[1]), ylab="Aposteriorní hustota", main="Hmotnost A", col=rainbow_hcl(1, start=80, end=80)) lines(x[[1]], fx2[[1]], col="darkgreen", lwd=3) hist(beta[, 2], prob=TRUE, xlim=c(100, 150), ylim=ylim.beta, xlab=expression(beta[2]), ylab="Aposteriorní hustota", main="Hmotnost B", col=rainbow_hcl(1, start=20, end=20)) lines(x[[2]], fx2[[2]], col="red4", lwd=3) hist(tau, prob=TRUE, xlim=c(0, 0.015), xlab=expression(tau), ylab="Aposteriorní hustota", main="Inverzní rozptyl", col=rainbow_hcl(1, start=230, end=230)) lines(xtau, fxtau, col="darkblue", lwd=3) hist(sigma, prob=TRUE, xlim=c(5, 32), ylim=range(fxsig), xlab=expression(sigma), ylab="Aposteriorní hustota", main="Směrodatná odchylka", col=rainbow_hcl(1, start=270, end=270)) lines(xsig, fxsig, col="blue", lwd=3) #dev.off() ### Odhady aposteriorních středních hodnot (Ebeta.approx <- apply(beta, 2, mean)) ### Odhady MC chyb (beta.MCchyba <- sqrt((1/M)*apply(beta, 2, var))) ### Odhady aposteriorních mediánů (Medbeta.approx <- apply(beta, 2, median)) ### MC odhad ET věrohodnostního intervalu pro beta1 a beta2 (dolni <- apply(beta, 2, quantile, prob=0.025)) (horni <- apply(beta, 2, quantile, prob=0.975)) ET.approx <- data.frame(dolni=dolni, horni=horni) rownames(ET.approx) <- c("beta1", "beta2") print(ET.approx) ### další s využitím coda balíčku vyber.beta <- mcmc(data.frame(beta1=beta[,1], beta2=beta[,2]))#, sigma=sigma, tau=tau)) summary(vyber.beta) HPDinterval(vyber.beta) vyber.tau <- mcmc(tau) summary(vyber.tau) HPDinterval(vyber.tau) vyber.sigma <- mcmc(sigma) summary(vyber.sigma) HPDinterval(vyber.sigma) ### jádrové odhady marginálních aposteriorních hustot, vše s defaultními volbami vyhlazovacích parametrů dbeta1 <- density(beta[,1]) dbeta2 <- density(beta[,2]) dtau <- density(tau) dsigma <- density(sigma) ### srovnání jádrových odhadů se skutečnými aposteriorními hustotami LWD <- 2 LTY2 <- 2 COL2 <- "grey30" ylim.beta <- c(0, 0.11) #PDF(paste("Vaha_lehkych_04-", M, sep = ""), width = 10, height = 7) par(mfrow=c(2, 2), bty="n", mar=c(4, 4, 4, 0)+0.1) plot(x[[1]], fx2[[1]], type="l", xlim=c(75, 125), ylim=ylim.beta, xlab=expression(beta[1]), ylab="Aposteriorní hustota", main="Hmotnost A", col=COL2, lwd=LWD, lty=LTY2) lines(dbeta1$x, dbeta1$y, col="darkgreen", lwd=LWD) plot(x[[2]], fx2[[2]], type="l", xlim=c(100, 150), ylim=ylim.beta, xlab=expression(beta[2]), ylab="Aposteriorní hustota", main="Hmotnost B", col=COL2, lwd=LWD, lty=LTY2) lines(dbeta2$x, dbeta2$y, col="red4", lwd=LWD) plot(xtau, fxtau, type="l", xlim=c(0, 0.015), ylim=c(0, 200), xlab=expression(tau), ylab="Aposteriorni hustota", main="Inverzní rozptyl", col=COL2, lwd=LWD, lty=LTY2) lines(dtau$x, dtau$y, col="darkblue", lwd=LWD) plot(xsig, fxsig, type="l", xlim=c(5, 32), ylim=range(fxsig), xlab=expression(sigma), ylab="Aposteriorní hustota", main="Inverzní rozptyl", col=COL2, lwd=LWD, lty=LTY2) lines(dsigma$x, dsigma$y, col="blue", lwd=LWD) #dev.off()