статистика — C ++ вычисляет стандартное отклонение

Я пытаюсь написать простую программу анализа данных, но у меня проблемы с вычислением стандартного отклонения массива из десяти элементов

Мы должны использовать квадратное приближение (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

Я перепробовал много методов, но мне никогда не удавалось добиться каких-то последовательных результатов.

Я прилагаю несколько строк данных (каждая строка представляет собой блок, из которого должны быть вычислены медиана и стандартное отклонение) и полный код в надежде на любую помощь

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;
}