Nonlinear fitting with lg(a+x)

[reposepuppy]

Hello,
I'm a college student who is attending a physical chemistry class now, and we're required to process the data obtained from experiments we've done to reach a plausible conclusion every week.
This time, I've got 2 columns of data, and one of them is the concentration of a solution [c], the other the surface tension of the solution [r]; the task is to draw the curve representing the function of F(c,r)=0, and use it to draw the curve representing the function of F(c, dr/dc)=0 to reach some conclusion。I found out on the Internet that there is a function between them stated as :
                             r=r0[1-a*lg(1+c/b)]
where r0 is the surface tension when there is no solute in the solution, or in other words, c=0; it could be measured during the experiment. The constants a and b could be determined by making a nonlinear fitting using the data I've got.
That's the question: I don't know how to transform this nonlinear function into a linear one; I could only make it (1-r/r0)=a*lg(b+c)-a*lg(b) and using the trial-and-error method(*See below) to find out what a and b are. But the error is too much for me(about 30% away from the theoretical one).

Does anyone know how to solve this fitting problem using solely Octave?
Thanks for help!

* The trial-and-error method I used is here, using Microsoft Excel (I could use Octave only to deal with some integration problem and linear fitting...):
First, input the data into Excel and draw a XY scattered diagram, where X=(1-r/r0) and Y=lg(b+c). A linear distribution of dots would appear if the value of b was right.
Second, calculate the R^2(the square of correlation coefficient) . It would be a value very close to 1 if the value of b was right.
Third, change the b value until I can't make the R^2 higher. First try from 0 to 10, and find out the optimal value, and minimize the scale by 10 fold and try again... I tried from 0 to 2, and found out the R^2 will be the biggest when b was 0.6; and I tried from 0.50 to 0.70 and I found out it's 0.59~0.62(the R^2 will not change any more during this range).  

[Jaroslav Hajek-2]

On Fri, Mar 13, 2009 at 9:20 AM, reposepuppy <reposepuppy@...> wrote:

>
> Hello,
> I'm a college student who is attending a physical chemistry class now, and
> we're required to process the data obtained from experiments we've done to
> reach a plausible conclusion every week.
> This time, I've got 2 columns of data, and one of them is the concentration
> of a solution [c], the other the surface tension of the solution [r]; the
> task is to draw the curve representing the function of F(c,r)=0, and use it
> to draw the curve representing the function of F(c, dr/dc)=0 to reach some
> conclusion。I found out on the Internet that there is a function between them
> stated as :
>                             r=r0[1-a*lg(1+c/b)]
> where r0 is the surface tension when there is no solute in the solution, or
> in other words, c=0; it could be measured during the experiment. The
> constants a and b could be determined by making a nonlinear fitting using
> the data I've got.
> That's the question: I don't know how to transform this nonlinear function
> into a linear one; I could only make it (1-r/r0)=a*lg(b+c)-a*lg(b) and using
> the trial-and-error method(*See below) to find out what a and b are. But the
> error is too much for me(about 30% away from the theoretical one).
>
> Does anyone know how to solve this fitting problem using solely Octave?
> Thanks for help!
>

fsolve is now capable of solving over-determined nonlinear systems -
in other words, nonlinear fitting (with moderate residuals).

function r = model (r0, a, b, c)
  r = r0*(1 - a*log(1 + c/b));
endfunction

... get c, r

model_error = @(par) r - model (par(1), par(2), par(3), c);
par0 = [... initial guesses for r0, a, b]

par = fsolve (model_error, par0, ...options)

It will be available in 3.2, so your options are
1. wait for 3.2 release
2. compile development sources
3. send me your data and I can make the fit for you

> * The trial-and-error method I used is here, using Microsoft Excel (I could
> use Octave only to deal with some integration problem and linear
> fitting...):
> First, input the data into Excel and draw a XY scattered diagram, where
> X=(1-r/r0) and Y=lg(b+c). A linear distribution of dots would appear if the
> value of b was right.
> Second, calculate the R^2(the square of correlation coefficient) . It would
> be a value very close to 1 if the value of b was right.
> Third, change the b value until I can't make the R^2 higher. First try from
> 0 to 10, and find out the optimal value, and minimize the scale by 10 fold
> and try again... I tried from 0 to 2, and found out the R^2 will be the
> biggest when b was 0.6; and I tried from 0.50 to 0.70 and I found out it's
> 0.59~0.62(the R^2 will not change any more during this range).

These seem to be fairly good guesses. I think fsolve should not have
problems refining them.
For a more automatized method of obtaining initial guesses, you can
sample expected ranges of
parameters and pick the best match.

cheers

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Michael Creel]

reposepuppy wrote:

Hello,
I'm a college student who is attending a physical chemistry class now, and we're required to process the data obtained from experiments we've done to reach a plausible conclusion every week.
This time, I've got 2 columns of data, and one of them is the concentration of a solution [c], the other the surface tension of the solution [r]; the task is to draw the curve representing the function of F(c,r)=0, and use it to draw the curve representing the function of F(c, dr/dc)=0 to reach some conclusion。I found out on the Internet that there is a function between them stated as :
                             r=r0[1-a*lg(1+c/b)]
where r0 is the surface tension when there is no solute in the solution, or in other words, c=0; it could be measured during the experiment. The constants a and b could be determined by making a nonlinear fitting using the data I've got.
That's the question: I don't know how to transform this nonlinear function into a linear one; I could only make it (1-r/r0)=a*lg(b+c)-a*lg(b) and using the trial-and-error method(*See below) to find out what a and b are. But the error is too much for me(about 30% away from the theoretical one).

Does anyone know how to solve this fitting problem using solely Octave?
Thanks for help!

* The trial-and-error method I used is here, using Microsoft Excel (I could use Octave only to deal with some integration problem and linear fitting...):
First, input the data into Excel and draw a XY scattered diagram, where X=(1-r/r0) and Y=lg(b+c). A linear distribution of dots would appear if the value of b was right.
Second, calculate the R^2(the square of correlation coefficient) . It would be a value very close to 1 if the value of b was right.
Third, change the b value until I can't make the R^2 higher. First try from 0 to 10, and find out the optimal value, and minimize the scale by 10 fold and try again... I tried from 0 to 2, and found out the R^2 will be the biggest when b was 0.6; and I tried from 0.50 to 0.70 and I found out it's 0.59~0.62(the R^2 will not change any more during this range).  

You don't need to linearize anything, just use nonlinear least squares to minimize the function

     sum_over_observations [ r-r0[1-a*lg(1+c/b)] ]^2  with respect to the unknown constants.

If you linearize and then estimate, the resulting estimator is biased and inconsistent. All you need to do is write Octave code to calculate the objective function. Then you can use one of the available minimizers to get the estimates. Octave provides sqp, which can handle this problem.

[reposepuppy]

Sorry for late reply due to lack of time...
Thanks for the replies!

Jaroslav Hajek-2 wrote:

On Fri, Mar 13, 2009 at 9:20 AM, reposepuppy <[e-mail]> wrote:
>
> Hello,
> I'm a college student who is attending a physical chemistry class now, and
> we're required to process the data obtained from experiments we've done to
> reach a plausible conclusion every week.
> This time, I've got 2 columns of data, and one of them is the concentration
> of a solution [c], the other the surface tension of the solution [r]; the
> task is to draw the curve representing the function of F(c,r)=0, and use it
> to draw the curve representing the function of F(c, dr/dc)=0 to reach some
> conclusion。I found out on the Internet that there is a function between them
> stated as :
>                             r=r0[1-a*lg(1+c/b)]
> where r0 is the surface tension when there is no solute in the solution, or
> in other words, c=0; it could be measured during the experiment. The
> constants a and b could be determined by making a nonlinear fitting using
> the data I've got.
> That's the question: I don't know how to transform this nonlinear function
> into a linear one; I could only make it (1-r/r0)=a*lg(b+c)-a*lg(b) and using
> the trial-and-error method(*See below) to find out what a and b are. But the
> error is too much for me(about 30% away from the theoretical one).
>
> Does anyone know how to solve this fitting problem using solely Octave?
> Thanks for help!
>

fsolve is now capable of solving over-determined nonlinear systems -
in other words, nonlinear fitting (with moderate residuals).

function r = model (r0, a, b, c)
  r = r0*(1 - a*log(1 + c/b));
endfunction

... get c, r

model_error = @(par) r - model (par(1), par(2), par(3), c);
par0 = [... initial guesses for r0, a, b]

par = fsolve (model_error, par0, ...options)

It will be available in 3.2, so your options are
1. wait for 3.2 release
2. compile development sources
3. send me your data and I can make the fit for you


> * The trial-and-error method I used is here, using Microsoft Excel (I could
> use Octave only to deal with some integration problem and linear
> fitting...):
> First, input the data into Excel and draw a XY scattered diagram, where
> X=(1-r/r0) and Y=lg(b+c). A linear distribution of dots would appear if the
> value of b was right.
> Second, calculate the R^2(the square of correlation coefficient) . It would
> be a value very close to 1 if the value of b was right.
> Third, change the b value until I can't make the R^2 higher. First try from
> 0 to 10, and find out the optimal value, and minimize the scale by 10 fold
> and try again... I tried from 0 to 2, and found out the R^2 will be the
> biggest when b was 0.6; and I tried from 0.50 to 0.70 and I found out it's
> 0.59~0.62(the R^2 will not change any more during this range).

These seem to be fairly good guesses. I think fsolve should not have
problems refining them.
For a more automatized method of obtaining initial guesses, you can
sample expected ranges of
parameters and pick the best match.

cheers

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
_______________________________________________
Help-octave mailing list
[e-mail]
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

Thanks for the suggestions, but I guess I can't wait for the 3.2 version and I can't even understand the codes you've written(because I'm really a newbie in Octave...), while the function produced from fitting is flexible(polynomial ones are allowed). So I made a polynomial one instead.
Still, I want to find out what the codes represents, and I can't find what "@par" and the codes below mean:
par0 = [... initial guesses for r0, a, b]

par = fsolve (model_error, par0, ...options)
Could you make a brief explanation for me?

Michael Creel wrote:

You don't need to linearize anything, just use nonlinear least squares to minimize the function

     sum_over_observations [ r-r0[1-a*lg(1+c/b)] ]^2  with respect to the unknown constants.

If you linearize and then estimate, the resulting estimator is biased and inconsistent. All you need to do is write Octave code to calculate the objective function. Then you can use one of the available minimizers to get the estimates. Octave provides sqp, which can handle this problem.

I can't find the function sum_over_observations... Is it a function I have to make up by myself, or some packages submitted somewhere other than octave-forge?
And..."All you need to do is write Octave code to calculate the objective function."
It might sound silly, but I really don't know (and I don't know where I can find the tutorial) how to calculate it, and how to use the minimizers. I've read the octave.pdf file and the function reference on octave-forge, and I could only found one minimizer called polyfit(x,y,n), which is used to make polynomial fittings...
Could you give me a detailed demonstration/code to work it out? The data has been collected in the attachment.data+c+and+r.txt

[Michael Creel]

reposepuppy wrote:

Michael Creel wrote:

You don't need to linearize anything, just use nonlinear least squares to minimize the function

     sum_over_observations [ r-r0[1-a*lg(1+c/b)] ]^2  with respect to the unknown constants.

If you linearize and then estimate, the resulting estimator is biased and inconsistent. All you need to do is write Octave code to calculate the objective function. Then you can use one of the available minimizers to get the estimates. Octave provides sqp, which can handle this problem.

I can't find the function sum_over_observations... Is it a function I have to make up by myself, or some packages submitted somewhere other than octave-forge?
And..."All you need to do is write Octave code to calculate the objective function."
It might sound silly, but I really don't know (and I don't know where I can find the tutorial) how to calculate it, and how to use the minimizers. I've read the octave.pdf file and the function reference on octave-forge, and I could only found one minimizer called polyfit(x,y,n), which is used to make polynomial fittings...
Could you give me a detailed demonstration/code to work it out? The data has been collected in the attachment.data+c+and+r.txt

For an example of minimization of a function, you could install the package "optim", and then run "bfgsmin_example". http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/optim/inst/bfgsmin_example.m?revision=HEAD&content-type=text/plain

For your problem, you  would just need to replace the objective function in that example with the one for your problem, Doing that is not hard if you know a little about Octave.

[Bertrand Roessli]

Hello,

I am happy to see that there is some discussion about data fitting.

With the octave-forge package, bfgsmin_example works fine.
However, if I try to fit an exponential function (from leasqrexample.m),
bfgsmin is unable to converge (which is easily done with the other
available functions (simplex, samin,...).

The main problem seems to originate from the numerical derivatives.
The output reads

octave:1> leasqrexamp
chisq =  0.86656
chisq =  0.86656
chisq =  0.86656
chisq =  0.86657
chisq =  0.86657
chisq =  0.86656
error: __bfgsmin_gradient: gradient contains NaNs or Inf
error: feval: function `chi2' not found
error: octave_base_value::double_value (): wrong type argument `<unknown
type>'
warning: __bfgsmin_obj: objective function could not be evaluated -
setting to DBL_MAX
error: feval: function `chi2' not found

If somebody has an idea I would appreciate.

This is the code
----------------

1;
function [obj_value]=chi2(p,x,y)

f=leasqrfunc(x,p);
v=length(y);
wt=1./sqrt(y);
chisq=sum(((y.-real(f)).*wt).^2 )/v

obj_value=chisq;

endfunction


t = [1:100];
p=[10; 0.1];
data=leasqrfunc(t,p);

%fprintf(1,'\nAdding Random Noise.\n data=data+0.05*rand(1,100);\n');
data=data+0.05*rand(1,100);
plot(t,data,"@12");
fflush(stdout);
pause;

control = {-1;1};  # maxiters, verbosity, conv. reg., arg_to_min

tic;
param=[11;1];
bfgsmin("chi2",{param,t,data},control)
toc

 
function y = leasqrfunc(x,p)
% leasqrfunc : this is a function of 'x' and 'p' parameters for
leasqrexamp.

% Description: leasqr example fit function

% example of function.
y=p(1)*exp(-p(2)*x);

endfunction


On Mon, 2009-03-16 at 04:27 -0700, Michael Creel wrote:

>
>
> reposepuppy wrote:
> >
> >
> >
> > Michael Creel wrote:
> >>
> >> You don't need to linearize anything, just use nonlinear least squares to
> >> minimize the function
> >>
> >>      sum_over_observations [ r-r0[1-a*lg(1+c/b)] ]^2  with respect to the
> >> unknown constants.
> >>
> >> If you linearize and then estimate, the resulting estimator is biased and
> >> inconsistent. All you need to do is write Octave code to calculate the
> >> objective function. Then you can use one of the available minimizers to
> >> get the estimates. Octave provides sqp, which can handle this problem.
> >>
> >
> > I can't find the function sum_over_observations... Is it a function I have
> > to make up by myself, or some packages submitted somewhere other than
> > octave-forge?
> > And..."All you need to do is write Octave code to calculate the objective
> > function."
> > It might sound silly, but I really don't know (and I don't know where I
> > can find the tutorial) how to calculate it, and how to use the minimizers.
> > I've read the octave.pdf file and the function reference on octave-forge,
> > and I could only found one minimizer called polyfit(x,y,n), which is used
> > to make polynomial fittings...
> > Could you give me a detailed demonstration/code to work it out? The data
> > has been collected in the attachment.
> > http://www.nabble.com/file/p22535164/data%2Bc%2Band%2Br.txt
> > data+c+and+r.txt
> >
>
> For an example of minimization of a function, you could install the package
> "optim", and then run "bfgsmin_example".
> http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/optim/inst/bfgsmin_example.m?revision=HEAD&content-type=text/plain
>
> For your problem, you  would just need to replace the objective function in
> that example with the one for your problem, Doing that is not hard if you
> know a little about Octave.
>

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

[Jaroslav Hajek-2]

On Mon, Mar 16, 2009 at 11:21 AM, reposepuppy <reposepuppy@...> wrote:

> Thanks for the suggestions, but I guess I can't wait for the 3.2 version and
> I can't even understand the codes you've written(because I'm really a newbie
> in Octave...), while the function produced from fitting is
> flexible(polynomial ones are allowed). So I made a polynomial one instead.
> Still, I want to find out what the codes represents, and I can't find what
> "@par" and the codes below mean:
> par0 = [... initial guesses for r0, a, b]
>
> par = fsolve (model_error, par0, ...options)
> Could you make a brief explanation for me?
>

It's simple - fsolve can't work with the parameters as individual
variables, but needs them packed into a vector, because it uses linear
algebra.

The statement
model_error = @(par) r - model (par(1), par(2), par(3), c);

defines a new function with a single parameter, par, that is expected
to be a vector of [r0, a, b] and calculates teh residual vector by
calling the model function with the parameters par(1), par(2), par(3)
and independent variables c and subtracts the dependent observations r
to get a vector of the residuals.

fsolve then attempts to find a vector par_opt such that
norm(model_error(par_opt)) is as small as possible, which is what you
need.

par0 = [... initial guesses for r0, a, b]

means that you need to supply initial guesses for the parameters. If
you have absolutely no idea, just use zeros, but judging by your post
this was not the case. The better the initial guess, the better
results, obviously.

Another code for nonlinear least squares is leasqr in the optimization package.

You can also use a general nonlinear minimization function, as
suggested by M. Creel. However, that usually means not exploiting the
special structure of the minimization problem (sum of squares).
The general rule is that Newton-like solvers (fsolve and leasqr) are
usually better when the expected residual is small (good fit) and
general minimization solvers (like bfgsmin) usually do equally good or
better job when the residual is large (bad fit).

cheers

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Michael Creel]

Bertrand Roessli wrote:

Hello,

I am happy to see that there is some discussion about data fitting.

With the octave-forge package, bfgsmin_example works fine.
However, if I try to fit an exponential function (from leasqrexample.m),
bfgsmin is unable to converge (which is easily done with the other
available functions (simplex, samin,...).

The main problem seems to originate from the numerical derivatives.
The output reads

octave:1> leasqrexamp
chisq =  0.86656
chisq =  0.86656
chisq =  0.86656
chisq =  0.86657
chisq =  0.86657
chisq =  0.86656
error: __bfgsmin_gradient: gradient contains NaNs or Inf
error: feval: function `chi2' not found
error: octave_base_value::double_value (): wrong type argument `<unknown
type>'
warning: __bfgsmin_obj: objective function could not be evaluated -
setting to DBL_MAX
error: feval: function `chi2' not found

If somebody has an idea I would appreciate.

This is the code
----------------

1;
function [obj_value]=chi2(p,x,y)

f=leasqrfunc(x,p);
v=length(y);
wt=1./sqrt(y);
chisq=sum(((y.-real(f)).*wt).^2 )/v

obj_value=chisq;

endfunction


t = [1:100];
p=[10; 0.1];
data=leasqrfunc(t,p);

%fprintf(1,'\nAdding Random Noise.\n data=data+0.05*rand(1,100);\n');
data=data+0.05*rand(1,100);
plot(t,data,"@12");
fflush(stdout);
pause;

control = {-1;1};  # maxiters, verbosity, conv. reg., arg_to_min

tic;
param=[11;1];
bfgsmin("chi2",{param,t,data},control)
toc

 
function y = leasqrfunc(x,p)
% leasqrfunc : this is a function of 'x' and 'p' parameters for
leasqrexamp.

% Description: leasqr example fit function

% example of function.
y=p(1)*exp(-p(2)*x);

endfunction

Hmm, I ran your code and it works for me:
chisq =  0.0084524                            
chisq =  0.0084524                            
chisq =  0.0084524                            
chisq =  0.0084524                            
chisq =  0.0084524                            
chisq =  0.0084524                            
chisq =  0.0084524                            
------------------------------------------------
bfgsmin iteration 15  convergence (f g p): 1 1 1

function value: 0.00845243  stepsize: 0.1375  

          param    gradient (n)          change
        9.69388        -0.00000         0.00000
        0.09558         0.00001        -0.00000
------------------------------------------------
bfgsmin final results: 15 iterations

function value: 0.00845243

STRONG CONVERGENCE
Function conv 1  Param conv 1  Gradient conv 1

          param    gradient (n)          change
        9.69388        -0.00000         0.00000
        0.09558         0.00001        -0.00000
ans =

   9.693878
   0.095576

Elapsed time is 0.084275 seconds.
octave:3>

Maybe you have an old version of bfgsmin? Every once in a while I find a bug and fix it. I believe that it has been in pretty good shape for a while now.
Michael

[Bertrand Roessli]

Hello,

Thanks for trying the code.

I have the latest package from octave sourceforge (optim 1.04) and
octave 3.1.54.

The thing is that on my laptop the result depends on the random
numbers...

e.g. If I run leasqrexample 2 times I get

first run:
---------
octave:1> leasqrexamp
------------------------------------------------
bfgsmin final results: 14 iterations

function value: 0.00775013

STRONG CONVERGENCE
Function conv 1  Param conv 1  Gradient conv 1

          param    gradient (n)          change
        9.73062         0.00000         0.00000
        0.09587        -0.00000         0.00000
ans =

   9.730625
   0.095865

Elapsed time is 0.15632 seconds.

Second run:
-----------
octave:2> leasqrexamp
error: __bfgsmin_gradient: gradient contains NaNs or Inf
error: feval: function `chi2' not found
error: octave_base_value::double_value (): wrong type argument `<unknown
type>'

Any idea?

Bertrand
 
On Mon, 2009-03-16 at 05:43 -0700, Michael Creel wrote:

>
>
> Bertrand Roessli wrote:
> >
> > Hello,
> >
> > I am happy to see that there is some discussion about data fitting.
> >
> > With the octave-forge package, bfgsmin_example works fine.
> > However, if I try to fit an exponential function (from leasqrexample.m),
> > bfgsmin is unable to converge (which is easily done with the other
> > available functions (simplex, samin,...).
> >
> > The main problem seems to originate from the numerical derivatives.
> > The output reads
> >
> > octave:1> leasqrexamp
> > chisq =  0.86656
> > chisq =  0.86656
> > chisq =  0.86656
> > chisq =  0.86657
> > chisq =  0.86657
> > chisq =  0.86656
> > error: __bfgsmin_gradient: gradient contains NaNs or Inf
> > error: feval: function `chi2' not found
> > error: octave_base_value::double_value (): wrong type argument `<unknown
> > type>'
> > warning: __bfgsmin_obj: objective function could not be evaluated -
> > setting to DBL_MAX
> > error: feval: function `chi2' not found
> >
> > If somebody has an idea I would appreciate.
> >
> > This is the code
> > ----------------
> >
> > 1;
> > function [obj_value]=chi2(p,x,y)
> >
> > f=leasqrfunc(x,p);
> > v=length(y);
> > wt=1./sqrt(y);
> > chisq=sum(((y.-real(f)).*wt).^2 )/v
> >
> > obj_value=chisq;
> >
> > endfunction
> >
> >
> > t = [1:100];
> > p=[10; 0.1];
> > data=leasqrfunc(t,p);
> >
> > %fprintf(1,'\nAdding Random Noise.\n data=data+0.05*rand(1,100);\n');
> > data=data+0.05*rand(1,100);
> > plot(t,data,"@12");
> > fflush(stdout);
> > pause;
> >
> > control = {-1;1};  # maxiters, verbosity, conv. reg., arg_to_min
> >
> > tic;
> > param=[11;1];
> > bfgsmin("chi2",{param,t,data},control)
> > toc
> >
> >  
> > function y = leasqrfunc(x,p)
> > % leasqrfunc : this is a function of 'x' and 'p' parameters for
> > leasqrexamp.
> >
> > % Description: leasqr example fit function
> >
> > % example of function.
> > y=p(1)*exp(-p(2)*x);
> >
> > endfunction
> >
> >
> >
>
> Hmm, I ran your code and it works for me:
> chisq =  0.0084524                            
> chisq =  0.0084524                            
> chisq =  0.0084524                            
> chisq =  0.0084524                            
> chisq =  0.0084524                            
> chisq =  0.0084524                            
> chisq =  0.0084524                            
> ------------------------------------------------
> bfgsmin iteration 15  convergence (f g p): 1 1 1
>
> function value: 0.00845243  stepsize: 0.1375  
>
>           param    gradient (n)          change
>         9.69388        -0.00000         0.00000
>         0.09558         0.00001        -0.00000
> ------------------------------------------------
> bfgsmin final results: 15 iterations
>
> function value: 0.00845243
>
> STRONG CONVERGENCE
> Function conv 1  Param conv 1  Gradient conv 1
>
>           param    gradient (n)          change
>         9.69388        -0.00000         0.00000
>         0.09558         0.00001        -0.00000
> ans =
>
>    9.693878
>    0.095576
>
> Elapsed time is 0.084275 seconds.
> octave:3>
>
> Maybe you have an old version of bfgsmin? Every once in a while I find a bug
> and fix it. I believe that it has been in pretty good shape for a while now.
> Michael
>

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

[Michael Creel]

Bertrand Roessli wrote:

Hello,

Thanks for trying the code.

I have the latest package from octave sourceforge (optim 1.04) and
octave 3.1.54.

It seems that the optim package was last released in August of 2008, while the last change to __bfgsmin.cc was in October 2008. I just ran your code 10 times without problems, so I think that the problem is a bug in the old version of __bfgsmin.cc. Just download that file to your working directory, run mkoctfile __bfgsmin.cc, and I think you'll be set.
Michael

p.s., I have not tested it with Octave 3.1.x,

[Bertrand Roessli]

Sorry I do not get it,

"Just download that file to your working directory,"

which file?

If you mean from optim-1.04, it is what I did before.

Bertrand

On Mon, 2009-03-16 at 07:00 -0700, Michael Creel wrote:

>
> Bertrand Roessli wrote:
> >
> > Hello,
> >
> > Thanks for trying the code.
> >
> > I have the latest package from octave sourceforge (optim 1.04) and
> > octave 3.1.54.
> >
> >
>
> It seems that the optim package was last released in August of 2008, while
> the last change to __bfgsmin.cc was in October 2008. I just ran your code 10
> times without problems, so I think that the problem is a bug in the old
> version of __bfgsmin.cc. Just download that file to your working directory,
> run mkoctfile __bfgsmin.cc, and I think you'll be set.
> Michael
>
> p.s., I have not tested it with Octave 3.1.x,
>

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

[Michael Creel]

Bertrand Roessli wrote:

Sorry I do not get it,

"Just download that file to your working directory,"

which file?

If you mean from optim-1.04, it is what I did before.

Bertrand

On Mon, 2009-03-16 at 07:00 -0700, Michael Creel wrote:
>
> Bertrand Roessli wrote:
> >
> > Hello,
> >
> > Thanks for trying the code.
> >
> > I have the latest package from octave sourceforge (optim 1.04) and
> > octave 3.1.54.
> >
> >
>
> It seems that the optim package was last released in August of 2008, while
> the last change to __bfgsmin.cc was in October 2008. I just ran your code 10
> times without problems, so I think that the problem is a bug in the old
> version of __bfgsmin.cc. Just download that file to your working directory,
> run mkoctfile __bfgsmin.cc, and I think you'll be set.
> Michael
>
> p.s., I have not tested it with Octave 3.1.x,
>

_______________________________________________
Help-octave mailing list
[e-mail]
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

You can browse the SVN repository to get individual files. __bfgsmin.cc is here: http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/optim/src/__bfgsmin.cc?view=log

Just click on that, then download the latest version. You can browse around the SNV repo to see the stucture, if you're curious. Because __bfgsmin.cc is written in C++, it needs to be compiled before it can be used by Octave. That's why you need the mkoctfile step I mentioned before.
Cheers, Michael

[Bertrand Roessli]

Hello,


I compiled everything (mkoctifle) from the SVN repository but it is not
better. The bfgsmin_example.m works perfectly, but my function not.

The funny thing is that bfgsmin makes a few cycles, chi2 is not modified
as the parameters are not varied and then crashes as soon as there is a
change in the parameters.


octave:1> leasqrexamp
chisq =  1.8953
p =

   11
    1

chisq =  1.8953
p =

   11
    1

chisq =  1.8953
p =

   11.0000
    1.0000

chisq =  1.8953
p =

   11.0000
    1.0000

chisq =  1.8953
p =

   11.0000
    1.0000

chisq =  1.8953
p =

   11.00000
    0.99999

error: __bfgsmin_gradient: gradient contains NaNs or Inf
error: feval: function `chi2' not found
error: octave_base_value::double_value (): wrong type argument `<unknown
type>'
warning: __bfgsmin_obj: objective function could not be evaluated -
setting to DBL_MAX
error: feval: function `chi2' not found

Thanks,

Bertrand

On Mon, 2009-03-16 at 10:02 -0700, Michael Creel wrote:

>
>
> Bertrand Roessli wrote:
> >
> > Sorry I do not get it,
> >
> > "Just download that file to your working directory,"
> >
> > which file?
> >
> > If you mean from optim-1.04, it is what I did before.
> >
> > Bertrand
> >
> > On Mon, 2009-03-16 at 07:00 -0700, Michael Creel wrote:
> >>
> >> Bertrand Roessli wrote:
> >> >
> >> > Hello,
> >> >
> >> > Thanks for trying the code.
> >> >
> >> > I have the latest package from octave sourceforge (optim 1.04) and
> >> > octave 3.1.54.
> >> >
> >> >
> >>
> >> It seems that the optim package was last released in August of 2008,
> >> while
> >> the last change to __bfgsmin.cc was in October 2008. I just ran your code
> >> 10
> >> times without problems, so I think that the problem is a bug in the old
> >> version of __bfgsmin.cc. Just download that file to your working
> >> directory,
> >> run mkoctfile __bfgsmin.cc, and I think you'll be set.
> >> Michael
> >>
> >> p.s., I have not tested it with Octave 3.1.x,
> >>
> >
> > _______________________________________________
> > Help-octave mailing list
> > Help-octave@...
> > https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
> >
> >
>
> You can browse the SVN repository to get individual files. __bfgsmin.cc is
> here:
> http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/optim/src/__bfgsmin.cc?view=log
>
> Just click on that, then download the latest version. You can browse around
> the SNV repo to see the stucture, if you're curious. Because __bfgsmin.cc is
> written in C++, it needs to be compiled before it can be used by Octave.
> That's why you need the mkoctfile step I mentioned before.
> Cheers, Michael
>

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

[Michael Creel]

Bertrand Roessli wrote:

Hello,


I compiled everything (mkoctifle) from the SVN repository but it is not
better. The bfgsmin_example.m works perfectly, but my function not.

The funny thing is that bfgsmin makes a few cycles, chi2 is not modified
as the parameters are not varied and then crashes as soon as there is a
change in the parameters.


octave:1> leasqrexamp
chisq =  1.8953
p =

   11
    1

chisq =  1.8953
p =

   11
    1

chisq =  1.8953
p =

   11.0000
    1.0000

chisq =  1.8953
p =

   11.0000
    1.0000

chisq =  1.8953
p =

   11.0000
    1.0000

chisq =  1.8953
p =

   11.00000
    0.99999

error: __bfgsmin_gradient: gradient contains NaNs or Inf
error: feval: function `chi2' not found
error: octave_base_value::double_value (): wrong type argument `<unknown
type>'
warning: __bfgsmin_obj: objective function could not be evaluated -
setting to DBL_MAX
error: feval: function `chi2' not found

Thanks,

Bertrand

On Mon, 2009-03-16 at 10:02 -0700, Michael Creel wrote:
>
>
> Bertrand Roessli wrote:
> >
> > Sorry I do not get it,
> >
> > "Just download that file to your working directory,"
> >
> > which file?
> >
> > If you mean from optim-1.04, it is what I did before.
> >
> > Bertrand
> >
> > On Mon, 2009-03-16 at 07:00 -0700, Michael Creel wrote:
> >>
> >> Bertrand Roessli wrote:
> >> >
> >> > Hello,
> >> >
> >> > Thanks for trying the code.
> >> >
> >> > I have the latest package from octave sourceforge (optim 1.04) and
> >> > octave 3.1.54.
> >> >
> >> >
> >>
> >> It seems that the optim package was last released in August of 2008,
> >> while
> >> the last change to __bfgsmin.cc was in October 2008. I just ran your code
> >> 10
> >> times without problems, so I think that the problem is a bug in the old
> >> version of __bfgsmin.cc. Just download that file to your working
> >> directory,
> >> run mkoctfile __bfgsmin.cc, and I think you'll be set.
> >> Michael
> >>
> >> p.s., I have not tested it with Octave 3.1.x,
> >>
> >
> > _______________________________________________
> > Help-octave mailing list
> > [e-mail]
> > https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
> >
> >
>
> You can browse the SVN repository to get individual files. __bfgsmin.cc is
> here:
> http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/optim/src/__bfgsmin.cc?view=log
>
> Just click on that, then download the latest version. You can browse around
> the SNV repo to see the stucture, if you're curious. Because __bfgsmin.cc is
> written in C++, it needs to be compiled before it can be used by Octave.
> That's why you need the mkoctfile step I mentioned before.
> Cheers, Michael
>

_______________________________________________
Help-octave mailing list
[e-mail]
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

I have run the code again several times, and it alway works for me  running Octave 3.0.3 on 64 bit Kubuntu 8.10. However, I just tried it on PelicanHPC 64 bit version (essentially, Debian Lenny, with Octave 3.0.1) and it crashes, even though __bfgsmin is identical on both machines. Sooo, something is going on here.  

The gradient at the initial parameter values has elements that differ by 3 orders of magnitude. That is never a good thing for a gradient-based minimizer. Nevertheless, it should not cause a crash, just a failure to converge at worst. I'll look into it. I'm not much of a C++ programmer, though, so if any experts would like to look at __bfgsmin.cc, help is welcome.

Michael

[Michael Creel]

Michael Creel wrote:

I have run the code again several times, and it alway works for me  running Octave 3.0.3 on 64 bit Kubuntu 8.10. However, I just tried it on PelicanHPC 64 bit version (essentially, Debian Lenny, with Octave 3.0.1) and it crashes, even though __bfgsmin is identical on both machines. Sooo, something is going on here.  

The gradient at the initial parameter values has elements that differ by 3 orders of magnitude. That is never a good thing for a gradient-based minimizer. Nevertheless, it should not cause a crash, just a failure to converge at worst. I'll look into it. I'm not much of a C++ programmer, though, so if any experts would like to look at __bfgsmin.cc, help is welcome.

Michael

Oops (well, actually, good), the problem on PelicanHPC does not exist. Your code runs fine on PelicanHPC (Debian Lenny) using the latest __bfgsmin.cc, as well as on Kubuntu 8.10, Perhaps __bfgsmin does not work properly with Octave 3.1.x. I will not have time to look into that possibility for the foreseeable future. If you are running Octave on Windows, I suggest trying it on Linux.
Michael