Leasqr vs. Weibull

help
,
#1

Dear all, I am having a reoccuring issue with leasqr, I am trying to interpolate data with a weibull function [something like y = e^((a*t)^k)].

This is my data

A =
0.00000 0.00000 0.00000
22.50000 0.90672 1.10453
47.00000 1.02958 1.26305
72.16667 1.11205 1.46428
95.50000 1.18885 1.57736
166.50000 1.34828 1.88047
192.25000 1.39940 1.94198
216.00000 1.45110 2.01961
262.50000 1.51369 2.14295
334.50000 1.60073 2.32114
385.00000 1.65601 2.42147
431.00000 1.69672 2.50011
504.00000 1.74113 2.58324
552.00000 1.90475 2.65936
603.00000 1.79520 2.71121
674.00000 1.84116 2.79580
720.00000 1.86893 2.66708
863.00000 1.91990 2.93466
934.66666 1.94194 2.98087
1008.00000 1.96009 3.01555

And this is my code

A = dlmread(“data.txt”);
##A = A([2:end],:);
xi = A(:,1);
yi = A(:,2);
qrfunc = @(xi, p) yi(end) - yi(end)*exp(p(1).*xi).^p(2);
p = [1, 0.4];
wt = linspace(1,1,length(yi))’;
dp = [.001; .001]; % bidirectional numeric gradient stepsize
dFdp = “dfdp”; % function for gradient (numerical)
[f, pint, kvg, iter, corp, covp, covr, stdresid, Z, r] = leasqr(xi, yi, p, qrfunc, 0.001, 200, wt, dp, dFdp, [-Inf, Inf; 0.74, 0.76]);

There are a few things that puzzle me.

First, I am getting the error

error: svd: cannot take SVD of matrix containing Inf or NaN values
error: called from
lm_svd at line 207 column 17
leasqr at line 660 column 25
data at line 13 column 56

I do not understand. There are no Inf or NaN values, I thought? Ok, I can work around by scaling xi and yi such that their maximum is 1 (backscale after the interpolation):

a = max(xi);
b = max(yi);
xi = xi/a;
yi = yi/b;

and

xinterp = linspace(xi(1), xi(end),40)’;
yinterp = yi(end) - yi(end)exp(pint(1).xinterp).^pint(2);
plot(axi, byi,‘bo’, axinterp, byinterp)
legend(‘Data’, ‘Interpolation’, ‘location’, ‘southeast’)

Then the code runs through, but the interpolation is unsatisfying. Why does the code require the maximum of xi and yi to be 1 (or any other number perhaps)? And why is the interpolation not so good? I mean, the chosen function can interpolate the data much better in terms of smaller r.

The second thing I noticed is that my bounds are not respected.

I would appreciate any input with regards to the “scaling requirement” and what can be changed such that the bounds are respected.

Please let me know if any information is missing.

data.m (682 Bytes)
data.txt (578 Bytes)

#2

Hello JSP

your definition of the function might not be what you really want

Try
A = dlmread(“data.txt”);
A = A([2:end],:);
xi = A(:,1);
yi = A(:,2);

qrfunc = @(xi, p) yi(end) - yi(end)*exp(p(1)*xi.^p(2));
p = [1, 0.4];
wt = ones(size(yi));
dp = [.001; .001]; % bidirectional numeric gradient stepsize
dFdp = “dfdp”; % function for gradient (numerical)

[f, pint, kvg, iter, corp, covp, covr, stdresid, Z, r] = …
leasqr(xi, yi, p, qrfunc, 0.000001, 200, wt, dp, dFdp, [-Inf, Inf; 0.74, 0.76]);

xinterp = linspace(xi(1), xi(end),40)’;
yinterp = qrfunc(xinterp,pint);

figure(1)
plot(xi, yi,‘bo’, xinterp, yinterp)
legend(‘Data’, ‘Interpolation’, ‘location’, ‘southeast’)

I hope it helps

Andreas

2 Likes
#3

Dear Andreas,

Thanks for your input!!! Unfortunately, changing the dot operator and the “)” from

qrfunc = @(xi, p) yi(end) - yi(end)*exp(p(1).*xi).^p(2);

to

qrfunc = @(xi, p) yi(end) - yi(end)*exp(p(1)*xi.^p(2));

did not improve the interpolation quality, but the code runs through - no scaling requires. This is already good progress!!!

However, the boundaries are still not respected - I am not saying that the boundaries will lead to a “good” interpolation, I would just like to see how this works.

#4

Dear JPS

With leasqr you are searching for a nonlinear least square approximation of a function
y(end)(1-exp(p1x^p2)) to your given data. The value of kvg = 1 indicates that leasqr() converged, i.e. it found the best possible solution.

I do not see your problem with the interpolation quality?
What do you mean by “quality”?

What do you mean by “boundaries respected”?
The algorithm leasqr() can not deal with constraints, it is an unconstraint least square
algorithm, based on Levenberg–Marquardt.

Your first approch used exp(p1x)^p2, with equals exp(p2p1*x), thus you did not have two independnet parameters. This lead to the NaN and Inf results. The approach had to fail.

Hopefully this helps

Andreas

1 Like
#5

For those exponential type of fits, it is important to know (or, at least, have a good idea) of the saturated
value at X = infinity. Here you assume it is yi(end), which is obviously not correct. If you change yi(end) in your function to something like 2.7 you will get a better looking fit.

1 Like
#6

Hello again
dasergatskov’s remark is spot on.
Try

A = dlmread("data.txt");  A = A([2:end],:);
xi = A(:,1);   yi = A(:,2);

qrfunc = @(xi, p) p(3)*(1-exp(p(1)*xi.^p(2)));
p = [0, 0.5,yi(end)];
pkg load optim
[f, pint, kvg, iter, corp, covp, covr, stdresid, Z, r] = ...
   leasqr(xi, yi, p, qrfunc);

xinterp = linspace(xi(1), xi(end),40)';
yinterp = qrfunc(xinterp,pint);

figure(1)
plot(xi, yi,'bo', xinterp, yinterp)
legend('Data', 'Interpolation', 'location', 'southeast')
1 Like
#7

Dear Andreas,

Thanks, that clarifies the situation a bit! I got confused with the exponent to e, thanks for the clarification. Still, I am irritated. If leasqr cannot deal with constraints, why does it offer the possibility to do so? Ok but maybe this goes a bit far.

#8

Oh yes!!! Thank you both for your input, the underlying issue is solved now. I earlier hat a scaling parameter (just as Andreas with p(3)) but had dropped it during my error search.

I only struggle to understand the constraint topic, it would be nice to learn how the contradiction between “leasqr can’t work with constraints” and “leasqr offers constraints” can be resolved.

#9

Hello again

You are correct in that the current version of leasqr() in the Octave package optim can deal with constraints on the parameters, thanks to Olaf Till. I had the old version in my mind.
The constraints have to be used by setting an option, see the code below.
Hopefully this helps

Andreas

A = dlmread("data.txt"); A = A([2:end],:);
xi = A(:,1); yi = A(:,2);

qrfunc = @(xi, p) p(3)*(1 -exp(p(1)*xi.^p(2)));
p = [0, 0.4,yi(end)]; wt = ones(size(yi));
dp = [0.001; 0.001; 0.001]; % bidirectional numeric gradient stepsize

[f, pint, kvg, iter, corp, covp, covr, stdresid, Z, r] =  leasqr(xi, yi, p, qrfunc);
pint

xinterp = linspace(xi(1), xi(end),40)';  yinterp = qrfunc(xinterp,pint);
figure(1)
plot(xi, yi,'bo', xinterp, yinterp)
legend('Data', 'Interpolation', 'location', 'southeast')

%% now try with a constraint on p(2)
options.bounds = [-Inf, Inf; 0, 0.1;-Inf, Inf];
pint(2) = 0;

[f, pint, kvg, iter, corp, covp, covr, stdresid, Z, r] = ...
   leasqr(xi, yi, p, qrfunc, 0.000001, 200, wt, dp, "dfdp", options);
pint
yinterp = qrfunc(xinterp,pint);
figure(2)
plot(xi, yi,'bo', xinterp, yinterp)
legend('Data', 'Interpolation', 'location', 'southeast')
1 Like
#10

Hello Andreas,

Yes, this works well now. Thanks for the example!

One last question: I understand that “dfdp” is a representation of qrfunc’s derivative. Is it in the form of a function which can conceptually be evaluated similar to

yinterp = qrfunc(xinterp,pint);

E. g.

yinterp_der = dfdp(xinterp, pint);

or how can I obtain the values so that I can print the derivative? Calling dfdp yields an error (“func” is undefined, called from dfdp).

Thanks for your time and input
JPS

#11

Hello JPS

the function dfdp() uses a finite difference method to aproximate the partial derivatives of the given function (your qrfunc()). leasqr() will internally call dfdp() with the correct arguments. Using dfdp() directly requires a few more arguments.
Try
pkg load optim; help dfdp
to acces some information.
For leasqr() you can provide a function to return the partial derivatives and then
leasqr() will use your derivatives and not call dfdp().

Andreas

1 Like
#12

Hi Andreas,

When I entered

help dfdp

I found another derivative function, which is supposedly “more general”:

function jac = dfxpdp (xinterp, pint, qrfunc);

The above did not yield results, the prompt looked like an infinite loop or very very long/slow calulation. So slow that either it does not work this way or I interrupted too early.
However, I thought the derivative is part of leasqr’s optional arguments as in

[f, pint, kvg, iter, corp, covp, covr, stdresid, Z, r] = …
leasqr(xi, yi, p, qrfunc, 0.000001, 200, wt, dp, “dfdp”, options);

I was looking for a simple way to access the “dfdp” information of the leasqr command. I understand the approach you suggested and which I tried (unsuccessfully, so far) is a completely new computation?

#13

I succeeded using dfdp
see below

pkg load optim
ff = @(p,x)sin(px)
x0 = pi/4; p0 = 1;
f = ff(p0,x0)
jacobian = dfdp(x0,f,p0,[1e-6],ff)
results_analytic = cos(p0x0)*x0

1 Like
#14

Hi Andreas,

thanks for your above code, it runs. However, I struggle to transfer it to my problem, as I am getting confused with the description given in the octave function reference. Which of the x’s must be a scalar, which an array?

ff = @ (p,x) sin (p*x)

here, x should be an array? What is p for, intended as scaling parameter? Required by the structure of dfdp and hence set to 1?

x0 = pi/4; p0 = 1;

both are scalars, is this some kind of initial condition? Or some arbitrary value for testing? Should these be arrays as in y = sin(x) with x being an array and hence is y?

f = ff(p0,x0)

is this the function ff, evaluated at the initial condition or at the arrays as in y = sin(x)? Is it then not identical to f?

jacobian = dfdp(x0,f,p0,[1e-6],ff)

x0 I take is the independent variable, an array. f is an array which requires initialisation. Ok, but is this not identical to p0 (dependent variable?) then? [1e-6] is the fractional increment, ff is the first function above?

I find this rather confusing, which is why I was looking for the “dfdp” functionality of leasqr in the first place.

So far I have consolidated a code snippen which throws an error (fixed(2): out of bound 1):

%qrfunc = @(xinterp, pint) pint(3)*(1-exp(pint(1)xinterp.^pint(2))); %function from far above
ff = @(xinterp, pint) pint(3)(1-exp(pint(1)*xinterp.^pint(2)));
x0 = xinterp; p0=pint;
f = qrfunc5(xinterp,pint);
jacobian = dfdp(xinterp,f,yinterp,[1e-6],qrfunc)