Chi Square Minimization using initial guesses – C PROGRAM

Method # 1: Chi-square minimization

Fitting a dataset (xi,yi,si) using a curve y=f(x)=a*exp(bx), where values of the parameters a and b are unknown, but defined within a set of range.

Dataset (xi,yi,si) is the following:

Define 𝜒2(𝑎,𝑏)= ∑[(𝑦𝑖−𝑓(𝑥𝑖))/𝜎𝑖]2

Prob1: Tabulate and plot 𝜒2(𝑎=2.101,𝑏) vs b (in steps of 0.1).
(i) Find the value of b for which 𝜒2(𝑎=2.101,𝑏) is minimum 𝜒𝑚𝑖𝑛2.
(ii) Find the values of b for which 𝜒2 value is 𝜒𝑚𝑖𝑛2+1
(iii) Find the values of b for which 𝜒2 value is 𝜒𝑚𝑖𝑛2+4
(iv) Find the values of b for which 𝜒2 value is 𝜒𝑚𝑖𝑛2+9

Prob2: Tabulate and plot 𝜒2(𝑎,𝑏=0.498) vs a (in steps of 0.1).
(i) Find the value of a for which 𝜒2(𝑎,𝑏=0.498) is minimum 𝜒𝑚𝑖𝑛2.
(ii) Find the values of a for which 𝜒2 value is 𝜒𝑚𝑖𝑛2+1
(iii) Find the values of a for which 𝜒2 value is 𝜒𝑚𝑖𝑛2+4
(iv) Find the values of a for which 𝜒2 value is 𝜒𝑚𝑖𝑛2+9

Prob3: Tabulate and plot 𝜒2(𝑎,𝑏) vs a, b (both in steps of 0.1).
(i) Find the values of (𝑎,𝑏) for which 𝜒2(𝑎,𝑏) is minimum 𝜒𝑚𝑖𝑛2. Plot a vs b.
(ii) Find the values of (𝑎,𝑏) for which 𝜒2 value is 𝜒𝑚𝑖𝑛2+1. Plot a vs b.
(iii) Find the values of (𝑎,𝑏) for which 𝜒2 value is 𝜒𝑚𝑖𝑛2+4. Plot a vs b.
(iv) Find the values of (𝑎,𝑏) for which 𝜒2 value is 𝜒𝑚𝑖𝑛2+9. Plot a vs b.

CODE:

/******************************************************
*************Chi-square fitting**************
******************************************************/
#include<stdio.h>
#include<math.h>
double f(double a, double b, double x){
	return a*exp(b*x);
}
double chi(double a, double b, int n, double x[n], double y[n], double sig[n], double f(double a, double b, double x)){
	double sum=0;
	int i;
	for(i=0;i<n;i++){
		sum=sum+pow((y[i]-f(a,b,x[i]))/sig[i],2);
	}
	return sum;
}

main(){
	int n=7;
	double x[7]={1,2,3,4,5,6,7};
	double y[7]={4,5,8,16,30,38,70};
	double sig[7]={2,2,3,3,4,5,5};
	double a=2.101, b;
	double h=0.000001;
	int N=(2-0.1)/h;
	FILE *fp=NULL;
	fp=fopen("chi2.txt","w");
	int i=0;
	for(b=0.1;b<=2;b=b+h){
		fprintf(fp,"%lf\t%lf\n",b,chi(a,b,n,x,y,sig,f));
		i++;
	}
	
	FILE *fp1=NULL;
	fp1=fopen("chi2.txt","r");
	double bmin, chi2min, amin;
	fscanf(fp1,"%lf\t%lf\n",&bmin,&chi2min);
	int minIndex=0;
	for(i=0;i<N;i++){
		double btemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&btemp,&chitemp);
		if(chitemp<chi2min){
			chi2min=chitemp;
			bmin=btemp;
			minIndex=i;
		}
	}
	printf("The min. value of Chi2 is %lf for b=%lf\n\n",chi2min,bmin);
	fp1=fopen("chi2.txt","r");
	for(i=0;i<N;i++){
		double btemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&btemp,&chitemp);
		if(fabs(chitemp-chi2min-1)<=0.0005&&i<minIndex){
			printf("\n%lf\t%lf",btemp,chitemp);
			break;
		}
	}
	fp1=fopen("chi2.txt","r");
	for(i=0;i<N;i++){
		double btemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&btemp,&chitemp);
		if(fabs(chitemp-chi2min-1)<=0.0005&&i>minIndex){
			printf("\n%lf\t%lf",btemp,chitemp);
			break;
		}
	}
	fp1=fopen("chi2.txt","r");
	for(i=0;i<N;i++){
		double btemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&btemp,&chitemp);
		if(fabs(chitemp-chi2min-4)<=0.0005&&i<minIndex){
			printf("\n%lf\t%lf",btemp,chitemp);
			break;
		}
	}	
	fp1=fopen("chi2.txt","r");
	for(i=0;i<N;i++){
		double btemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&btemp,&chitemp);
		if(fabs(chitemp-chi2min-4)<=0.0005&&i>minIndex){
			printf("\n%lf\t%lf",btemp,chitemp);
			break;
		}
	}
	fp1=fopen("chi2.txt","r");
	for(i=0;i<N;i++){
		double btemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&btemp,&chitemp);
		if(fabs(chitemp-chi2min-9)<=0.0005&&i<minIndex){
			printf("\n%lf\t%lf",btemp,chitemp);
			break;
		}
	}
	fp1=fopen("chi2.txt","r");
	for(i=0;i<N;i++){
		double btemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&btemp,&chitemp);
		if(fabs(chitemp-chi2min-9)<=0.0005&&i>minIndex){
			printf("\n%lf\t%lf",btemp,chitemp);
			break;
		}
	}
	/************************************************************
	Part 2
	************************************************************/
	fp=fopen("chi3.txt","w");
	b=0.498;
	i=0;
	for(a=1;a<=4;a=a+h){
		fprintf(fp,"%lf\t%lf\n",a,chi(a,b,n,x,y,sig,f));
		i++;
	}
	fp1=fopen("chi3.txt","r");
	fscanf(fp1,"%lf\t%lf\n",&amin,&chi2min);
	minIndex=0;
	for(i=0;i<N;i++){
		double atemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&atemp,&chitemp);
		if(chitemp<chi2min){
			chi2min=chitemp;
			amin=atemp;
			minIndex=i;
		}
	}
	printf("\n\nThe min. value of Chi2 is %lf for a=%lf\n\n",chi2min,amin);
	fp1=fopen("chi3.txt","r");
	for(i=0;i<N;i++){
		double atemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&atemp,&chitemp);
		if(fabs(chitemp-chi2min-1)<=0.0005&&i<minIndex){
			printf("\n%lf\t%lf",atemp,chitemp);
			break;
		}
	}
	fp1=fopen("chi3.txt","r");
	for(i=0;i<N;i++){
		double atemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&atemp,&chitemp);
		if(fabs(chitemp-chi2min-1)<=0.0005&&i>minIndex){
			printf("\n%lf\t%lf",atemp,chitemp);
			break;
		}
	}
	fp1=fopen("chi3.txt","r");
	for(i=0;i<N;i++){
		double atemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&atemp,&chitemp);
		if(fabs(chitemp-chi2min-4)<=0.0005&&i<minIndex){
			printf("\n%lf\t%lf",atemp,chitemp);
			break;
		}
	}	
	fp1=fopen("chi3.txt","r");
	for(i=0;i<N;i++){
		double atemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&atemp,&chitemp);
		if(fabs(chitemp-chi2min-4)<=0.0005&&i>minIndex){
			printf("\n%lf\t%lf",atemp,chitemp);
			break;
		}
	}
	fp1=fopen("chi3.txt","r");
	for(i=0;i<N;i++){
		double atemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&atemp,&chitemp);
		if(fabs(chitemp-chi2min-9)<=0.0005&&i<minIndex){
			printf("\n%lf\t%lf",atemp,chitemp);
			break;
		}
	}
	fp1=fopen("chi3.txt","r");
	for(i=0;i<N;i++){
		double atemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\n",&atemp,&chitemp);
		if(fabs(chitemp-chi2min-9)<=0.0005&&i>minIndex){
			printf("\n%lf\t%lf",atemp,chitemp);
			break;
		}
	}	
	/*************************************************************
	Part 3
	*************************************************************/
	fp=fopen("chi4.txt","w");
	h=0.001;
	i=0;
	N=(2-0.1)/h*(4-1)/h;
	for(a=1;a<=4;a=a+h){
		for(b=0.1;b<=2;b=b+h){
			fprintf(fp,"%lf\t%lf\t%lf\n",a,b,chi(a,b,n,x,y,sig,f));
			i++;	
		}
	}
	fp1=fopen("chi4.txt","r");
	fscanf(fp1,"%lf\t%lf\t%lf\n",&amin,&bmin,&chi2min);
	minIndex=0;
	for(i=0;i<N;i++){
		double atemp;
		double btemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\t%lf\n",&atemp,&btemp,&chitemp);
		if(chitemp<chi2min){
			chi2min=chitemp;
			amin=atemp;
			bmin=btemp;
			minIndex=i;
		}
	}
	printf("\n\nThe min. value of Chi2 is %lf for a=%lf and b=%lf\n\n",chi2min,amin,bmin);
	FILE *fp2=fopen("chi2+1.txt","w");
	fp1=fopen("chi4.txt","r");
	for(i=0;i<N;i++){
		double atemp;
		double btemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\t%lf\n",&atemp,&btemp,&chitemp);
		if(fabs(chitemp-chi2min-1)<=0.08){
			fprintf(fp2,"\n%lf\t%lf\t%lf",atemp,btemp,chitemp);
			//break;
		}
	}
	fp2=fopen("chi2+4.txt","w");
	fp1=fopen("chi4.txt","r");
	for(i=0;i<N;i++){
		double atemp;
		double btemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\t%lf\n",&atemp,&btemp,&chitemp);
		if(fabs(chitemp-chi2min-4)<=0.01){
			fprintf(fp2,"\n%lf\t%lf\t%lf",atemp,btemp,chitemp);
			//break;
		}
	}
	fp2=fopen("chi2+9.txt","w");
	fp1=fopen("chi4.txt","r");
	for(i=0;i<N;i++){
		double atemp;
		double btemp;
		double chitemp;
		fscanf(fp1,"%lf\t%lf\t%lf\n",&atemp,&btemp,&chitemp);
		if(fabs(chitemp-chi2min-9)<=0.01){
			fprintf(fp2,"\n%lf\t%lf\t%lf",atemp,btemp,chitemp);
			//break;
		}
	}
}

OUTPUT:


I'm a physicist specializing in theoretical, computational and experimental condensed matter physics. I like to develop Physics related apps and softwares from time to time. Can code in most of the popular languages. Like to share my knowledge in Physics and applications using this Blog and a YouTube channel.



Leave a Reply

Your email address will not be published. Required fields are marked *