Last active
January 15, 2018 16:46
-
-
Save ScratchyCode/c153207b389adf0f4eda14e05e639e1d to your computer and use it in GitHub Desktop.
Generate N events that follow the gaussian distribution and print statistics.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// Coded by ScratchyCode | |
// Compile in gcc with option -lm | |
#include <math.h> | |
#include <stdio.h> | |
#include <stdlib.h> | |
#include <time.h> | |
#define SUBSAMPLES 10 | |
int control(void); | |
double input(double min, double max, int divisor); | |
double realRand(double min, double max); | |
double randGaussInit(void); | |
double randGauss(double mu, double sigma); | |
void average(double N, double *mu, double array[]); | |
void sigma(double N, double *mu, double *sigma, double array[]); | |
int main(){ | |
int N, i, j; | |
double mu_dist, sigma_dist, *sample; | |
// control that gnuplot is present as program to plot data | |
int gnupl = control(); | |
if(gnupl == 1){ | |
printf("\nYou need gnuplot to graph the results."); | |
printf("\nInstall it with: sudo apt-get install gnuplot\n\n"); | |
exit(2); | |
} | |
srand48(time(NULL)); | |
printf("\nEnter events's number to be simulated: "); | |
N = input(10,1E12,SUBSAMPLES); | |
printf("\nEnter the distribution's average: "); | |
mu_dist = input(-1E12,1E12,1); | |
printf("\nEnter the distribution's standard deviation: "); | |
sigma_dist = input(0,1E12,1); | |
FILE *output = fopen("gauss.dat","w"); | |
if(output == NULL){ | |
perror("\nError"); | |
exit(1); | |
} | |
sample = calloc(N,sizeof(double)); | |
if(sample == NULL){ | |
perror("\nError"); | |
exit(1); | |
} | |
// events generation | |
for(i=0; i<N; i++){ | |
sample[i] = randGauss(mu_dist,sigma_dist); | |
fprintf(output,"%lf\n",sample[i]); | |
} | |
fflush(output); | |
// evaluate events's average and standard deviation | |
double mu_emp, sigma_emp; | |
average(N,&mu_emp,sample); | |
sigma(N,&mu_emp,&sigma_emp,sample); | |
printf("\n"); | |
printf("Events's average: %lf\n",mu_emp); | |
printf("Events's standard deviation: %lf\n",sigma_emp); | |
// evaluete SUBSAMPLES | |
double *subsample = calloc(N/SUBSAMPLES,sizeof(double)); | |
if(subsample == NULL){ | |
perror("\nError"); | |
exit(1); | |
} | |
double mu_sample[SUBSAMPLES], sigma_sample[SUBSAMPLES]; | |
for(i=0; i<SUBSAMPLES; i++){ | |
for(j=0; j<N/SUBSAMPLES; j++){ | |
subsample[j] = sample[i*SUBSAMPLES+j]; | |
} | |
average(N/SUBSAMPLES,mu_sample+i,subsample); | |
sigma(N/SUBSAMPLES,mu_sample+i,sigma_sample+i,subsample); | |
} | |
double averageAverages, sigmaDev; | |
average(SUBSAMPLES,&averageAverages,mu_sample); | |
sigma(SUBSAMPLES,&averageAverages,&sigmaDev,mu_sample); | |
printf("The average of averages of %d subsamples is %lf\n",SUBSAMPLES,averageAverages); | |
printf("The sigma of averages of %d subsamples is: %lf\n",SUBSAMPLES,sigmaDev); | |
fflush(output); | |
fclose(output); | |
return 0; | |
} | |
double input(double min, double max, int divisor){ | |
double x; | |
do{ | |
scanf("%lf",&x); | |
if(!(min<x && x<max && (int)(x)%divisor==0)){ | |
if(divisor == 1){ | |
printf("\nThe number must be between %.1lf and %.1lf.\n",min,max); | |
}else{ | |
printf("\nThe number must be between %.1lf and %.1lf and be multiple of %d\n",min,max,divisor); | |
} | |
printf("Try again: "); | |
} | |
}while(!(min<x && x<max && (int)(x)%divisor==0)); | |
return x; | |
} | |
void average(double N, double *mu, double array[]){ | |
int i; | |
double sum = 0; | |
for(i=0; i<N; i++){ | |
sum += array[i]; | |
} | |
*mu = sum/N; | |
return ; | |
} | |
void sigma(double N, double *mu, double *sigma, double array[]){ | |
int i; | |
double sum = 0; | |
for(i=0; i<N; i++){ | |
sum += pow((array[i] - *mu),2); | |
} | |
*sigma = sqrt(sum/(N-1)); | |
return ; | |
} | |
double realRand(double min, double max){ | |
double range = (max - min); | |
double divisore = RAND_MAX / range; | |
return min + (lrand48() / divisore); | |
} | |
double randGaussInit(void){ | |
double g = 0; | |
int i; | |
for(i=0; i<12; i++){ | |
g += realRand(-0.5,0.5); | |
} | |
return g; | |
} | |
double randGauss(double mu, double sigma){ | |
return randGaussInit() * sigma + mu; | |
} | |
int control(void){ | |
char path[] = "/usr/bin/gnuplot"; | |
FILE *pf = fopen(path,"r"); | |
if(pf == NULL){ | |
fclose(pf); | |
return 1; | |
}else{ | |
fclose(pf); | |
return 0; | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment