Curve Fitting with erfc

Question

Russell Evans 2020-3-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/513218-curve-fitting-with-erfc

评论： Walter Roberson 2020-3-27

采纳的回答： Walter Roberson

Data.mat

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.

2 个评论
显示无隐藏无

John D'Errico 2020-3-26

编辑：John D'Errico 2020-3-26

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 2020-3-26

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 2020-3-26

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/513218-curve-fitting-with-erfc#answer_422185

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.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Alex Sha 2020-3-27

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 2020-3-27

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);