
library(doBy)

#address folder
setwd("C:/Users/Elena/Desktop/ILTER/EjemploModel")

###############################################################################
##############################INPUTS###################################

rm(list=ls(all=TRUE))

InventoryData<-read.csv("InventoryDataFCP.csv")

###############################################################################
############################## PARAMETERS MODELS ###################################

####### Error in meassurement of tree diameter
MeassurementeErrorDBH<-0.5

####### Errors in biomass regresion
#Model: exp(a+b*ln(DBH))
#parameter "a"
a<-(-0.46213)

#parameter "b"
b<-(1.65874)

#MSE of the model
MSE<-(0.2137)

#sample size (of sample used to estimate the allometric model)
n<-(46)

#Mean x (of DBH of sample used to estimate the allometric model)
MeanX<-(3.7016)

#sum of squared deviations of x
SdX<-(80.2768)

###############################################################################
############################## SIMULATION OF BIOMASS###################################

#Number of simulations
k<-100

# A vector to save biomass simulations
IterationVector<-as.data.frame(c(1,2,3,4,5,6,7,8,9,10))
colnames(IterationVector)<-"Plot"

####### Simulations of biomass at plot level
for (i in 1:k){
  attach(InventoryData)
  for (i in 1:length(Tree.diameter)){
  # Tree diameter with meassurement error
  InventoryData$DBHplusMeassuremntError[i]<-rnorm(1,Tree.diameter[i],MeassurementeErrorDBH)
  }
  detach(InventoryData)

  # Natural Logaritm of "DBH + mesurement error"
  InventoryData$lnDBHerror<-log(InventoryData$DBHplusMeassuremntError)

  attach(InventoryData)
  for (i in 1:length(Tree.diameter)){
  # Error in prediction of the mean
  InventoryData$ErrorPredictionMean[i]<-rnorm(1,0,sqrt(MSE*(1/n+
                                      ((InventoryData$lnDBHerror[i]-MeanX)^2/(SdX)))))

  # Error in prediction of an indivual
  InventoryData$ErrorPredictionIndividual[i]<-rnorm(1,0,sqrt(MSE*(1+1/n+
                                      ((InventoryData$lnDBHerror[i]-MeanX)^2/(SdX)))))
  }
  detach(InventoryData)

  # Error in both sources
  InventoryData$BothSources<-(InventoryData$ErrorPredictionMean + 
                              InventoryData$ErrorPredictionIndividual)

  # Natural Logaritm of "Biomass estimation with two sources of error"
  InventoryData$Biomas2sourcesErrors<-(a+b*InventoryData$lnDBHerror+InventoryData$BothSources)

  # "Biomass estimation with two sources of error" (kg)
  InventoryData$Biomas2sourcesErrors<-exp(InventoryData$Biomas2sourcesErrors)

  #Estimation of total biomass at Plot level
    Fun_Sum <- function(x){c(Sum=sum(x[!is.na(x)]))}
    IterationTable<-0
    IterationTable<-summaryBy(InventoryData$Biomas2sourcesErrors/0.04/1000 ~ 
                    InventoryData$Plot, data=InventoryData,
                    FUN=Fun_Sum,keep.names=TRUE, var.names=c("BiomassSimulation_TonHa"))

  #Merge of simulations
  IterationVector<-merge(IterationVector,IterationTable,by.x="Plot",by.y="Plot")

}

####### Estimation of uncertainty of biomass at plot levels
#Estimation of quantiles 2.5% and 97.5%
for (j in 1:10){
IterationVector$MeanTonHectare[j]<-sum(IterationVector[j,2:(k+1)])/k
IterationVector$Quantile2_5[j]<-quantile(IterationVector[j,2:(k+1)],0.025, na.rm=TRUE)
IterationVector$Quantile97_5[j]<-quantile(IterationVector[j,2:(k+1)],0.975, na.rm=TRUE)
}
#Estimation of uncertinty at plot level
IterationVector$UinfPercentage<-abs(((IterationVector$MeanTonHectare)-IterationVector$Quantile2_5)/
                      (IterationVector$MeanTonHectare))*100
IterationVector$UsupPercentage<-abs(((IterationVector$MeanTonHectare)-IterationVector$Quantile97_5)/
                      (IterationVector$MeanTonHectare))*100

####### Resutls and Saving
#Selection of columns
TableIterationVector<-IterationVector[,c(1,(k+2),(k+3),(k+4),(k+5),(k+6))]
TableIterationVector

#Saving results
write.csv(TableIterationVector, file = "TableIterationVector.csv")

####### END



