Use leasqr.m in octave to fit data with error

[higgs ()]

Dear Octave users,

I have a data file looks like,

data = [
 10 1684.023253675213 142923.5184171736
 11 1154.101457435897 79000.44034169338
 12 793.8254991880336 44430.77751038589
 13 548.7675124786326 24452.17758038106
...]

where the first column are variable values, the second one are data values and the third one are corresponding errors of data on second column.

I've used function leasqr.m to fit these data to some analytic function. It seems to work well but I'm still wondering that I've done right things with inputs. My question is:
Is it correct to set variable "wt" of the function equal to the inverse of error: wt = 1./data(:,3)

Many thanks,

For reference:

[f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]= ...
      leasqr(x,y,pin,F,stol,niter,wt,dp,dFdp,options)

% wt = column vector (dim=length(x)) of statistical weights.  These
%   should be set to be proportional to (sqrt of var(y))^-1; (That is,
%   the covariance matrix of the data is assumed to be proportional to
%   diagonal with diagonal equal to (wt.^2)^-1.  The constant of
%   proportionality will be estimated.); default = ones(length(y),1).

[Francesco Potortì]

>I have a data file looks like,
>
>data = [
> 10 1684.023253675213 142923.5184171736
> 11 1154.101457435897 79000.44034169338
> 12 793.8254991880336 44430.77751038589
> 13 548.7675124786326 24452.17758038106
>...]
>
>where the first column are variable values, the second one are data values
>and the third one are corresponding errors of data on second column.
>
>I've used function leasqr.m to fit these data to some analytic function. It
>seems to work well but I'm still wondering that I've done right things with
>inputs. My question is:
>Is it correct to set variable "wt" of the function equal to the inverse of
>error: wt = 1./data(:,3)

Yes, it is, as far as I know.

--
Francesco Potortì (ricercatore)        Voice: [phone] (op.2111)
ISTI - Area della ricerca CNR          Fax:   [phone]
via G. Moruzzi 1, I-56124 Pisa         Email: Potorti@...
(entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[higgs]

Thank you very much for your confirmation.
If it is familiar to you, could you please tell how can I extract "Chi-squared/ Number of Degree of freedom" from the program. Is there some option or I have to code it my own?

Best Regards,

Francesco Potorti`-2 wrote:

>I have a data file looks like,
>
>data = [
> 10 1684.023253675213 142923.5184171736
> 11 1154.101457435897 79000.44034169338
> 12 793.8254991880336 44430.77751038589
> 13 548.7675124786326 24452.17758038106
>...]
>
>where the first column are variable values, the second one are data values
>and the third one are corresponding errors of data on second column.
>
>I've used function leasqr.m to fit these data to some analytic function. It
>seems to work well but I'm still wondering that I've done right things with
>inputs. My question is:
>Is it correct to set variable "wt" of the function equal to the inverse of
>error: wt = 1./data(:,3)

Yes, it is, as far as I know.

--
Francesco Potortì (ricercatore)        Voice: [phone] (op.2111)
ISTI - Area della ricerca CNR          Fax:   [phone]
via G. Moruzzi 1, I-56124 Pisa         Email: [e-mail]
(entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/
_______________________________________________
Help-octave mailing list
[e-mail]
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Francesco Potortì]

>If it is familiar to you, could you please tell how can I extract
>"Chi-squared/ Number of Degree of freedom" from the program. Is there some
>option or I have to code it my own?

I do not know.  Anyway, have a look at the help and the source code:
there are various return variables that are not dispayed or returned by
default.  They should all printed out if you set

global verbose=1

before calling leasqr.  Maybe you can get what you want starting from
those.

--
Francesco Potortì (ricercatore)        Voice: [phone] (op.2111)
ISTI - Area della ricerca CNR          Fax:   [phone]
via G. Moruzzi 1, I-56124 Pisa         Email: Potorti@...
(entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[higgs]

Thank you very much. I read the code and found that chi2 out put has been removed. Thus I have to add some piece of code to get what I want.

Many thanks,

--- On Thu, 4/9/09, Francesco Potorti` <Potorti@...> wrote:

From: Francesco Potorti` <Potorti@...>
Subject: Re: Use leasqr.m in octave to fit data with error
To: "higgs" <nguyenoh@...>
Cc: help-octave@...
Date: Thursday, April 9, 2009, 4:59 AM

>If it is familiar to you, could you please tell how can I extract
>"Chi-squared/ Number of Degree of freedom" from the program. Is there some
>option or I have to code it my own?

I do not know.  Anyway, have a look at the help and the source code:
there are various return variables that are not dispayed or returned by
default.  They should all printed out if you set

global verbose=1

before calling leasqr.  Maybe you can get what you want starting from
those.

--
Francesco Potortì (ricercatore)        Voice: [phone] (op.2111)
ISTI - Area della ricerca CNR          Fax:   [phone]
via G. Moruzzi 1, I-56124 Pisa         Email: Potorti@...
(entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Francesco Potortì]

>Thank you very much. I read the code and found that chi2 out put has
>been removed. Thus I have to add some piece of code to get what I want.

Once you manage to get what you want, possibly in a generically useful
way, would you consider submitting a patch, for boh the code and the
documentation?

--
Francesco Potortì (ricercatore)        Voice: [phone] (op.2111)
ISTI - Area della ricerca CNR          Fax:   [phone]
via G. Moruzzi 1, I-56124 Pisa         Email: Potorti@...
(entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[higgs]

As I understood after skimming throught the code, it is very easy to extract Chi-squared since it was already there.. You can simply add the variable "sbest" (or "ss") to the output list than you have total chi-squared sum. After that if you want "Chi-squared / No. of degree of freedom", you just do a divison. I am not 100% sure for this issue because I have not had enough time to read the code carefully.

 

I just changed one line in the code from:

...

function [f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2]= ...
      leasqr(x,y,pin,F,stol,niter,wt,dp,dFdp,options)
....

to:

 

function [f,p,kvg,iter,corp,covp,covr,stdresid,Z,r2,sbest]= ...
      leasqr(x,y,pin,F,stol,niter,wt,dp,dFdp,options)

Best Regards,

--- On Thu, 4/9/09, Francesco Potorti` <Potorti@...> wrote:

From: Francesco Potorti` <Potorti@...>
Subject: Re: Use leasqr.m in octave to fit data with error
To: "Oanh Nguyen" <nguyenoh@...>
Cc: help-octave@...
Date: Thursday, April 9, 2009, 6:18 PM

>Thank you very much. I read the code and found that chi2 out put has
>been removed. Thus I have to add some piece of code to get what I want.

Once you manage to get what you want, possibly in a generically useful
way, would you consider submitting a patch, for boh the code and the
documentation?

--
Francesco Potortì (ricercatore)        Voice: [phone] (op.2111)
ISTI - Area della ricerca CNR          Fax:   [phone]
via G. Moruzzi 1, I-56124 Pisa         Email: Potorti@...
(entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave