####################All the functions used for computing the complete monotone LSE


#-----------------------------------------------------------
# This grid is a discretization of the interval [0,1].
# It should be named grid but I just thought it would be OK
# to call it "interval". Maybe we do not need a step
# as small as 10^{-5}. This can be given as an option to the user.
#-----------------------------------------------------------

interval = seq(0,1,10^{-5})

#--------------------------------------------------------------
#This function empirical distribution based on 
# observed data "x". It returns a vector of length equal 
# to max(x) +1 because we deal here with data of the form 
# 0,1,..., max(x). If some value is not observed, the function
#returns 0.
#--------------------------------------------------------------

Emp <- function(x) {  
  # compute the empirical frequence of x
  n = length(x)
  supp <- seq(0,max(x),by=1)
  count <- rep(0,length(supp))
  for(i in 1:length(count)) {
    count[i] <- sum(x==supp[i])
  }

  return(list(supp=supp,freq=count/n))


}

#----------------------------------------------------------
#For given observed data "x" and given candidate "alpha" in
#[0,1] (or rather the grid above),  the following function
# gives the corresponding optimal weight: the function
#returns "\hat{pi}" such that the completely monotone
# function "\hat{\pi} (1-\alpha) \alpha^i, i \in \NN", is the
#the LS projection of the empirical estimator obtained by
# "Emp" 
#---------------------------------------------------------

InitializeC = function(x, alpha){
V  = 0:max(x)
T = sum(Emp(x)$freq*alpha^V)

return((1+alpha) * T)

}


##########This intermediate/auxiliary function is used for computing the directional
##########derivative below


Interm = function(x,alpha){
V=0:max(x)
return(sum(Emp(x)$freq*alpha^V))
}


######Now, given a current candidate for the Least Squares estimtion problem
#########with support vector "S" and weights vector "C", the following function computes the directional
#########derivative of the LS criterion in the direction of a geoemtric pmf with prob "\alpha"
############Note that this function calls this "Interm" function to make things less cumbersome.

DirDeriv = function(alpha, x, S, C){
# C is the vector of weights  of the current iterate
# S is the support vector of the current iterate

T1 = sum(C* (1-S)/(1-S*alpha))
T2 = Interm(x=x, alpha=alpha)

return((1-alpha)* (T1 - T2))

}

}


##########We need to find the minimum in "\alpha"  of the directional derivative
##########defined above. That's what this function is doing. I have decide to look
##########for the minimum using my hands and not giving this task to some optimizing  
#############function

FindMinfun = function(x, S, C){

V = 0:max(x)
Out1 = rep(0, length(interval))
Out2 = rep(0, length(interval))

for(j in 1:length(S)){

Out1 = Out1 + C[j]*(1-S[j])/(1-S[j]*interval)
}


for(j in 1:length(V)){

Out2 = Out2 + Emp(x)$freq[j]*interval^{V[j]}
}

Out = (1-interval)*(Out1 - Out2)
Ind =match(min(Out), Out)

return(list(minimum=interval[Ind],objective=min(Out)))

}


##########The following function solves the LS problem in the "unconstrained" set, that is, 
##########for a given support vector "S", it finds the optimal weight vector "C" so that the
########## elements of "C" are allowed to be negative. The system to solve is linear: M C = B
##########where "M" is a matric which can be very ill-conditionned. Here, WE NEED HIGH PRECISION.
###########The method I use to solve the system can be changed of course. The only things that should
###########REMAIN are "M" and "B".


SolveUncons  = function(x, S){

M = matrix(0, length(S), length(S))
for(i in 1:length(S)){

 for(j in 1:length(S)){
 M[i,j]  = (1-S[j])/(1-S[i]*S[j])
  }

}
#M1=rep(1,length(S))%*%t(1-S)
#M2 = 1-S%*%t(S)
#M =M1/M2

B= rep(0, length(S))
for(i in 1:length(S)){
B[i]  = Interm(x=x, alpha=S[i])
}
svdM = svd(M)
uM = svdM$u
vM = svdM$v
diagM =diag(1/svdM$d)
Sol1 = solve(uM, B)
Sol2 = diagM%*%Sol1
Sol  = vM%*%Sol2

#Sol = solve(M, B)
Sol= as.vector(Sol)
List = list(Sol = Sol, Check1=as.vector(M%*%Sol)-B)

return(List$Sol)
}



################The following function makes the support reduction step works. When the unconstrained
##############solution is not permissible, that is, some elements of "C" are negative, the function
############## tells us which support point to remove from "S" by computing a special convex combination
############## of the previous permissible iterate and the new non-permissible one. 
##################The way this function works is described in comments I've put in the function itself

ReduceSuppfun = function(S1, C1, S2, C2){

S.m = c(S1, S2)
S.m = unique(sort(S.m))

# S1 will always be the support vector of the previous permissible iterate
#  S2 will be the current support vector. An example:

# S1=  c(0.2, 0.3, 0.4)
# S2 = c(0.1, 0.2, 0.3, 0.4)  (this means that the point 0.1 was added to S1
# Then S.m = c(0.1, 0.2, 0.3, 0.4) 

S2.m = rep(10, length(S.m)) 

# S2.m is a vector filled with the value 10. This choice is really arbitrary.
 
C1.m = rep(0, length(S.m))  # we will store weights in this vector

C2.m = rep(0, length(S.m))  # idem

ind1 = pmatch(S1, S.m)
ind2 = pmatch(S2, S.m)

# The vector ind1 identifies the indices where S.m has the same values as in S1
# The same thing for ind1
# In the example given above we have: ind1 = c(2,3,4) and ind2 = c(1,2,3,4)
# 
 


C1.m[ind1]  = C1   # we fill the right positions with the weights from C1
C2.m[ind2]  = C2   # idem

S2.m[ind2]  = S2    # we fill the right positions with the support points from S2

Lambda = C1.m/(C1.m-C2.m)  # we compute this ratio (componentwise) as the mimimal value of the positive
                           # elements of Lambda gives the right convex combination (which will delete
                           # one of the support point of S2

index = match(min(Lambda[Lambda > 0]), Lambda)  # We look for the position of the minimum

S2new  = S2.m[-index]   # We remove that support points from the augmented support vector corresponding to S2
S2new = S2new[S2new < 10]  # We get rid of the values 10 to get finally the reduced suppport vector 

return(S2new)
}


####This function is the main one which should give the completely monotone LSE


#-------------------------------------------------------------------------
#The Support reduction algorithm for computing the completely monotone LSE
# if epsilon > 10^{-9}, then one does not get stuck. So, typically one can
# choose epsilon no smaller than 10^{-9} but no bigger than 10^{-7}
#alpha0 can be given also to the user, the default is 0.1 
#-------------------------------------------------------------------------

CompMonLSE= function(x,alpha0,epsilon){

S0=alpha0   # any intial value of choice

C0 = InitializeC(x=x, alpha=alpha0)  # the initial optimal weight
#print(C0)

Resmin = FindMinfun(x=x, S=S0,C=C0)  # Looking for the optimal Dirac direction in which we will go 
thetamin=Resmin$minimum              # This is the new support point that will be added to "\alpha0" 
valmin = Resmin$objective

count =0
Count = 0
while(abs(valmin) > epsilon ){   # a stopping condition
  count <- count + 1	
	#cat("Main loup numb = ",count,"\n")  
      Snew=c(S0,thetamin)  # The new suport
            Snew = sort(Snew)
      
            
      Cnew=SolveUncons(x=x,S=Snew)  # Solve the unconstrained LS problem
      
      min.C <- min(Cnew) 
	#Count <- 0
	   while(min.C < 0){  # condition indicating that the unconstrained solution is not permissible
                             # and that we should go into the reduction step

		Count <- Count+1
		#cat("Sub loup numb = ",Count," of the main loop numb=", count, "\n")
            Snew = ReduceSuppfun(S1=S0, C1=C0, S2=Snew, C2=Cnew)
             
            if(length(Snew) > 1){  # two expressions for the solutions depending on whether we are left with 
                                   # more than 1 point or just one
            Cnew = SolveUncons(x=x, S=Snew)
             }
            else{
            Cnew=InitializeC(x=x, alpha=as.numeric(Snew))
            
            }  #  the program will go out of the inner loop (reduction loop) if 
               # and only if the solution is permissible

            min.C = min(Cnew)
             
            
        }#min.C
  
     S0= Snew  # uppdate the support vector
     C0= Cnew  # updfate the weight vector
    
     Resmin = FindMinfun(x=x, S=S0,C=C0)  # compute the new minimizer of the directional derivative
                                          # to build up a new support or to stop if the stopping 
                                          #crierion is met. 
    
     thetamin=Resmin$minimum
     valmin = Resmin$objective
     #print(valmin)
}#valmin

     return(list(S0=S0,C0=C0, valmin=valmin))

}


#--------------------------------------------------------------------------------------------------------


library(changepoint)   #### call the library changepoint

################  This function uses the other parametrization for a geomtric distribution
 
myrgeom = function(n, p){
rgeom(n=n,prob=1-p)
}



################################### The complete monotone estimator of "N" (given an eta)


EstimateInver1minp0.2 = function(S, C, eta=0.1){
List =list(S=S, C=C)
estim=  (1/eta-1) * sum(List$C[List$S <= eta]) +  sum(List$C[List$S > eta]*(1/List$S[List$S > eta]  - 1))

return(estim  +1)
}


###################################  Estimating the parameter eta

Selectioneta2 = function(Data, B){

Grid = seq(0, 0.5, 0.01)

n=length(Data)

Matx =  matrix(0, nrow=length(Grid), ncol=B)

for(b in 1:B){

#cat("b = ", b, "\n")

Rescm.b = CompMonLSE(x=sample(x=Data,size=n, prob=rep(1/n, n), replace=T),alpha0=0.1, epsilon=10^{-8})
S0cm.b = Rescm.b$S0
C0cm.b = Rescm.b$C0

Stat.b =  apply(X=matrix(Grid, ncol=1), FUN = cdf, MARGIN=1, S=S0cm.b,  C=C0cm.b) 


Matx[, b]  =  sqrt(n)*Stat.b

}

hateta = Grid[as.numeric(cpt.meanvar(apply(X=Matx, FUN=mean, MARGIN=1), class=FALSE)[1])]
#print(apply(X=Matx, FUN=mean, MARGIN=1))

return(hateta)
}
##################################  Coding the estimator of "N" based on the $k$-monotone model as done in Giguelay (2017a)

Nhatk = function(Obs, k){

Estim = Emp(Obs)$freq
Estim = Estim[-1]


if(length(Estim) > k){
Estim =Estim[1:k]
}
if(length(Estim) < k){

Estim = c(Estim, rep(0, k-length(Estim)))
}

EstimInv =  1- sum((-1)^{1:k}*choose(k, 1:k)*Estim) 
EstimNk = floor(length(Obs)*EstimInv)  
return(EstimNk)

}

########################################This function computes the statistics of the complete monotone and Giguelay's
#########################################estimator for k=2, 4, 10  (reported in Table 2 and 3) and also for k=20 (not
#########################################reported) 

SimSpeciesRichness = function(N=1000, M=200, Model="I"){
Nvec = NULL
#NvecEM = NULL
Nveck2 = NULL
Nveck6 = NULL
Nveck10 = NULL
Nveck20 = NULL

j=0

if(Model =="I"){
while(j < M){
j = j+1
cat("j out of 200 = ", j, "\n") 
#pb = 1-(0.7*(1-0.8) + 0.3*(1-0.1))
#S=rbinom(n=1, size=N,prob=pb)
Z = rbinom(n=N, size=1, prob=0.7)
X = Z*rgeom(n=N, prob=0.2)+(1-Z)*rgeom(n=N, prob=0.9)
X = X[X > 0]
hateta = try(Selectioneta2(Data=X-1, B=100))

while(inherits(hateta, "try-error"))
  {
    j = j-1
    #error handling code, maybe just skip this iteration using
    Z = rbinom(n=N, size=1, prob=0.7)
    X = Z*rgeom(n=N, prob=0.2)+(1-Z)*rgeom(n=N, prob=0.9)
    X = X[X > 0]
    hateta = try(Selectioneta2(Data=X-1, B=100))

  }

LSEout = CompMonLSE(x=X-1,alpha0=0.1,epsilon=10^{-8})
LSEsupp =  LSEout$S0
LSEweights= LSEout$C0
Nhat = floor(EstimateInver1minp0.2(S=LSEsupp, C=LSEweights, eta=hateta)*length(X))
Nvec = c(Nvec, Nhat)
Nveck2  =c(Nveck2, Nhatk(Obs=X, k=2))
Nveck6  =c(Nveck6, Nhatk(Obs=X, k=6))
Nveck10  =c(Nveck10, Nhatk(Obs=X, k=10))
Nveck20  =c(Nveck20, Nhatk(Obs=X, k=20))
#NvecEM = c(NvecEM, NhatEM(Obs=X))

}#M
 
}#Model "I"

if(Model =="II"){
while(j < M){
j = j+1
cat("j out of 200 = ", j, "\n") 
#pb = 1-(0.7*(1-0.8) + 0.3*(1-0.1))
#S=rbinom(n=1, size=N,prob=pb)
Z = rbinom(n=N, size=1, prob=0.5)
X = Z*rgeom(n=N, prob=0.15)+(1-Z)*rgeom(n=N, prob=0.85)
X = X[X > 0]
hateta = try(Selectioneta2(Data=X-1, B=100))
while(inherits(hateta, "try-error"))
  {
    j = j-1
    #error handling code, maybe just skip this iteration using
    Z = rbinom(n=N, size=1, prob=0.5)
    X = Z*rgeom(n=N, prob=0.15)+(1-Z)*rgeom(n=N, prob=0.85)
    X = X[X > 0]


    hateta = try(Selectioneta2(Data=X-1, B=100))

  }
LSEout = CompMonLSE(x=X-1,alpha0=0.1,epsilon=10^{-8})
LSEsupp =  LSEout$S0
LSEweights= LSEout$C0
Nhat = floor(EstimateInver1minp0.2(S=LSEsupp, C=LSEweights, eta=hateta)*length(X))

Nvec = c(Nvec, Nhat)
Nveck2  =c(Nveck2, Nhatk(Obs=X, k=2))
Nveck6  =c(Nveck6, Nhatk(Obs=X, k=6))
Nveck10  =c(Nveck10, Nhatk(Obs=X, k=10))
Nveck20  =c(Nveck20, Nhatk(Obs=X, k=20))
#NvecEM = c(NvecEM, NhatEM(Obs=X))



}#M
 
}#Model "II"

if(Model =="III"){
while(j < M){
j = j+1
cat("j out of 200 = ", j, "\n") 
#pb = 1-(0.7*(1-0.8) + 0.3*(1-0.1))
#S=rbinom(n=1, size=N,prob=pb)
Z = rbinom(n=N, size=1, prob=0.8)
X = Z*rpois(n=N, lambda=0.2)+(1-Z)*rpois(n=N, lambda=1.3)
X = X[X > 0]
hateta = Selectioneta2(Data=X-1, B=100)
while(inherits(hateta, "try-error"))
  {
    j = j-1
    #error handling code, maybe just skip this iteration using
    Z = rbinom(n=N, size=1, prob=0.8)
    X = Z*rpois(n=N, lambda=0.2)+(1-Z)*rpois(n=N, lambda=1.3)
    X = X[X > 0]
    hateta = try(Selectioneta2(Data=X-1, B=100))

  }

LSEout = CompMonLSE(x=X-1,alpha0=0.1,epsilon=10^{-8})
LSEsupp =  LSEout$S0
LSEweights= LSEout$C0
Nhat = floor(EstimateInver1minp0.2(S=LSEsupp, C=LSEweights, eta=hateta)*length(X))
Nvec = c(Nvec, Nhat)
Nveck2  =c(Nveck2, Nhatk(Obs=X, k=2))
Nveck6  =c(Nveck6, Nhatk(Obs=X, k=6))
Nveck10  =c(Nveck10, Nhatk(Obs=X, k=10))
Nveck20  =c(Nveck20, Nhatk(Obs=X, k=20))
#NvecEM = c(NvecEM, NhatEM(Obs=X))


}#M
 
}#Model "III"

if(Model =="IV"){
while(j < M){
j = j +1
cat("j out of 200 = ", j, "\n") 
#pb = 1-(0.7*(1-0.8) + 0.3*(1-0.1))
#S=rbinom(n=1, size=N,prob=pb)
Z = rbinom(n=N, size=1, prob=0.7)
X = Z*rnbinom(n=N, size=0.5, mu=0.8)+(1-Z)*rnbinom(n=N, size=2, mu=2)
X = X[X > 0]

hateta = try(Selectioneta2(Data=X-1, B=100))
while(inherits(hateta, "try-error"))
  {
    j = j-1
    #error handling code, maybe just skip this iteration using
    Z = rbinom(n=N, size=1, prob=0.7)
    X = Z*rnbinom(n=N, size=0.5, mu=0.8)+(1-Z)*rnbinom(n=N, size=2, mu=2)
    X = X[X > 0]

    hateta = try(Selectioneta2(Data=X-1, B=100))

  }

LSEout = CompMonLSE(x=X-1,alpha0=0.1,epsilon=10^{-8})
LSEsupp =  LSEout$S0
LSEweights= LSEout$C0
Nhat = floor(EstimateInver1minp0.2(S=LSEsupp, C=LSEweights, eta=hateta)*length(X))
Nvec = c(Nvec, Nhat)
Nveck2  =c(Nveck2, Nhatk(Obs=X, k=2))
Nveck6  =c(Nveck6, Nhatk(Obs=X, k=6))
Nveck10  =c(Nveck10, Nhatk(Obs=X, k=10))
Nveck20  =c(Nveck20, Nhatk(Obs=X, k=20))
#NvecEM = c(NvecEM, NhatEM(Obs=X))


}#M
 
}#Model "IV"


if(Model =="V"){
while(j < M){
j = j +1
cat("j out of 200 = ", j, "\n") 
#pb = 1-(0.7*(1-0.8) + 0.3*(1-0.1))
#S=rbinom(n=1, size=N,prob=pb)
X = SimuMixGeomn(n=N, alpha=c(0.2, 0.4, 0.8), prob=c(0.25, 0.5, 0.25))
X = X[X > 0]

hateta = try(Selectioneta2(Data=X-1, B=100))
while(inherits(hateta, "try-error"))
  {
    j = j-1
    #error handling code, maybe just skip this iteration using
    X = SimuMixGeomn(n=N, alpha=c(0.2, 0.4, 0.8), prob=c(0.25, 0.5, 0.25))
    X = X[X > 0]

    hateta = try(Selectioneta2(Data=X-1, B=100))

  }

LSEout = CompMonLSE(x=X-1,alpha0=0.1,epsilon=10^{-8})
LSEsupp =  LSEout$S0
LSEweights= LSEout$C0
Nhat = floor(EstimateInver1minp0.2(S=LSEsupp, C=LSEweights, eta=hateta)*length(X))
Nvec = c(Nvec, Nhat)
Nveck2  =c(Nveck2, Nhatk(Obs=X, k=2))
Nveck6  =c(Nveck6, Nhatk(Obs=X, k=6))
Nveck10  =c(Nveck10, Nhatk(Obs=X, k=10))
Nveck20  =c(Nveck20, Nhatk(Obs=X, k=20))


}#M
 
}#Model "V"

if(Model =="VI"){
while(j < M){
j = j +1
cat("j out of 200 = ", j, "\n") 
#0.4 Geom(0.25)  + 0.1 Geom(0.45)  + 0.3 Geom(0.6) + 0.2 Geom(0.9)

X = SimuMixGeomn(n=N, alpha=c(0.25, 0.45, 0.6, 0.9), prob=c(0.4, 0.1, 0.3, 0.2))
X = X[X > 0]

hateta = try(Selectioneta2(Data=X-1, B=100))
while(inherits(hateta, "try-error"))
  {
    j = j-1
    #error handling code, maybe just skip this iteration using
    X = SimuMixGeomn(n=N, alpha=c(0.25, 0.45, 0.6, 0.9), prob=c(0.4, 0.1, 0.3, 0.2))
    X = X[X > 0]


    hateta = try(Selectioneta2(Data=X-1, B=100))

  }

LSEout = CompMonLSE(x=X-1,alpha0=0.1,epsilon=10^{-8})
LSEsupp =  LSEout$S0
LSEweights= LSEout$C0
Nhat = floor(EstimateInver1minp0.2(S=LSEsupp, C=LSEweights, eta=hateta)*length(X))
Nvec = c(Nvec, Nhat)
Nveck2  =c(Nveck2, Nhatk(Obs=X, k=2))
Nveck6  =c(Nveck6, Nhatk(Obs=X, k=6))
Nveck10  =c(Nveck10, Nhatk(Obs=X, k=10))
Nveck20  =c(Nveck20, Nhatk(Obs=X, k=20))
#NvecEM = c(NvecEM, NhatEM(Obs=X))


}#M
 
}#Model "VI"

if(Model =="VII"){
while(j < M){
j = j +1
cat("j out of 200 = ", j, "\n") 
# Poisson of rate smaller than 2 - sqrt(2)

X = rpois(n=N, lambda=0.5)
X = X[X > 0]

hateta = try(Selectioneta2(Data=X-1, B=100))
while(inherits(hateta, "try-error"))
  {
    j = j-1
    #error handling code, maybe just skip this iteration using
    X = rpois(n=N, lambda=0.5)
    X = X[X > 0]




    hateta = try(Selectioneta2(Data=X-1, B=100))

  }

LSEout = CompMonLSE(x=X-1,alpha0=0.1,epsilon=10^{-8})
LSEsupp =  LSEout$S0
LSEweights= LSEout$C0
Nhat = floor(EstimateInver1minp0.2(S=LSEsupp, C=LSEweights, eta=hateta)*length(X))
Nvec = c(Nvec, Nhat)
Nveck2  =c(Nveck2, Nhatk(Obs=X, k=2))
Nveck6  =c(Nveck6, Nhatk(Obs=X, k=6))
Nveck10  =c(Nveck10, Nhatk(Obs=X, k=10))
Nveck20  =c(Nveck20, Nhatk(Obs=X, k=20))


}#M
 
}#Model "VII"

if(Model =="VIII"){
while(j < M){
j = j +1
cat("j out of 200 = ", j, "\n") 
# Poisson of rate smaller than 2 - sqrt(2)

X = rSpline(n=N, k=10, supp=10)
X = X[X > 0]

hateta = try(Selectioneta2(Data=X-1, B=100))
while(inherits(hateta, "try-error"))
  {
    j = j-1
    #error handling code, maybe just skip this iteration using
    X = rSpline(n=N, k=10, supp=10)
    X = X[X > 0]
   hateta = try(Selectioneta2(Data=X-1, B=100))

  }

LSEout = CompMonLSE(x=X-1,alpha0=0.1,epsilon=10^{-8})
LSEsupp =  LSEout$S0
LSEweights= LSEout$C0
Nhat = floor(EstimateInver1minp0.2(S=LSEsupp, C=LSEweights, eta=hateta)*length(X))
Nvec = c(Nvec, Nhat)
Nveck2  =c(Nveck2, Nhatk(Obs=X, k=2))
Nveck6  =c(Nveck6, Nhatk(Obs=X, k=6))
Nveck10  =c(Nveck10, Nhatk(Obs=X, k=10))
Nveck20  =c(Nveck20, Nhatk(Obs=X, k=20))


}#M
 
}#Model "VIII"

return(rbind(Nvec, Nveck2,Nveck6,Nveck10,Nveck20))

}


SimSpeciesRichnessModelI=SimSpeciesRichness(N=1000, M=200, Model ="I")
cat("Mean of cm =", mean(SimSpeciesRichnessModelI[1,]), "\n")
cat("Mean of 2-mon =", mean(SimSpeciesRichnessModelI[2,]), "\n")
cat("Mean of 6-mon =", mean(SimSpeciesRichnessModelI[3,]), "\n")
cat("Mean of 10-mon =", mean(SimSpeciesRichnessModelI[4,]), "\n")
cat("Mean of 20-mon =", mean(SimSpeciesRichnessModelI[5,]), "\n")



cat("Median of cm =", median(SimSpeciesRichnessModelI[1,]), "\n")
cat("Median of 2-mon =", median(SimSpeciesRichnessModelI[2,]), "\n")
cat("Median of 6-mon =", median(SimSpeciesRichnessModelI[3,]), "\n")
cat("Median of 10-mon =", median(SimSpeciesRichnessModelI[4,]), "\n")
cat("Median of 20-mon =", median(SimSpeciesRichnessModelI[5,]), "\n")

cat("MSE of cm=", (mean((SimSpeciesRichnessModelI[1,]-1000)^2))^{1/2}, "\n")
cat("MSE of 2-mon=", (mean((SimSpeciesRichnessModelI[2,]-1000)^2))^{1/2}, "\n")
cat("MSE of 6-mon=", (mean((SimSpeciesRichnessModelI[3,]-1000)^2))^{1/2}, "\n")
cat("MSE of 10-mon=", (mean((SimSpeciesRichnessModelI[4,]-1000)^2))^{1/2}, "\n")
cat("MSE of 20-mon=", (mean((SimSpeciesRichnessModelI[5,]-1000)^2))^{1/2}, "\n")


cat("MAE of cm=", mean(abs(SimSpeciesRichnessModelI[1,]-1000)), "\n")
cat("MAE of 2-mon=", mean(abs(SimSpeciesRichnessModelI[2,]-1000)), "\n")
cat("MAE of 6-mon=", mean(abs(SimSpeciesRichnessModelI[3,]-1000)), "\n")
cat("MAE of 10-mon=", mean(abs(SimSpeciesRichnessModelI[4,]-1000)), "\n")
cat("MAE of 20-mon=", mean(abs(SimSpeciesRichnessModelI[5,]-1000)), "\n")

cat("sd of cm=", sd(SimSpeciesRichnessModelI[1,]), "\n")
cat("sd of 2-mon=", sd(SimSpeciesRichnessModelI[2,]), "\n")
cat("sd of 6-mon=", sd(SimSpeciesRichnessModelI[3,]), "\n")
cat("sd of 10-mon=", sd(SimSpeciesRichnessModelI[4,]), "\n")
cat("sd of 20-mon=", sd(SimSpeciesRichnessModelI[5,]), "\n")

cat("conf int of cm=", quantile(SimSpeciesRichnessModelI[1,], prob=c(0.025,0.975)), "\n")
cat("conf int of 2-mon=", quantile(SimSpeciesRichnessModelI[2,], prob=c(0.025,0.975)), "\n")
cat("conf int of 6-mon=", quantile(SimSpeciesRichnessModelI[3,], prob=c(0.025,0.975)), "\n")
cat("conf int of 10-mon=", quantile(SimSpeciesRichnessModelI[4,], prob=c(0.025,0.975)), "\n")
cat("conf int of 20-mon=", quantile(SimSpeciesRichnessModelI[5,], prob=c(0.025,0.975)), "\n")




SimSpeciesRichnessModelII=SimSpeciesRichness(N=1000, M=200, Model ="II")

cat("Mean of cm =", mean(SimSpeciesRichnessModelII[1,]), "\n")
cat("Mean of 2-mon =", mean(SimSpeciesRichnessModelII[2,]), "\n")
cat("Mean of 6-mon =", mean(SimSpeciesRichnessModelII[3,]), "\n")
cat("Mean of 10-mon =", mean(SimSpeciesRichnessModelII[4,]), "\n")
cat("Mean of 20-mon =", mean(SimSpeciesRichnessModelII[5,]), "\n")



cat("Median of cm =", median(SimSpeciesRichnessModelII[1,]), "\n")
cat("Median of 2-mon =", median(SimSpeciesRichnessModelII[2,]), "\n")
cat("Median of 6-mon =", median(SimSpeciesRichnessModelII[3,]), "\n")
cat("Median of 10-mon =", median(SimSpeciesRichnessModelII[4,]), "\n")
cat("Median of 20-mon =", median(SimSpeciesRichnessModelII[5,]), "\n")

cat("MSE of cm=", (mean((SimSpeciesRichnessModelII[1,]-1000)^2))^{1/2}, "\n")
cat("MSE of 2-mon=", (mean((SimSpeciesRichnessModelII[2,]-1000)^2))^{1/2}, "\n")
cat("MSE of 6-mon=", (mean((SimSpeciesRichnessModelII[3,]-1000)^2))^{1/2}, "\n")
cat("MSE of 10-mon=", (mean((SimSpeciesRichnessModelII[4,]-1000)^2))^{1/2}, "\n")
cat("MSE of 20-mon=", (mean((SimSpeciesRichnessModelII[5,]-1000)^2))^{1/2}, "\n")


cat("MAE of cm=", mean(abs(SimSpeciesRichnessModelII[1,]-1000)), "\n")
cat("MAE of 2-mon=", mean(abs(SimSpeciesRichnessModelII[2,]-1000)), "\n")
cat("MAE of 6-mon=", mean(abs(SimSpeciesRichnessModelII[3,]-1000)), "\n")
cat("MAE of 10-mon=", mean(abs(SimSpeciesRichnessModelII[4,]-1000)), "\n")
cat("MAE of 20-mon=", mean(abs(SimSpeciesRichnessModelII[5,]-1000)), "\n")

cat("sd of cm=", sd(SimSpeciesRichnessModelII[1,]), "\n")
cat("sd of 2-mon=", sd(SimSpeciesRichnessModelII[2,]), "\n")
cat("sd of 6-mon=", sd(SimSpeciesRichnessModelII[3,]), "\n")
cat("sd of 10-mon=", sd(SimSpeciesRichnessModelII[4,]), "\n")
cat("sd of 20-mon=", sd(SimSpeciesRichnessModelII[5,]), "\n")

cat("conf int of cm=", quantile(SimSpeciesRichnessModelII[1,], prob=c(0.025,0.975)), "\n")
cat("conf int of 2-mon=", quantile(SimSpeciesRichnessModelII[2,], prob=c(0.025,0.975)), "\n")
cat("conf int of 6-mon=", quantile(SimSpeciesRichnessModelII[3,], prob=c(0.025,0.975)), "\n")
cat("conf int of 10-mon=", quantile(SimSpeciesRichnessModelII[4,], prob=c(0.025,0.975)), "\n")
cat("conf int of 20-mon=", quantile(SimSpeciesRichnessModelII[5,], prob=c(0.025,0.975)), "\n")




SimSpeciesRichnessModelIII=SimSpeciesRichness(N=1000, M=200, Model ="III")

cat("Mean of cm =", mean(SimSpeciesRichnessModelIII[1,]), "\n")
cat("Mean of 2-mon =", mean(SimSpeciesRichnessModelIII[2,]), "\n")
cat("Mean of 6-mon =", mean(SimSpeciesRichnessModelIII[3,]), "\n")
cat("Mean of 10-mon =", mean(SimSpeciesRichnessModelIII[4,]), "\n")
cat("Mean of 20-mon =", mean(SimSpeciesRichnessModelIII[5,]), "\n")



cat("Median of cm =", median(SimSpeciesRichnessModelIII[1,]), "\n")
cat("Median of 2-mon =", median(SimSpeciesRichnessModelIII[2,]), "\n")
cat("Median of 6-mon =", median(SimSpeciesRichnessModelIII[3,]), "\n")
cat("Median of 10-mon =", median(SimSpeciesRichnessModelIII[4,]), "\n")
cat("Median of 20-mon =", median(SimSpeciesRichnessModelIII[5,]), "\n")

cat("MSE of cm=", (mean((SimSpeciesRichnessModelIII[1,]-1000)^2))^{1/2}, "\n")
cat("MSE of 2-mon=", (mean((SimSpeciesRichnessModelIII[2,]-1000)^2))^{1/2}, "\n")
cat("MSE of 6-mon=", (mean((SimSpeciesRichnessModelIII[3,]-1000)^2))^{1/2}, "\n")
cat("MSE of 10-mon=", (mean((SimSpeciesRichnessModelIII[4,]-1000)^2))^{1/2}, "\n")
cat("MSE of 20-mon=", (mean((SimSpeciesRichnessModelIII[5,]-1000)^2))^{1/2}, "\n")


cat("MAE of cm=", mean(abs(SimSpeciesRichnessModelIII[1,]-1000)), "\n")
cat("MAE of 2-mon=", mean(abs(SimSpeciesRichnessModelIII[2,]-1000)), "\n")
cat("MAE of 6-mon=", mean(abs(SimSpeciesRichnessModelIII[3,]-1000)), "\n")
cat("MAE of 10-mon=", mean(abs(SimSpeciesRichnessModelIII[4,]-1000)), "\n")
cat("MAE of 20-mon=", mean(abs(SimSpeciesRichnessModelIII[5,]-1000)), "\n")

cat("sd of cm=", sd(SimSpeciesRichnessModelIII[1,]), "\n")
cat("sd of 2-mon=", sd(SimSpeciesRichnessModelIII[2,]), "\n")
cat("sd of 6-mon=", sd(SimSpeciesRichnessModelIII[3,]), "\n")
cat("sd of 10-mon=", sd(SimSpeciesRichnessModelIII[4,]), "\n")
cat("sd of 20-mon=", sd(SimSpeciesRichnessModelIII[5,]), "\n")

cat("conf int of cm=", quantile(SimSpeciesRichnessModelIII[1,], prob=c(0.025,0.975)), "\n")
cat("conf int of 2-mon=", quantile(SimSpeciesRichnessModelIII[2,], prob=c(0.025,0.975)), "\n")
cat("conf int of 6-mon=", quantile(SimSpeciesRichnessModelIII[3,], prob=c(0.025,0.975)), "\n")
cat("conf int of 10-mon=", quantile(SimSpeciesRichnessModelIII[4,], prob=c(0.025,0.975)), "\n")
cat("conf int of 20-mon=", quantile(SimSpeciesRichnessModelIII[5,], prob=c(0.025,0.975)), "\n")



SimSpeciesRichnessModelIV=SimSpeciesRichness(N=1000, M=200, Model ="IV")


cat("Mean of cm =", mean(SimSpeciesRichnessModelIV[1,]), "\n")
cat("Mean of 2-mon =", mean(SimSpeciesRichnessModelIV[2,]), "\n")
cat("Mean of 6-mon =", mean(SimSpeciesRichnessModelIV[3,]), "\n")
cat("Mean of 10-mon =", mean(SimSpeciesRichnessModelIV[4,]), "\n")
cat("Mean of 20-mon =", mean(SimSpeciesRichnessModelIV[5,]), "\n")



cat("Median of cm =", median(SimSpeciesRichnessModelIV[1,]), "\n")
cat("Median of 2-mon =", median(SimSpeciesRichnessModelIV[2,]), "\n")
cat("Median of 6-mon =", median(SimSpeciesRichnessModelIV[3,]), "\n")
cat("Median of 10-mon =", median(SimSpeciesRichnessModelIV[4,]), "\n")
cat("Median of 20-mon =", median(SimSpeciesRichnessModelIV[5,]), "\n")

cat("MSE of cm=", (mean((SimSpeciesRichnessModelIV[1,]-1000)^2))^{1/2}, "\n")
cat("MSE of 2-mon=", (mean((SimSpeciesRichnessModelIV[2,]-1000)^2))^{1/2}, "\n")
cat("MSE of 6-mon=", (mean((SimSpeciesRichnessModelIV[3,]-1000)^2))^{1/2}, "\n")
cat("MSE of 10-mon=", (mean((SimSpeciesRichnessModelIV[4,]-1000)^2))^{1/2}, "\n")
cat("MSE of 20-mon=", (mean((SimSpeciesRichnessModelIV[5,]-1000)^2))^{1/2}, "\n")


cat("MAE of cm=", mean(abs(SimSpeciesRichnessModelIV[1,]-1000)), "\n")
cat("MAE of 2-mon=", mean(abs(SimSpeciesRichnessModelIV[2,]-1000)), "\n")
cat("MAE of 6-mon=", mean(abs(SimSpeciesRichnessModelIV[3,]-1000)), "\n")
cat("MAE of 10-mon=", mean(abs(SimSpeciesRichnessModelIV[4,]-1000)), "\n")
cat("MAE of 20-mon=", mean(abs(SimSpeciesRichnessModelIV[5,]-1000)), "\n")

cat("sd of cm=", sd(SimSpeciesRichnessModelIV[1,]), "\n")
cat("sd of 2-mon=", sd(SimSpeciesRichnessModelIV[2,]), "\n")
cat("sd of 6-mon=", sd(SimSpeciesRichnessModelIV[3,]), "\n")
cat("sd of 10-mon=", sd(SimSpeciesRichnessModelIV[4,]), "\n")
cat("sd of 20-mon=", sd(SimSpeciesRichnessModelIV[5,]), "\n")

cat("conf int of cm=", quantile(SimSpeciesRichnessModelIV[1,], prob=c(0.025,0.975)), "\n")
cat("conf int of 2-mon=", quantile(SimSpeciesRichnessModelIV[2,], prob=c(0.025,0.975)), "\n")
cat("conf int of 6-mon=", quantile(SimSpeciesRichnessModelIV[3,], prob=c(0.025,0.975)), "\n")
cat("conf int of 10-mon=", quantile(SimSpeciesRichnessModelIV[4,], prob=c(0.025,0.975)), "\n")
cat("conf int of 20-mon=", quantile(SimSpeciesRichnessModelIV[5,], prob=c(0.025,0.975)), "\n")



SimSpeciesRichnessModelV=SimSpeciesRichness(N=1000, M=200, Model ="V")


cat("Mean of cm =", mean(SimSpeciesRichnessModelV[1,]), "\n")
cat("Mean of 2-mon =", mean(SimSpeciesRichnessModelV[2,]), "\n")
cat("Mean of 6-mon =", mean(SimSpeciesRichnessModelV[3,]), "\n")
cat("Mean of 10-mon =", mean(SimSpeciesRichnessModelV[4,]), "\n")
cat("Mean of 20-mon =", mean(SimSpeciesRichnessModelV[5,]), "\n")



cat("Median of cm =", median(SimSpeciesRichnessModelV[1,]), "\n")
cat("Median of 2-mon =", median(SimSpeciesRichnessModelV[2,]), "\n")
cat("Median of 6-mon =", median(SimSpeciesRichnessModelV[3,]), "\n")
cat("Median of 10-mon =", median(SimSpeciesRichnessModelV[4,]), "\n")
cat("Median of 20-mon =", median(SimSpeciesRichnessModelV[5,]), "\n")

cat("MSE of cm=", (mean((SimSpeciesRichnessModelV[1,]-1000)^2))^{1/2}, "\n")
cat("MSE of 2-mon=", (mean((SimSpeciesRichnessModelV[2,]-1000)^2))^{1/2}, "\n")
cat("MSE of 6-mon=", (mean((SimSpeciesRichnessModelV[3,]-1000)^2))^{1/2}, "\n")
cat("MSE of 10-mon=", (mean((SimSpeciesRichnessModelV[4,]-1000)^2))^{1/2}, "\n")
cat("MSE of 20-mon=", (mean((SimSpeciesRichnessModelV[5,]-1000)^2))^{1/2}, "\n")


cat("MAE of cm=", mean(abs(SimSpeciesRichnessModelV[1,]-1000)), "\n")
cat("MAE of 2-mon=", mean(abs(SimSpeciesRichnessModelV[2,]-1000)), "\n")
cat("MAE of 6-mon=", mean(abs(SimSpeciesRichnessModelV[3,]-1000)), "\n")
cat("MAE of 10-mon=", mean(abs(SimSpeciesRichnessModelV[4,]-1000)), "\n")
cat("MAE of 20-mon=", mean(abs(SimSpeciesRichnessModelV[5,]-1000)), "\n")

cat("sd of cm=", sd(SimSpeciesRichnessModelV[1,]), "\n")
cat("sd of 2-mon=", sd(SimSpeciesRichnessModelV[2,]), "\n")
cat("sd of 6-mon=", sd(SimSpeciesRichnessModelV[3,]), "\n")
cat("sd of 10-mon=", sd(SimSpeciesRichnessModelV[4,]), "\n")
cat("sd of 20-mon=", sd(SimSpeciesRichnessModelV[5,]), "\n")

cat("conf int of cm=", quantile(SimSpeciesRichnessModelV[1,], prob=c(0.025,0.975)), "\n")
cat("conf int of 2-mon=", quantile(SimSpeciesRichnessModelV[2,], prob=c(0.025,0.975)), "\n")
cat("conf int of 6-mon=", quantile(SimSpeciesRichnessModelV[3,], prob=c(0.025,0.975)), "\n")
cat("conf int of 10-mon=", quantile(SimSpeciesRichnessModelV[4,], prob=c(0.025,0.975)), "\n")
cat("conf int of 20-mon=", quantile(SimSpeciesRichnessModelV[5,], prob=c(0.025,0.975)), "\n")



SimSpeciesRichnessModelVI=SimSpeciesRichness(N=1000, M=200, Model ="VI")


cat("Mean of cm =", mean(SimSpeciesRichnessModelVI[1,]), "\n")
cat("Mean of 2-mon =", mean(SimSpeciesRichnessModelVI[2,]), "\n")
cat("Mean of 6-mon =", mean(SimSpeciesRichnessModelVI[3,]), "\n")
cat("Mean of 10-mon =", mean(SimSpeciesRichnessModelVI[4,]), "\n")
cat("Mean of 20-mon =", mean(SimSpeciesRichnessModelVI[5,]), "\n")



cat("Median of cm =", median(SimSpeciesRichnessModelVI[1,]), "\n")
cat("Median of 2-mon =", median(SimSpeciesRichnessModelVI[2,]), "\n")
cat("Median of 6-mon =", median(SimSpeciesRichnessModelVI[3,]), "\n")
cat("Median of 10-mon =", median(SimSpeciesRichnessModelVI[4,]), "\n")
cat("Median of 20-mon =", median(SimSpeciesRichnessModelVI[5,]), "\n")

cat("MSE of cm=", (mean((SimSpeciesRichnessModelVI[1,]-1000)^2))^{1/2}, "\n")
cat("MSE of 2-mon=", (mean((SimSpeciesRichnessModelVI[2,]-1000)^2))^{1/2}, "\n")
cat("MSE of 6-mon=", (mean((SimSpeciesRichnessModelVI[3,]-1000)^2))^{1/2}, "\n")
cat("MSE of 10-mon=", (mean((SimSpeciesRichnessModelVI[4,]-1000)^2))^{1/2}, "\n")
cat("MSE of 20-mon=", (mean((SimSpeciesRichnessModelVI[5,]-1000)^2))^{1/2}, "\n")


cat("MAE of cm=", mean(abs(SimSpeciesRichnessModelVI[1,]-1000)), "\n")
cat("MAE of 2-mon=", mean(abs(SimSpeciesRichnessModelVI[2,]-1000)), "\n")
cat("MAE of 6-mon=", mean(abs(SimSpeciesRichnessModelVI[3,]-1000)), "\n")
cat("MAE of 10-mon=", mean(abs(SimSpeciesRichnessModelVI[4,]-1000)), "\n")
cat("MAE of 20-mon=", mean(abs(SimSpeciesRichnessModelVI[5,]-1000)), "\n")

cat("sd of cm=", sd(SimSpeciesRichnessModelVI[1,]), "\n")
cat("sd of 2-mon=", sd(SimSpeciesRichnessModelVI[2,]), "\n")
cat("sd of 6-mon=", sd(SimSpeciesRichnessModelVI[3,]), "\n")
cat("sd of 10-mon=", sd(SimSpeciesRichnessModelVI[4,]), "\n")
cat("sd of 20-mon=", sd(SimSpeciesRichnessModelVI[5,]), "\n")

cat("conf int of cm=", quantile(SimSpeciesRichnessModelVI[1,], prob=c(0.025,0.975)), "\n")
cat("conf int of 2-mon=", quantile(SimSpeciesRichnessModelVI[2,], prob=c(0.025,0.975)), "\n")
cat("conf int of 6-mon=", quantile(SimSpeciesRichnessModelVI[3,], prob=c(0.025,0.975)), "\n")
cat("conf int of 10-mon=", quantile(SimSpeciesRichnessModelVI[4,], prob=c(0.025,0.975)), "\n")
cat("conf int of 20-mon=", quantile(SimSpeciesRichnessModelVI[5,], prob=c(0.025,0.975)), "\n")



SimSpeciesRichnessModelVII=SimSpeciesRichness(N=1000, M=200, Model ="VII")


cat("Mean of cm =", mean(SimSpeciesRichnessModelVII[1,]), "\n")
cat("Mean of 2-mon =", mean(SimSpeciesRichnessModelVII[2,]), "\n")
cat("Mean of 6-mon =", mean(SimSpeciesRichnessModelVII[3,]), "\n")
cat("Mean of 10-mon =", mean(SimSpeciesRichnessModelVII[4,]), "\n")
cat("Mean of 20-mon =", mean(SimSpeciesRichnessModelVII[5,]), "\n")



cat("Median of cm =", median(SimSpeciesRichnessModelVII[1,]), "\n")
cat("Median of 2-mon =", median(SimSpeciesRichnessModelVII[2,]), "\n")
cat("Median of 6-mon =", median(SimSpeciesRichnessModelVII[3,]), "\n")
cat("Median of 10-mon =", median(SimSpeciesRichnessModelVII[4,]), "\n")
cat("Median of 20-mon =", median(SimSpeciesRichnessModelVII[5,]), "\n")

cat("MSE of cm=", (mean((SimSpeciesRichnessModelVII[1,]-1000)^2))^{1/2}, "\n")
cat("MSE of 2-mon=", (mean((SimSpeciesRichnessModelVII[2,]-1000)^2))^{1/2}, "\n")
cat("MSE of 6-mon=", (mean((SimSpeciesRichnessModelVII[3,]-1000)^2))^{1/2}, "\n")
cat("MSE of 10-mon=", (mean((SimSpeciesRichnessModelVII[4,]-1000)^2))^{1/2}, "\n")
cat("MSE of 20-mon=", (mean((SimSpeciesRichnessModelVII[5,]-1000)^2))^{1/2}, "\n")


cat("MAE of cm=", mean(abs(SimSpeciesRichnessModelVII[1,]-1000)), "\n")
cat("MAE of 2-mon=", mean(abs(SimSpeciesRichnessModelVII[2,]-1000)), "\n")
cat("MAE of 6-mon=", mean(abs(SimSpeciesRichnessModelVII[3,]-1000)), "\n")
cat("MAE of 10-mon=", mean(abs(SimSpeciesRichnessModelVII[4,]-1000)), "\n")
cat("MAE of 20-mon=", mean(abs(SimSpeciesRichnessModelVII[5,]-1000)), "\n")

cat("sd of cm=", sd(SimSpeciesRichnessModelVII[1,]), "\n")
cat("sd of 2-mon=", sd(SimSpeciesRichnessModelVII[2,]), "\n")
cat("sd of 6-mon=", sd(SimSpeciesRichnessModelVII[3,]), "\n")
cat("sd of 10-mon=", sd(SimSpeciesRichnessModelVII[4,]), "\n")
cat("sd of 20-mon=", sd(SimSpeciesRichnessModelVII[5,]), "\n")

cat("conf int of cm=", quantile(SimSpeciesRichnessModelVII[1,], prob=c(0.025,0.975)), "\n")
cat("conf int of 2-mon=", quantile(SimSpeciesRichnessModelVII[2,], prob=c(0.025,0.975)), "\n")
cat("conf int of 6-mon=", quantile(SimSpeciesRichnessModelVII[3,], prob=c(0.025,0.975)), "\n")
cat("conf int of 10-mon=", quantile(SimSpeciesRichnessModelVII[4,], prob=c(0.025,0.975)), "\n")
cat("conf int of 20-mon=", quantile(SimSpeciesRichnessModelVII[5,], prob=c(0.025,0.975)), "\n")



SimSpeciesRichnessModelVIII=SimSpeciesRichness(N=1000, M=200, Model ="VIII")


cat("Mean of cm =", mean(SimSpeciesRichnessModelVIII[1,]), "\n")
cat("Mean of 2-mon =", mean(SimSpeciesRichnessModelVIII[2,]), "\n")
cat("Mean of 6-mon =", mean(SimSpeciesRichnessModelVIII[3,]), "\n")
cat("Mean of 10-mon =", mean(SimSpeciesRichnessModelVIII[4,]), "\n")
cat("Mean of 20-mon =", mean(SimSpeciesRichnessModelVIII[5,]), "\n")



cat("Median of cm =", median(SimSpeciesRichnessModelVIII[1,]), "\n")
cat("Median of 2-mon =", median(SimSpeciesRichnessModelVIII[2,]), "\n")
cat("Median of 6-mon =", median(SimSpeciesRichnessModelVIII[3,]), "\n")
cat("Median of 10-mon =", median(SimSpeciesRichnessModelVIII[4,]), "\n")
cat("Median of 20-mon =", median(SimSpeciesRichnessModelVIII[5,]), "\n")

cat("MSE of cm=", (mean((SimSpeciesRichnessModelVIII[1,]-1000)^2))^{1/2}, "\n")
cat("MSE of 2-mon=", (mean((SimSpeciesRichnessModelVIII[2,]-1000)^2))^{1/2}, "\n")
cat("MSE of 6-mon=", (mean((SimSpeciesRichnessModelVIII[3,]-1000)^2))^{1/2}, "\n")
cat("MSE of 10-mon=", (mean((SimSpeciesRichnessModelVIII[4,]-1000)^2))^{1/2}, "\n")
cat("MSE of 20-mon=", (mean((SimSpeciesRichnessModelVIII[5,]-1000)^2))^{1/2}, "\n")


cat("MAE of cm=", mean(abs(SimSpeciesRichnessModelVIII[1,]-1000)), "\n")
cat("MAE of 2-mon=", mean(abs(SimSpeciesRichnessModelVIII[2,]-1000)), "\n")
cat("MAE of 6-mon=", mean(abs(SimSpeciesRichnessModelVIII[3,]-1000)), "\n")
cat("MAE of 10-mon=", mean(abs(SimSpeciesRichnessModelVIII[4,]-1000)), "\n")
cat("MAE of 20-mon=", mean(abs(SimSpeciesRichnessModelVIII[5,]-1000)), "\n")
I
cat("sd of cm=", sd(SimSpeciesRichnessModelVII[1,]), "\n")
cat("sd of 2-mon=", sd(SimSpeciesRichnessModelVIII[2,]), "\n")
cat("sd of 6-mon=", sd(SimSpeciesRichnessModelVIII[3,]), "\n")
cat("sd of 10-mon=", sd(SimSpeciesRichnessModelVIII[4,]), "\n")
cat("sd of 20-mon=", sd(SimSpeciesRichnessModelVIII[5,]), "\n")

cat("conf int of cm=", quantile(SimSpeciesRichnessModelVIII[1,], prob=c(0.025,0.975)), "\n")
cat("conf int of 2-mon=", quantile(SimSpeciesRichnessModelVIII[2,], prob=c(0.025,0.975)), "\n")
cat("conf int of 6-mon=", quantile(SimSpeciesRichnessModelVIII[3,], prob=c(0.025,0.975)), "\n")
cat("conf int of 10-mon=", quantile(SimSpeciesRichnessModelVIII[4,], prob=c(0.025,0.975)), "\n")
cat("conf int of 20-mon=", quantile(SimSpeciesRichnessModelVIII[5,], prob=c(0.025,0.975)), "\n")