MATLAB Answers

Goodness of fit , residual plot for fminbnd fitting

조회 수: 3(최근 30일)
Anand Rathnam
Anand Rathnam 2021년 7월 27일
댓글: Anand Rathnam 2021년 7월 29일
I have the data fitting using fminbnd solver. I am wondering how to evaluate the goodness of the fit and obtain residual plots. I dont see any options to display the residuals for fminbnd solver like other solvers. Below is my code:
t1 = [0:300:28800]'; % input X data
% input Y data
y_obs = [
0
0.0350666
0.170773
0.298962
0.400482
0.481344
0.541061
0.588307
0.626498
0.657928
0.684406
0.705545
0.721963
0.738828
0.753222
0.765903
0.776001
0.786196
0.795698
0.804062
0.81206
0.820732
0.825598
0.832848
0.837778
0.8436
0.848495
0.852999
0.858091
0.863251
0.86657
0.870919
0.875362
0.879617
0.882049
0.884957
0.887106
0.889922
0.894813
0.896395
0.900105
0.903234
0.905787
0.907843
0.909099
0.913799
0.914104
0.916195
0.920424
0.922772
0.923837
0.922742
0.924935
0.927408
0.92851
0.930684
0.930988
0.933917
0.935012
0.938209
0.940926
0.942448
0.943642
0.942436
0.94564
0.946308
0.949709
0.950971
0.951911
0.954338
0.955225
0.958114
0.958801
0.962341
0.963808
0.965617
0.965214
0.966752
0.971954
0.971949
0.973827
0.977233
0.977157
0.980893
0.979747
0.981409
0.984914
0.986015
0.986951
0.990709
0.990882
0.991937
0.992701
0.996347
0.998733
0.999351
1
];
%************
% Guessing the initial assumption d0 by finding the minimum using min function
%Implementing fminbnd instead of lsqnonlin
pfun = @(d) norm( ypred(d, t1) - y_obs);
dsamps=linspace(0,15e-10,50);
[~,imin]=min( arrayfun(pfun,dsamps) );
[best_d,fval,exitflag,output]=fminbnd(pfun,dsamps(imin-1), dsamps(imin+1),...
optimset('TolX',1e-14))
best_d =
6.545737727815822e-10
fval =
4.578375975388594e-01
exitflag =
1
output = struct with fields:
iterations: 8 funcCount: 9 algorithm: 'golden section search, parabolic interpolation' message: 'Optimization terminated:↵ the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-14 ↵'
exitflag,
exitflag =
1
best_d,
best_d =
6.545737727815822e-10
%****************
%************
t2 = [0:5/60:8]';
predicted_y = ypred(best_d, t1);
figure(2)
plot(t2, y_obs, 'ro', t2, predicted_y, 'b-');
grid on
set(gca,'XLim',[0 8])
set(gca,'XTick',(0:0.5:8))
ylim([ 0 1.4])
ylabel('At/Ainf')
xlabel('Time in h')
legend({'observed', 'predicted'})
title('fitting upto full 8 hours; thickness = 2.63')
%***************
%Plotting the verification of minimum
figure(1)
fplot(pfun,[0,4e-10])
Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.
hold on; plot(best_d*[1,1],ylim,'--rx');
title('Verification of whether minimum was correctly found by fminbnd solver')
hold off
%***** Fitting model equation
function y_pred = ypred(d, t1)
a=0.00263;
gama = 0.01005;
L2 = zeros(14,1);
L3 = zeros(100,1);
L4 = zeros(100,1);
L5 = zeros(100,1);
S= zeros(97,1);
y_pred = zeros(97,1);
% t = 0;
L1 = ((8*gama)/((pi*(1-exp(-2*gama*a)))));
format longE
for t = t1(:).'
for n=0:1:100
L2(n+1) = exp((((2*n + 1)^2)*-d*pi*pi*t)/(4*a*a));
L3(n+1) = (((-1)^n)*2*gama)+(((2*n+1)*pi)*exp(-2*gama*a))/(2*a);
L4(n+1)= ((2*n)+1)*((4*gama*gama)+((((2*n)+1)*pi)/(2*a))^2);
L5(n+1) = ((L2(n+1)*L3(n+1))/L4(n+1));
end
S((t/300) +1) = sum(L5);
y_pred((t/300)+1)= 1 -(L1*S((t/300) +1)); % predicted data
end
end
  댓글 수: 2
Adam Danz
Adam Danz 2021년 7월 27일
The prediction-observation plot does not look like a good fit to me.

댓글을 달려면 로그인하십시오.

채택된 답변

Adam Danz
Adam Danz 2021년 7월 27일
편집: Adam Danz 2021년 7월 29일
See matlab's documentation on goodness of fit. There are lots of options and you've got all the variables you need to move forward.
The residuals are just the difference between the predicted y-values and the actual y values. Here's how to plot them:
figure()
tiledlayout(2,1)
nexttile
stem(t2, predicted_y-y_obs)
xlabel('t2')
ylabel('residuals')
nexttile()
histogram(predicted_y-y_obs)
xlabel('residual')
ylabel('frequency')

추가 답변(0개)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by