5 * Created by Sarah Westcott on 1/7/09.
6 * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
12 /***********************************************************************/
14 EstOutput Ace::getValues(SAbundVector* rank) {
17 double ace, acelci, acehci;
23 double Cace, term1, gamace;
26 double maxRank = (double)rank->getMaxRank();
28 for(int i=1;i<=maxRank;i++){
30 srare += rank->get(i);
31 nrare += i*rank->get(i);
32 numsum += (i-1)*i*rank->get(i);
34 else if(i>abund) {sabund += rank->get(i);}
36 double sobs = srare + sabund;
38 if (nrare == 0){ Cace = 0.0000; }
39 else { Cace = 1.0000 -(double)rank->get(1)/(double)nrare; }
41 double denom = Cace * (double)(nrare * (nrare-1));
43 if(denom <= 0.0){ term1=0.0000; } else { term1 = (double)(srare * numsum)/(double)denom - 1.0; }
44 if(term1 >= 0.0){ gamace = term1; } else { gamace = 0.0; }
47 gamace = gamace * (1 + (nrare * (1 - Cace) * numsum) / denom);
48 if(gamace<0){ gamace = 0; }
54 ace = (double)sabund+((double)srare+(double)rank->get(1)*gamace)/Cace;//ace
58 The following code was obtained from Anne Chao for calculating the SE for her ACE estimator
59 My modification was to reset the frequencies so that a singleton is found in rank[1] insted
62 I have also added the forumlae to calculate the 95% confidence intervals.
65 double j,D_s=0,nn=0,ww=0;
66 int Max_Index=rank->getMaxRank()+1;
67 double pp, temp1, temp2;
68 vector<double> Part_N_Part_F(Max_Index+1,0.0);
70 for (j=1; j<Max_Index; j++) if(j<=abund) D_s += rank->get(j);
71 for (j=1; j<Max_Index; j++){
73 nn += rank->get(j) * j;
74 ww += rank->get(j) * j * ( j - 1);
77 double C_hat = 1.-rank->get(1)/double(nn);
78 double Gamma = ( D_s * ww) / ( C_hat * nn * ( nn - 1.)) - 1.;
79 temp1 = double(nn - rank->get(1));
80 temp2 = double(nn - 1.);
83 Part_N_Part_F[1] = ( D_s + nn) * ( 1. + rank->get(1) * ww / temp1 / temp2) / temp1 + nn * D_s * ww * ( temp1 - 1.) /
84 ( temp1 * temp1 * temp2 * temp2) - ( nn + rank->get(1)) / temp1;
85 for ( j=2; j<=Max_Index; j++){
87 Part_N_Part_F[j] = ( nn * temp1 - j * rank->get(1) * D_s) / temp1 / temp1 * ( 1. + rank->get(1) * ww / temp1 / temp2)
88 + j * rank->get(1) * D_s * nn * ( ( j - 1.) * temp1 * temp2 - ww * ( temp1 + temp2))
89 / temp1 / temp1 / temp1 / temp2 / temp2 + j * rank->get(1) * rank->get(1) / temp1 / temp1;
94 Part_N_Part_F[1] = ( nn + D_s ) / temp1;
95 for ( j=2; j<=Max_Index; j++){
97 Part_N_Part_F[j-1] = ( nn * temp1 - j * rank->get(1) * D_s ) / temp1 / temp1;
102 for ( j=abund+1; j<=Max_Index; j++){
103 Part_N_Part_F[j-1] = 1.;
106 for ( temp1=0., temp2=0., j=0; j<Max_Index; j++) {
107 pp = Part_N_Part_F[j];
108 temp1 += pp * rank->get(j);
109 temp2 += pp * pp * rank->get(j);
112 double se = temp2 - temp1 * temp1 / ace;
114 if(toString(se) == "nan"){
123 double ci = 1.96*pow(se,0.5);
124 acelci = ace-ci; //ace lci
125 acehci = ace+ci; //ace hci
127 double denom = pow(ace-sobs,2);
128 double c = exp(1.96*pow((log(1+se/denom)),0.5));
129 acelci = sobs+(ace-sobs)/c; //ace lci
130 acehci = sobs+(ace-sobs)*c; //ace hci
137 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
138 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
139 if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
143 catch(exception& e) {
144 m->errorOut(e, "Ace", "getValues");
149 /***********************************************************************/