### Prognosebaender für Regression ### Streuung der Daten mit Normalverteilung N(., sigma) ### Aufruf: prognoseRegression(x,y,xlow,xhigh,ylow,yhigh, ### formel,parameter,kette, NBand) ### Autor: D. Baettig ## formel ist Regression mit x als unabhaengiger Variable: ## z.B. expression(a+b*x+c*x^2) ## parameter sind Parameter Regression: z.B. c("a","b","c") ## kette ist MCMC-Kette aus STAN-Simulation ## NBand: Anzahl Bänder ## Fenster: [xlow,xhigh] X [ylow,yhigh] ## Daten x (horizontal) und y (vertikal) prognoseRegression <- function(x,y,xlow,xhigh,ylow,yhigh, formel,parameter,kette,NBand){ plot(x,y,xlim=c(xlow,xhigh),ylim=c(ylow,yhigh)) grid() sigma <- extract(kette)$sigma N <- length(sigma) xNeu <- seq(xlow,xhigh,length.out=NBand) yNeu <- matrix(nrow=NBand, ncol =N) p <- extract(kette) m <- length(parameter) listeA <- list(p[[parameter[1]]]) for ( i in 2:m) { listeP <- list(p[[parameter[i]]]) listeA <- c(listeA,listeP) } for (i in 1:NBand) { liste <- list(x=xNeu[i]) liste <- c(liste,listeA) names(liste) <- c("x",parameter) yNeu[i,] <- rnorm(N, eval(formel,liste), sigma) } q975 <- apply(yNeu, 1, quantile, 0.975) q025 <- apply(yNeu, 1, quantile, 0.025) segments(xNeu, q025, xNeu, q975, lwd=7, col=gray(0.85)) q75 <- apply(yNeu, 1, quantile, 0.75) q25 <- apply(yNeu, 1, quantile, 0.25) segments(xNeu, q25, xNeu, q75, lwd=7, col=gray(0.65)) lines(xNeu, (q75+q25)/2, lwd=3) points(x, y, pch=16) }