###################################################
# 'online' alogrithm for the mean
# reading example file 
# 1981-2010 monthly mean tavg anomalies
###################################################

# clear memory
rm(list=ls())

###################################################
# User's program control parameters
###################################################
# station
stationname<-"USW00014735"
# variable
varname<-"tavg"
# first and last year for averaging
yr1<-1981
yr2<-2010

mymean<-function(x,m=NA,n=0){
  # input is a single sample value
  # m is a previous estimate of the
  # sample mean
  # n is the counter of the sample
  # if it is called the first time
  # without a previous estimate
  # it returns x as a first estimate
  if(n>0) {
    n<-n+1
    m<-m+1/n*(x-m)
  } else {
    n<-1
    m<-x
  }
  res<-list(m=m,n=n)
  return(res)
}







###################################################
# Main program 
###################################################
source("scripts/myfunctions.R")
infile<-paste("data/",stationname,"_",varname,"_mon_mean_ano.csv",sep='')

buffer<-read.csv(infile)
# data are in column 2
time<-buffer[,2]
x<-buffer[,3]
inum<-length(x)
test<-array(NA,dim=c(inum))
# call first time 
res<-mymean(x[1])
print(res)
test[1]<-res$m
for ( n in 2:inum) {
  old<-res
  res<-mymean(x[n],m=old$m,n=old$n)
  print(paste("n=,",n,"mean=",round(res$m,4),"after n-samples",res$n))
  test[n]<-res$m
}

plot(time,x,typ='l',xlab='time',ylab='tavg ano [C]')
lines(time,test,typ='b',col=2)