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) {
16 cout << "abund = " << abund << "\n";
18 // vector<double> aceData(3,0);
19 double ace, acelci, acehci;
25 double Cace, term1, gamace;
28 double maxRank = (double)rank->getMaxRank();
30 for(int i=1;i<=maxRank;i++){
32 srare += rank->get(i);
33 nrare += i*rank->get(i);
34 numsum += (i-1)*i*rank->get(i);
36 else if(i>abund) {sabund += rank->get(i);}
38 int sobs = srare + sabund;
40 if (nrare == 0){ Cace = 0.0000; }
41 else { Cace = 1.0000 -(double)rank->get(1)/(double)nrare; }
43 double denom = Cace * (double)(nrare * (nrare-1));
45 if(denom <= 0.0){ term1=0.0000; } else { term1 = (double)(srare * numsum)/(double)denom - 1.0; }
46 if(term1 >= 0.0){ gamace = term1; } else { gamace = 0.0; }
49 gamace = gamace * (1 + (nrare * (1 - Cace) * numsum) / denom);
50 if(gamace<0){ gamace = 0; }
56 ace = (double)sabund+((double)srare+(double)rank->get(1)*gamace)/Cace;//ace
60 The following code was obtained from Anne Chao for calculating the SE for her ACE estimator
61 My modification was to reset the frequencies so that a singleton is found in rank[1] insted
64 I have also added the forumlae to calculate the 95% confidence intervals.
67 int j,D_s=0,nn=0,ww=0,Max_Index=rank->getMaxRank()+1;
68 double pp, temp1, temp2;
69 vector<double> Part_N_Part_F(Max_Index+1,0.0);
71 for (j=1; j<Max_Index; j++) if(j<=abund) D_s += rank->get(j);
72 for (j=1; j<Max_Index; j++){
74 nn += rank->get(j) * j;
75 ww += rank->get(j) * j * ( j - 1);
78 double C_hat = 1.-rank->get(1)/double(nn);
79 double Gamma = ( D_s * ww) / ( C_hat * nn * ( nn - 1.)) - 1.;
80 temp1 = double(nn - rank->get(1));
81 temp2 = double(nn - 1.);
84 Part_N_Part_F[1] = ( D_s + nn) * ( 1. + rank->get(1) * ww / temp1 / temp2) / temp1 + nn * D_s * ww * ( temp1 - 1.) /
85 ( temp1 * temp1 * temp2 * temp2) - ( nn + rank->get(1)) / temp1;
86 for ( j=2; j<=Max_Index; j++){
88 Part_N_Part_F[j] = ( nn * temp1 - j * rank->get(1) * D_s) / temp1 / temp1 * ( 1. + rank->get(1) * ww / temp1 / temp2)
89 + j * rank->get(1) * D_s * nn * ( ( j - 1.) * temp1 * temp2 - ww * ( temp1 + temp2))
90 / temp1 / temp1 / temp1 / temp2 / temp2 + j * rank->get(1) * rank->get(1) / temp1 / temp1;
95 Part_N_Part_F[1] = ( nn + D_s ) / temp1;
96 for ( j=2; j<=Max_Index; j++){
98 Part_N_Part_F[j-1] = ( nn * temp1 - j * rank->get(1) * D_s ) / temp1 / temp1;
103 for ( j=abund+1; j<=Max_Index; j++){
104 Part_N_Part_F[j-1] = 1.;
107 for ( temp1=0., temp2=0., j=0; j<Max_Index; j++) {
108 pp = Part_N_Part_F[j];
109 temp1 += pp * rank->get(j);
110 temp2 += pp * pp * rank->get(j);
113 double se = temp2 - temp1 * temp1 / ace;
115 if(toString(se) == "nan"){
124 double ci = 1.96*pow(se,0.5);
125 acelci = ace-ci; //ace lci
126 acehci = ace+ci; //ace hci
128 double denom = pow(ace-sobs,2);
129 double c = exp(1.96*pow((log(1+se/denom)),0.5));
130 acelci = sobs+(ace-sobs)/c; //ace lci
131 acehci = sobs+(ace-sobs)*c; //ace hci
138 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
139 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
140 if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
144 catch(exception& e) {
145 cout << "Standard Error: " << e.what() << " has occurred in the Ace class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
149 cout << "An unknown error has occurred in the Ace class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
154 /***********************************************************************/