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[1:31];    summary01 <- summary(mo01)
  mo02 <- t[32:59];   summary02 <- summary(mo02)
  mo03 <- t[60:90];   summary03 <- summary(mo03)
  mo04 <- t[91:120];  summary04 <- summary(mo04)
  mo05 <- t[121:151]; summary05 <- summary(mo05)
  mo06 <- t[152:181]; summary06 <- summary(mo06)
  mo07 <- t[182:212]; summary07 <- summary(mo07)
  mo08 <- t[213:243]; summary08 <- summary(mo08)
  mo09 <- t[244:273]; summary09 <- summary(mo09)
  mo10 <- t[274:304]; summary10 <- summary(mo10)
  mo11 <- t[305:334]; summary11 <- summary(mo11)
  mo12 <- t[335:365]; summary12 <- summary(mo12)
  month.means <- round(c(summary01[4],summary02[4],summary03[4],
                         summary04[4],summary05[4],summary06[4],
                         summary07[4],summary08[4],summary09[4],
                         summary10[4],summary11[4],summary12[4]),1)
  month.mins <- c(summary01[1],summary02[1],summary03[1],
                  summary04[1],summary05[1],summary06[1],
                  summary07[1],summary08[1],summary09[1],
                  summary10[1],summary11[1],summary12[1])
  month.maxs <- c(summary01[6],summary02[6],summary03[6],
                  summary04[6],summary05[6],summary06[6],
                  summary07[6],summary08[6],summary09[6],
                  summary10[6],summary11[6],summary12[6])
  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 [°C]")
    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[1:12],month.means, type="s", lwd=2, lend="square",
      col="darkred")
    lines(x=c(x1.at[12:13]),y=c(month.means[12],month.means[12]),
      type="l", lwd=2, lend="square", col="darkred")
    points(x1.at[1:12],month.maxs, type="s", col="red", lwd=1,
      lty = "dotted")
    lines(x=c(x1.at[12:13]),y=c(month.maxs[12],month.maxs[12]),
      type="l", lwd=1, lty = "dotted", lend="square", col="red")
    points(x1.at[1],month.maxs[1], type="h", lwd=1, lty = "dotted",
      lend="square", col="red")
    points(x1.at[13],month.maxs[12], type="h", lwd=1, lty = "dotted",
      lend="square", col="red")
    points(x1.at[1:12],month.mins, type="s", lwd=1, lty = "dotted",
      lend="square", col="blue")
    lines(x=c(x1.at[12:13]),y=c(month.mins[12],month.mins[12]),
      type="l", lwd=1, lty = "dotted", lend="square", col="blue")
    points(x1.at[1],month.mins[1], type="h", lend="square", lwd=1,
      lty = "dotted", col="blue")
    points(x1.at[13],month.mins[12], 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