//
// verifica della relazione esistente tra
// una variabile aleatoria X e una variabile aleatoria y=f(x)
//
#include "TRandom.h"
const double PI =3.141592653589;
// funzione per calcolare la superficie dei cerchi
double f(double x) {
    return PI*x*x;
}
void verifica1(){
// media della distribuzione dei raggi
    double mX=5;
// Sigma della distribuzione dei raggi (errore di misura statistico)
    double sX=1;
// numero di esperimenti
    int N=100000;
    TRandom R(123456);
    double mediaX=0,mediaY=0;
    double mediaX2=0,mediaY2=0;
// loop principale per simulare gli N esperimenti
    for (int k=0;k<N;k++) {
// raggio del cerchio dell'esperimento k-esimo
        double x=R.Gaus(mX,sX);
// supoerficie del cerchio dell'esperimento k-esimo
        double y=f(x);
// per il calcolo di medie e varianze....
        mediaX+=x;
        mediaY+=y;
        mediaX2+=x*x;
        mediaY2+=y*y;
    }
// medie e varianze delle due grandezze
    mediaX=mediaX/N;
    mediaY=mediaY/N;
    mediaX2=mediaX2/N;
    mediaY2=mediaY2/N;
    double varX=(mediaX2-mediaX*mediaX);
    double varY=(mediaY2-mediaY*mediaY);
    double errX=sqrt(varX)/sqrt(N);  // perchè questa formula????
    double errY=sqrt(varY)/sqrt(N);  // perchè questa formula????
// calcolo teorico di media e errore sulla media della Y usando X
    double varYth=pow((2.*PI*mediaX),2)*errX*errX;
    double Yth=f(mediaX);
// termine correttivo
    double correct=0.5*(2*PI)*varX;
// output
    printf("\n\n\n media di X e var(X) %g %g\n",mediaX,varX);
    printf("Y=PI*X^2\n");
    printf("varY=(2*PI*x)^2*varX\n\n");
    printf("Risultati attesi\n superficie  media=%g errore=%g \n\n\n",Yth, sqrt(varYth));
    printf("Risultati ottenuti con la simulazione\n superficie media e errore %g %g \n",mediaY,errY);
    printf("Termine correttivo %g \n",correct);
}