####################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")