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:
Ph.D. researcher at Friedrich-Schiller University Jena, Germany. I’m a physicist specializing in computational material science. I write efficient codes for simulating light-matter interactions at atomic scales. I like to develop Physics, DFT, and Machine Learning related apps and software from time to time. Can code in most of the popular languages. I like to share my knowledge in Physics and applications using this Blog and a YouTube channel.