Frustrating. You mention a constrained least squares in your question. But no place did you EVER suggest what the constraint was. I had to dig way down into the comments to find anything.
Your model is :
Y = a*X.^b
The simple answer is to first plot EVERYTHING. And that model implies an important suggestion, if the model is even remotely valid. I.e., you probably want to log the model. That is, if we log things to get...
log(Y) = log(a) + log(X)*b
then we should see a roughly linear relationship, between log(X) and log(Y).
loglog(X,Y,'o')
grid on
Hmm. That does suggest a roughly linear model could be appropriate in log space. Of course if we try that, we will see:
BA = polyfit(log(X),log(Y),1)
BA =
-0.443558076798707 3.77850311606457
Admittedly, this is not a nonlinear regression, but it should give us some feeling for what is happening. And, we would need to exponentiate the constant term to get back to your original parameters. That would yield roughly 44.
So then I started reading through your comments. As well, I plotted the data, and the model you posed as "correct". The model did not come close to even passing through the data. Finally, I saw the only comment where you stated the blasted objective of the problem!!!!!!!!!!
You want to solve a least squares problem, subject to the constraint that Yhat(i) >= a*X(i)^b. So your constraint is non-negativity of the residuals. Could you maybe have told us that in the first place?
The simplest way to solve this problem in a nonlinear least squares sense is to use fmincon. But we can use the log trick from above to convert the problem to a linear one, then use lsqlin. (Both of these tools are in the optimization toolbox. I think you do not have that toolbox from another of your comments.) It does help us if you tell us information like that up front.
N = numel(X);
Amat = [log(X'),ones(N,1)];
BA = lsqlin(Amat,log(Y'),Amat,log(Y'))
Minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.
<stopping criteria details>
BA =
-0.381649116037055
2.70805020110219
exp(BA(2))
ans =
14.9999999999997
So, your model would be estimated as:
Yhat = 14.9999999999997 * X.^(-0.381649116037055)
Again, that uses a linear regression, and it uses a toolfrom the optimization toolbox, and you don't have that, so let me find a reasonable compromise.
We can always use a penalty constrained fminsearch. So, something like this, as an objective function:
function R = fmsconobj(ab,X,Y)
R = Y - ab(1) * X.^ab(2);
if any(R < 0)
R = inf;
else
R = norm(R);
end
Now, call fminsearch:
ab = fminsearch(@(ab) fmsconobj(ab,X,Y),[10, -.5])
ab =
14.9999818205206 -0.381653090395052
That does what you want, but you do need to be careful. Start from the wrong place, and fminsearch will fail.
