Comments and Ratings

Sorry about that. Only the binomfactors code uses consolidator, and only nchoosek calls binomfactors. I never caught it in testing, since consolidator is on my search path.

Anyway, and would be more efficient to use a call to unique there, then accumarray anyway, since no tolerance is employed.

I've submitted a new release. It will be there on Monday morning.

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

I might point out that this is identical to what the already existing function consolidator does. In addition, consolidator allows a tolerance for floating point data.

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

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.

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

Derek - it points out that even though the function mtimes is provided, this package only claims to support ELEMENTWISE operations.

The quadratic behavior tells me that this code did not do a true matrix multiply. Did you check that A64*B64 was correct in your test? Element wise multiplication will take O(N^2) flops, whereas a matrix multiply (as * is supposed to generate) is O(N^3).

I would have preferred that an error be generated if you use * on a pair of matrices, when that operation will not return the proper result. This would normally make me downrate this code. Since I cannot test it without compiling it, I won't give any rating at all.

For schur, I'd like to help. But schurly, then I would have to call the solution the schurshuffle.

http://www.youtube.com/watch?v=0V-VgRqsEcg

Seriously, it seems the same methodology would work. And with a name like the schurshuffle, I'd hate to pass up the opportunity.

Sorry, but I fail to see any value here. This does not teach anything about Simpson's rule. It uses only a 3 point rule, not a composite over multiple intervals, so this will be of no use to anyone who actually needs to do numerical integration.

On top of this, there are several other Simpson's rule submissions on the FEX.

By the way, why use the symbolic toolbox at all?

The code is even misleading. The input statement (a terrible way to build an interface, by the way) tells you that you need to supply f(x) = 0. Sorry, but that has no meaning in the context of a numerical integration. ANY function may be integrated. So either the author misunderstands integration, or they misunderstand what a function is.

This is nothing more than a homework assignment. As such, the only place for it is on your parents refrigerator, and then only if you are actually proud of it.

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

Nifty

Hi Gwendolyn,

I was going to submit a new version today again. The simplest thing to do is to save that file in a form that is compatible with older releases that go back to version 7.

John

John - I think the problem is shown by your path. Are you using an older release of MATLAB? Perhaps release 11? If so, I'll not be at all surprised at seeing a problem. I would not expect vpi to run on older versions, and I have not tested it there.

If you actually have a more recent release, then please contact me via mail and I will help to resolve the problem.

I see some changes have been made. These extra comments really help the person reading the code. My thanks there.

Personally, I'd still hope that beginning PE solvers try the problem first, before looking at this. Only afterwards should you read through this submission to compare this code to yours. But this is a reasonable approach to the problem.

I'll still abstain on an actual rating, but if I were to give this a numeric rating, my rating would be a decent one.

The problem with this code is it is potentially inaccurate. It uses an identity that people are oh so impressed with, but squares the numbers unnecessarily. I'll talk about this problem eventually, but first, let me discuss another major flaw with this code.

It does not actually compute the Euclidean distance. It computes the square of that distance. As such, it has the wrong units. Don't forget that if your data has units miles, then all "distances" produced by this function will have units of miles squared.

Note that this code claims to produce the same numbers as does the pdist function. But it does not. pdist produces TRUE distances, not the square of the distance. There is a difference.

A = randn(3,5);
B = randn(3,4);

distance(A,B)
ans =
       2.6333 7.7299 4.5349 4.6107
       3.1854 4.636 9.1202 1.0953
        2.277 9.2126 3.9594 9.9298
       5.9439 1.416 6.8336 2.748
       3.1344 0.96642 7.4918 1.6271

The true interpoint distances, as one will normally see by any correct tool available to compute distance is closer to this:

sqrt(distance(A,B))
ans =
       1.6228 2.7803 2.1295 2.1472
       1.7848 2.1531 3.02 1.0466
        1.509 3.0352 1.9898 3.1512
        2.438 1.19 2.6141 1.6577
       1.7704 0.98306 2.7371 1.2756

I imagine the author will claim that by not computing the square root, it saves time. I suppose so, but to then claim that the result is a distance as most expect to see would be misleading. And surely to claim that it produces the same as other codes is very misleading.

A nice property for a distance measure to have is linearity. We would like to see that

distance(k*A,k*B) = k*distance(A,B)

It is something that pdist will give you. But this is not a property one will gain from distance.

On to the other problem, the one that I see as most damaging. Recall the results of the above experiment for distance(A,B). Now, try adding any fixed offset to both A and B. I'll add a large one to the matrices to show that complete trash is generated.

distance(A+1e8,B+1e8)
ans =
     8 8 8 8
     4 -4 4 4
     4 4 4 12
     8 0 8 0
     0 0 8 0

See that NO significant digits remain in the result.

A higher quality tool (like pdist for example) will survive even this extreme test quite nicely. As it turns out, this is not an extreme test for pdist. (Since pdist does only interpoint distances for a single matrix, these numbers will not be comparable. See only that adding 1e8 did not change any digits printed out.)

pdist(A'+1e8)
ans =
  Columns 1 through 8
       1.9891 1.8594 3.139 2.8911 2.7152 2.6316 1.5634 3.9224
  Columns 9 through 10
       3.1472 1.7458

pdist(A')
ans =
  Columns 1 through 8
       1.9891 1.8594 3.139 2.8911 2.7152 2.6316 1.5634 3.9224
  Columns 9 through 10
       3.1472 1.7458

I will point out that all of the above mentioned flaws could have been repaired with a few simple changes to the code. Had the code been written as follows, it would take only about 20% longer to execute in my quick test, but it would have worked properly and robustly.

if nargin == 1
    B = A;
    mu = mean(A,2);
else
    mu = (mean(A,2) + mean(B,2))/2;
end
A = bsxfun(@minus,A,mu);
B = bsxfun(@minus,B,mu);
D = full(-2*A'*B);
D = bsxfun(@plus,D,full(sum(B.^2)));
D = sqrt(bsxfun(@plus,D,full(sum(A.^2)')));

Finally, the help for this code is itself wrong. Here is what it says when you do "help distance".

Compute distances between all sample pairs
  X: d x n data matrix
  D: n x n pairwise distance matrix
  Written by Mo Chen (mochen@ie.cuhk.edu.hk). March 2009.

See that it tells you that X is a d by n matrix of data. Does it tell you that the function actually accepts TWO arguments? That if there are two arguments, then it computes distances between all pairs of columns of the two matrices?

No error checks to verify the data actually conforms for distance computation. I'll be generous and give this a 2 star rating.

Benn,

My guess is that you have written your own factorial function, probably as a homework assignment, since the error messages show a recursive implementation of factorial. This causes a failure in derivest, since your version is not properly vectorized. By vectorized, I mean that when called with a list of numbers, the Mathworks version of factorial computes the factorial of each number.

You can learn which factorial version is called by typing this at the command line:

which factorial -all

MATLAB will show which version of the factorial function will be used. You should see the version that you wrote listed at the very top.

At first, I thought that possibly you simply are using an old version of MATLAB, but this is unlikely, since the Mathworks would never have written the factorial function to have a recursive definition as I see. (Recursive factorial definitions are a nice way to teach a student recursive coding styles, but they are a terrible coding style in practice because of the extraneous, wasted function call overhead just to multiply a list of consecutive numbers together.)

The solution is to remove the bad (homework assignment) version of factorial from your search path. It will be inefficient anyway compared to the supplied version of factorial that MATLAB already comes with, and your version is not properly vectorized. At the very least, move it to another place far down on your search path so that matlab will use the correct version of factorial instead.

John

Ok. I had to look here. The problem is, the formulas are poor. Her is a paste of part of the code:

===========================================

% Solve for the regression coefficients using ordinary least-squares
% ------------------------------------------------------------------

   if (INTCPT == 1)
     X = [ ones(n,1) X ] ;
   end
 
    XTXI = inv(X' * X);
    Coefficients = XTXI * X' * Y ;

===========================================

See that the normal equations are used, as well as (shudder) inv.

This is the WRONG way to solve this problem. It is a mistake that many novices to regression make. It is a mistake that many people who have done regression for years make, simply because they have never been told they are doing it the wrong way. It is a mistake that many authors of statistics texts make, simply because nobody has ever told them either.

The normal equations, are

   C = inv(X' * X) * X' * Y

are ONE way to solve the over-determined problem

   X*C = Y

They are a POOR solution. First of all, this formula uses inv. Somewhat better would be to write it as

   C = (X' * X) \ (X' * Y)

The use of backslash here avoids the use of inv. But it still fails because the normal equations are still involved, just in another form.

Look at a better solution.

   C = X \ Y;

This uses the built-in backslash to solve the problem. How would backslash solve it internally? (I'll assume that no pivoting is done to make things simple.) Suppose we choose to form a QR factorization of X.

   X = Q*R

Here, Q is an orthogonal matrix, R is upper triangular. (Again, I'll ignore any pivoting considerations.) Our problem is now to solve

  Q*R*C = Y

Since Q is orthogonal, then we can left multiply by the transpose of Q, to reduce the problem to

  R*C = Q' * Y

Remember, Q was orthogonal, so that (Q' * Q) was an identity matrix.

We now solve for C as

  C = R \ (Q' * Y)

Backslash is smart here. It knows that R is upper triangular, so backslash will do the solution as a back-substitution, an O(n^2) operation. No explicit inverse is required at all.

Again, backslash does all of this internally. But it is a far better solution than the normal equations. What happens in the normal equations to make that solution poor? The problem is the condition number of the matrix. If your matrix is at all poorly conditioned (something that is very common in multiple regression problems) then when you form (X'*X), the condition number of the matrix will be SQUARED. So the normal equations are a terribly poor way to solve the system posed.

We can get the condition number of a matrix from the cond function. It is defined as the ratio of the largest to the smallest singular value of the matrix. When this number is large, then you should expect to loose accuracy in the solution of your equations. Look at a simple, random matrix. (By the way, the condition number will be small here. But most regression problems have a much larger condition number.)

X = rand(10,5);
svd(X)
ans =
          3.70360600537412
          1.16305761722647
         0.819599668934314
          0.77149157109882
         0.359840732709897

See that the largest and smallest singular values would have a ratio of a wee bit over 10.

cond(X)
ans =
          10.2923478881418

What happens when I form X'*X?

cond(X'*X)
ans =
          105.932425050538

Yes. As expected, the condition number was squared. Try this again with a more typical regression problem. Here, just a simple 10th order polynomial in one variable.

x = rand(20,1);
X = bsxfun(@power,x,[10 9 8 7 6 5 4 3 2 1 0]);

cond(X)
ans =
          32218698.6067888

cond(X'*X)
ans =
      1.03682639151296e+15

See that the condition number of (X' * X) is roughly 1e15. In fact, if we use the normal equations to solve this problem, the coefficients will have few correct digits in them, because (X' * X) is nearly a numerically singular matrix.

Now try an 11th order polynomial. Make some random data too. Since the actual numbers don't matter, just see if there is a difference in the solution.

X = bsxfun(@power,x,[11 10 9 8 7 6 5 4 3 2 1 0]);
Y = rand(20,1);

X\Y
ans =
          213864.116016404
         -1151194.24570325
          2686678.25926054
         -3564335.59383269
          2961378.38376311
         -1600304.51887482
          565978.359556758
          -128793.21499217
          18153.6837583098
         -1482.77560777251
          58.7870441102109
        0.0478989431356096

Now compare that to the normal equations. See that inv complains here, because the system is now effectively singular.

inv(X'*X)*X'*Y
Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 8.102118e-18.
ans =
          223829.921190437
         -1205573.39200913
          2815035.41413421
         -3736079.41837729
          3104808.97033064
         -1677906.49043969
          593313.714534699
         -134931.578398688
          18988.6654796978
         -1545.18491246855
           60.905800953616
        0.0268631240254516

See that many of the resulting coefficients do not even agree in a single digit of the result. The problem is when we formed (X' * X).

Well done.

Agustin found his problem - this is always a path problem. Something got put where MATLAB does not look.

Perhaps fig may be available on the internet, but where? Give a link. Is it really that difficult to be friendly to those who might try to download your work? Is it really necessary for others to be forced to do these searches to find out how to use your code?

IF your code uses something that will not be available in MATLAB proper, then tell people about it. If it uses a toolbox function, then show that toolbox as a dependency. If it uses another FEX submission, then tell people about it. Put a link in your code that shows where to find it. Do NOT post code that has dependencies without that information, as this then is non-working code. Non-working code is by definition poor.

Macs will return either of MACI or MAC, depending upon whether the system is an older or new inter Mac system. The help for computer.m states this fact.

So hostname still fails to work on older MAC systems, although it will probably now run on intel based macs.

Also, this submission too has the flaw that it has a blank comment line BEFORE what serves as an H1 line. Since that disables lookfor, it is a flaw, and a silly one. I'm not at all sure why the author insists on doing this. It is unfriendly to disable a very useful feature of MATLAB for absolutely no good reason.

I see that the code has been updated to provide vectorized operation. As always, I should consider raising my rating.

First, I looked for an H1 line. No. Still a blank line as the very first line of the help. This is utterly silly, since the second line of the help is a perfectly serviceable H1 line. This disables lookfor. Interestingly, the help facility that automatically builds a contents file for directories using the first lines of the help blocks finds roundl. But lookfor is itself disabled here. So this problem has not been repaired.

How bout the vectorized question?

x = rand(1,5)*100
x =
    74.555 73.627 56.186 18.419 59.721

roundl(x)
ans =
   100 100 50 20 50

By default, roundl uses a rounding strategy, i.e., a rounding mode of 'n'. But it can round up or down, more like a ceiling or floor operation then. All of these operations work nicely, and a default appears to be supplied for the rounding mode if necessary. Or is it?

roundl(x,'n')
ans =
   100 100 50 20 50

roundl(x,'u')
ans =
   100 100 100 20 100

roundl(x,'d')
ans =
    50 50 50 10 50

Now try the code by adding the third argument. The third argument to roundl allows you to control the levels rounded to.

roundl(x,'n',7)
ans =
    80 80 60 20 60

No problem so far. But what happens if you don't wish to change the default for the rounding mode, only the default for n? The paradigm in MATLAB when you have variables that will take on a default is if the variable is left empty ([]), then use the default. Does this work? No. This fails. I tried these things:

roundl(x,[],7)
??? Error using ==> roundl at 120
Invalid rounding mode:

roundl(x,'',7)
??? Error using ==> roundl at 120
Invalid rounding mode:

I even tried this:

roundl(x,7)
??? Error using ==> roundl at 120
Invalid rounding mode: 

As I said, this fails to follow the standard MATLAB paradigm for default arguments, which makes the code not very friendly. If you wish to change the rounding level, you must also provide a default for the other parameter. It is a truly silly flaw, because is it really that difficult to write your code as:

if (nargin < 2) || isempty(n), mode= 'n'; end

versus the line as it was written?

if nargin < 2, mode= 'n'; end

Sorry, but it is easy to write friendly code, and it is easy to make the code compatible with lookfor. Given this newly found problem, a 3 rating still applies here.

I do appreciate the review by Dr. Feldman. However, it shows a few basic misunderstandings of MATLAB. They are logical mistakes, but they reflect the point of view of someone who still thinks they are using perhaps C or JAVA or FORTRAN or PYTHON, pick your favorite language other than MATLAB.

To start with, he suggests that this code should use integers, to make it faster. The first problem (a truly fundamental one) with that suggestion is that support for operations between integers is limited in matlab. For example, in R7007b of MATLAB (not all users use the latest release, so it makes sense to do this test in a slightly older release anyway)

>> N = int64(101);
>> P = int64(23);
>> mod(N,P)
??? Undefined function or method 'mod' for input arguments of type 'int64'.

>> N = int64(101);
>> P = int64(23);
>> N*P
??? Undefined function or method 'mtimes' for input arguments of type 'int64'.

See that even the most trivial operations are not supported for int64. Likewise, uint64 is also lacking in support for these operations in this release.

>> P = uint64(23);
>> N = uint64(101);
>> mod(N,P)
??? Undefined function or method 'mod' for input arguments of type 'uint64'.

So it might be a good idea to suggest use of integers. That might possibly have been a good suggestion were we writing in one of those other languages rather than MATLAB. But we are using MATLAB here. Dr. Feldman's suggestion is a practical impossibility when arithmetic for that variable type does not exist. While I admit that that support is present for the INT32 variable type:

>> N = int32(101);
>> P = int32(23);
>> N*P
ans =
        2323
>> mod(N,P)
ans =
           9
 
Does it even speed things up? Use of int32 instead of double would be a tremendously silly thing to do, since int32 operations are not even faster than operations on doubles! In fact, use of int32 would slow down this code!

>> N = int32(1:10000);
>> P = int32(23);
>> timeit(@() mod(N,P))
ans =
    0.0036231

>> N = double(1:10000);
>> P = double(23);
>> timeit(@() mod(N,P))
ans =
    0.0026903

In fact, arithmetic using doubles is FASTER than int32 arithmetic. int64 or uint64 arithmetic is not even defined in the target release, so the speed is incomparable there. ;-)

NEXTPRIME does work for my vpi class, in fact, it was originally designed for those tools. However, they are slower to use than the same direct operations on double variables. So there is no reason to force all operations into the vpi class. Worse, this would restrict the utility of this function, forcing the user to have my vpi toolbox installed even to do something as trivial as finding the next prime number above 20. The reason I posted this submission separately from my vpi toolbox is because it has a broader range of utility than for just those people who would work with truly large integers.

Finally, the author suggests that NEXTPRIME should support doubles beyond 2^32. The main problem here is with MATLAB again.

>> isprime(2^32+1)
??? Error using ==> isprime at 22
The maximum value of X allowed is 2^32.

The isprime function does not work beyond 2^32, and nextprime uses isprime to make the final determination of whether a candidate number is prime, having weeded out most alternatives with a partial number sieve first. So this operation will fail. One solution here is to convert to my vpi class if the user supplies a number as large or larger than 2^32, testing first to see if the vpi tools are installed. A different solution would be to replace the call to isprime with a custom isprime call. I'll consider these alternatives as options.

Despite the claim, this fails to run on all operating systems. Try it on a PPC Mac. No go.

unix('hostname') works fairly well if you are on a PPC mac. I've not verified it is correct for an intel mac.