Comments and Ratings

John,

The demo_vpi.m needs your consolidator function:

http://www.mathworks.com/matlabcentral/fileexchange/8354-consolidator

Great package.

Derek O'Connor

As I explained the first time, you need to put polyfitn into a directory that is on your search path.

help addpath
help savepath
help pathtool

Your particular test problem gives me the results:

uv = rand(100,2);
w = sin(sum(uv,2));
p = polyfitn(uv,w,'u, v');

polyn2sympoly(p)
ans =
    0.76152*u + 0.73572*v

Nothing special needs be done except putting it in a directory where matlab can find it.

If you are trying to place downloaded files into the Matlab/toolbox directories, this is generally wrong. Those directories are cached by matlab when it starts up, so any files that you put in those directories will not be seen until you restart matlab. Also, anytime you make any changes to files that you put in those directories, the changes will not be recognized by matlab.

Better is to create your own new (and separate) directory for downloaded files. Add that directory on your search path, using the pathtool function, or using addpath and then savepath.

James,

You don't need every possible combination of (x,y,z) to be known to estimate a model, certainly not a first order model like this. To estimate the 3 unknown parameters in this model, in theory, you need at least 3 data points at a minimum. In practice, you won't get very good estimates of those parameters with only 3 data points.

In the past, I often had clients ask me how much data they needed to obtain. After all, data is often scarce, expensive, hard to obtain. It costs money, and time. You can't have too much data, at least if it is good data.

In practice, if you end up with 10 or 20 points to estimate these three parameters, that will (probably) be entirely adequate. More data will give you higher quality estimates. Even 6 points, carefully chosen, if the noise is low enough, will often be adequate. Only you know the quality of your data, and how much you have, as well as how accurately you need the model to be estimated.

If you have further questions, I can only suggest you send me an e-mail, with your data, and any other information you have on the problem.

John

Very nice John, I will cite you in our papers.

I also got this warning:
Warning: Options LargeScale = 'off' and Algorithm = 'trust-region-reflective' conflict.
Ignoring Algorithm and running active-set method. To run trust-region-reflective, set
LargeScale = 'on'. To run active-set without this warning, use Algorithm = 'active-set'.

and added this line after line numer 158 in your slmengine.m file;
fminconoptions.Algorithm = 'active-set';

now I don't get the warning anymore.

The warning message that fmincon returns is new, due apparently to a change in the optimization toolbox. It is only a warning though, that does not hurt the operation of the code itself.

I'll fix the problem.

Very nice and thorough compilation of tips and tricks.

N-I-C-E.
I cannot understand why something similar is not implemented in matlab yet.
It would be nice if one could add several selections by pressing e. g. "Ctrl".

Dear John,

I wonder if it is possible with your functions to fit a spline to a set of (x,y) coordinates over time (each data point also has an associated weight). I want to constrain the velocity of the fitted spline trajectory, which means I can't just fit (x,t) and (y,t) separately. If this is not possible with your SLM tools, do you have any suggestions how to approach this problem?

Thanks, Fabian

Simplest is to replace the nans by interpolation. Use my inpaint_nans (also on the file exchange) to do this. Then apply this code to compute the moving standard deviation. High variance sections of the curve will still be high variance after this process. The virtue of this scheme is that both parts are fully vectorized codes.

The alternative is to use a loop, computing a windowed standard deviation on varying numbers of points by excluding the nan elements. This scheme may be tripped up if you have enough nan elements. Of course, with very many nan elements, no scheme will be perfect.

Hi John,

As pointed out by Eric, I guess, for newer versions, you need an update.

Shaun

function stop = optimplot(x, optimValues, state)
% plots the current point of a 2-d otimization
stop = false;
hold on;
plot(x(1),x(2),'.');
drawnow

Ryan,

This is really not a bound constrained problem. Your constraint is a general one, either a linear or a nonlinear inequality constraint.

The simple solution is to use a code that allows explicit linear or nonlinear constraints. There are a few on the FEX. In fact my fminsearchcon on the FEX does that, by applying a penalty to the problem when a constraint is violated. You might also look at optimize, by Rody Oldenhuis. This code allows explicit constraints of all forms.

However, there is a simple way to solve this class of problem with fminsearchbnd. Use another class of transformation. For example, suppose one wishes to minimize the function f(x,y), subject to the "bound" constraint that x <= erf(y).

Transform your problem such that

u = x + erf(y)
v = x - erf(y)

Clearly, v must be bounded above by 0. So use fminsearchbnd to optimize over the two dimensional domain (u,v). Inside your objective function, for any value of (u,v) you will compute the parameters (x,y) as

x = (u + v)/2;
y = erfinv(u - x);

Now you can finish evaluating the function f(x,y).

The only problem here happens if you also had other fixed bounds on x and y. But many other transformations would also have worked. For example, this transformation:

u = x
v = x - erf(y)

will still allow simple bound constraints on the parameter x, as well as allowing the nonlinear constraint x <= erf(y) as a bound constraint.

John

Pete - A goal to be finished as soon as I can is to write a gui that will wrap SLM inside. Note that the computational tool in SLM is called SLMENGINE. This was purposeful, since my expectation is that most users would use the gui form when I am able to offer it. So I have definitely been planning for a gui wrapper.

Even at that though, SLM would not have allowed you to just draw a curve through your data freehand, and then return the coefficients of that curve. A freehand drawn curve may be arbitrarily complex, or not even a single-valued function. SLM is best when you fit the curve, then apply your own special knowledge of a system to be modeled, in the form of constraints on the overall shape of the underlying functional form. But a drawn curve to follow is too broad of a constraint to follow. So SLM might not be capable in general of returning such a functional form in that eventuality.

As well, this becomes a unique problem. Is the problem then to find the curve that fits the original data, and has the general shape of the freehand drawn curve? This involves two separate problems of approximation. It seems one must use a tool to smooth and approximate the drawn curve. Then take that same tool and try to fit the drawn curve through the data? This task would take some serious amount of effort, and would never be possible to do so where it would be transparent to the user. Worse, suppose the curve in the end had some lack of fit to the data (as perceived by the user?) Where did the lack of fit arise? Is it lack of fit to the freehand drawn curve? Or is it lack of fit to the data itself?

What is the advantage of using gridfit instead of griddata followed by smoothing?

I thought so, that it couldn't be any harder than eigenvalues, but then I looked here
http://www.math.uu.nl/publications/preprints/1180.ps

Actually the code there seems to work pretty well. So schurly, no need to reinvent schurshuffle

Is it possible to write the same functionality but for schur decomposition? Schur function in matlab returns the eigenvalues on diagonal in same order as eig function..which is not always nice. Thanks.

Hallo John,
Nice work!
I am a newbie in matlab and still got some questions about 3d fit.
I have a 3d dataset F(128,128,16) and would like to get the polynomial function F(x,y,z), how should I do it?

thanks!

Hi Michelle,

The problem is, when you define myfun as you have, it is not a scalar valued function. It returns a vector argument. Thus we can do this:

myfun=inline('[exp(x);exp(2*x)]');

myfun(2)
ans =
       7.3891
       54.598

See that there are two numbers returned for a scalar input.

(By the way, I'd strongly recommend the use of function handles and anonymous functions instead of inline functions. An inline function is comparatively very slow to evaluate.)

The derivest code computes the derivative of a scalar valued function. Change your function to return a scalar value, as do myfun1 and myfun2 below, and derivest works as designed.

myfun1 = inline('exp(x)');
[d,err]=derivest(myfun1,1,'Vectorized','no')

d =
       2.7183
err =
   9.3127e-15

myfun2 = inline('exp(2*x)');
[d,err]=derivest(myfun2,1,'Vectorized','no')
d =
       14.778
err =
   1.8401e-13

Now, how might one define the derivative of a vector valued function? This is defined mathematically as the Jacobian. In fact, in the derivest suite, I have provided the jacobianest function to compute what you wish. It works on a vector valued function, returning a jacobian matrix. Here, since the parameter of these functions is a scalar, the Jacobian matrix is just a vector of length 2.

[d,err] = jacobianest(myfun,1)
d =
       2.7183
       14.778
err =
   3.6841e-14
   1.4736e-13
 
John