1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
| parfor ii=1:5
% debut=ii;
% fin=ii*10;
for ipopp=ii:ii*10
pfitm1=fitm1;
pteta00=teta00;
pfitness=fitness;
if individus(ipopp)==0;
gene_outputs=gene_outputsbis{ipopp};
x=nnn(ipopp);
m=mmm(ipopp);
N=num_data_points;
matreg=zeros(1196,m+x+1);
% optimization loop for Steiglitz method
loop=1;
l=0;
while loop
l=l+1;
if l~=1,
tt=pteta00{ipopp};
pfitm1(ipopp)=bestfit;
den=1+gene_outputs(:,m+2:end)*tt(m+2:end);
else
pfitm1(ipopp)=Inf;
den=ones(N,1);
bestfit(ipopp)=Inf;
end
% prepare LS matrix
vectreg=y./den;
vectreg=vectreg.*ponder;
% clear matreg;
for k=1:N
matreg(k,1:m+1) = gene_outputs(k,1:m+1)/den(k);
matreg(k,m+2:m+x+1) = -gene_outputs(k,m+2:end)*vectreg(k);
end
for i=1:m+1, matreg(:,i)=matreg(:,i).*ponder; end
prj=matreg'*matreg;
% calculate coefs using SVD least squares on full training data set
try
warning off
teta1=pinv(prj)*(matreg'*vectreg);
teta2=matreg\vectreg;
job=qr([matreg vectreg]);
teta3=triu(job(1:m+1+x,1:m+1+x))\job(1:m+1+x,end);
warning on
catch
pfitness(ipopp)=Inf;
%return;
keyboard
end
% assign poor fitness if any NaN or Inf
if any(isinf(teta1)) || any(isnan(teta1))
pfitness(ipopp)=Inf;
% return;
keyboard
end
% compute prediction of full training data set using the estimated weights
warning off
ypre1=(gene_outputs(:,1:m+1)*teta1(1:m+1))./(1+gene_outputs(:,m+2:end)*teta1(m+2:end));
ypre2=(gene_outputs(:,1:m+1)*teta2(1:m+1))./(1+gene_outputs(:,m+2:end)*teta2(m+2:end));
ypre3=(gene_outputs(:,1:m+1)*teta3(1:m+1))./(1+gene_outputs(:,m+2:end)*teta3(m+2:end));
warning on
% unscale before reporting fitness, if required
% calculate RMS prediction error (fitness)
fitness1=sqrt(mean(((y-ypre1).*ponder).^2));
if fitness1<bestfit, bestfit=fitness1; pteta00{ipopp}=teta1; end
fitness2=sqrt(mean(((y-ypre2).*ponder).^2));
if fitness2<bestfit, bestfit=fitness2; pteta00{ipopp}=teta2; end
fitness3=sqrt(mean(((y-ypre3).*ponder).^2));
if fitness3<bestfit, bestfit=fitness3; pteta00{ipopp}=teta3; end
if l==5, loop=0;
elseif l~=1 & abs(fitm1(ipopp)-bestfit) < 0.01*bestfit, loop=0;
end
end % Steiglitz loop
end
end %for ipop... (boucle loop)
end |
