473,320 Members | 1,845 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,320 software developers and data experts.

trying to minimise the chi-squared statistic

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
Jul 18 '06 #1
2 2615
Banfa
9,065 Expert Mod 8TB
In theory you ought to be able to calculate K to give minimum chi-squared.

I note that

1. chi-squared is dependent on chi
2. chi is dependent on a[]
3. a[] is dependent on b[] squared and is therefore positive
4. b[] is dependent on I[+1] and input value and i[], as i[] -> I[+1] so chi will tend to zero

So to minimise chi-squared you need to minimise

(I[+1] - Imax * exp((lambda[j] - average)/X) + K) squared

At this point I think it is time to resort to a good maths book and some differentiation and may be the newton method (an iterative where you start with a guess and refine it).


BTW I notice that you never fill in lamba and I for the same index value j and that you aways refer to I[+1] and that where the data is declared you have some comments suggesting that the indexing is not working correctly. This is indicative of a programming error. I assume you want wavelegth intencity pairs but you input routine does not do that because it's logic is wrong.

You need something more like

Expand|Select|Wrap|Line Numbers
  1. bool intesityFirst;
  2. bool firstRead = true;
  3.  
  4. for (j = 0; j < m; j++) 
  5. {
  6.   read >> L;
  7.  
  8.   if (firstRead)
  9.   {
  10.     if (L - 1000 > 0) //this is to ensure only wavelengths, which have value > 1000, are equated to lambda[j]
  11.     {
  12.       intesityFirst = false;
  13.     }
  14.     else
  15.     {
  16.       intesityFirst = true;
  17.     }
  18.  
  19.     firstRead = false;
  20.   }
  21.  
  22.   if(intesityFirst)
  23.   {
  24.     I[j] = L;
  25.     cout << "I[" << j+1 << "] = " << I[j] << endl;
  26.   }
  27.   else
  28.   {
  29.     lambda[j] = L;
  30.     cout << "lambda[" << j+1 << "] = " << lambda[j] << endl;
  31.   }
  32.  
  33.   read >> L;
  34.   if(!intesityFirst)
  35.   {
  36.     I[j] = L;
  37.     cout << "I[" << j+1 << "] = " << I[j] << endl;
  38.   }
  39.   else
  40.   {
  41.     lambda[j] = L;
  42.     cout << "lambda[" << j+1 << "] = " << lambda[j] << endl;
  43.   }
  44. }
  45.  
This isn't terribly elegant, it would help it you knew if the file contained wavelengths or intensities first but it will result in the date for 1 entry all being in the same index of all the variables and enable you to stop allocating twice as much memory as required. In fact I would almost be inclined to use an array of structure to hold the data with all the data in 1 structure being associated with a single entry.
Jul 18 '06 #2
thanks for the suggestions Banfa, didn't think of using boolean expressions until you mentioned it
Pixie
Jul 19 '06 #3

Sign in to post your reply or Sign up for a free account.

Similar topics

5
by: jdh2358 | last post by:
I have a python file that is trying to read raw data from a raw partition on a dying dist, eg f = file('/dev/sda') f.seek(SOMEWHERE) s = f.read(SOMEBYTES) On some blocks, the read succeeds,...
1
by: Tom | last post by:
Hi friends. How do I override minimise event ? so that it does not minimise and simply hide my form ? the reason is I have the problem with the form after its been minimised calling activate...
2
by: A.M | last post by:
Hi, How can I make a window application that when I minimise it, it minimise to tray? Thanks, Alan
1
by: Paul [Paradise Solutions] | last post by:
Hi all Problem I want to minimise my app to the sytem tray ensuring that: 1) Application is not listed in the taskbar. 2) Application window does not minimise and hover as a small block above...
1
by: c language | last post by:
Hello Everybody, I am looking for a C code that I can use it to get Alpha of a Chi-Square value. Actually I want to give the Chi Squre value and Degree of Freedom (df) and get the alpha value...
5
by: smarty | last post by:
I have an example of code that allows and application to minise to the notification tray but how can I override the minimise and close buttons to ensure they always go back to the notify tray...
1
by: Mohsen | last post by:
Hello everyone, I have set values and also the degree of freedom of those values in every round of my simulation. I need to know the Chi-Square probabilities of those values in every round of...
5
by: DeanL | last post by:
Hi all, I'm trying to set up a query that runs from a command button on a form (simple enough so far), what I want the query to do is take values from the fields on the form (seven fields in...
1
by: Tradeorganizer | last post by:
Hi i want to use an excel formula in asp , is it possible the formula chidist in excel gives the value of chi - square distribution over the degree of frequency. but when i am using the same...
20
by: ziycon | last post by:
I have this code and when i minimise the application should disappear from the taskbar and only be visible in the tray and when i double click on the tray icon it show maximise the app. but neither...
0
by: DolphinDB | last post by:
The formulas of 101 quantitative trading alphas used by WorldQuant were presented in the paper 101 Formulaic Alphas. However, some formulas are complex, leading to challenges in calculation. Take...
0
by: DolphinDB | last post by:
Tired of spending countless mintues downsampling your data? Look no further! In this article, you’ll learn how to efficiently downsample 6.48 billion high-frequency records to 61 million...
0
by: ryjfgjl | last post by:
ExcelToDatabase: batch import excel into database automatically...
0
by: jfyes | last post by:
As a hardware engineer, after seeing that CEIWEI recently released a new tool for Modbus RTU Over TCP/UDP filtering and monitoring, I actively went to its official website to take a look. It turned...
0
by: ArrayDB | last post by:
The error message I've encountered is; ERROR:root:Error generating model response: exception: access violation writing 0x0000000000005140, which seems to be indicative of an access violation...
0
by: Defcon1945 | last post by:
I'm trying to learn Python using Pycharm but import shutil doesn't work
1
by: Shællîpôpï 09 | last post by:
If u are using a keypad phone, how do u turn on JavaScript, to access features like WhatsApp, Facebook, Instagram....
0
by: af34tf | last post by:
Hi Guys, I have a domain whose name is BytesLimited.com, and I want to sell it. Does anyone know about platforms that allow me to list my domain in auction for free. Thank you
0
by: Faith0G | last post by:
I am starting a new it consulting business and it's been a while since I setup a new website. Is wordpress still the best web based software for hosting a 5 page website? The webpages will be...

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.