Experiment and Model not fitting well

filename = 'dotun2.csv';

%delimiterIn = ',';

headerlinesIn = 1;

A = importdata(filename,delimiterIn,headerlinesIn);

ndata = A.data(:,2);

ydata = A.data(:,1);

plot(ydata,ndata,'o')

xlabel 'n'

ylabel 'b'

f = optimvar('f');

g = optimvar('g');

fun = (ndata*18)/(0.1*2.8*10^7)+g*18.*exp(f*log(ndata/(0.1*284)));

obj = sum((fun - ydata).^2);

lsqproblem = optimproblem("Objective",obj);

x0.f = 0.45;

x0.g = 0.184;

show(lsqproblem)

[sol,fval] = solve(lsqproblem,x0);

disp(sol)

figure

responsedata = evaluate(fun,sol);

plot(ydata,ndata,'bo',responsedata,ndata,'o')

plot(ydata,ndata,responsedata,ndata)

%text(0.1,20, '0.1in 100C')

set(gca,'FontSize',30)

legend('Expt','Model', 'text','FontSize', 18)

xlabel ('\it d','FontSize',16,'FontWeight','bold')

ylabel ('\it F','FontSize',16,'FontWeight','bold')

set(gca,'FontSize',20)

set(gca, 'FontName', 'Cambria')

%title("Fitted Response")

axis([0 0.16 0 300])

[Torsten]

Are you sure that your "fun" is adequate ? Shouldn't it asymtotically approach a fixed value as d grows ? Your "fun" grows without limit.

Answer 1

[John D'Errico]

We don't actually have the data. so it is difficult to be more clear about this. But you seem to think what the model must be able to fit better. Let me offer a simple analogy.

Suppose you asked to model a cicle, using a square. Yes, the two shapes have similarities. They are both closed. But the similarities end almost there. You need to see that not all models are appropriate for all problems.

The model you show here is:

fun = (ndata*18)/(0.1*2.8*10^7)+g*18.*exp(f*log(ndata/(0.1*284)));

We can throw away most of those constants as being irrelevant to the basic shape of the curve. I'll write it more cleanly like this, so that we can look at the basic shape of this model.

syms x g f k1 k2
y(x) = x*k1 + g*exp(f*log(x*k2))

y(x) = 

Effectively, your model has a linear term in it, where the parameter k1 (and also k2) is FIXED by you in advance. Essentially, that imposes a linear trend on the curve.

The second term LOOKS like an exponential thing, but it still lacks the correct shape to fit your data. We can write it in a simpler form if we remember two identities.

exp(a*b) = exp(b)^a

exp(log(a)) = a

that means we can write the model in a far simpler way, as

y(x) = x*k1 + g*(x*k2)^f

y(x) = 

Again, k1 and k2 are known variables, and g and f are to be estimated. So your model is just a polynomial like thing, partly linear with a fixed slope, partly just a basic power-law model. It still will not have the shape of your data.

Honestly, you are trying to fit a square peg in a round hole.

[Adeb]

Thank you so much for your detailed response, I have in time expressed the model just as in the pure linear and exponentail for you expressed it in the last equation. Here is the data to provide more clarity.

b n

2.17E-05 3.3722

4.33E-05 6.7443

7.58E-05 11.803

0.000119 18.547

0.000163 25.291

0.000206 32.035

0.000249 38.78

0.000293 45.524

0.000336 52.268

0.000379 59.013

0.000426 65.757

0.00048 72.501

0.000537 79.246

0.000597 85.99

0.000661 92.734

0.000728 99.478

0.0008 106.22

0.000878 112.97

0.000961 119.71

0.00105 126.46

0.001147 133.2

0.001253 139.94

0.001368 146.69

0.001494 153.43

0.001634 160.18

0.00179 166.92

0.001964 173.67

0.00216 180.41

0.002383 187.15

0.002638 193.9

0.002934 200.64

0.003282 207.39

0.003696 214.13

0.004197 220.88

0.004815 227.62

0.005598 234.36

0.006622 241.11

0.008016 247.85

0.010025 254.6

0.013171 261.34

0.018794 268.09

0.03166 274.83

0.091768 281.57

0.15281 288.32

0.15281 288.32

0.15281 288.32

0.15281 288.32

[John D'Errico]

Perhaps the answer is to just play with your data. My rule number 1 is to plot EVERYTHING. And then, when you get bored and can't think of what else to do, plot something else. Every once in a while something makes sense, something works. Can't hurt, right?

bn = [ 2.17E-05  3.3722 
    4.33E-05  6.7443 
    7.58E-05  11.803 
    0.000119  18.547 
    0.000163  25.291 
    0.000206  32.035 
    0.000249  38.78 
    0.000293  45.524 
    0.000336  52.268 
    0.000379  59.013 
    0.000426  65.757 
    0.00048  72.501 
    0.000537  79.246 
    0.000597  85.99 
    0.000661  92.734 
    0.000728  99.478 
    0.0008  106.22 
    0.000878  112.97 
    0.000961  119.71 
    0.00105  126.46 
    0.001147  133.2 
    0.001253  139.94 
    0.001368  146.69 
    0.001494  153.43 
    0.001634  160.18 
    0.00179  166.92 
    0.001964  173.67 
    0.00216  180.41 
    0.002383  187.15 
    0.002638  193.9 
    0.002934  200.64 
    0.003282  207.39 
    0.003696  214.13 
    0.004197  220.88 
    0.004815  227.62 
    0.005598  234.36 
    0.006622  241.11 
    0.008016  247.85 
    0.010025  254.6 
    0.013171  261.34 
    0.018794  268.09 
    0.03166  274.83 
    0.091768  281.57 
    0.15281  288.32 
    0.15281  288.32 
    0.15281  288.32 
    0.15281  288.32];
b = bn(:,1);
n = bn(:,2);
plot(b,n,'o')

Yep. Not very useful there. One simple rule to consider is that when you have data that spans multiple orders of magnitude, that you need to consider a log transformation. In this case, both variables seems to be like that. So we might try this:

semilogx(b,n,'o')

Or this. See if something pops out.

loglog(b,n,'o')

In either case, it suggests a potentially viable model. For example, I could easily consider a model for the first curve using an erf, or a logit shape to fit the semi-log x curve. Try this, for example:

mdl = fittype('a1 + a2*erf((x - a3)/a4)','indep','x');
abcd0 = [150 150 mean(log(b)) std(log(b))]

abcd0 = 1×4

150.0000 150.0000 -6.2681 2.1241

Now fit that to n, as a function of log(b).

fittedmdl = fit(log(b),n,mdl,'start',abcd0)

fittedmdl = 
     General model:
     fittedmdl(x) = a1 + a2*erf((x - a3)/a4)
     Coefficients (with 95% confidence bounds):
       a1 =       142.9  (141.5, 144.4)
       a2 =       141.8  (140, 143.5)
       a3 =       -6.63  (-6.657, -6.603)
       a4 =       2.171  (2.12, 2.222)

plot(fittedmdl,log(b),n,'o')

And that is not too bad. It misses a bit in the upper tail. I might try a simple logistic form instead.

mdl2 = fittype('a1 + a2./(1 + exp(-(x - a3)/a4))','indep','x');
abcd0 = [0 300 mean(log(b)) std(log(b))]

abcd0 = 1×4

0 300.0000 -6.2681 2.1241

fittedmdl2 = fit(log(b),n,mdl2,'start',abcd0)

fittedmdl2 = 
     General model:
     fittedmdl2(x) = a1 + a2./(1 + exp(-(x - a3)/a4))
     Coefficients (with 95% confidence bounds):
       a1 =      -3.423  (-5.409, -1.437)
       a2 =         291  (288.3, 293.7)
       a3 =      -6.645  (-6.663, -6.626)
       a4 =      0.9609  (0.9434, 0.9783)

plot(fittedmdl2,log(b),n,'o')

That seems to fit quite nicely. Now look at the original set of axes.

plot(b,n,'bo',b,fittedmdl2(log(b)),'r-')

And that looks reasonably ok to me. Far better than what you had. Is there any sane reason to believe that what I wrote is the correct model for this process? OF COURSE NOT! Be serious. I just tried a couple of standard model forms, picked out of a hat, based on experience. One worked reasonably well.

The only issue might be that your data never goes below zero on the y axis, but this last model would predict that, if you look at the parameter a1.

If you want to force the curve to not predict below zero, then you could modeify the fit slightly.

mdl3 = fittype('a2./(1 + exp(-(x - a3)/a4))','indep','x');
bcd0 = [300 mean(log(b)) std(log(b))]

bcd0 = 1×3

300.0000 -6.2681 2.1241

fittedmdl3 = fit(log(b),n,mdl3,'start',bcd0)

fittedmdl3 = 
     General model:
     fittedmdl3(x) = a2./(1 + exp(-(x - a3)/a4))
     Coefficients (with 95% confidence bounds):
       a2 =       286.8  (285.4, 288.3)
       a3 =      -6.623  (-6.638, -6.608)
       a4 =       0.939  (0.9258, 0.9522)

plot(fittedmdl3,log(b),n,'o')

And that still looks reasonably good.

plot(b,n,'bo',b,fittedmdl3(log(b)),'r-')

That also seems reasonable. It never predicts a value below zero. If I wanted to look more I could probably find a model for the loglog case too that would fit reasonably well, but why bother.

[Adeb]

Thanks for this effort.