+++ /dev/null
-/*
- * ace.cpp
- * Dotur
- *
- * Created by Sarah Westcott on 1/7/09.
- * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
- *
- */
-
-#include "ace.h"
-
-/***********************************************************************/
-
-EstOutput Ace::getValues(SAbundVector* rank) {
- try {
- data.resize(3,0);
- double ace, acelci, acehci;
-
- double nrare = 0;
- double srare = 0;
- double sabund = 0;
-
- double Cace, term1, gamace;
- double numsum = 0;
-
- double maxRank = (double)rank->getMaxRank();
-
- for(int i=1;i<=maxRank;i++){
- if(i<=abund){
- srare += rank->get(i);
- nrare += i*rank->get(i);
- numsum += (i-1)*i*rank->get(i);
- }
- else if(i>abund) {sabund += rank->get(i);}
- }
- double sobs = srare + sabund;
-
- if (nrare == 0){ Cace = 0.0000; }
- else { Cace = 1.0000 -(double)rank->get(1)/(double)nrare; }
-
- double denom = Cace * (double)(nrare * (nrare-1));
-
- if(denom <= 0.0){ term1=0.0000; } else { term1 = (double)(srare * numsum)/(double)denom - 1.0; }
- if(term1 >= 0.0){ gamace = term1; } else { gamace = 0.0; }
-
- if(gamace >= 0.64){
- gamace = gamace * (1 + (nrare * (1 - Cace) * numsum) / denom);
- if(gamace<0){ gamace = 0; }
- }
-
- if(Cace == 0.0){
- ace = 0.00;}//ace
- else{
- ace = (double)sabund+((double)srare+(double)rank->get(1)*gamace)/Cace;//ace
- }
-
- /*
- The following code was obtained from Anne Chao for calculating the SE for her ACE estimator
- My modification was to reset the frequencies so that a singleton is found in rank[1] insted
- of in rank[0], etc.
-
- I have also added the forumlae to calculate the 95% confidence intervals.
- */
-
- double j,D_s=0,nn=0,ww=0;
- int Max_Index=rank->getMaxRank()+1;
- double pp, temp1, temp2;
- vector<double> Part_N_Part_F(Max_Index+1,0.0);
-
- for (j=1; j<Max_Index; j++) if(j<=abund) D_s += rank->get(j);
- for (j=1; j<Max_Index; j++){
- if(j<=abund){
- nn += rank->get(j) * j;
- ww += rank->get(j) * j * ( j - 1);
- }
- }
- double C_hat = 1.-rank->get(1)/double(nn);
- double Gamma = ( D_s * ww) / ( C_hat * nn * ( nn - 1.)) - 1.;
- temp1 = double(nn - rank->get(1));
- temp2 = double(nn - 1.);
-
- if ( Gamma > 0.){
- Part_N_Part_F[1] = ( D_s + nn) * ( 1. + rank->get(1) * ww / temp1 / temp2) / temp1 + nn * D_s * ww * ( temp1 - 1.) /
- ( temp1 * temp1 * temp2 * temp2) - ( nn + rank->get(1)) / temp1;
- for ( j=2; j<=Max_Index; j++){
- if(j<=abund){
- Part_N_Part_F[j] = ( nn * temp1 - j * rank->get(1) * D_s) / temp1 / temp1 * ( 1. + rank->get(1) * ww / temp1 / temp2)
- + j * rank->get(1) * D_s * nn * ( ( j - 1.) * temp1 * temp2 - ww * ( temp1 + temp2))
- / temp1 / temp1 / temp1 / temp2 / temp2 + j * rank->get(1) * rank->get(1) / temp1 / temp1;
- }
- }
- }
- else{
- Part_N_Part_F[1] = ( nn + D_s ) / temp1;
- for ( j=2; j<=Max_Index; j++){
- if(j<=abund){
- Part_N_Part_F[j-1] = ( nn * temp1 - j * rank->get(1) * D_s ) / temp1 / temp1;
- }
- }
- }
- if(Max_Index>abund){
- for ( j=abund+1; j<=Max_Index; j++){
- Part_N_Part_F[j-1] = 1.;
- }
- }
- for ( temp1=0., temp2=0., j=0; j<Max_Index; j++) {
- pp = Part_N_Part_F[j];
- temp1 += pp * rank->get(j);
- temp2 += pp * pp * rank->get(j);
- }
-
- double se = temp2 - temp1 * temp1 / ace;
-
- if(toString(se) == "nan"){
- acelci = ace;
- acehci = ace;
- }
- else if(ace==0.000){
- acelci = ace;
- acehci = ace;
- }
- else if(ace==sobs){
- double ci = 1.96*pow(se,0.5);
- acelci = ace-ci; //ace lci
- acehci = ace+ci; //ace hci
- }else{
- double denom = pow(ace-sobs,2);
- double c = exp(1.96*pow((log(1+se/denom)),0.5));
- acelci = sobs+(ace-sobs)/c; //ace lci
- acehci = sobs+(ace-sobs)*c; //ace hci
- }
-
- data[0] = ace;
- data[1] = acelci;
- data[2] = acehci;
-
- if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
- if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
- if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
-
- return data;
- }
- catch(exception& e) {
- m->errorOut(e, "Ace", "getValues");
- exit(1);
- }
-}
-
-/***********************************************************************/
projects / mothur.git / blobdiff