Curve Fitting with erfc

Greetings,

I am trying to fit the attached data (Ydata vs time) with this function attached as a picture.

The coefficients for fitting are B and D and must be positive. L is 1000, Q is -0.0038

I am trying to use cftool with the custom equation

((D*-0.0038/1000^2)*B) .*exp((B.^2) .*(((D.*x))./1000.^2)) .*erfc(B.*(sqrt(((D.*x))./1000^2)))

but it just says my imputs are not real. It would be helpful to see how to fit the data.

[John D'Errico]

Why does your custom equation not look like the model you claim to want to use?

Here, for example, we see

exp((B.^2) .*((D.*x)./1000.^2))

what appears to be a term like D*t/L^2.

But then in the next term, if x is supposed to be t, then all you have is sqrt(t).

As far as knowing why your inputs are not real, we cannot conjecture what other wrong things you did, since you already have make a last one significant mistake. That suggests your skills at coding are not incredibly great. While there is nothing overtly wrong with that, it does suggest that you have surely made at least a few other mistakes in your code, especially given the vague error you have described. If you want help, the easiest way to get help is to show what you tried. Then someone can just correct your mistakes.

As well, when you have an error, rather than just trying to describe what you think is pertinent about the error message, show the complete error text, so everything in red.

[Russell Evans]

I edited the question to reflect

((D*-0.0038/1000^2)*B) .*exp((B.^2) .*(((D.*x))./1000.^2)) .*erfc(B.*(sqrt(((D.*x))./1000^2)))

This inital typo was my fault. "x" is supposed to be "t"

As you mention, my skills at coding are not incredibly great. Thanks for pointing that out.

No red error shows. I do get a yellow warning message

Error in fittype expression ==> ((D.*-0.0038./1000.^2).*B) .*exp((B.^2) .*(((D.*x))./1000.^2)) .*erfc(B.*sqrt(x))

??? Input must be real and full.

I have made sure that my x data is real even when imput into erfc(1/sqrt(x)) outside of cftool

Answer 1

[Walter Roberson]

I was able to track this down.

What is happening is that cftool likes to draw the smooth line, from before the start of the x data, to after the end of the x data. It uses a 5% margin, starting the projection at min(x) - (max(x)-min(x))*0.05 . That happens to be negative for your data, and the sqrt(x) goes complex and erfc complains.

Unfortunately there are no controls for where the line is drawn. When I changed the sqrt(x) to sqrt(max(x,0)) to prevent square root of negatives, the smooth line doesn't appear and the fit is not good.

The easiest solution, if you have the symbolic toolbox, is build the function for sum-of-squared residues, and minimize the function:

syms D B
x = time; y = ypap;
residue = sum((((D*-0.0038/1000^2)*B) .*exp((B.^2) .*((D.*x)./1000.^2)) .*erfc(B.*sqrt(x)) - y).^2);
F = matlabFunction(residue, 'vars', {[B D]});
goodBD = fminsearch(F, [1 2]);

goodB = goodBD(1);
goodD = goodBD(2);

The fminsearch is fast, and the solution it gets you is pretty much within round-off noise of the best solution I have been able to find with more careful analysis.

[Russell Evans]

Thank you. This has gotten me closest so far. I did make the mistake of not typing the function correctly

((D*-0.0038/1000^2)*B) .*exp((B.^2) .*(((D.*x))./1000.^2)) .*erfc(B.*(sqrt(((D.*x))./1000^2)))

is the updated answer.

I am attempting to verify data in a peer reviewed journal article where that function above is the red fit to the data i attached. They transform the x axis from time to the inverse root of time when they plot it, but thats the only difference here. They claim using cftool in matlab to do this fit.

[Russell Evans]

Using the sqrt(max(x,0)) I was able to get here,so I think I am close. Maybe just missing one thing.

[Alex Sha]

Hi, refer the result below (with the function: y=((D*(-0.0038)/1000^2)*B)*exp((B^2)*(((D*x))/1000^2))*erfc(B*(sqrt(((D*x))/1000^2)))), the result looks good but be careful in the end of chart lines (with the blue circle), it seems to have a problem with that fitting function.

Root of Mean Square Error (RMSE): 1.09533665147734E-8
Sum of Squared Residual: 4.07919209223658E-15
Correlation Coef. (R): 0.99976230572428
R-Square: 0.999524667947129

Parameter	Best Estimate
----------	-------------
b	122.725546738343
d	5.47654549122705

[Walter Roberson]

Why did you set minimum for B and D to be 1, when earlier you just said "positive" ?

syms D B real
x = time; y = ypap;

ERFC = @(G)simplify(arrayfun(@(H)piecewise(imag(H)~=0,1e10,erfc(H)),G));
I = @(X) ((D*-0.0038/1000^2)*B) .*exp((B.^2) .*(((D.*X))./1000.^2)) .*ERFC(B.*(sqrt(((D.*X))./1000^2)));
parts = arrayfun(I, x)-y;
residue = simplify(sum(simplify(parts.^2)));
F = matlabFunction(residue, 'vars', {[B D]}, 'file', 'residue_fun.m');
goodBD = fminsearch(F, [100 1]);

goodB = goodBD(1);
goodD = goodBD(2);

The starting point was important for fast convergence.