A serious problem with this is what I would suggest looks like almost two functions merged into one. Worse, the two "curves" at the bottom end are not consistent in either case with what we see at the top end, so schizophrenic data. And worse yet, you have a least two large outliers. Large residual nonlinear least square problems are difficult to solve, because one data point can really swings things around. It has an effectively huge residual. What else? Oh yeah, your data varies by an order of magnitude. The result of that is the points at one end of the curve impart effectively more influence on the result than they should.
First what happens to your model if we log y? This actually solves one of the issues I mention above, since we no longer have data that varies by powers of 10 when the fit is performed. Effectively, I am suggesting we think of your data in a semi-logy form.
semilogy(x,y,'o')
grid on
xlabel X
ylabel Y
The strange behavior of your data should be clear here. Thus we see the outliers (be careful not to miss the outlier on the very bottom.) We see data that on the right half seems to be nearly a perfectly straight line. The left half we see two pieces that clearly represent two completely different and inconsistent processes. They are also inconsistent with the data on the right. Two different things are happening there, neither of which is consistent with the right half of your data.
What will happen when we try to estimate a model, is your curve will not represent all three of those behaviors.
Before I get into a fit, first, log your model. That is, when we form the semilogy plot, we are changing your model from this:
y = a*exp(-b/(x^p))
into this:
log(y) = log(a) - b*x^(-p)
Dividing by x^p is just a multiply by x^(-p). And log(a) is still just a constant that we don't know. Effectively, we can think of the model in a logged form as almost a polynomial. Most seem to call that a power-law model. Again though, these curves really do seem to be purely linear when we look at the semi-logy plot. That implies p==-1 is about as close as my eye wants to fit this. Remember, you have that -p in there, so two negatives do make a positive in this case. As such, the proper model for your curve prbably has p as very close to -1.
Time to do a first cut at a model. I'll just use the basic fitting tools from the plot figure itself to create that model.
plot(x,log(y),'o')
grid on
On the right half, again, that data really is nearly perfectly linear, implying p==-1. But on the left half, which of those curves should it follow? A curve fit cannot follow both curves. And if it chases one, then it misses the other by much more. And worse, if you do go chase one of those lower personalities, then the curve will miss the upper part by more yet.
The funny thing is, those two nasty looking outlier points almost seem to be the least of your problems. Even so, I would still delete them from the analysis before proceeding further. They are still influencing your estimates.
What would I do?
- Delete the two outliers.
- Determine why it is your data is so schizophrenic. What happens that drives the two behaviors? Why is there a split at roughly x~150? Choose the process you wish to accept as truth. Only when you excise the split personalities in your data can you build a model that works.
- Finally, either split the problem into two pieces, performing a broken-stick regression with the break between the two halves of your curve, or try to build once global model that won't fit either half very well.
So first, in your data set (ignoring the NaNs) we see
Right half of the curve: points # 10:47
The two outliers: # 1, 49
Left curve (lower branch: # 2:9
Left curve, upper bransh: # 48:60
So if I had to guess, you are doing something strange to generate the data, that creates different behaviors that are inconsistent with each other. A model is not going to bring it all together in some magic way.
I might split the data like this:
ind1 = 10:47;
ind2 = 2:9;
ind3 = 48:60;
ind4 = [1 49];
If we just consider the piece at the right, we would see:
p1 = polyfit(x(ind1),log(y(ind1)),1)
p1 =
0.0185627813572888 -27.2729761320311
Nothing more sophisticated than polyfit was needed. So a very good model for that curve segment will be simply
y = exp(-27.2729761320311) * exp(0.0185627813572888 * x)
Note the fact that x in that model did not even require a negative exponent.
Plot the data now, with this model. As you will see, the fit really is very good. I'll even view it in an unlogged form, to show I was sucessful in posing a reasonable model for the entire process, even though the model I would pose is far simpler than what you have intended.
plot(x(ind1),y(ind1),'bo',x(ind2),y(ind2),'go',x(ind3),y(ind3),'ko',x(ind4),y(ind4),'r*')
grid on
hold on
xhat = 40:230;
yhat = exp(p1(2))*exp(p1(1)*xhat);
plot(xhat,yhat,'b-')
In fact, the model fits splendidly well for the upper half. There simply is no good way to resolve the split personality in the lower half though. At least not so until you understand why the split exists, and how to properly deal with it. My suggestion at this point would be to get better data.
Finally, can I reconcile your claim that p should be on the order of 0.25 to 0.5? That seems to be wildly inconsistent with the data you have provided.
