Fit of a function with a weighted sum

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)

[Mathieu NOE]

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

Answer 1

[arushi]

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, the chi2gof 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,the coeff 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

Answer 2

[Star Strider]

@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.

.

[Camilla Bandinelli]

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]

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.

.