RMSE and trend fit relationship - https://github.com/USGS-R/mda.lakes/issues/91

RMSE and trend fit relationship

# work-through of https://github.com/USGS-R/mda.lakes/issues/91

set.seed(443)

trend = 0.04 #in units per year, kinda like wtemp
freq = 26 #sample freq
year = (1:(30*freq)/freq)

yobs = 10*sin(year*2*pi)+ 10 + trend*year + rnorm(length(year), sd=0)
ymod = 10*sin(year*2*pi) + 10 + (trend*max(year)/2)

plot(year, yobs)
fit <- stl(ts(yobs, frequency = 26), "periodic")
trend_fit <- fit$time.series[,"trend"]
lines(year, trend_fit, col = "red")
slope <- (trend_fit[length(trend_fit)] - trend_fit[1]) / year[length(year)]


plot(year, yobs, main = round(sqrt(mean((yobs - ymod)^2)), 2))
lines(year, ymod, col = "red")

#####
yobs2 <- yobs + rnorm(length(year), sd = 0.6)

plot(year, yobs2, main = round(sqrt(mean((yobs2 - ymod)^2)), 2))
lines(year, ymod, col = "red")

#####
par(mfrow = c(2, 1))
yobs3 <- yobs + rnorm(length(year))
title <- paste0("RMSE=", round(sqrt(mean((yobs3 - ymod)^2)), 2),
                " ME=", round(mean(yobs3 - ymod), 3),
                " MAE=", round(mean(abs(yobs3 - ymod)), 3),
                " MAPE=", round(mean(abs((yobs3 - ymod) / yobs3)), 3),
                " RMSE w trend=", round(sqrt(mean((trend_fit - ymod)^2)), 2))
plot(year, yobs3, main = title)
lines(year, ymod, col = "red")

# hist(yobs3 - ymod)

#####

ymod2 = 10*sin(year*2*pi) + 10 + trend*year
title <- paste0("RMSE=", round(sqrt(mean((yobs3 - ymod2)^2)), 2),
                " ME=", round(mean(yobs3 - ymod2), 3),
                " MAE=", round(mean(abs(yobs3 - ymod2)), 3),
                " MAPE=", round(mean(abs((yobs3 - ymod2) / yobs3)), 3),
                " RMSE w trend=", round(sqrt(mean((trend_fit - ymod2)^2)), 2))
plot(year, yobs3, main = title)
lines(year, ymod2, col = "red")

# hist(yobs3 - ymod2)

#####