Hi,
I have a progarm which models some data I'm working with using a gaussian profile and then calculates chi-squared in order to compare the gaussian values to the real thing. I want to minimise chi-squared using the arbitrary parameter I have edfined as 'K', problem is I have no idea how to proceed and the book I have doesn't explain it very well. For clarity the data I am working with consists of 2 columns of numbers representing wavelengths and the intensity of electromagentic radiation at that wavelength form a word document, each column containing 2048 numbers (in total the document contain 4096 numbers). What I have so far follows;
--------------------------------------------------------------------------------------------------------------------
#include <iostream.h>
#include <fstream.h>
#include <math.h>
int main()
{
int m, j, x;
float lambda[9096], L, I[9096]; // note: lambda[9096] is a lie, 4096 not used as it result in calculation error
cout << "the value of m is "; //use 4096 to calculate average over all wavelengths
cin >> m;
x = m/2; //we only want every second term, but must read all terms. Thus m != number of terms for calculating average.
ifstream read("C:\\Documents and Settings\\Administrator\\My Documents\\My Honours\\qso.doc");
ofstream prt("C:\\Documents and Settings\\Administrator\\My documents\\My Honours\\average wavelength, standard deviation and maximum intenesity value for qso data.doc", ios::trunc);
if(!read)
{
cout << "input file failed to open";
system ("pause");
exit(1);
}
if(!prt)
{
cout << "output file 1 failed to open";
system ("pause");
exit(1);
}
for (j = 0; j < m; j++)
{
read >> L;
if(L - 1000 > 0) //this is to ensure only wavelengths, which have value > 1000, are equated to lambda[j]
{
lambda[j] = L;
cout << "lambda[" << j/2 + 1 << "] = " << lambda[j] << endl; //(j/2 + 1) as intensity values need to be ignored so zero isn't used for output
}
if(L - 1000 < 0) //this is to ensure only intensities, which have value , 1000, are equated to I[j]
{
I[j] = L;
cout << "I[" << j/2 +1 << "] = " << I[j] << endl; //(j/2 + 1) as wavelength values need to be ignored
}
}
read.close();
float average; //calculating average
average = 0;
for (j = 0; j < m; j++)
{
if(lambda[j] - 1000 > 0)
average = average + (lambda[j]/x);
else if(lambda[j] - 1000 < 0) //to avert intensity values being included in calculation
average = average;
}
float Imax; //determining maximum intensity value
Imax = I[0];
for (j = 0; j < m; j++)
{
if(I[j] - Imax > 0)
Imax = I[0] + I[j];
else if(I[j] - Imax <= 0)
Imax = Imax;
}
cout << "average wavelength = " << average << endl;
prt << "average wavelength = " << average << endl;
cout << "maximum intensity value = " << Imax << endl;
prt << "maximum intensity value = " << Imax << endl;
ofstream prt2("C:\\Documents and Settings\\Administrator\\My documents\\My Honours\\gaussian.doc", ios::trunc);
if(!prt2)
{
cout << "output file 2 failed to open";
system ("pause");
exit(1);
}
float sigma, S, X; //calculating standard deviation
sigma = 0;
for (j = 0; j < m; j++)
{
if(lambda[j] - 1000 > 0)
sigma = sigma + (((average - lambda[j])*(average - lambda[j]))/x);
else if(lambda[j] - 1000 < 0) //to avert intensity values being included in calculation
sigma = sigma;
//cout << "cumulative sigma = " << sigma << endl;
}
S = sqrt(sigma);
X = 2 * S * S * -1; //needed for gaussian form later
cout << "standard deviation = " << S << "\n" << "X = " << X << endl;
prt << "standard deviation = " << S << endl;
float i[9096], a[9096], b[9096], K; //using the gaussian formula to calculate intensity
for (j = 0; j < m; j++)
{
i[j] = Imax * exp((lambda[j] - average)/X) + K; //K an arbitrary parameter introduced to minimise chi squared
b[j] = I[j+1] - i[j]; //we require (lambda[j] - 1000 > 0) for i[j], but this causes I[j] = wavelength, not intensity
a[j] = pow(b[j] /S, 2);
//prt2 << "I[" << j/2 +1 << "] = " << I[j] << endl;
if(lambda[j] - 1000 > 0)
{
cout << "i[" << j/2 +1 << "] = " << i[j] << " ";
prt2 << "i[" << j/2 +1 << "] = " << i[j] << " ";
cout << "b[" << j/2 +1 << "] = " << b[j] << " ";
prt2 << "b[" << j/2 +1 << "] = " << b[j] << " ";
cout << "a[" << j/2 +1 << "] = " << a[j] << endl;
prt2 << "a[" << j/2 +1 << "] = " << a[j] << endl;
}
}
float chi; //calculating chi squared
chi = 0;
for (j = 0; j < m; j++)
{
if(lambda[j] - 1000 > 0)
chi = chi + a[j];
else if(lambda[j] - 1000 < 0)
chi = chi;
}
cout << "chi squared = " << chi << endl;
prt << "chi squared = " << chi << endl;
prt.close();
prt2.close();
system ("pause");
}
--------------------------------------------------------------------------------------------------------------------
as you may see, I've calculated chi-squared. What I need to do is minimise it by varying K somehow.
Pixie