library(tidyverse)
getErr <- function(true,measured){
paste0(round(mean(abs(measured-true)/true) * 100,digits = 1),
"%") %>% cat("\n")
}
indata <- tribble(
~Year,~Woodbasedpanels,~Sawnwood ,~Paperandpaperboard,
1960, NA ,454000 ,6200,
1961, 2500, 430000, 6200,
1962, 1400, 475000,12200,
1963, 700 ,510000 ,21100,
1964, 700 ,510000, 20200,
1965, 1700, 474000 ,21900)
vec_of_mean_errs <- c(0.05,0.1,0.2)
standard_deviations <- tribble(
~ Woodbasedpanels,~ Sawnwood, ~Paperandpaperboard,
1,100,10,
2,200,20,
5,5000,100
)
permute_vec <- function(vec,meanerr,sd){
newvec_p <- vec + rnorm(n=seq_along(vec),mean = mean(vec,na.rm=TRUE)*meanerr,sd=sd)
newvec_n <- vec - rnorm(n=seq_along(vec),mean = mean(vec,na.rm=TRUE)*meanerr,sd=sd)
# getErr(vec,newvec_p)
# getErr(vec,newvec_n)
pindex <- sample(x = c(TRUE,FALSE),
size = length(vec),
replace = TRUE)
newvec <- newvec_p
newvec[!pindex] <- newvec_n[!pindex]
# getErr(vec,newvec)
newvec
}
innerloop <- function(err){
suppressMessages({
map_dfc(names(standard_deviations),
~ permute_vec(indata[[.]],err,standard_deviations[[.]])) %>% setNames(names(standard_deviations)) %>%
mutate(err=err) %>% bind_cols(Year=indata$Year)
})
}
outerloop <- function(x){
cat(".")
if(x%%50==0){
cat("\n",x,"\n")
}
map_dfr(vec_of_mean_errs,innerloop)
}
results <- map_dfr(1:500
,outerloop)
indata_tidy <- pivot_longer(indata,
cols = -Year )
results_tidy <- pivot_longer(results,
cols = c(-Year,-err ))
ggplot(data=indata_tidy,
mapping = aes(x=Year,y=value))+geom_point(size=4,color="blue") + facet_wrap(vars(name,err),
nrow=3,
ncol=3,
scales="free")+
geom_point(data=results_tidy,color="red",alpha=0.01,size=.5)