# simulation

#constants
stdev<- array(NA, dim=c(20,2))
colnames(stdev)<-c("mean","median")
zvalue <- array(NA, dim=c(712*712,20))
average <- array(NA, dim=c((712),(712),1))
median<- array(NA, dim=c((712),(712),1))

#random numbers
for (i in 1:20) {
rnorm(712*712, mean=0, sd=1)->zvalue[,i]
}
dim(zvalue)<-c(712,712,20)

#simulation
for (i in 1:20) {
for (c_row in 3:(712-3)) {
for (c_col in 3:(712-3)) {
b<-(zvalue[(c_row-2):(c_row+2),(c_col-2):(c_col+2),i])
bna<-b[!is.na(b)]
if(length(bna)>(5*5*0.75)) {
bmed<-(median(bna, na.rm=T))
bave<-(mean(bna, na.rm=T))
median[c_row, c_col, 1]<- bmed
average[c_row, c_col, 1]<- bave
} } }
stdev[i,1] <- sd(average, na.rm=T)
stdev[i,2]<- sd(median, na.rm=T)
}

#output
cat("average=", mean(stdev[,2]),"\n")
stdev


# 01 Sep. 2006
# Tomokazu Konishi