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 |
Partager