Clear Filters
Clear Filters

Fit of a function with a weighted sum

4 views (last 30 days)
Camilla Bandinelli
Camilla Bandinelli on 20 Dec 2023
Commented: Star Strider on 2 Jan 2024
Hi, I'm trying to fit an ACF of data (that you can find in the attached document) with this weighted sum where the coefficients are . I implemented a while cicle,to consider each time in the sum the previous fit coefficient, but is not working. Also I noticed that the intial values in the fit have a relevant role, so I tried to use the coefficients of an exponential fit as starting point in each trial.
I tried also to implement a chi-squared test to stop the fit.
load('wind_value_restricted.mat');
wind_speed=new_magnitude;
maxLag=100;
acf=autocorr(training_ws,'NumLag', maxLag);
lag=(0:maxLag)';
%evaluate the start point with an exponential fit
exp0=fit(lag,acf,'exp1');
coeff0=[coeffvalues(exp0),0];
y=0;
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
i=1;
% residual=1;
max_iterations=5;
threshold=1e-2;
iter=1;
h=1;
coeff=[0,0,0];
while h==1 && iter <= max_iterations
ft=fittype(@(w,a,b,x)w*exp(-a*x).*cos(b*x)+coeff(1)*exp(-coeff(2)*x).*cos(coeff(3)));
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = coeff0;
f=fit(lag,acf,ft,opts);
coeff_component(i,:)=coeffvalues(f);
% coeff0=coeffvalues(f);
coeff=sum(coeff_component,1);
y=feval(f,lag);
[h, p, stats] = chi2gof(y, 'expected', acf);
i=i+1;
iter=iter+1;
end
figure
plot(lag,feval(f,lag),lag,acf)
  1 Comment
Mathieu NOE
Mathieu NOE on 20 Dec 2023
hello
I am not sure to help you here but a few comments / suggestions
your data are not zero mean , so I was first looking at what is the mean signal here and I ended up founding a exponential decaying (with pulsation = 0 ) signal
so that should be the first term of your serie :
C = 2.27e+00 * e^{-2.218e-05*t}
to get to this result I needed to remove the first 7000 samples , otherwise the fit is quite wrong.
NB you have not supplied any time interval info / sampling rate so I assumed dt = 1 (unit ??)
once we have retrieved that first term,the residual signal can be analysed with fft to pick the dominant remaining frequencies. Those values could be used as IC for a (more robust) fit (I dont have the curve fitting toolbox, so that portion of work I leave it to you)
code :
load('wind_value_restricted.mat');
wind_speed=new_magnitude(6e3:end);
wind_speedS = smoothdata(wind_speed,'gaussian',1e4); % smooth heavily the data to get the DC trend (background)
% exp fit of wind_speedS (DC trend (background))
% model : y = b*exp(m*x);
t = (0:numel(wind_speed)-1)';
P = polyfit(t, log(wind_speed), 1);
m = P(1);
b = exp(P(2));
yfit = b*exp(t*m);
figure(1)
plot(t,wind_speed,t,wind_speedS, t, yfit, '--k')
TE = sprintf('C = %0.2e * e^{%0.3e*t}',b, m);
tt = text(t(round(0.75*end)), yfit(round(0.75*end))*5, TE);
set (tt,'FontSize',15)
legend('raw data','smoothed ','fit');
figure(2)
wind_speed2 = wind_speed-wind_speedS; % remove exponential DC term background
plot(t,wind_speed2)
% fft of the wind_speed2 signal
L = length(wind_speed2); % Length of signal
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Spectrum = fft(wind_speed2, NFFT)/L;
fs = 1; % as we have no time spacing value
f = fs/2*linspace(0,1,NFFT/2+1);
Spectrum = Spectrum(1:NFFT/2+1);
figure(3);
plot(f, abs(Spectrum));
title('FFT of wind speed2')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
hold on
xlim([0 0.005*fs]);
ylim([0 max(abs(Spectrum))*1.5]);
% add findpeaks to display dominant frequencies
NP = 5;
[PKS,LOCS] = findpeaks(abs(Spectrum),f,'SortStr','descend','NPeaks',NP);
plot(LOCS,PKS,'dr');
for ck = 1:NP
ht = text(LOCS(ck),PKS(ck)*1.1,[sprintf('%.2e',LOCS(ck)) 'Hz']);
set (ht,'Rotation',90)
end
hold off

Sign in to comment.

Answers (2)

arushi
arushi on 27 Dec 2023
Hi Camilla,
I understand that you want to fit an ACF of data and also implement a chi-squared test to stop the fit. To fit an autocorrelation function (ACF) in MATLAB, you can use the built-in autocorr function to calculate the ACF of your data and then use the fit function from the Curve Fitting Toolbox to fit a model to the ACF data. Here are a few suggestions for the code you provided:
  1. Use of chi2gof: Assuming the data to be continuous in nature, thechi2gof function is meant for categorical data, not continuous data like an ACF. This function is not appropriate for the type of data you are working with.
  2. The fittype function inside the while loop references coeff, which is supposed to be updated in each iteration, but the way it's written, it will always use the initial value of coeff.Please make sure the values are getting updated correctly after each iteration.Also,thecoeff array is initialized with zeros, which may not be a suitable starting point for the fit.
Link to the documentation of autocorr function - https://www.mathworks.com/help/econ/autocorr.html
Hope this helps.
Thank you

Star Strider
Star Strider on 31 Dec 2023
@Camilla Bandinelli — Responding to your email, first thank you for clarifying my concerns, especially with respect to the ‘training_ws’ vector. Regarding the results in the paper (that I have not read because I do not know what it is) you are attempting to reproduce, it appears that they are essentially doing something similar to the genetic algorithm, however using the goodness-of-fit to determine the fitness of a particular set of parameter estimates to the data actually determines if it is normally-distributed, not that those parameter estimates fit the data. (You can code a test to do that, however the one provided as chi2gof does not do that, at least as I read the documentation.) You can certainly use a statistical test to determine the best fir to the data, however I fail to understand the reasoning behind it. The ‘typical‘ fitness metric for a genetic algorithm (the Global Optimization Toolbox ga function) is the norm, essentially the square root of the sum of squares in this application of it. I doubt that any statistical test could improve on it with respect to determining the fittest individual. Although it would be relatively straightforward to use a statistical test as a fitness metric, it would likely slow the algorithm down considerably. So I am using the norm here.
I ran the ga code 40 times in a loop here (that limit is imposed by the 55-second time limit here, so you can expand it if you wish), and eventually came up with the fittest individual of those runs.
Try this —
load('wind_value_restricted.mat');
% whos
wind_speed=new_magnitude;
training_ws = wind_speed; % Detail Supplied By Camilla Bandinelli
maxLag=100;
acf=autocorr(training_ws,'NumLag', maxLag);
lag=(0:maxLag)';
%evaluate the start point with an exponential fit
exp0=fit(lag,acf,'exp1')
exp0 =
General model Exp1: exp0(x) = a*exp(b*x) Coefficients (with 95% confidence bounds): a = 0.8205 (0.8094, 0.8315) b = -0.00178 (-0.002024, -0.001537)
coeff0=[coeffvalues(exp0),0];
y=0;
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
i=1;
% residual=1;
max_iterations=5;
threshold=1e-2;
iter=1;
h=1;
coeff=[0,0,0];
while h==1 && iter <= max_iterations
ft=fittype(@(w,a,b,x)w*exp(-a*x).*cos(b*x)+coeff(1)*exp(-coeff(2)*x).*cos(coeff(3)));
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = coeff0;
f=fit(lag,acf,ft,opts)
coeff_component(i,:)=coeffvalues(f);
% coeff0=coeffvalues(f);
coeff=sum(coeff_component,1);
y=feval(f,lag);
[h, p, stats] = chi2gof(y, 'expected', acf);
i=i+1;
iter=iter+1;
end
f =
General model: f(x) = w*exp(-a*x).*cos(b*x)+coeff(1)*exp(-coeff(2)*x).*cos(coeff(3)) Coefficients (with 95% confidence bounds): w = 0.8443 (0.8266, 0.862) a = 0.002225 (0.001207, 0.003244) b = 2.889e-05 (-0.351, 0.351)
f =
General model: f(x) = w*exp(-a*x).*cos(b*x)+coeff(1)*exp(-coeff(2)*x).*cos(coeff(3)) Coefficients (with 95% confidence bounds): w = 0.1557 (0.1062, 0.2052) a = 27.31 b = 18.75
f =
General model: f(x) = w*exp(-a*x).*cos(b*x)+coeff(1)*exp(-coeff(2)*x).*cos(coeff(3)) Coefficients (with 95% confidence bounds): w = 0.7729 (0.725, 0.8209) a = 0.0005166 (-0.002353, 0.003386) b = 0.003048 (-0.005935, 0.01203)
f =
General model: f(x) = w*exp(-a*x).*cos(b*x)+coeff(1)*exp(-coeff(2)*x).*cos(coeff(3)) Coefficients (with 95% confidence bounds): w = 0.7151 (0.63, 0.8002) a = -0.002319 (-0.007267, 0.002629) b = 0.007126 (0.001381, 0.01287)
f =
General model: f(x) = w*exp(-a*x).*cos(b*x)+coeff(1)*exp(-coeff(2)*x).*cos(coeff(3)) Coefficients (with 95% confidence bounds): w = 0.6695 (0.553, 0.786) a = -0.004477 (-0.01117, 0.002215) b = 0.008802 (0.003164, 0.01444)
figure
plot(lag,feval(f,lag),lag,acf)
% ---------- GENETIC ALGORITHM SOLUTIONS ----------
objfcn = @(w,a,b,x)w*exp(-a*x).*cos(b*x)+coeff(1)*exp(-coeff(2)*x).*cos(coeff(3));
ftns = @(b) norm(acf - objfcn(b(1),b(2),b(3),lag));
PopSz = 500;
Parms = 3;
optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10); % Options Structure For 'Answers' Problems
tic
for kga = 1:40
t0 = clock;
fprintf('\n----------\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)
GAdata(kga,:) = [fval theta(:).'];
ParamChr = ["w" "a" "b"];
fprintf(1,'\tParameters:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\t%s = %8.5f\n', ParamChr(k1), theta(k1))
end
end
---------- Start Time: 2023-12-31 14:25:55.7734
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:25:57.4655
GA_Time = 1.6922
Elapsed Time: 1.692154000000002E+00 00:00:01.692
Fitness value at convergence = 2.9098 Generations = 263
Parameters:
w = 0.72052 a = 0.00000 b = 6.28318
---------- Start Time: 2023-12-31 14:25:57.4820
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:25:58.5811
GA_Time = 1.0991
Elapsed Time: 1.099117999999997E+00 00:00:01.099
Fitness value at convergence = 2.9098 Generations = 230
Parameters:
w = 0.72060 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:25:58.5894
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:25:59.6697
GA_Time = 1.0803
Elapsed Time: 1.080308000000002E+00 00:00:01.080
Fitness value at convergence = 2.9098 Generations = 227
Parameters:
w = 0.72056 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:25:59.6776
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:00.6108
GA_Time = 0.9332
Elapsed Time: 9.331519999999998E-01 00:00:00.933
Fitness value at convergence = 2.9099 Generations = 196
Parameters:
w = 0.72066 a = 0.00001 b = 6.28319
---------- Start Time: 2023-12-31 14:26:00.6203
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:02.2622
GA_Time = 1.6419
Elapsed Time: 1.641901000000000E+00 00:00:01.641
Fitness value at convergence = 2.9099 Generations = 269
Parameters:
w = 0.72079 a = 0.00001 b = 6.28319
---------- Start Time: 2023-12-31 14:26:02.2724
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:03.4270
GA_Time = 1.1546
Elapsed Time: 1.154634000000000E+00 00:00:01.154
Fitness value at convergence = 2.9099 Generations = 245
Parameters:
w = 0.72072 a = 0.00001 b = 6.28319
---------- Start Time: 2023-12-31 14:26:03.4351
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:05.2050
GA_Time = 1.7700
Elapsed Time: 1.769961000000000E+00 00:00:01.769
Fitness value at convergence = 2.9098 Generations = 385
Parameters:
w = 0.72048 a = 0.00000 b = 0.00000
---------- Start Time: 2023-12-31 14:26:05.2137
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:06.1972
GA_Time = 0.9835
Elapsed Time: 9.835009999999995E-01 00:00:00.983
Fitness value at convergence = 2.9098 Generations = 205
Parameters:
w = 0.72051 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:06.2060
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:07.5222
GA_Time = 1.3162
Elapsed Time: 1.316179999999999E+00 00:00:01.316
Fitness value at convergence = 2.9098 Generations = 275
Parameters:
w = 0.72050 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:07.5314
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:08.5327
GA_Time = 1.0013
Elapsed Time: 1.001320000000000E+00 00:00:01.001
Fitness value at convergence = 2.9099 Generations = 213
Parameters:
w = 0.72094 a = 0.00001 b = 6.28319
---------- Start Time: 2023-12-31 14:26:08.5424
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:09.8923
GA_Time = 1.3499
Elapsed Time: 1.349884000000001E+00 00:00:01.349
Fitness value at convergence = 2.9099 Generations = 284
Parameters:
w = 0.72068 a = 0.00001 b = 6.28319
---------- Start Time: 2023-12-31 14:26:09.9020
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:11.0644
GA_Time = 1.1624
Elapsed Time: 1.162423000000000E+00 00:00:01.162
Fitness value at convergence = 2.9098 Generations = 243
Parameters:
w = 0.72053 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:11.0745
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:12.4028
GA_Time = 1.3283
Elapsed Time: 1.328329999999999E+00 00:00:01.328
Fitness value at convergence = 2.9098 Generations = 272
Parameters:
w = 0.72052 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:12.4187
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:13.3915
GA_Time = 0.9728
Elapsed Time: 9.727519999999998E-01 00:00:00.972
Fitness value at convergence = 2.9098 Generations = 192
Parameters:
w = 0.72054 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:13.4031
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:14.4954
GA_Time = 1.0923
Elapsed Time: 1.092252000000000E+00 00:00:01.092
Fitness value at convergence = 2.9099 Generations = 228
Parameters:
w = 0.72084 a = 0.00001 b = 6.28319
---------- Start Time: 2023-12-31 14:26:14.5066
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:15.5074
GA_Time = 1.0008
Elapsed Time: 1.000772000000000E+00 00:00:01.000
Fitness value at convergence = 2.9098 Generations = 209
Parameters:
w = 0.72059 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:15.5189
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:16.6926
GA_Time = 1.1737
Elapsed Time: 1.173739000000001E+00 00:00:01.173
Fitness value at convergence = 2.9098 Generations = 247
Parameters:
w = 0.72052 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:16.7053
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:18.2200
GA_Time = 1.5147
Elapsed Time: 1.514737000000000E+00 00:00:01.514
Fitness value at convergence = 2.9098 Generations = 322
Parameters:
w = 0.72053 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:18.2318
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:19.5148
GA_Time = 1.2830
Elapsed Time: 1.282973999999999E+00 00:00:01.282
Fitness value at convergence = 2.9099 Generations = 273
Parameters:
w = 0.72086 a = 0.00001 b = 6.28319
---------- Start Time: 2023-12-31 14:26:19.5273
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:20.9990
GA_Time = 1.4716
Elapsed Time: 1.471627000000002E+00 00:00:01.471
Fitness value at convergence = 2.9098 Generations = 310
Parameters:
w = 0.72048 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:21.0120
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:22.6157
GA_Time = 1.6038
Elapsed Time: 1.603783999999997E+00 00:00:01.603
Fitness value at convergence = 2.9098 Generations = 339
Parameters:
w = 0.72049 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:22.6284
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:23.7898
GA_Time = 1.1614
Elapsed Time: 1.161424000000000E+00 00:00:01.161
Fitness value at convergence = 2.9098 Generations = 262
Parameters:
w = 0.72049 a = 0.00000 b = 0.00000
---------- Start Time: 2023-12-31 14:26:23.8026
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:25.1569
GA_Time = 1.3543
Elapsed Time: 1.354309000000001E+00 00:00:01.354
Fitness value at convergence = 2.9098 Generations = 288
Parameters:
w = 0.72051 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:25.1700
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:26.6721
GA_Time = 1.5021
Elapsed Time: 1.502136000000000E+00 00:00:01.502
Fitness value at convergence = 2.9098 Generations = 318
Parameters:
w = 0.72059 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:26.6859
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:27.6855
GA_Time = 0.9997
Elapsed Time: 9.996550000000006E-01 00:00:00.999
Fitness value at convergence = 2.9099 Generations = 212
Parameters:
w = 0.72080 a = 0.00001 b = 6.28319
---------- Start Time: 2023-12-31 14:26:27.6994
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:28.6990
GA_Time = 0.9996
Elapsed Time: 9.995829999999977E-01 00:00:00.999
Fitness value at convergence = 2.9098 Generations = 212
Parameters:
w = 0.72057 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:28.7133
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:29.7831
GA_Time = 1.0698
Elapsed Time: 1.069754000000000E+00 00:00:01.069
Fitness value at convergence = 2.9098 Generations = 226
Parameters:
w = 0.72058 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:29.7977
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:31.3042
GA_Time = 1.5065
Elapsed Time: 1.506498000000001E+00 00:00:01.506
Fitness value at convergence = 2.9098 Generations = 323
Parameters:
w = 0.72053 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:31.3198
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:32.2299
GA_Time = 0.9101
Elapsed Time: 9.100859999999962E-01 00:00:00.910
Fitness value at convergence = 2.9099 Generations = 191
Parameters:
w = 0.72072 a = 0.00001 b = 6.28319
---------- Start Time: 2023-12-31 14:26:32.2448
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:33.6549
GA_Time = 1.4101
Elapsed Time: 1.410072000000000E+00 00:00:01.410
Fitness value at convergence = 2.9098 Generations = 302
Parameters:
w = 0.72055 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:33.6711
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:35.0149
GA_Time = 1.3438
Elapsed Time: 1.343785999999994E+00 00:00:01.343
Fitness value at convergence = 2.9098 Generations = 284
Parameters:
w = 0.72049 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:35.0311
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:36.6171
GA_Time = 1.5861
Elapsed Time: 1.586091000000003E+00 00:00:01.586
Fitness value at convergence = 2.9098 Generations = 337
Parameters:
w = 0.72049 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:36.6327
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:37.5239
GA_Time = 0.8912
Elapsed Time: 8.912260000000032E-01 00:00:00.891
Fitness value at convergence = 2.9098 Generations = 187
Parameters:
w = 0.72062 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:37.5395
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:38.3660
GA_Time = 0.8265
Elapsed Time: 8.265200000000021E-01 00:00:00.826
Fitness value at convergence = 2.9098 Generations = 174
Parameters:
w = 0.72052 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:38.3843
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:39.2827
GA_Time = 0.8984
Elapsed Time: 8.984390000000033E-01 00:00:00.898
Fitness value at convergence = 2.9098 Generations = 188
Parameters:
w = 0.72049 a = 0.00000 b = 6.28318
---------- Start Time: 2023-12-31 14:26:39.3012
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:40.5822
GA_Time = 1.2810
Elapsed Time: 1.280996000000002E+00 00:00:01.280
Fitness value at convergence = 2.9098 Generations = 273
Parameters:
w = 0.72054 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:40.5991
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:42.3464
GA_Time = 1.7473
Elapsed Time: 1.747278999999999E+00 00:00:01.747
Fitness value at convergence = 2.9099 Generations = 375
Parameters:
w = 0.72079 a = 0.00001 b = 6.28319
---------- Start Time: 2023-12-31 14:26:42.3634
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:43.4641
GA_Time = 1.1007
Elapsed Time: 1.100693000000000E+00 00:00:01.100
Fitness value at convergence = 2.9098 Generations = 236
Parameters:
w = 0.72055 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:43.4813
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:44.8134
GA_Time = 1.3321
Elapsed Time: 1.332079000000000E+00 00:00:01.332
Fitness value at convergence = 2.9098 Generations = 284
Parameters:
w = 0.72051 a = 0.00000 b = 6.28319
---------- Start Time: 2023-12-31 14:26:44.8332
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2023-12-31 14:26:46.0816
GA_Time = 1.2484
Elapsed Time: 1.248381999999999E+00 00:00:01.248
Fitness value at convergence = 2.9098 Generations = 264
Parameters:
w = 0.72055 a = 0.00000 b = 6.28319
toc
Elapsed time is 50.332348 seconds.
format long
GAdata
GAdata = 40×4
2.909816779054780 0.720521355271339 0.000001100420952 6.283184939384460 2.909841922361221 0.720604048013687 0.000003346204758 6.283185431122780 2.909826948488437 0.720555702447891 0.000002009630203 6.283185371518135 2.909861215382331 0.720663586854935 0.000005064368248 6.283189453125000 2.909894837503769 0.720789332032204 0.000008047938347 6.283192091941833 2.909876062026678 0.720715050339699 0.000006383895874 6.283186399698257 2.909804496722381 0.720482776165009 0.000000000715256 0.000000008821487 2.909811882022058 0.720506573319435 0.000000662088394 6.283188470959663 2.909807694298695 0.720495337843895 0.000000287175179 6.283186059832573 2.909948798620759 0.720944041609764 0.000012812137604 6.283185073494911
[MaxFtns,idx] = min(GAdata(:,1)) % Best Fitness & Index
MaxFtns =
2.909804496722381
idx =
7
theta = GAdata(idx, 2:end) % Paramters: 'w', 'a', 'b'
theta = 1×3
0.720482776165009 0.000000000715256 0.000000008821487
t = lag;
c = acf;
tv = linspace(min(t), max(t));
Cfit = objfcn(theta(1),theta(2),theta(3), tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
This is the best I can do with these data. The fittest individual (as defined by the norm of the rtesiduals) has a fitness value of about 2.909804496722381.
.
  2 Comments
Camilla Bandinelli
Camilla Bandinelli on 2 Jan 2024
Thank you for your help.
It's the first time that I see the genetic algorithm.
Maybe I didn't explain myself well but in your script you just find only one set of optimal coefficient for the fit. What I need is a sum of function with different coefficents which togheter constituites the best fit for the acf.
The paper that I'm tryng to reproduce http://faraday1.ucd.ie/archive/papers/windauto.pdf . You can find the ACf fit in section 4.1
Star Strider
Star Strider on 2 Jan 2024
All the parameters are adjusted individually, each on its own trajectory. I initially did not pick up on the fact that my code did not actually estimate all the necessary parameters. The objective function I wrote to use with ga needed to be changed to incorporate (and estimate) the ‘coeff’ parameters as well.
Then it works —
load('wind_value_restricted.mat');
% whos
wind_speed=new_magnitude;
training_ws = wind_speed; % Detail Supplied By Camilla Bandinelli
maxLag=100;
acf=autocorr(training_ws,'NumLag', maxLag);
lag=(0:maxLag)';
%evaluate the start point with an exponential fit
exp0=fit(lag,acf,'exp1')
exp0 =
General model Exp1: exp0(x) = a*exp(b*x) Coefficients (with 95% confidence bounds): a = 0.8205 (0.8094, 0.8315) b = -0.00178 (-0.002024, -0.001537)
coeff0=[coeffvalues(exp0),0];
y=0;
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
i=1;
% residual=1;
max_iterations=5;
threshold=1e-2;
iter=1;
h=1;
coeff=[0,0,0];
% while h==1 && iter <= max_iterations
% ft=fittype(@(w,a,b,x)w*exp(-a*x).*cos(b*x)+coeff(1)*exp(-coeff(2)*x).*cos(coeff(3)));
% opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
% opts.Display = 'Off';
% opts.StartPoint = coeff0;
% f=fit(lag,acf,ft,opts)
% coeff_component(i,:)=coeffvalues(f);
% % coeff0=coeffvalues(f);
% coeff=sum(coeff_component,1);
% y=feval(f,lag);
% [h, p, stats] = chi2gof(y, 'expected', acf);
% i=i+1;
% iter=iter+1;
% end
%
% figure
% plot(lag,feval(f,lag),lag,acf)
% ---------- GENETIC ALGORITHM SOLUTIONS ----------
objfcn = @(w,a,b,coeff1,coeff2,coeff3,x) w*exp(-a*x).*cos(b*x)+coeff1*exp(-coeff2*x).*cos(coeff3);
ftns = @(b) norm(acf - objfcn(b(1),b(2),b(3),b(4),b(5),b(6),lag));
PopSz = 500;
Parms = 6;
optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10); % Options Structure For 'Answers' Problems
tic
for kga = 1:9
t0 = clock;
fprintf('\n----------\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)
GAdata(kga,:) = [fval theta(:).'];
ParamChr = ["w" "a" "b" "coeff(1)", "coeff(2)", "coeff(3)"];
fprintf(1,'\tParameters:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\t%-8s = %8.5f\n', ParamChr(k1), theta(k1))
end
end
---------- Start Time: 2024-01-02 19:05:57.9561
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2024-01-02 19:06:03.2327
GA_Time = 5.2766
Elapsed Time: 5.276624999999996E+00 00:00:05.276
Fitness value at convergence = 0.0601 Generations = 918
Parameters:
w = 0.18562 a = 0.18669 b = 6.28319 coeff(1) = 4.36904 coeff(2) = 0.00100 coeff(3) = 1.39168
---------- Start Time: 2024-01-02 19:06:03.2637
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2024-01-02 19:06:07.4812
GA_Time = 4.2175
Elapsed Time: 4.217514000000000E+00 00:00:04.217
Fitness value at convergence = 0.0601 Generations = 795
Parameters:
w = 0.18553 a = 0.18763 b = 6.28320 coeff(1) = 6.11055 coeff(2) = 0.00100 coeff(3) = 1.44302
---------- Start Time: 2024-01-02 19:06:07.4927
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2024-01-02 19:06:16.5702
GA_Time = 9.0775
Elapsed Time: 9.077470000000000E+00 00:00:09.077
Fitness value at convergence = 0.0602 Generations = 1784
Parameters:
w = 0.77954 a = 0.00102 b = 6.28319 coeff(1) = 2.97657 coeff(2) = 0.19077 coeff(3) = 4.77469
---------- Start Time: 2024-01-02 19:06:16.5823
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2024-01-02 19:06:24.4028
GA_Time = 7.8204
Elapsed Time: 7.820421999999997E+00 00:00:07.820
Fitness value at convergence = 0.0602 Generations = 1520
Parameters:
w = 0.18533 a = 0.19057 b = 0.00001 coeff(1) = 2.31002 coeff(2) = 0.00102 coeff(3) = 5.05658
---------- Start Time: 2024-01-02 19:06:24.4147
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2024-01-02 19:06:27.7778
GA_Time = 3.3631
Elapsed Time: 3.363118999999998E+00 00:00:03.363
Fitness value at convergence = 0.0601 Generations = 707
Parameters:
w = 0.18563 a = 0.18676 b = 0.00002 coeff(1) = 5.87659 coeff(2) = 0.00100 coeff(3) = 1.43795
---------- Start Time: 2024-01-02 19:06:27.7944
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2024-01-02 19:06:33.3864
GA_Time = 5.5920
Elapsed Time: 5.592003999999996E+00 00:00:05.592
Fitness value at convergence = 0.0601 Generations = 973
Parameters:
w = 0.18553 a = 0.18776 b = 6.28319 coeff(1) = 6.25435 coeff(2) = 0.00100 coeff(3) = 1.44597
---------- Start Time: 2024-01-02 19:06:33.3970
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2024-01-02 19:06:36.1812
GA_Time = 2.7843
Elapsed Time: 2.784284999999997E+00 00:00:02.784
Fitness value at convergence = 0.0601 Generations = 603
Parameters:
w = 0.18551 a = 0.18794 b = 0.00001 coeff(1) = 7.38458 coeff(2) = 0.00100 coeff(3) = 7.74833
---------- Start Time: 2024-01-02 19:06:36.1908
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2024-01-02 19:06:42.4817
GA_Time = 6.2909
Elapsed Time: 6.290863000000002E+00 00:00:06.290
Fitness value at convergence = 0.0601 Generations = 1310
Parameters:
w = 0.18550 a = 0.18819 b = 6.28319 coeff(1) = 5.35260 coeff(2) = 0.00101 coeff(3) = 4.85841
---------- Start Time: 2024-01-02 19:06:42.4913
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.
Stop Time: 2024-01-02 19:06:49.6758
GA_Time = 7.1844
Elapsed Time: 7.184440000000002E+00 00:00:07.184
Fitness value at convergence = 0.0604 Generations = 1426
Parameters:
w = 0.18517 a = 0.19262 b = 6.28319 coeff(1) = 3.80063 coeff(2) = 0.00103 coeff(3) = 4.91910
toc
Elapsed time is 51.736998 seconds.
format long
GAdata
GAdata = 9×7
0.060098662512042 0.185622520327568 0.186686070322990 6.283189035177231 4.369037237286568 0.000996715664864 1.391682169318199 0.060115087946716 0.185532517313957 0.187632030963898 6.283198549270629 6.110545777916909 0.001002467513084 1.443020674586296 0.060239242473639 0.779542485594749 0.001020976185799 6.283185716509819 2.976571087837219 0.190768869161606 4.774685595393181 0.060228070996826 0.185329077243805 0.190565874695778 0.000010191440582 2.310021248579025 0.001019797563553 5.056580905079842 0.060099325649120 0.185626501560211 0.186757689476013 0.000017189145088 5.876594552040101 0.000997043490410 1.437948901772499 0.060117861569760 0.185530356645584 0.187762442827225 6.283186613082886 6.254346811413765 0.001003182768822 1.445968616008759 0.060122505315540 0.185513157606125 0.187941393971443 0.000005973815918 7.384582336425781 0.001004297256470 7.748329963326454 0.060129330473981 0.185495641946793 0.188193074584007 6.283186741471290 5.352595896601677 0.001005787372589 4.858409527897835 0.060364289794760 0.185173566222191 0.192624524831772 6.283186549544334 3.800628127813339 0.001031921744347 4.919103551864624
[MaxFtns,idx] = min(GAdata(:,1)) % Best Fitness & Index
MaxFtns =
0.060098662512042
idx =
1
theta1 = GAdata(idx, 2:end) % Paramters: 'w', 'a', 'b'
theta1 = 1×6
0.185622520327568 0.186686070322990 6.283189035177231 4.369037237286568 0.000996715664864 1.391682169318199
ACFmodel = fitnlm(lag, acf, @(b,x)objfcn(b(1),b(2),b(3),b(4),b(5),b(6),x), theta) % Calculate & Return Statistics, Tweak Fit
Warning: The Jacobian at the solution is ill-conditioned, and some model parameters may not be estimated well (they are not identifiable). Use caution in making predictions.
ACFmodel =
Nonlinear regression model: y ~ objfcn(b1,b2,b3,b4,b5,b6,x) Estimated Coefficients: Estimate SE tStat pValue ___________________ ____________________ __________________ _____________________ b1 0.185173566222207 0.00616277657451726 30.0471003586094 1.27234979988776e-50 b2 0.192624524831774 0.0158186304285436 12.1770671425635 3.51211348505923e-21 b3 6.28318577015851 92.8752658899986 0.0676518738325908 0.946203463501781 b4 3.80062812781334 3.91458533183139e-05 97088.9073973835 0 b5 0.00103192174430148 5.00150533492884e-05 20.6322232047797 4.63566969048404e-37 b6 4.91910355186462 0.000709744237260698 6930.81154254977 2.20484873960777e-275 Number of observations: 101, Error degrees of freedom: 96 Root Mean Squared Error: 0.00616 R-Squared: 0.984, Adjusted R-Squared 0.983 F-statistic vs. zero model: 3.02e+05, p-value = 9.64e-200
theta2 = ACFmodel.Coefficients.Estimate
theta2 = 6×1
0.185173566222207 0.192624524831774 6.283185770158513 3.800628127813339 0.001031921744301 4.919103551864624
t = lag;
c = acf;
tv = linspace(min(t), max(t));
Cfit1 = objfcn(theta1(1),theta1(2),theta1(3),theta1(4),theta1(5),theta1(6), tv);
Cfit2 = objfcn(theta2(1),theta2(2),theta2(3),theta2(4),theta2(5),theta2(6), tv);
figure
hd = plot(t, c, 'p', 'DisplayName','Data');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
plot(tv, Cfit2, '-r', 'DisplayName','GA Tweaked');
hlp = plot(tv, Cfit1, 'DisplayName','GA');
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
legend('Location','best')
I used the fitnlm function here to provide the statistics of the parameters, and also to tweak the parameters ga estimated, using a gradient descent approach. (This is commonly done, and ga actually provides a 'HybridFcn' option to do just that as part of the ga parameter estimation process. I did it separately here because I wanted the statistics.) The tweaking did not change the parameters or the fit noticably.
Except for ‘b(3)’ corresponding to ‘b’ in the original objective function, all the parameters are highly statistically significant, as is the fit.
I kept the number of iterations in the loop to 9 because of the 55-second time limit here. Increase it to get more ga runs and perhaps better parameter estimates.
.

Sign in to comment.

Products


Release

R2022b

Community Treasure Hunt

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

Start Hunting!