Datei:Jahresmittel Simulation 3.svg

d          <- seq(1:365)
nsim       <- 1E4
bias       <- 0
bias.lo    <- 0
bias.hi    <- 0
plot.image <- FALSE
for(i in 1:nsim){
  t <- 15*sin(2*pi*d/365-2) + 15
  noise <- rnorm(365, mean=0, sd=2.5)
  t <- round(t + noise,1)
  annual.mean1 <- round(mean(t),1)
  mo01 <- t;    summary01 <- summary(mo01)
  mo02 <- t;   summary02 <- summary(mo02)
  mo03 <- t;   summary03 <- summary(mo03)
  mo04 <- t;  summary04 <- summary(mo04)
  mo05 <- t; summary05 <- summary(mo05)
  mo06 <- t; summary06 <- summary(mo06)
  mo07 <- t; summary07 <- summary(mo07)
  mo08 <- t; summary08 <- summary(mo08)
  mo09 <- t; summary09 <- summary(mo09)
  mo10 <- t; summary10 <- summary(mo10)
  mo11 <- t; summary11 <- summary(mo11)
  mo12 <- t; summary12 <- summary(mo12)
  month.means <- round(c(summary01,summary02,summary03,
                         summary04,summary05,summary06,
                         summary07,summary08,summary09,
                         summary10,summary11,summary12),1)
  month.mins <- c(summary01,summary02,summary03,
                  summary04,summary05,summary06,
                  summary07,summary08,summary09,
                  summary10,summary11,summary12)
  month.maxs <- c(summary01,summary02,summary03,
                  summary04,summary05,summary06,
                  summary07,summary08,summary09,
                  summary10,summary11,summary12)
  annual.mean2 <- round(mean(month.means),1)
  bias.actual <- 100*(annual.mean2-annual.mean1)/annual.mean1
  bias <- bias + bias.actual
  if(bias.actual < 0){bias.lo <- bias.lo + 1}
    else{bias.hi <- bias.hi + 1}
  if(plot.image){
    par(bg = "transparent")
    plot(d,t, type="n", xlim=c(1,365), ylim=c(-10,40), axes = FALSE,
    frame.plot = FALSE, xlab = "day #", ylab = "t ")
    x1.at <- c(1,31,59,90,120,151,181,212,243,273,304,334,365)
    axis(1, at = x1.at, lwd = 2, lwd.ticks = 1,
    labels = formatC(x1.at, format="fg"))
    x2.at <- c(16,46,76,106,137,167,198,229,259,290,320,351)
    axis(3, at = x1.at, labels = FALSE, lwd = 2, lwd.ticks = 1)
    axis(3, at = x2.at, tick = FALSE,
    labels = c("Jan","Feb","Mar", "Apr","May","Jun",
      "Jul","Aug","Sep","Oct","Nov","Dec"))
    axis(2, lwd = 2, lwd.ticks = 1)
    axis(4, lwd = 2, lwd.ticks = 1)
    lines(x=c(1,365), y=c(0,0))
    lines(d,t, type="s", lend="square", col="darkorange")
    points(x1.at,month.means, type="s", lwd=2, lend="square",
      col="darkred")
    lines(x=c(x1.at),y=c(month.means,month.means),
      type="l", lwd=2, lend="square", col="darkred")
    points(x1.at,month.maxs, type="s", col="red", lwd=1,
      lty = "dotted")
    lines(x=c(x1.at),y=c(month.maxs,month.maxs),
      type="l", lwd=1, lty = "dotted", lend="square", col="red")
    points(x1.at,month.maxs, type="h", lwd=1, lty = "dotted",
      lend="square", col="red")
    points(x1.at,month.maxs, type="h", lwd=1, lty = "dotted",
      lend="square", col="red")
    points(x1.at,month.mins, type="s", lwd=1, lty = "dotted",
      lend="square", col="blue")
    lines(x=c(x1.at),y=c(month.mins,month.mins),
      type="l", lwd=1, lty = "dotted", lend="square", col="blue")
    points(x1.at,month.mins, type="h", lend="square", lwd=1,
      lty = "dotted", col="blue")
    points(x1.at,month.mins, type="h", lend="square", lwd=1,
      lty = "dotted", col="blue")
    lines(x=c(1,365),y=c(annual.mean1,annual.mean1), col="darkgreen")
    lines(x=c(1,365),y=c(annual.mean2,annual.mean2), col="red")
    legend(x="topleft",
      c("green: mean of daily values or weighted mean of monthly means",
      "red: mean of monthly means"), bty = "n", cex=0.9)
    }
  }
cat("\n", nsim, "simulated datasets\n",
  "Bias = ", round(bias/nsim,4), "%\n",
  "Bias (negative) = ", round(100*bias.lo/nsim,4), "% of datasets\n",
  "Bias (positive) = ", round(100*bias.hi/nsim,4), "% of datasets\n")

nsim <- 1

plot.image = TRUE