C++ 计算标准差

Question

我正在尝试编写一个简单的数据分析程序，但我在计算十元素数组的标准差时遇到了问题

我们应该使用平方近似 (stddev = sqrt(meanofsquares/90-squareofmeans900)) 而不是标准方法，但我得到的数字几乎是随机的。

为了计算总和，我正在使用以下函数

void rms(int *mc, float &Sum, float &Sum2){  //restituisce
    Sum = 0.; Sum2 = 0.;
    for(int i = 0; i<10;i++){
    Sum += (float) mc[i];
    Sum2+=(float)mc[i]*mc[i];}
}

如下调用

   flux[idep]=((float) mc[4]+ mc[5])*.5; //anzichè dividere per 2, moltiplico per .5

    rms(mc, Sum, Sum2);

    errflux[idep] = 1.177 * sqrt(Sum2*div1 - (Sum1*Sum1*div2)); //notice corrective term 1.177 it's not an error

我尝试了很多方法，但我从来没有设法得到一些连贯的结果。

我附上几行数据(每行是一个 block ，应该计算中位数和标准差)和完整代码，希望对您有所帮助

Dumand.dat

 1023001  998765 1109975  865876 1407325 1498650 1065880  999623 1300568 1400421
  178433  438761  165001  234121  765999  999650  876500  190001  206443  180065
   88951  501003   13760   21880  197550  199902   46909   54331   76008   70913
   10569   15900   29761   20769   22331   20117   21555   26700   45600   37654
    4535    4289    4059    4099    4300    4401    4221    4390    4101    4203
    1402    1436    1420    1450    1398    1590    1360    1570    1531    1466
 1693243 1653128 1837200 1433174 2329366 2480524 1764215 1654548 2152664 2317938
  295337  726225  273105  387511 1267860 1654593 1450758  314484  341699  298035
  147229  829246   22775   36215  326979  330872   77644   89920  125800  117400
   17490   26317   49260   34376   36960   33297   35677   44193   75475   62323
    7500    7099    6718    6785    7118    7284    6987    7280    6790    6961
    2321    2370    2350    2398    2314    2630    2251    2600    2534    2420

完整代码是

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iomanip>
#include <string>
#include <algorithm>      //algoritmi di ordinamento
#include <unistd.h>    //Useremo lo standard posix per le api per correggere delle carenze in gestione i/o del c++
using namespace std;

bool myfunct(double i, double j){return j<i;}
//funzione per calcolarci i contributi all'errore

void rms(int *mc, float &Sum, float &Sum2){  //restituisce
    Sum = 0.; Sum2 = 0.;
    for(int i = 0; i<10;i++){
    Sum += (float) mc[i];
    Sum2+=(float)mc[i]*mc[i];}

}



//ocio: usando il namespace std le funzioni del namespace ios danno errore. Calma e sangue freddo, ogni volta che ne usiamo una accanto ci sta il metodo alternativo senza usarle
int main (){
//input data
int mc[10];  //le dieci misure inerenti ad una profondità
float C[2];  // le due costanti di calibrazione
float depth[4] = {1500., 2500., 3500., 4500.}; //profondità delle misure (servono per il best fit)
float Sum, Sum2, Sum1;
//output data
float flux[4], errflux[4]; //flussi e errori alle varie profondità
float fluxb, errfluxb;     //flusso e errore sul fondo dell'oceano
float bkg, errbkg;        //flusso e errore sul rumore (liv del mare)


//calcoliamo ora le divisioni in modo da trasformarle in prodotti. Su programmi più pesanti questo trucco accelera un botto il tempo di esecuzione
float div1= 1./90.;
float div2 = 1./900.;

//lettura delle costanti di calibrazione

cout << "type calibration constants C1, C2 ... " << endl;
cin >> C[0] >> C[1];
    float invc[2] = {1./C[0], 1./C[1]};
// external files
/* ====== */
ifstream in("DUMAND.dat");//, ios::nocreate);
if(!in){cerr << "The input file DUMAND.dat does not exist. Check path. " << endl; return -1;} //se il file non esiste, avverti e killa.
/* ====== */

//ALTERNATIVA PER L'INPUT USANDO IL POSIX
/*
if(access("DUMAND.dat", 0)==-1){ //regola i permessi di accesso a un file, la useremo per controllare l'esistenza di un file
    cerr << "The input file DUMAND.dat does not exist. Check path. " << endl; return -1;} //se il file non esiste, avverti e killa.}
//altrimenti SOLO A QUESTO PUNTO dichiaro l'ifstream ecc ecc
ifstream input("DUMAND.dat");
*/

//if(access("results.dat", 0)==0){cerr << "The file results.dat already exists. Aborting operation"; return -1;}
ofstream out ("results.dat");//, ios::noreplace);


//error strings
string erreof = "Unexpected EoF at rec # ";
string errmisc = "Unexpected error reading rec # ";


//loop over PMs
for (int ipm = 0; ipm < 2; ipm++){
    int nrec;
//loop over depths
    for (int idep = 0; idep < 4; idep++){

        nrec = 6*ipm + idep+1;
        //loop for reading from file
        for (int i=0; i < 10; i++){

                in >> mc[i];

                if (in.eof()){cerr << erreof << nrec << endl;
                              return -1;} 
                if (in.bad()){
                               cerr << errmisc << nrec << endl;
                               return -1;}

                } //close acquisition loop
//compute median. Il secondo argomento è l'idirizzo dell'ultimo elemento di mc +1, il terzo il criterio //(ascendente o discendente)

        sort(mc,mc+10, myfunct); 

        flux[idep]=((float) mc[4]+ mc[5])*.5; //anzichè dividere per 2, moltiplico per .5

        rms(mc, Sum, Sum2);
//        Sum2=0.;
//        Sum=.1*Sum;
//        for(int i =0; i<10;i++){
//          Sum2+=(mc[i]-Sum)*(mc[i]-Sum);
//        }
        errflux[idep] = 1.177 * sqrt(Sum2*div1 - (Sum1*Sum1*div2)); //notice corrective term 1.177
       cout << flux[idep] << " " << errflux[idep] << endl;
        system("PAUSE");
        } // close loop over depth
    //bottom level analysis
    ++nrec; 
            for(int i = 0; i<10; i++){
            in >> mc[i];
            if (in.eof()){cerr << erreof << nrec << endl;
                          return -1;} 
            if (in.bad()){
                           cerr << errmisc << nrec << endl;
                           return -1;}
            } //close acquisition loop
    //per fondo e lab calcoliamo il val medio e non la mediana perchè le dist sono a code piccole
    rms(mc, Sum, Sum2);
    errfluxb = sqrt(Sum2*div1 - Sum*Sum * div2); //std dev bottom level
    fluxb = Sum*.1; //average bottom level   

//lab level analysis (noise analysis)
    ++nrec; 
            for(int i = 0; i<10; i++){
            in >> mc[i];
            if (in.eof()){cerr << erreof << nrec << endl;
                          return -1;} 
            if (in.bad()){
                           cerr << errmisc << nrec << endl;
                           return -1;}
            } //close acquisition loop
    //per fondo e lab calcoliamo il val medio e non la mediana perchè le dist sono a code piccole
    rms(mc, Sum, Sum2);
    errbkg = sqrt(Sum2*div1 - Sum*Sum * div2); //std dev laboratory level
    bkg = Sum*.1; //average laboratory level   

    //compute true fluxes removing noise and compensating calibration
    for(int idep=0; idep <4; idep++ ){
        flux[idep]-=bkg;
        flux[idep]*=invc[ipm]; //rinormalizzazione con la cost di calibrazione
        errflux[idep]=(sqrt(errflux[idep]*errflux[idep] + errbkg*errbkg)) *invc[ipm];
    }    
    fluxb -= bkg;
    fluxb *= invc[ipm];
    errflux[idep]=(sqrt(errfluxb*errfluxb + errbkg*errbkg)) *invc[ipm];


    //best fit
    float S1 = 0, SX = 0, SY = 0, SXX = 0, SXY = 0;


    for (int idep =0; idep < 4; idep++){
        float y = log(flux[idep]);
        float erryinv = flux[idep]/errflux[idep]; //since we need 1/error, we compute it directly
        erryinv = erryinv * erryinv;
        S1+=erryinv;
        SX += depth[idep]*erryinv;
        SY += y * erryinv;
        SXX = depth[idep]*depth[idep] * erryinv;
        SXY = depth[idep]*y * erryinv;

    }//close depth loop
    float invD = 1./((S1*SXX)-(SX*SX)); //again, we define 1/D because we need it and not the original D   
    cout << SXX << " "<<S1 << " " << SX << "  " << invD << endl << endl;

    float a = (SY*SXX - SX*SXY)*invD;
    float b = (S1*SXY - SX*SY)*invD;
    float erra = sqrt((SXX*invD));
    float errb = sqrt((S1*invD));

    //true fit parameters
    a = exp(a);
    b = -1/b;
    erra = a * erra;
    errb = errb * b * b;

    //table of results
    string t24 = "\t \t \t " , title = "results for PM( ", tit1 = "Flux and related errors for each depth", tit2 = "flux and error at bottom level", tit3 = "fit parameters and errors";

    out << t24 << title << ipm+1 << " )" << endl << endl;
    out << '\t' << tit1 << endl;
    out.flags(ios::scientific|ios::uppercase); //per salvare i numeri in notazione scientifica con la E maiuscola
    for (int l =0; l < 4; l++){
    out<<setw(11) << setprecision(4) << flux[l] << " +- " << errflux[l] << endl; //imposto il numero di caratteri da usare nella notazione scientifica (setw) e la precisione(setprecision). oss: l'argomento di setw è pari a 7 + l'arg di setprecision
    }    
out << endl;
    out<< '\t' << tit2 << endl;
    out<<setw(11) << setprecision(4) << fluxb << " +- " << errfluxb << endl;   
out << endl;
    out<< '\t' << tit3 << endl;    
    out<<setw(11) << setprecision(4) << "a = " <<  a << " +- " << erra << endl;   
    out<<setw(11) << setprecision(4) << "b = " <<  b << " +- " << errb << endl;
    }// close pm loop   

in.close();
out.close();
return 1;
}

Answer 1

好吧，你的代码中没有计算方法

应该除以 n。

void rms(int *mc, int n, float &Sum, float &Sum2) {  //restituisce
    Sum = 0.;
    Sum2 = 0.;
    for(int i = 0; i != n; i++) {
        float d = float(mc[i]);
        Sum  += d;
        Sum2 += d*d;
    }
    Sum  /= float(n); // !!!
    Sum2 /= float(n); // !!!
}