RResourcesEtc/kre_adaptive_AMSE_minimization.R
In GerritEichner/kader: Kernel Adaptive Density Estimation and Regression
# ### ~/Beratung/StuteW/KernelAdjustedNonparRegression/AMSE_Minimierung.R
# ### Schaetzung und Minimierung der AMSE-Funktion fuer die Simulationen in "Ker-
# ### nel Adjusted Nonparametric Regression" von Eichner & Stute
# ### R 2.12.2, 17./18./22.3./5./6./7./8./14./18./26./27./29.4./27.5./3.6.2011
# ###********************************************************************************
#
# ## Hilfsfunktionen
# ##*****************
# # R-Paket fuer "Kernel Regression Smoothing with Local Plug-in Bandwidth"
# #*************************************************************************
# library( lokern)
#
# # Grafikfunktion
# #****************
# TitleAndTopAxis <- function( at, digits = 4, ...) {
# op <- par( ...); at0 = par( "usr")[ 1]; lh <- 1.1 * par( "cex.axis")
# B <- signif( Biases, digits); V <- signif( Variances, digits)
# M <- signif( MSEs, digits)
# axis( 3, at = c( at0, at), tick = FALSE, hadj = 0, line = 0,
# labels = c( expression( paste( widehat(plain(MSE)), ":", sep = "")),
# bquote( .(M[ 1])), bquote( .(M[ 2])), bquote( .(M[ 3])) ) )
# axis( 3, at = c( at0, at), tick = FALSE, hadj = 0, line = 1 * lh,
# labels = c( expression( paste( widehat(plain(Var)), ":", sep = "")),
# bquote( .(V[ 1])), bquote( .(V[ 2])), bquote( .(V[ 3])) ) )
# axis( 3, at = c( at0, at), tick = FALSE, hadj = 0, line = 2 * lh,
# labels = c( expression( paste( widehat(plain(Bias)), ":", sep = "")),
# bquote( .(B[ 1])), bquote( .(B[ 2])), bquote( .(B[ 3])) ) )
# axis( 3, at = at, tick = FALSE, hadj = 0, line = 2.9 * lh,
# labels = c( expression( hat(m)[n](x[0])), expression( m[n](x[0])),
# expression( plain(lokerns)(x[0])) ) )
#
# title( bquote( list( h[adj] == {n^{-1/5} == .(signif( h, digits))},
# ~~h[ideal] == .(signif( h_opt, digits)),
# ~~bar(h[lokerns]) == .(signif( lox0.meanh, digits)),
# ~~~x[0] == .(x0) )),
# line = 4.2 * lh, outer = TRUE)
# title( bquote( list( "AMSE-optimal, K = Gauss kernel", ~~S == .(MCsamplesize),
# ~~n == .(n), ~~sigma(paste(Y,"|",X==x)) == .(sd.eps) )),
# line = 5.4 * lh, outer = TRUE)
#
# invisible( par( op))
# }
#
#
# # Potenzielle Kern- und damit assoziierte Funktionen
# #****************************************************
# Kern <- function( x) dnorm( x) # Gausskern
#
# Integral_Kquadrat <- integrate( function( u) { k <- Kern( u); k*k },
# lower = -Inf, upper = Inf)$value
# # = 1 / (2 * sqrt( pi)), falls Kern Gausskern
# Integral_xquadratK <- integrate( function( u) u*u * Kern( u),
# lower = -Inf, upper = Inf)$value
# # = 1, falls Kern Gausskern
#
#
# # Potenzielle Regressionsfunktionen
# #***********************************
#
# #m.poly3 <- function( x, a0 = 1, a1 = -3, a2 = 1, a3 = 1)
# # a0 + a1 * x + a2 * x*x + a3 * x*x*x # m'(x) = a1 + 2 a2 x + 3 a3 x^2
# # # m''(x) = 2 a2 + 6 a3 x
# #m.poly4 <- function( x, a0 = 1, a1 = -3, a2 = 1, a3 = 1, a4 = -1)
# # a0 + a1 * x + a2 * x*x + a3 * x*x*x + a4 * x*x*x*x
# # # m'(x) = a1 + 2 a2 x + 3 a3 x^2 + 4 a4 x^3
# # # m''(x) = 2 a2 + 6 a3 x + 12 a4 x^2
# # # Effizienter: a0 + (a1 + (a2 + (a3 + a4 * x) * x) * x) * x
#
# # Gewuenscht ist ein Polynom 4. Grades m mit einem
# # Tief-/Hochpunkt, wo also m'(x) = 0 und m''(x) <> 0, einem
# # Wendepunkt, wo also m''(x) = 0 und m'(x) <> 0 sowie einem
# # Sattelpunkt, wo also m'(x) = 0 = m''(x).
# # Memo: Fuer f(x) = a * (x - x1) * (x - x2)^3 + b ist der
# # TP/HP an (3 * x1 + x2) / 4, der WP an (x1 + x2) / 2 und der SP an x2.
# m.poly4 <- function( x, x1 = 0, x2 = 8, a = 0.01, b = 0)
# a * (x - x1) * (x - x2)^3 + b
#
#
# #m.sin <- function( x, a0 = 0, a1 = 1, a2 = 1) # m'(x) = a1 a2 cos( a2 x)
# # a0 + a1 * sin( a2 * x) # m''(x) = -a1 a2^2 sin( a2 x)
#
#
# # Nadaraya-Watson-Schaetzer ("NW-Schaetzer", m_n im Paper)
# #**********************************************************
# NadWat <- function( x, X, Y, h, K = Kern) {
# # x: Vektor (x_1, ..., x_r) oder Skalar
# # X, Y: Vektoren (X_1, ..., X_n), (Y_1, ..., Y_n)
# # h: Skalar
# # K: Fkt. mit vektorieller Ein- & Ausgabe
# M <- K( outer( x/h, X/h, "-")) # r x n
# drop( (M %*% Y) / rowSums( M)) # r x 1
# }
#
#
# # Gewichte
# #**********
# Wn <- function( sigma, h, xXh, thetaXh, K = Kern) {
# # sigma: Vektor (sigma_1, ..., sigma_s) oder Skalar
# # h: Skalar
# # xXh: Vektor ((x - X_1)/h, ..., (x - X_n)/h)
# # thetaXh: Vektor ((\theta - X_1)/h, ..., (\theta - X_n)/h)
# # K: Fkt. mit vektorieller Ein- & Ausgabe
#
# A <- outer( sigma/h, xXh) # s x n
# A <- outer( A, thetaXh, "+") # s x n x n
# A <- rowSums( K( A), dims = 2) # s x n
# A / rowSums( A) # s x n: (W_{ni}( x, \sigma_r))_{1<=r<=s, 1<=i<=n}
# }
#
#
# # Zwei Versionen des Bias-Schaetzers: BiasA() mit m_n, BiasB() mit \hat m_n
# #****************************************************************************
# BiasA <- function( sigma, h, xXh, thetaXh, K = Kern, mmDiff) {
# # sigma, X, h, xXh, thetaXh, K: Siehe Wn()!
# # mmDiff: Vektor (m_n(X_1) - m_n(x), ..., m_n(X_n) - m_n(x))
#
# drop( Wn( sigma, h, xXh, thetaXh, K) %*% mmDiff) # s x 1
# }
#
#
# BiasB <- function( sigma, Y, h, xXh, XXh, thetaXh, K = Kern) {
# # sigma, h, xXh, thetaXh, K: Siehe Wn()!
# # Y: Vektor (Y_1, ..., Y_n)
# # XXh: Matrix ((X_j - X_i)/h)_{1<=j<=n, 1<=i<=n}
#
# wx <- Wn( sigma, h, xXh, thetaXh, K) # s x n
#
# wX <- apply( XXh, 2, function( XXih) Wn( sigma, h, XXih, thetaXh, K))
# # (1<=r<=s)*(1<=j<=n) x (1<=i<=n)
# wX <- array( wX, dim = c( length( sigma), dim( XXh)))
# # (1<=r<=s) x (1<=j<=n) x (1<=i<=n)
# wX <- aperm( wX - as.vector( wx), c( 2, 1, 3))
# # (1<=j<=n) x (1<=r<=s) x (1<=i<=n)
# rowSums( wx * colSums( wX * Y)) # s x 1
# # (1<=r<=s) x (1<=i<=n)
# }
#
#
# # Zwei Versionen des Varianzschaetzer: VarA() mit m_n, VarB() mit \hat m_n
# #***************************************************************************
# VarA <- function( sigma, h, xXh, thetaXh, K = Kern, YmDiff2) {
# # sigma, h, xXh, thetaXh, K: Siehe Wn()!
# # YmDiff2: Vektor ( (Y_1 - m_n(x))^2, ..., (Y_n - m_n(x))^2 )
#
# wx <- Wn( sigma, h, xXh, thetaXh, K) # s x n
# drop( (wx * wx) %*% YmDiff2) # s x 1
# }
#
#
# VarB <- function( sigma, Y, h, xXh, XXh, thetaXh, K = Kern) {
# # sigma, h, xXh, thetaXh, K: Siehe Wn()!
# # Y: Vektor (Y_1, ..., Y_n)
# # XXh: Matrix ((X_j - X_i)/h)_{1<=j<=n, 1<=i<=n}
#
# wx <- Wn( sigma, h, xXh, thetaXh, K) # s x n
#
# wX <- apply( XXh, 2, function( XXih) Wn( sigma, h, XXih, thetaXh, K))
# # (1<=r<=s)*(1<=j<=n) x (1<=i<=n)
# wX <- aperm( array( wX, dim = c( length( sigma), dim( XXh))), c( 2, 1, 3))
# # (1<=j<=n) x (1<=r<=s) x (1<=i<=n)
# Y_twXY <- Y - t( colSums( wX * Y)) # (1<=i<=n) x (1<=r<=s)
#
# rowSums( (wx * wx) * t( Y_twXY * Y_twXY)) # s x 1
# }
#
#
# # Zwei Versionen des geschaetzten "Asymptotic Mean Squared Errors"
# # (\widehat AMSE): AMSE_A() mit m_n, AMSE_B() mit \hat m_n
# #*****************************************************************
# AMSE_A <- function( sigma, h, xXh, thetaXh, K = Kern, mmDiff, YmDiff2) {
# # sigma, h, xXh, thetaXh, K: Siehe Wn()!
# # mmDiff: Vektor (m_n(X_1) - m_n(x), ..., m_n(X_n) - m_n(x))
# # YmDiff2: Vektor ( (Y_1 - m_n(x))^2, ..., (Y_n - m_n(x))^2 )
# W <- Wn( sigma, h, xXh, thetaXh, K) # s x n
# Bn <- W %*% mmDiff # = BiasA( sigma, h, xXh, thetaXh, K, mmDiff)
# Vn2 <- (W * W) %*% YmDiff2 # = VarA( sigma, h, xXh, thetaXh, K, YmDiff2)
# drop( Vn2 + Bn * Bn) # s x 1
# }
#
#
# AMSE_B <- function( sigma, Y, h, xXh, XXh, thetaXh, K = Kern) {
# # sigma, h, xXh, thetaXh, K: Siehe Wn()!
# # Y: Vektor (Y_1, ..., Y_n)
# # XXh: Matrix ((X_j - X_i)/h)_{1<=j<=n, 1<=i<=n}
#
# wx <- Wn( sigma, h, xXh, thetaXh, K) # s x n
#
# wX <- apply( XXh, 2, function( XXih) Wn( sigma, h, XXih, thetaXh, K))
# # (1<=r<=s)*(1<=j<=n) x (1<=i<=n)
# wX <- array( wX, dim = c( length( sigma), dim( XXh)))
# # (1<=r<=s) x (1<=j<=n) x (1<=i<=n)
#
# wX_wx <- aperm( wX - as.vector( wx), c( 2, 1, 3))
# # (1<=j<=n) x (1<=r<=s) x (1<=i<=n)
# Bn <- rowSums( wx * colSums( wX_wx * Y))
# # = BiasB( sigma, Y, h, xXh, XXh, thetaXh, K), s x 1
#
# Y_twXY <- Y - t( colSums( aperm( wX, c( 2, 1, 3)) * Y)) # (1<=i<=n) x (1<=r<=s)
# Vn2 <- rowSums( (wx * wx) * t( Y_twXY * Y_twXY))
# # = VarB( sigma, Y, h, xXh, XXh, thetaXh, K), s x 1
#
# drop( Vn2 + Bn * Bn) # s x 1
# }
#
#
#
# ## Die Simulations-Szenarien
# ##***************************
#
# ##Version <- "A" # Varianz- und Bias-Schaetzung mit NW-Schaetzer m_n fuer m
# Version <- "B" # Varianz- und Bias-Schaetzung mit \hat m_n fuer m
#
# # Simulationsparameter
# x00 <- c( 2, 4, 5, 8) # Punkte, an denen die AMSE-optim. "kernel adjusted
# # nonparametric regressions" von m stattfinden sollen.
# # Fuer m.poly4()'s Voreinstellung ist TP an 2, WP an 4
# # und SP an 8; ein "beliebiger" Punkt liegt z. B. bei 5.
# nn <- c( 25, 50, 100) # Stichprobenumfaenge.
# sdd <- c( 0.5, 1, 2) # Std.-abwn. der Fehler \epsilon_i: \sqrt(Var(Y|X=x)).
#
# Resultate <- array( NA, dim = c( length( x00), length( nn), length( sdd), 3, 3),
# dimnames = list( x0 = paste( "x0", x00, sep = "="),
# n = paste( "n", nn, sep = "="),
# sd = paste( "sd", sdd, sep = "="),
# Schaetzer = c( "Bias", "Var", "MSE"),
# Methode = c( "adjusted", "idealistic",
# "classic")) )
#
# mm <- m.poly4 # Regressionsfunktion und ihre erste sowie zweite Ableitung
# mm_prime <- function() {}; body( mm_prime) <- D( body( mm), "x")
# mm_second <- function() {}; body( mm_second) <- D( body( mm_prime), "x")
# formals( mm_second) <- formals( mm_prime) <- formals( mm)
#
# # Dichte der Normalverteilung (der X_i) und ihre erste Ableitung
# f <- dnorm; f_prime <- function( x) -x * dnorm( x)
#
#
# for( x0 in x00) { # Punkt x_0
#
# design.mean <- x0 + 0.5 # "Zentrum" und "Skala" des Designs = Erwartungswert
# design.sd <- 1 # bzw. Standardabweichung der Verteilung der X_i.
# xx <- design.mean + seq( -3.3 * design.sd, 3.3 * design.sd, length = 101)
# # x-Gitter, worauf m und der NW-Schaetzer m_n ausgewertet & gezeichnet werden.
#
# sig.range <- c( 0.01, 1.5 * diff( range( xx))) # Zu betrachtender Wertebereich
# sig.grid <- seq( sig.range[ 1], # samt Gitter fuer den Skalenparame-
# sig.range[ 2], length = 501)# ter \sigma des adaptiv. Verfahrens.
#
# for( n in nn) { # Stichprobenumfang
# for( sd.eps in sdd) { # Standardabweichung der \epsilon_i
#
# set.seed( 20110426) # Startwert des Pseudo-Zufallszahlengenerators
#
# # Berechnung der fuer x_0 optimalen, "idealistischen" Fensterbreite h_opt
# # zur Berechnung des NW-Schaetzers an x_0
# numerator <- sd.eps^2 * Integral_Kquadrat * f( x0)
# denominator <- n * Integral_xquadratK^2 *
# (f( x0) * mm_second( x0) + 2 * f_prime( x0) * mm_prime( x0))^2
# h_opt <- (numerator / denominator)^(1/5)
#
# # Wahl von h fuer das adaptive Verfahren
# h <- n^(-1/5)
#
# # Vorbereitung des Grafik-Ausgabe
# now <- format( Sys.time(), "%Y%m%d_%H%M%S")
# pdf( paste( "Plots/mBVA_x", x0, "_n", formatC( n, width = 3, flag = 0), "_s",
# formatC( 10*sd.eps, width = 2, flag = 0), "_D", now, ".pdf",
# sep = ""),
# paper = "a4", width = 0, height = 0)
# op <- par( mfrow = c( 3,1), mar = c( 3,3,3,0.1)+0.1, mgp = c( 1.5,0.5,0),
# tcl = -0.3, cex.main = 2)
#
# ## Ein paar beispielhafte Simulationslaeufe mit Grafikgenerierung
# ##****************************************************************
# for( bsp in 1:4) {
# # Datengenerierung Y_i = m( X_i) + \epsilon_i, i = 1, ..., n
# X <- rnorm( n, mean = design.mean, sd = design.sd)
# Y <- mm( X) + rnorm( length( X), sd = sd.eps)
#
# # NW-Schaetzer auf einem x-Gitter, um ihn zur Anschauung zu zeichnen
# mnxx <- NadWat( xx, X, Y, h_opt)
#
# # Kern-Schaetzer mit lokal-adaptiver "plug-in"-Bandbreite (aus dem R-Paket
# # "lokern"), auf selbem x-Gitter, um auch ihn zur Anschauung zu zeichnen
# loxx <- lokerns( x = X, y = Y, x.out = xx)$est
#
# # Vorratsberechnung von Groessen, die nicht von \sigma abhaengen, aber
# # i. F. wiederholt gebraucht werden; Algor. wird so effizienter.
# x0Xh <- (x0 - X) / h; thetaXh <- (mean( X) - X) / h
#
# if( Version == "A") { # Var.- & Bias-Schaetzung mit m_n
#
# # NW-Schaetzer m_n(x0) und m_n(X_1), ..., m_n(X_n) mit Fensterbreite
# # h = n^(-1/5) sowie daraus gebildete Groessen zur Schaetzung von
# # Var_x0(\sigma) und Bias_x0(\sigma)
# mnx0 <- NadWat( x0, X, Y, h) # m_n(x0)
# mnX <- NadWat( X, X, Y, h) # m_n(X_i) fuer i = 1, ..., n
# mnX_mnx0 <- mnX - mnx0 # m_n(X_i) - m_n(x0) fuer i = 1, ..., n
# Y_mnX2 <- (Y - mnX)^2 # (Y_i - m_n(X_i))^2 fuer i = 1, ..., n
#
# # Schaetzer von Var_x0(\sigma) und Bias_x0(\sigma) auf dem \sigma-Gitter,
# # um sie zur Anschauung zu zeichnen
# Bn <- BiasA( sigma = sig.grid, h = h, xXh = x0Xh, thetaXh = thetaXh,
# mmDiff = mnX_mnx0)
# Vn2 <- VarA( sigma = sig.grid, h = h, xXh = x0Xh, thetaXh = thetaXh,
# YmDiff2 = Y_mnX2)
#
# # # Versuch 1: Minimizer und Minimum des Schaetzers von AMSE_x0(\sigma):
# # opt <- optimize( f = AMSE_A, interval = sig.range, h = h, xXh = x0Xh,
# # thetaXh = thetaXh, mmDiff = mnX_mnx0, YmDiff2 = Y_mnX2)
# # sig_opt_x0 <- opt$minimum; amse_min <- opt$objective
# #
# # # Versuch 2: Minimizer und Minimum des Schaetzers von AMSE_x0(\sigma):
# # opt <- optimize( f = AMSE_A, interval = sig.range/c(1,3), h = h, xXh = x0Xh,
# # thetaXh = thetaXh, mmDiff = mnX_mnx0, YmDiff2 = Y_mnX2)
# # sig_opt_x0_2 <- opt$minimum; amse_min_2 <- opt$objective
#
# } else if( Version == "B") { # Var.- & Bias-Schaetzung mit \hat m_n
#
# # Weitere Vorratsgroesse, die nicht von \sigma abhaengt, aber i. F.
# # wiederholt gebraucht wird
# XjXih <- outer( X/h, X/h, "-")
#
# # Schaetzer von Var_x0(\sigma) und Bias_x0(\sigma) auf dem \sigma-Gitter,
# # um sie zur Anschauung zu zeichnen
# Bn <- BiasB( sigma = sig.grid, Y = Y, h = h, xXh = x0Xh, XXh = XjXih,
# thetaXh = thetaXh)
# Vn2 <- VarB( sigma = sig.grid, Y = Y, h = h, xXh = x0Xh, XXh = XjXih,
# thetaXh = thetaXh)
#
# # # Versuch 1: Minimizer und Minimum des Schaetzers von AMSE_x0(\sigma):
# # opt <- optimize( f = AMSE_B, interval = sig.range, Y = Y, h = h,
# # xXh = x0Xh, XXh = XjXih, thetaXh = thetaXh)
# # sig_opt_x0 <- opt$minimum; amse_min <- opt$objective
# #
# # # Versuch 2: Minimizer und Minimum des Schaetzers von AMSE_x0(\sigma):
# # opt <- optimize( f = AMSE_B, interval = sig.range/c(1,3), Y = Y, h = h,
# # xXh = x0Xh, XXh = XjXih, thetaXh = thetaXh)
# # sig_opt_x0_2 <- opt$minimum; amse_min_2 <- opt$objective
#
# } else stop( "Version '", Version, "' gibt es nicht!")
#
# # # Versuch 2 "besser" als Versuch 1?
# # if( amse_min_2 < amse_min) {
# # sig_opt_x0 <- sig_opt_x0_2; amse_min <- amse_min_2
# # }
#
# # "Komposition" von AMSE_x0(\sigma) aus den Schaetzern von Var_x0(\sigma)
# # und Bias_x0(\sigma) auf dem \sigma-Gitter
# An <- Vn2 + Bn * Bn
#
# # Minimizer und Minimum des Schaetzers von AMSE_x0(\sigma)
# # durch "grid search":
# min_idx <- which.min( An)
# sig_opt_x0 <- sig.grid[ min_idx]; amse_min <- An[ min_idx]
#
# # Der AMSE-optimale "kernel adjusted nonparametric regression"-Schaetzwert
# # an x_0, also \hat m_n( x_0)
# W <- Wn( sigma = sig_opt_x0, h = h, xXh = x0Xh, thetaXh = thetaXh)
# mnhatx0 <- W %*% Y
#
# # Grafische Darstellung fuer den aktuellen Datensatz
# plot( Y ~ X, xlim = range( X, xx),
# ylim = range( Y, mnxx, mnhatx0, loxx, na.rm = TRUE),# mm(xx)
# main = bquote( n == .(n)), xlab = "x", ylab = "y") # Das Streudiagramm
# lines( xx, mm( xx), lty = 2) # Die "wahre" Regressionsfunktion m
# lines( xx, mnxx) # Der NW-Schaetzer m_n von m
# lines( xx, loxx, col = "blue") # Der lokal-adaptive Schaetzer von m
# points( x0, mnhatx0, col = "red", pch = 8, cex = 2) # (x_0, \hat m_n(x_0))
# legend( "bottomleft", lty = c( 2, 1, 1), col = c( "black", "black", "blue"),
# bty = "n", cex = 1.5,
# legend = c( bquote( paste( "m, ", sigma(paste( Y,"|", X==x))
# == .(sd.eps)) ),
# expression( m[n]), "lokern()"))
# legend( "top", pch = 8, col = "red", bty = "n", cex = 1.5,
# legend = bquote( paste( hat(m)[n](x[0]), ", ", x[0] == .(x0))))
#
# matplot( sig.grid, cbind( Bn*Bn, Vn2), type = "l", lty = 1:2,
# col = c( "black", "red"), xlab = expression( sigma), ylab = "")
# legend( "top", lty = 1:2, col = c( "black", "red"), bty = "n", cex = 1.5,
# legend = c( expression( paste( widehat(plain(Bias))[x[0]]^2, (sigma))),
# expression( widehat(plain(Var))[x[0]](sigma))))
#
# plot( sig.grid, An, type = "l", xlab = expression( sigma), ylab = "",
# main = expression( widehat(plain(AMSE))[x[0]](sigma)))
# points( sig_opt_x0, amse_min, pch = 13, col = "red", cex = 2)
# legend( "top", pch = 13, col = "red", bty = "n", cex = 1.5,
# legend = substitute( group( "(", list( plain( Minimizer),
# plain( Minimum)), ")")
# == group( "(", list( x, y), ")"),
# list( x = signif( amse_min, 4),
# y = signif( sig_opt_x0, 4))))
# } # bsp-loop
# par( op)
# dev.off() ## Ende der beispielhaften Simulationslaeufe
#
#
# ## Die eigentliche Simulation (Kommentierung siehe oben!)
# ##****************************
# MCsamplesize <- 5000
#
# now <- format( Sys.time(), "%Y%m%d_%H%M%S")
# set.seed( 20110426) # Derselbe Startwert des Pseudo-RNGs wie oben.
# Sim <- replicate( MCsamplesize, { # Schleife ueber Monte-Carlo-Samples
# X <- rnorm( n, mean = design.mean, sd = design.sd)
# Y <- mm( X) + rnorm( length( X), sd = sd.eps)
#
# if( n > 25 || sd.eps == 1) {
# c( mnhatx0 = NA, mnx0_hopt = NA, lox0 = NA, lox0.h = NA)
# } else {
# mnx0_hopt <- NadWat( x0, X, Y, h_opt)
# lokerns.x0 <- lokerns( x = X, y = Y, x.out = x0)
# lox0 <- lokerns.x0$est; lox0.h <- lokerns.x0$bandwidth
#
# x0Xh <- (x0 - X) / h; thetaXh <- (mean( X) - X) / h
#
# if( Version == "A") {
# mnx0 <- NadWat( x0, X, Y, h); mnX <- NadWat( X, X, Y, h)
# mnX_mnx0 <- mnX - mnx0; Y_mnX2 <- (Y - mnX)^2
# # # Minimierungsversuch 1:
# # opt <- optimize( f = AMSE_A, interval = sig.range, h = h, xXh = x0Xh,
# # thetaXh = thetaXh, mmDiff = mnX_mnx0, YmDiff2 = Y_mnX2)
# # sig_opt_x0 <- opt$minimum; amse_min <- opt$objective
# # # Minimierungsversuch 2:
# # opt <- optimize( f = AMSE_A, interval = sig.range/c(1,3), h = h, xXh = x0Xh,
# # thetaXh = thetaXh, mmDiff = mnX_mnx0, YmDiff2 = Y_mnX2)
# # sig_opt_x0_2 <- opt$minimum; amse_min_2 <- opt$objective
#
# An <- AMSE_A( sig.grid, h = h, xXh = x0Xh, thetaXh = thetaXh,
# mmDiff = mnX_mnx0, YmDiff2 = Y_mnX2)
#
# } else if( Version == "B") {
# XjXih <- outer( X/h, X/h, "-")
# # # Minimierungsversuch 1:
# # opt <- optimize( f = AMSE_B, interval = sig.range, Y = Y, h = h,
# # xXh = x0Xh, XXh = XjXih, thetaXh = thetaXh)
# # sig_opt_x0 <- opt$minimum; amse_min <- opt$objective
# # # Minimierungsversuch 2:
# # opt <- optimize( f = AMSE_B, interval = sig.range/c(1, 3), Y = Y, h = h,
# # xXh = x0Xh, XXh = XjXih, thetaXh = thetaXh)
# # sig_opt_x0_2 <- opt$minimum; amse_min_2 <- opt$objective
#
# An <- AMSE_B( sig.grid, Y = Y, h = h, xXh = x0Xh, XXh = XjXih,
# thetaXh = thetaXh)
#
# } else stop( "Version '", Version, "' gibt es nicht!")
#
# # # Minimierungsversuch 2 "besser" als 1?
# # if( amse_min_2 < amse_min) sig_opt_x0 <- sig_opt_x0_2
#
# # Minimizer und Minimum des Schaetzers von AMSE_x0(\sigma)
# # durch "grid search":
# sig_opt_x0 <- sig.grid[ which.min( An)]
#
# W <- Wn( sigma = sig_opt_x0, h = h, xXh = x0Xh, thetaXh = thetaXh)
# mnhatx0 <- W %*% Y
#
# c( mnhatx0 = mnhatx0, mnx0_hopt = mnx0_hopt, lox0 = lox0, lox0.h = lox0.h)
# }
# } )
#
# lox0.meanh <- mean( Sim[ "lox0.h",]); Sim <- Sim[ -4,]
# mx0 <- mm( x0); Means <- rowMeans( Sim); Biases <- Means - mx0
# Variances <- apply( Sim, 1, var); MSEs <- Biases * Biases + Variances
#
# Resultate[ paste( "x0", x0, sep = "="), paste( "n", n, sep = "="),
# sd = paste( "sd", sd.eps, sep = "="),
# c( "Bias", "Var", "MSE"), ] <- rbind( Biases, Variances, MSEs)
#
# if( !(n > 25 || sd.eps == 1)) {
# ## Grafische Darstellung der Resultate (Boxplots)
# ##************************************************
# Bim1 <- as.list( as.data.frame( t( Sim)))
# Bim2 <- list( "'Idealistic' - Adjusted" = Sim[ "mnx0_hopt",] - Sim[ "mnhatx0",],
# "'Classic' - Adjusted" = Sim[ "lox0",] - Sim[ "mnhatx0",])
#
# pch <- c( 1, 4, 5); cols <- c( "red", "orange", "green"); cex <- 2
# pdf( paste( "Plots/AMSEopt_x", x0, "_n", formatC( n, width = 3, flag = 0),
# "_s", formatC( 10*sd.eps, width = 2, flag = 0), "_D", now,
# ".pdf", sep = ""),
# paper = "a4", width = 0, height = 0)
# par( pl <- list( oma = c( 0,0,8,0), mar = c( 3.5, 3.3, 1, 3)+.1,
# mgp = c( 2, 0.5, 0), tcl = -0.3, lab = c( 5, 10, 0)))
# boxplot( Bim1, names = c( "Adjusted ...", "'Idealistic' ...", "'Classic' ..."),
# ylab = expression( hat(m)[n](x[0])~~~~"bzw."~~~~m[n](x[0])
# ~~~~"bzw."~~~~lokerns(x[0])),
# cex.axis = 1.2, cex.lab = 1.2,
# xlab = expression( "kernel regression estimate values at"~x[0]))
#
# abline( h = mx0, col = "blue")
# axis( 4, at = mx0, tick = FALSE, labels = expression( m(x[0])), line = 0,
# las = 2, cex.axis = 1.2, col.axis = "blue")
#
# points( Means, col = cols, pch = pch, cex = cex, lwd = 3)
# legend( "top", pch = pch, col = cols, pt.lwd = 3, bty = "n", cex = 1.4,
# legend = c( expression( bar(hat(m)[n](x[0]))),
# expression( bar(m[n](x[0]))),
# expression( bar(lokerns(x[0]))) ), y.intersp = 1.7)
# TitleAndTopAxis( at = 1:3 - 0.18, cex.axis = 1.2)
#
#
# boxplot( Bim2, cex.axis = 1.1, cex.lab = 1.1,
# ylab = expression( m[n](x[0]) - hat(m)[n](x[0])~~~~~~"bzw."~~~~~~
# lokerns(x[0]) - hat(m)[n](x[0])),
# xlab = bquote( "Difference between kernel regression estimates at "~
# x[0] == .(x0)~"with"~
# sigma(paste(Y,"|",X==x)) == .(sd.eps)))
# abline( h = 0, col = "blue")
# points( c( Means[ "mnx0_hopt"] - Means[ "mnhatx0"],
# Means[ "lox0"] - Means[ "mnhatx0"]),
# col = c( "darkgreen", "brown"), pch = c( 8, 9), cex = cex, lwd = 2)
# legend( "topleft", pch = c( 8, 9, -1), col = c( "darkgreen", "brown", NA),
# pt.lwd = 2, bty = "n", cex = 1.4,
# legend = c( bquote( bar(m[n](x[0]) - hat(m)[n](x[0])) ==
# .(round( Means[ "mnx0_hopt"] - Means[ "mnhatx0"],
# 5))),
# bquote( bar(lokerns(x[0]) - hat(m)[n](x[0]))==
# .(round( Means[ "lox0"] - Means[ "mnhatx0"], 5))),
# expression() ), y.intersp = 1.7)
# TitleAndTopAxis( at = c( 0.833, 1.5, 2.165) - 0.12, cex.axis = 1.2)
# dev.off()
# }
# } # sd.eps-loop
# } # n-loop
# } # x0-loop
#
# # rm( list = ls()); dev.off()
#
#
# today <- format( Sys.time(), "%Y%m%d")
# sink( paste( "Resultate_", today, ".txt", sep = ""))
# print( ftable( Resultate))
# sink()
#
# #round( ftable( Resultate[ c( "x0=4", "x0=8", "x0=2"),
# # c( "n=25", "n=50", "n=100"),
# # c( "sd=0.5", "sd=2"), , ],
# # row.vars = c( "x0", "sd", "Schaetzer") ), 4)
#
# library( Hmisc)
# punkte <- c( "minimum" = "x0=2", "point of inflection" = "x0=4",
# "'arbitrary' point" = "x0=5", "saddle point" = "x0=8")
# Tab <- lapply( punkte, function( v)
# ftable( Resultate[ v, c( "n=50", "n=25", "n=100"),
# c( "sd=0.5", "sd=2"), , ],
# row.vars = c( "sd", "Schaetzer")) )
# names( Tab) <- paste( "Performance in", names( punkte),
# paste( "$(", gsub( "x0", "x_0", punkte, fixed = TRUE),
# ")$", sep = ""),
# "based on", MCsamplesize, "Monte Carlo samples" )
#
# CV <- attr( Tab[[1]], "col.vars")
# CV$n <- paste( "$", CV$n, "$", sep = "")
# CV$Methode2 <- c( "$\\hat m_n$", "$m_n$", "\\texttt{lokerns}")
#
# RV <- attr( Tab[[1]], "row.vars")
# RV$sd <- paste( "$", gsub( "sd", "\\sigma", RV$sd, fixed = TRUE), "$", sep = "")
# RV$Schaetzer <- paste( "$\\widehat{\\text{", RV$Schaetzer, "}}$", sep = "")
#
# w <- latex( object = Tab, file = paste( "Tab_", today, ".tex", sep = ""),
# title = "", where = "!htbp", dec = 3, first.hline.double = FALSE,
# cgroup = CV$n, colheads = rep( CV$Methode2, length( CV$n)),
# extracolheads = rep( CV$Methode, length( CV$n)),
# rgroup = RV$sd, rowname = rep( RV$Schaetzer, length( RV$sd))
# )
GerritEichner/kader documentation built on May 10, 2019, 1:14 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# ### ~/Beratung/StuteW/KernelAdjustedNonparRegression/AMSE_Minimierung.R
# ### Schaetzung und Minimierung der AMSE-Funktion fuer die Simulationen in "Ker-
# ### nel Adjusted Nonparametric Regression" von Eichner & Stute
# ### R 2.12.2, 17./18./22.3./5./6./7./8./14./18./26./27./29.4./27.5./3.6.2011
# ###********************************************************************************
#
# ## Hilfsfunktionen
# ##*****************
# # R-Paket fuer "Kernel Regression Smoothing with Local Plug-in Bandwidth"
# #*************************************************************************
# library( lokern)
#
# # Grafikfunktion
# #****************
# TitleAndTopAxis <- function( at, digits = 4, ...) {
# op <- par( ...); at0 = par( "usr")[ 1]; lh <- 1.1 * par( "cex.axis")
# B <- signif( Biases, digits); V <- signif( Variances, digits)
# M <- signif( MSEs, digits)
# axis( 3, at = c( at0, at), tick = FALSE, hadj = 0, line = 0,
# labels = c( expression( paste( widehat(plain(MSE)), ":", sep = "")),
# bquote( .(M[ 1])), bquote( .(M[ 2])), bquote( .(M[ 3])) ) )
# axis( 3, at = c( at0, at), tick = FALSE, hadj = 0, line = 1 * lh,
# labels = c( expression( paste( widehat(plain(Var)), ":", sep = "")),
# bquote( .(V[ 1])), bquote( .(V[ 2])), bquote( .(V[ 3])) ) )
# axis( 3, at = c( at0, at), tick = FALSE, hadj = 0, line = 2 * lh,
# labels = c( expression( paste( widehat(plain(Bias)), ":", sep = "")),
# bquote( .(B[ 1])), bquote( .(B[ 2])), bquote( .(B[ 3])) ) )
# axis( 3, at = at, tick = FALSE, hadj = 0, line = 2.9 * lh,
# labels = c( expression( hat(m)[n](x[0])), expression( m[n](x[0])),
# expression( plain(lokerns)(x[0])) ) )
#
# title( bquote( list( h[adj] == {n^{-1/5} == .(signif( h, digits))},
# ~~h[ideal] == .(signif( h_opt, digits)),
# ~~bar(h[lokerns]) == .(signif( lox0.meanh, digits)),
# ~~~x[0] == .(x0) )),
# line = 4.2 * lh, outer = TRUE)
# title( bquote( list( "AMSE-optimal, K = Gauss kernel", ~~S == .(MCsamplesize),
# ~~n == .(n), ~~sigma(paste(Y,"|",X==x)) == .(sd.eps) )),
# line = 5.4 * lh, outer = TRUE)
#
# invisible( par( op))
# }
#
#
# # Potenzielle Kern- und damit assoziierte Funktionen
# #****************************************************
# Kern <- function( x) dnorm( x) # Gausskern
#
# Integral_Kquadrat <- integrate( function( u) { k <- Kern( u); k*k },
# lower = -Inf, upper = Inf)$value
# # = 1 / (2 * sqrt( pi)), falls Kern Gausskern
# Integral_xquadratK <- integrate( function( u) u*u * Kern( u),
# lower = -Inf, upper = Inf)$value
# # = 1, falls Kern Gausskern
#
#
# # Potenzielle Regressionsfunktionen
# #***********************************
#
# #m.poly3 <- function( x, a0 = 1, a1 = -3, a2 = 1, a3 = 1)
# # a0 + a1 * x + a2 * x*x + a3 * x*x*x # m'(x) = a1 + 2 a2 x + 3 a3 x^2
# # # m''(x) = 2 a2 + 6 a3 x
# #m.poly4 <- function( x, a0 = 1, a1 = -3, a2 = 1, a3 = 1, a4 = -1)
# # a0 + a1 * x + a2 * x*x + a3 * x*x*x + a4 * x*x*x*x
# # # m'(x) = a1 + 2 a2 x + 3 a3 x^2 + 4 a4 x^3
# # # m''(x) = 2 a2 + 6 a3 x + 12 a4 x^2
# # # Effizienter: a0 + (a1 + (a2 + (a3 + a4 * x) * x) * x) * x
#
# # Gewuenscht ist ein Polynom 4. Grades m mit einem
# # Tief-/Hochpunkt, wo also m'(x) = 0 und m''(x) <> 0, einem
# # Wendepunkt, wo also m''(x) = 0 und m'(x) <> 0 sowie einem
# # Sattelpunkt, wo also m'(x) = 0 = m''(x).
# # Memo: Fuer f(x) = a * (x - x1) * (x - x2)^3 + b ist der
# # TP/HP an (3 * x1 + x2) / 4, der WP an (x1 + x2) / 2 und der SP an x2.
# m.poly4 <- function( x, x1 = 0, x2 = 8, a = 0.01, b = 0)
# a * (x - x1) * (x - x2)^3 + b
#
#
# #m.sin <- function( x, a0 = 0, a1 = 1, a2 = 1) # m'(x) = a1 a2 cos( a2 x)
# # a0 + a1 * sin( a2 * x) # m''(x) = -a1 a2^2 sin( a2 x)
#
#
# # Nadaraya-Watson-Schaetzer ("NW-Schaetzer", m_n im Paper)
# #**********************************************************
# NadWat <- function( x, X, Y, h, K = Kern) {
# # x: Vektor (x_1, ..., x_r) oder Skalar
# # X, Y: Vektoren (X_1, ..., X_n), (Y_1, ..., Y_n)
# # h: Skalar
# # K: Fkt. mit vektorieller Ein- & Ausgabe
# M <- K( outer( x/h, X/h, "-")) # r x n
# drop( (M %*% Y) / rowSums( M)) # r x 1
# }
#
#
# # Gewichte
# #**********
# Wn <- function( sigma, h, xXh, thetaXh, K = Kern) {
# # sigma: Vektor (sigma_1, ..., sigma_s) oder Skalar
# # h: Skalar
# # xXh: Vektor ((x - X_1)/h, ..., (x - X_n)/h)
# # thetaXh: Vektor ((\theta - X_1)/h, ..., (\theta - X_n)/h)
# # K: Fkt. mit vektorieller Ein- & Ausgabe
#
# A <- outer( sigma/h, xXh) # s x n
# A <- outer( A, thetaXh, "+") # s x n x n
# A <- rowSums( K( A), dims = 2) # s x n
# A / rowSums( A) # s x n: (W_{ni}( x, \sigma_r))_{1<=r<=s, 1<=i<=n}
# }
#
#
# # Zwei Versionen des Bias-Schaetzers: BiasA() mit m_n, BiasB() mit \hat m_n
# #****************************************************************************
# BiasA <- function( sigma, h, xXh, thetaXh, K = Kern, mmDiff) {
# # sigma, X, h, xXh, thetaXh, K: Siehe Wn()!
# # mmDiff: Vektor (m_n(X_1) - m_n(x), ..., m_n(X_n) - m_n(x))
#
# drop( Wn( sigma, h, xXh, thetaXh, K) %*% mmDiff) # s x 1
# }
#
#
# BiasB <- function( sigma, Y, h, xXh, XXh, thetaXh, K = Kern) {
# # sigma, h, xXh, thetaXh, K: Siehe Wn()!
# # Y: Vektor (Y_1, ..., Y_n)
# # XXh: Matrix ((X_j - X_i)/h)_{1<=j<=n, 1<=i<=n}
#
# wx <- Wn( sigma, h, xXh, thetaXh, K) # s x n
#
# wX <- apply( XXh, 2, function( XXih) Wn( sigma, h, XXih, thetaXh, K))
# # (1<=r<=s)*(1<=j<=n) x (1<=i<=n)
# wX <- array( wX, dim = c( length( sigma), dim( XXh)))
# # (1<=r<=s) x (1<=j<=n) x (1<=i<=n)
# wX <- aperm( wX - as.vector( wx), c( 2, 1, 3))
# # (1<=j<=n) x (1<=r<=s) x (1<=i<=n)
# rowSums( wx * colSums( wX * Y)) # s x 1
# # (1<=r<=s) x (1<=i<=n)
# }
#
#
# # Zwei Versionen des Varianzschaetzer: VarA() mit m_n, VarB() mit \hat m_n
# #***************************************************************************
# VarA <- function( sigma, h, xXh, thetaXh, K = Kern, YmDiff2) {
# # sigma, h, xXh, thetaXh, K: Siehe Wn()!
# # YmDiff2: Vektor ( (Y_1 - m_n(x))^2, ..., (Y_n - m_n(x))^2 )
#
# wx <- Wn( sigma, h, xXh, thetaXh, K) # s x n
# drop( (wx * wx) %*% YmDiff2) # s x 1
# }
#
#
# VarB <- function( sigma, Y, h, xXh, XXh, thetaXh, K = Kern) {
# # sigma, h, xXh, thetaXh, K: Siehe Wn()!
# # Y: Vektor (Y_1, ..., Y_n)
# # XXh: Matrix ((X_j - X_i)/h)_{1<=j<=n, 1<=i<=n}
#
# wx <- Wn( sigma, h, xXh, thetaXh, K) # s x n
#
# wX <- apply( XXh, 2, function( XXih) Wn( sigma, h, XXih, thetaXh, K))
# # (1<=r<=s)*(1<=j<=n) x (1<=i<=n)
# wX <- aperm( array( wX, dim = c( length( sigma), dim( XXh))), c( 2, 1, 3))
# # (1<=j<=n) x (1<=r<=s) x (1<=i<=n)
# Y_twXY <- Y - t( colSums( wX * Y)) # (1<=i<=n) x (1<=r<=s)
#
# rowSums( (wx * wx) * t( Y_twXY * Y_twXY)) # s x 1
# }
#
#
# # Zwei Versionen des geschaetzten "Asymptotic Mean Squared Errors"
# # (\widehat AMSE): AMSE_A() mit m_n, AMSE_B() mit \hat m_n
# #*****************************************************************
# AMSE_A <- function( sigma, h, xXh, thetaXh, K = Kern, mmDiff, YmDiff2) {
# # sigma, h, xXh, thetaXh, K: Siehe Wn()!
# # mmDiff: Vektor (m_n(X_1) - m_n(x), ..., m_n(X_n) - m_n(x))
# # YmDiff2: Vektor ( (Y_1 - m_n(x))^2, ..., (Y_n - m_n(x))^2 )
# W <- Wn( sigma, h, xXh, thetaXh, K) # s x n
# Bn <- W %*% mmDiff # = BiasA( sigma, h, xXh, thetaXh, K, mmDiff)
# Vn2 <- (W * W) %*% YmDiff2 # = VarA( sigma, h, xXh, thetaXh, K, YmDiff2)
# drop( Vn2 + Bn * Bn) # s x 1
# }
#
#
# AMSE_B <- function( sigma, Y, h, xXh, XXh, thetaXh, K = Kern) {
# # sigma, h, xXh, thetaXh, K: Siehe Wn()!
# # Y: Vektor (Y_1, ..., Y_n)
# # XXh: Matrix ((X_j - X_i)/h)_{1<=j<=n, 1<=i<=n}
#
# wx <- Wn( sigma, h, xXh, thetaXh, K) # s x n
#
# wX <- apply( XXh, 2, function( XXih) Wn( sigma, h, XXih, thetaXh, K))
# # (1<=r<=s)*(1<=j<=n) x (1<=i<=n)
# wX <- array( wX, dim = c( length( sigma), dim( XXh)))
# # (1<=r<=s) x (1<=j<=n) x (1<=i<=n)
#
# wX_wx <- aperm( wX - as.vector( wx), c( 2, 1, 3))
# # (1<=j<=n) x (1<=r<=s) x (1<=i<=n)
# Bn <- rowSums( wx * colSums( wX_wx * Y))
# # = BiasB( sigma, Y, h, xXh, XXh, thetaXh, K), s x 1
#
# Y_twXY <- Y - t( colSums( aperm( wX, c( 2, 1, 3)) * Y)) # (1<=i<=n) x (1<=r<=s)
# Vn2 <- rowSums( (wx * wx) * t( Y_twXY * Y_twXY))
# # = VarB( sigma, Y, h, xXh, XXh, thetaXh, K), s x 1
#
# drop( Vn2 + Bn * Bn) # s x 1
# }
#
#
#
# ## Die Simulations-Szenarien
# ##***************************
#
# ##Version <- "A" # Varianz- und Bias-Schaetzung mit NW-Schaetzer m_n fuer m
# Version <- "B" # Varianz- und Bias-Schaetzung mit \hat m_n fuer m
#
# # Simulationsparameter
# x00 <- c( 2, 4, 5, 8) # Punkte, an denen die AMSE-optim. "kernel adjusted
# # nonparametric regressions" von m stattfinden sollen.
# # Fuer m.poly4()'s Voreinstellung ist TP an 2, WP an 4
# # und SP an 8; ein "beliebiger" Punkt liegt z. B. bei 5.
# nn <- c( 25, 50, 100) # Stichprobenumfaenge.
# sdd <- c( 0.5, 1, 2) # Std.-abwn. der Fehler \epsilon_i: \sqrt(Var(Y|X=x)).
#
# Resultate <- array( NA, dim = c( length( x00), length( nn), length( sdd), 3, 3),
# dimnames = list( x0 = paste( "x0", x00, sep = "="),
# n = paste( "n", nn, sep = "="),
# sd = paste( "sd", sdd, sep = "="),
# Schaetzer = c( "Bias", "Var", "MSE"),
# Methode = c( "adjusted", "idealistic",
# "classic")) )
#
# mm <- m.poly4 # Regressionsfunktion und ihre erste sowie zweite Ableitung
# mm_prime <- function() {}; body( mm_prime) <- D( body( mm), "x")
# mm_second <- function() {}; body( mm_second) <- D( body( mm_prime), "x")
# formals( mm_second) <- formals( mm_prime) <- formals( mm)
#
# # Dichte der Normalverteilung (der X_i) und ihre erste Ableitung
# f <- dnorm; f_prime <- function( x) -x * dnorm( x)
#
#
# for( x0 in x00) { # Punkt x_0
#
# design.mean <- x0 + 0.5 # "Zentrum" und "Skala" des Designs = Erwartungswert
# design.sd <- 1 # bzw. Standardabweichung der Verteilung der X_i.
# xx <- design.mean + seq( -3.3 * design.sd, 3.3 * design.sd, length = 101)
# # x-Gitter, worauf m und der NW-Schaetzer m_n ausgewertet & gezeichnet werden.
#
# sig.range <- c( 0.01, 1.5 * diff( range( xx))) # Zu betrachtender Wertebereich
# sig.grid <- seq( sig.range[ 1], # samt Gitter fuer den Skalenparame-
# sig.range[ 2], length = 501)# ter \sigma des adaptiv. Verfahrens.
#
# for( n in nn) { # Stichprobenumfang
# for( sd.eps in sdd) { # Standardabweichung der \epsilon_i
#
# set.seed( 20110426) # Startwert des Pseudo-Zufallszahlengenerators
#
# # Berechnung der fuer x_0 optimalen, "idealistischen" Fensterbreite h_opt
# # zur Berechnung des NW-Schaetzers an x_0
# numerator <- sd.eps^2 * Integral_Kquadrat * f( x0)
# denominator <- n * Integral_xquadratK^2 *
# (f( x0) * mm_second( x0) + 2 * f_prime( x0) * mm_prime( x0))^2
# h_opt <- (numerator / denominator)^(1/5)
#
# # Wahl von h fuer das adaptive Verfahren
# h <- n^(-1/5)
#
# # Vorbereitung des Grafik-Ausgabe
# now <- format( Sys.time(), "%Y%m%d_%H%M%S")
# pdf( paste( "Plots/mBVA_x", x0, "_n", formatC( n, width = 3, flag = 0), "_s",
# formatC( 10*sd.eps, width = 2, flag = 0), "_D", now, ".pdf",
# sep = ""),
# paper = "a4", width = 0, height = 0)
# op <- par( mfrow = c( 3,1), mar = c( 3,3,3,0.1)+0.1, mgp = c( 1.5,0.5,0),
# tcl = -0.3, cex.main = 2)
#
# ## Ein paar beispielhafte Simulationslaeufe mit Grafikgenerierung
# ##****************************************************************
# for( bsp in 1:4) {
# # Datengenerierung Y_i = m( X_i) + \epsilon_i, i = 1, ..., n
# X <- rnorm( n, mean = design.mean, sd = design.sd)
# Y <- mm( X) + rnorm( length( X), sd = sd.eps)
#
# # NW-Schaetzer auf einem x-Gitter, um ihn zur Anschauung zu zeichnen
# mnxx <- NadWat( xx, X, Y, h_opt)
#
# # Kern-Schaetzer mit lokal-adaptiver "plug-in"-Bandbreite (aus dem R-Paket
# # "lokern"), auf selbem x-Gitter, um auch ihn zur Anschauung zu zeichnen
# loxx <- lokerns( x = X, y = Y, x.out = xx)$est
#
# # Vorratsberechnung von Groessen, die nicht von \sigma abhaengen, aber
# # i. F. wiederholt gebraucht werden; Algor. wird so effizienter.
# x0Xh <- (x0 - X) / h; thetaXh <- (mean( X) - X) / h
#
# if( Version == "A") { # Var.- & Bias-Schaetzung mit m_n
#
# # NW-Schaetzer m_n(x0) und m_n(X_1), ..., m_n(X_n) mit Fensterbreite
# # h = n^(-1/5) sowie daraus gebildete Groessen zur Schaetzung von
# # Var_x0(\sigma) und Bias_x0(\sigma)
# mnx0 <- NadWat( x0, X, Y, h) # m_n(x0)
# mnX <- NadWat( X, X, Y, h) # m_n(X_i) fuer i = 1, ..., n
# mnX_mnx0 <- mnX - mnx0 # m_n(X_i) - m_n(x0) fuer i = 1, ..., n
# Y_mnX2 <- (Y - mnX)^2 # (Y_i - m_n(X_i))^2 fuer i = 1, ..., n
#
# # Schaetzer von Var_x0(\sigma) und Bias_x0(\sigma) auf dem \sigma-Gitter,
# # um sie zur Anschauung zu zeichnen
# Bn <- BiasA( sigma = sig.grid, h = h, xXh = x0Xh, thetaXh = thetaXh,
# mmDiff = mnX_mnx0)
# Vn2 <- VarA( sigma = sig.grid, h = h, xXh = x0Xh, thetaXh = thetaXh,
# YmDiff2 = Y_mnX2)
#
# # # Versuch 1: Minimizer und Minimum des Schaetzers von AMSE_x0(\sigma):
# # opt <- optimize( f = AMSE_A, interval = sig.range, h = h, xXh = x0Xh,
# # thetaXh = thetaXh, mmDiff = mnX_mnx0, YmDiff2 = Y_mnX2)
# # sig_opt_x0 <- opt$minimum; amse_min <- opt$objective
# #
# # # Versuch 2: Minimizer und Minimum des Schaetzers von AMSE_x0(\sigma):
# # opt <- optimize( f = AMSE_A, interval = sig.range/c(1,3), h = h, xXh = x0Xh,
# # thetaXh = thetaXh, mmDiff = mnX_mnx0, YmDiff2 = Y_mnX2)
# # sig_opt_x0_2 <- opt$minimum; amse_min_2 <- opt$objective
#
# } else if( Version == "B") { # Var.- & Bias-Schaetzung mit \hat m_n
#
# # Weitere Vorratsgroesse, die nicht von \sigma abhaengt, aber i. F.
# # wiederholt gebraucht wird
# XjXih <- outer( X/h, X/h, "-")
#
# # Schaetzer von Var_x0(\sigma) und Bias_x0(\sigma) auf dem \sigma-Gitter,
# # um sie zur Anschauung zu zeichnen
# Bn <- BiasB( sigma = sig.grid, Y = Y, h = h, xXh = x0Xh, XXh = XjXih,
# thetaXh = thetaXh)
# Vn2 <- VarB( sigma = sig.grid, Y = Y, h = h, xXh = x0Xh, XXh = XjXih,
# thetaXh = thetaXh)
#
# # # Versuch 1: Minimizer und Minimum des Schaetzers von AMSE_x0(\sigma):
# # opt <- optimize( f = AMSE_B, interval = sig.range, Y = Y, h = h,
# # xXh = x0Xh, XXh = XjXih, thetaXh = thetaXh)
# # sig_opt_x0 <- opt$minimum; amse_min <- opt$objective
# #
# # # Versuch 2: Minimizer und Minimum des Schaetzers von AMSE_x0(\sigma):
# # opt <- optimize( f = AMSE_B, interval = sig.range/c(1,3), Y = Y, h = h,
# # xXh = x0Xh, XXh = XjXih, thetaXh = thetaXh)
# # sig_opt_x0_2 <- opt$minimum; amse_min_2 <- opt$objective
#
# } else stop( "Version '", Version, "' gibt es nicht!")
#
# # # Versuch 2 "besser" als Versuch 1?
# # if( amse_min_2 < amse_min) {
# # sig_opt_x0 <- sig_opt_x0_2; amse_min <- amse_min_2
# # }
#
# # "Komposition" von AMSE_x0(\sigma) aus den Schaetzern von Var_x0(\sigma)
# # und Bias_x0(\sigma) auf dem \sigma-Gitter
# An <- Vn2 + Bn * Bn
#
# # Minimizer und Minimum des Schaetzers von AMSE_x0(\sigma)
# # durch "grid search":
# min_idx <- which.min( An)
# sig_opt_x0 <- sig.grid[ min_idx]; amse_min <- An[ min_idx]
#
# # Der AMSE-optimale "kernel adjusted nonparametric regression"-Schaetzwert
# # an x_0, also \hat m_n( x_0)
# W <- Wn( sigma = sig_opt_x0, h = h, xXh = x0Xh, thetaXh = thetaXh)
# mnhatx0 <- W %*% Y
#
# # Grafische Darstellung fuer den aktuellen Datensatz
# plot( Y ~ X, xlim = range( X, xx),
# ylim = range( Y, mnxx, mnhatx0, loxx, na.rm = TRUE),# mm(xx)
# main = bquote( n == .(n)), xlab = "x", ylab = "y") # Das Streudiagramm
# lines( xx, mm( xx), lty = 2) # Die "wahre" Regressionsfunktion m
# lines( xx, mnxx) # Der NW-Schaetzer m_n von m
# lines( xx, loxx, col = "blue") # Der lokal-adaptive Schaetzer von m
# points( x0, mnhatx0, col = "red", pch = 8, cex = 2) # (x_0, \hat m_n(x_0))
# legend( "bottomleft", lty = c( 2, 1, 1), col = c( "black", "black", "blue"),
# bty = "n", cex = 1.5,
# legend = c( bquote( paste( "m, ", sigma(paste( Y,"|", X==x))
# == .(sd.eps)) ),
# expression( m[n]), "lokern()"))
# legend( "top", pch = 8, col = "red", bty = "n", cex = 1.5,
# legend = bquote( paste( hat(m)[n](x[0]), ", ", x[0] == .(x0))))
#
# matplot( sig.grid, cbind( Bn*Bn, Vn2), type = "l", lty = 1:2,
# col = c( "black", "red"), xlab = expression( sigma), ylab = "")
# legend( "top", lty = 1:2, col = c( "black", "red"), bty = "n", cex = 1.5,
# legend = c( expression( paste( widehat(plain(Bias))[x[0]]^2, (sigma))),
# expression( widehat(plain(Var))[x[0]](sigma))))
#
# plot( sig.grid, An, type = "l", xlab = expression( sigma), ylab = "",
# main = expression( widehat(plain(AMSE))[x[0]](sigma)))
# points( sig_opt_x0, amse_min, pch = 13, col = "red", cex = 2)
# legend( "top", pch = 13, col = "red", bty = "n", cex = 1.5,
# legend = substitute( group( "(", list( plain( Minimizer),
# plain( Minimum)), ")")
# == group( "(", list( x, y), ")"),
# list( x = signif( amse_min, 4),
# y = signif( sig_opt_x0, 4))))
# } # bsp-loop
# par( op)
# dev.off() ## Ende der beispielhaften Simulationslaeufe
#
#
# ## Die eigentliche Simulation (Kommentierung siehe oben!)
# ##****************************
# MCsamplesize <- 5000
#
# now <- format( Sys.time(), "%Y%m%d_%H%M%S")
# set.seed( 20110426) # Derselbe Startwert des Pseudo-RNGs wie oben.
# Sim <- replicate( MCsamplesize, { # Schleife ueber Monte-Carlo-Samples
# X <- rnorm( n, mean = design.mean, sd = design.sd)
# Y <- mm( X) + rnorm( length( X), sd = sd.eps)
#
# if( n > 25 || sd.eps == 1) {
# c( mnhatx0 = NA, mnx0_hopt = NA, lox0 = NA, lox0.h = NA)
# } else {
# mnx0_hopt <- NadWat( x0, X, Y, h_opt)
# lokerns.x0 <- lokerns( x = X, y = Y, x.out = x0)
# lox0 <- lokerns.x0$est; lox0.h <- lokerns.x0$bandwidth
#
# x0Xh <- (x0 - X) / h; thetaXh <- (mean( X) - X) / h
#
# if( Version == "A") {
# mnx0 <- NadWat( x0, X, Y, h); mnX <- NadWat( X, X, Y, h)
# mnX_mnx0 <- mnX - mnx0; Y_mnX2 <- (Y - mnX)^2
# # # Minimierungsversuch 1:
# # opt <- optimize( f = AMSE_A, interval = sig.range, h = h, xXh = x0Xh,
# # thetaXh = thetaXh, mmDiff = mnX_mnx0, YmDiff2 = Y_mnX2)
# # sig_opt_x0 <- opt$minimum; amse_min <- opt$objective
# # # Minimierungsversuch 2:
# # opt <- optimize( f = AMSE_A, interval = sig.range/c(1,3), h = h, xXh = x0Xh,
# # thetaXh = thetaXh, mmDiff = mnX_mnx0, YmDiff2 = Y_mnX2)
# # sig_opt_x0_2 <- opt$minimum; amse_min_2 <- opt$objective
#
# An <- AMSE_A( sig.grid, h = h, xXh = x0Xh, thetaXh = thetaXh,
# mmDiff = mnX_mnx0, YmDiff2 = Y_mnX2)
#
# } else if( Version == "B") {
# XjXih <- outer( X/h, X/h, "-")
# # # Minimierungsversuch 1:
# # opt <- optimize( f = AMSE_B, interval = sig.range, Y = Y, h = h,
# # xXh = x0Xh, XXh = XjXih, thetaXh = thetaXh)
# # sig_opt_x0 <- opt$minimum; amse_min <- opt$objective
# # # Minimierungsversuch 2:
# # opt <- optimize( f = AMSE_B, interval = sig.range/c(1, 3), Y = Y, h = h,
# # xXh = x0Xh, XXh = XjXih, thetaXh = thetaXh)
# # sig_opt_x0_2 <- opt$minimum; amse_min_2 <- opt$objective
#
# An <- AMSE_B( sig.grid, Y = Y, h = h, xXh = x0Xh, XXh = XjXih,
# thetaXh = thetaXh)
#
# } else stop( "Version '", Version, "' gibt es nicht!")
#
# # # Minimierungsversuch 2 "besser" als 1?
# # if( amse_min_2 < amse_min) sig_opt_x0 <- sig_opt_x0_2
#
# # Minimizer und Minimum des Schaetzers von AMSE_x0(\sigma)
# # durch "grid search":
# sig_opt_x0 <- sig.grid[ which.min( An)]
#
# W <- Wn( sigma = sig_opt_x0, h = h, xXh = x0Xh, thetaXh = thetaXh)
# mnhatx0 <- W %*% Y
#
# c( mnhatx0 = mnhatx0, mnx0_hopt = mnx0_hopt, lox0 = lox0, lox0.h = lox0.h)
# }
# } )
#
# lox0.meanh <- mean( Sim[ "lox0.h",]); Sim <- Sim[ -4,]
# mx0 <- mm( x0); Means <- rowMeans( Sim); Biases <- Means - mx0
# Variances <- apply( Sim, 1, var); MSEs <- Biases * Biases + Variances
#
# Resultate[ paste( "x0", x0, sep = "="), paste( "n", n, sep = "="),
# sd = paste( "sd", sd.eps, sep = "="),
# c( "Bias", "Var", "MSE"), ] <- rbind( Biases, Variances, MSEs)
#
# if( !(n > 25 || sd.eps == 1)) {
# ## Grafische Darstellung der Resultate (Boxplots)
# ##************************************************
# Bim1 <- as.list( as.data.frame( t( Sim)))
# Bim2 <- list( "'Idealistic' - Adjusted" = Sim[ "mnx0_hopt",] - Sim[ "mnhatx0",],
# "'Classic' - Adjusted" = Sim[ "lox0",] - Sim[ "mnhatx0",])
#
# pch <- c( 1, 4, 5); cols <- c( "red", "orange", "green"); cex <- 2
# pdf( paste( "Plots/AMSEopt_x", x0, "_n", formatC( n, width = 3, flag = 0),
# "_s", formatC( 10*sd.eps, width = 2, flag = 0), "_D", now,
# ".pdf", sep = ""),
# paper = "a4", width = 0, height = 0)
# par( pl <- list( oma = c( 0,0,8,0), mar = c( 3.5, 3.3, 1, 3)+.1,
# mgp = c( 2, 0.5, 0), tcl = -0.3, lab = c( 5, 10, 0)))
# boxplot( Bim1, names = c( "Adjusted ...", "'Idealistic' ...", "'Classic' ..."),
# ylab = expression( hat(m)[n](x[0])~~~~"bzw."~~~~m[n](x[0])
# ~~~~"bzw."~~~~lokerns(x[0])),
# cex.axis = 1.2, cex.lab = 1.2,
# xlab = expression( "kernel regression estimate values at"~x[0]))
#
# abline( h = mx0, col = "blue")
# axis( 4, at = mx0, tick = FALSE, labels = expression( m(x[0])), line = 0,
# las = 2, cex.axis = 1.2, col.axis = "blue")
#
# points( Means, col = cols, pch = pch, cex = cex, lwd = 3)
# legend( "top", pch = pch, col = cols, pt.lwd = 3, bty = "n", cex = 1.4,
# legend = c( expression( bar(hat(m)[n](x[0]))),
# expression( bar(m[n](x[0]))),
# expression( bar(lokerns(x[0]))) ), y.intersp = 1.7)
# TitleAndTopAxis( at = 1:3 - 0.18, cex.axis = 1.2)
#
#
# boxplot( Bim2, cex.axis = 1.1, cex.lab = 1.1,
# ylab = expression( m[n](x[0]) - hat(m)[n](x[0])~~~~~~"bzw."~~~~~~
# lokerns(x[0]) - hat(m)[n](x[0])),
# xlab = bquote( "Difference between kernel regression estimates at "~
# x[0] == .(x0)~"with"~
# sigma(paste(Y,"|",X==x)) == .(sd.eps)))
# abline( h = 0, col = "blue")
# points( c( Means[ "mnx0_hopt"] - Means[ "mnhatx0"],
# Means[ "lox0"] - Means[ "mnhatx0"]),
# col = c( "darkgreen", "brown"), pch = c( 8, 9), cex = cex, lwd = 2)
# legend( "topleft", pch = c( 8, 9, -1), col = c( "darkgreen", "brown", NA),
# pt.lwd = 2, bty = "n", cex = 1.4,
# legend = c( bquote( bar(m[n](x[0]) - hat(m)[n](x[0])) ==
# .(round( Means[ "mnx0_hopt"] - Means[ "mnhatx0"],
# 5))),
# bquote( bar(lokerns(x[0]) - hat(m)[n](x[0]))==
# .(round( Means[ "lox0"] - Means[ "mnhatx0"], 5))),
# expression() ), y.intersp = 1.7)
# TitleAndTopAxis( at = c( 0.833, 1.5, 2.165) - 0.12, cex.axis = 1.2)
# dev.off()
# }
# } # sd.eps-loop
# } # n-loop
# } # x0-loop
#
# # rm( list = ls()); dev.off()
#
#
# today <- format( Sys.time(), "%Y%m%d")
# sink( paste( "Resultate_", today, ".txt", sep = ""))
# print( ftable( Resultate))
# sink()
#
# #round( ftable( Resultate[ c( "x0=4", "x0=8", "x0=2"),
# # c( "n=25", "n=50", "n=100"),
# # c( "sd=0.5", "sd=2"), , ],
# # row.vars = c( "x0", "sd", "Schaetzer") ), 4)
#
# library( Hmisc)
# punkte <- c( "minimum" = "x0=2", "point of inflection" = "x0=4",
# "'arbitrary' point" = "x0=5", "saddle point" = "x0=8")
# Tab <- lapply( punkte, function( v)
# ftable( Resultate[ v, c( "n=50", "n=25", "n=100"),
# c( "sd=0.5", "sd=2"), , ],
# row.vars = c( "sd", "Schaetzer")) )
# names( Tab) <- paste( "Performance in", names( punkte),
# paste( "$(", gsub( "x0", "x_0", punkte, fixed = TRUE),
# ")$", sep = ""),
# "based on", MCsamplesize, "Monte Carlo samples" )
#
# CV <- attr( Tab[[1]], "col.vars")
# CV$n <- paste( "$", CV$n, "$", sep = "")
# CV$Methode2 <- c( "$\\hat m_n$", "$m_n$", "\\texttt{lokerns}")
#
# RV <- attr( Tab[[1]], "row.vars")
# RV$sd <- paste( "$", gsub( "sd", "\\sigma", RV$sd, fixed = TRUE), "$", sep = "")
# RV$Schaetzer <- paste( "$\\widehat{\\text{", RV$Schaetzer, "}}$", sep = "")
#
# w <- latex( object = Tab, file = paste( "Tab_", today, ".tex", sep = ""),
# title = "", where = "!htbp", dec = 3, first.hline.double = FALSE,
# cgroup = CV$n, colheads = rep( CV$Methode2, length( CV$n)),
# extracolheads = rep( CV$Methode, length( CV$n)),
# rgroup = RV$sd, rowname = rep( RV$Schaetzer, length( RV$sd))
# )
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.