somplsoptimise.m

有监督自组织映射-偏最小二乘算法(A supervised self-organising map–partial least squares algorithm),可以用语多变量数据的回归分析

function [OPT] = SOMPLSOptimise(Nval,Percentage,BestMLK,MLKP)

global p1 p2
sigma=[0.25 1.0 4.0 16.0 64.0 256];
omega=[0.0625 1 16 256];
MapMethod=MLKP.TrainingMethod;

MLKP.TrainingMethod='PLS';
if (upper(MLKP.CrossValDataSelection) == 'Y')
    MLKP.CreateDataSet='KStone';
else
    MLKP.CreateDataSet='Random';
end
MLKP.ScaleTypeX = 'Raw';   
MLKP.ScaleTypeY = 'Raw';
MLKP.FixedRandomData='N';
MLKP.BinnedData = 'N';
MLKP.SilentMode='Y';
    
NMLKP=MLKP;
NBestMLK=BestMLK;

Y=BestMLK.Ytr;
[Nobjects,NvarY]=size(Y);
ScaleFact=max(max(dist(BestMLK.Xtr,BestMLK.Xtr')));
XMap=BestMLK.XMap/ScaleFact;
X=BestMLK.Xtr/ScaleFact;
for ival=1:Nval
    [CV(ival).Xtr,CV(ival).Ytr,CV(ival).Xtest,CV(ival).Ytest] = MakeBalancedSet(X,Y,Percentage,MLKP);
    if (upper(MLKP.ProcMode) == 'KPLS')
        CV(ival).XMap=CV(ival).Xtr';
    end
end
if (upper(MLKP.ProcMode) == 'XMAP')
    [NvarX,MaxLV]=size(XMap);
else
    clear XMap;
    [NvarX,MaxLV]=size(CV(1).XMap);
end
if (MLKP.MaxLatentVar > 0)
    MaxLV=min(MLKP.MaxLatentVar,MaxLV);
end

if (upper(MLKP.ProblemType) == 'CLS')
    OldPCTest=0;
    MLXPercCorrTrain=BestMLK.PercCorrTrain;
    MLXPercCorrTest=BestMLK.PercCorrTest;
else
    OldPCTest=-realmax;
    MLXPercCorrTrain=BestMLK.RmseTrainAll;
    MLXPercCorrTest=BestMLK.RmseTestAll;
end

for sim=1:3
    b0max=1;
    bulentmax=2;
    omegamax=1;
    if (sim==1)
        DistTypeX = 'EUC';
    elseif (sim==2)
        DistTypeX = 'COR';
    elseif (sim==3)
        DistTypeX = 'PUK';
        b0max=0;
        bulentmax=length(sigma);
        omegamax=length(omega);
    end
    for bulent=bulentmax:-1:1
        for iomega=1:omegamax
            if (DistTypeX == 'PUK')
                p1=sigma(bulent);
                p2=omega(iomega);
            end
            for b0=0:b0max
                for ival=1:Nval
                    if (upper(MLKP.ProcMode) == 'KPLS')
                        XMap=CV(ival).XMap;
                    end
                    Xtr=CV(ival).Xtr;
                    Xtest=CV(ival).Xtest;
                    if (DistTypeX == 'COR')
                        PLS(ival).Xtr=(b0+Xtr*XMap).^bulent;
                        PLS(ival).Xtest=(b0+Xtest*XMap).^bulent;
                    elseif (DistTypeX == 'EUC')
                        PLS(ival).Xtr=(b0+dist(Xtr,XMap)).^bulent;
                        PLS(ival).Xtest=(b0+dist(Xtest,XMap)).^bulent;
                    elseif (DistTypeX == 'RBF')
                        PLS(ival).Xtr=exp(-(dist(Xtr,XMap).^2)/(2*sigma(bulent)^2));
                        PLS(ival).Xtest=exp(-(dist(Xtest,XMap).^2)/(2*sigma(bulent)^2));
                    elseif (DistTypeX == 'PUK')
                        [NobjTrain,NvarX]=size(Xtr);
                        [NobjTest,NvarX]=size(Xtest);
                        [NvarX,Nunit]=size(XMap);
                        for jobj=1:Nunit
                            for iobj=1:NobjTrain
                                PLS(ival).Xtr(iobj,jobj)=svkernel('puk',Xtr(iobj,:),XMap(:,jobj)');
                            end
                            for iobj=1:NobjTest
                                PLS(ival).Xtest(iobj,jobj)=svkernel('puk',Xtest(iobj,:),XMap(:,jobj)');
                            end
                        end
                    end
                end
                for lv=MaxLV:-1:1
                    MLKP.LatentVar=lv;
                    RmseTrain=zeros(1,NvarY);
                    RmseTest=zeros(1,NvarY);
                    CorrTrain=zeros(1,NvarY);
                    CorrTest=zeros(1,NvarY);
                    for ival=1:Nval
                        Xtr=PLS(ival).Xtr;
                        Xtest=PLS(ival).Xtest;
                        Ytr=CV(ival).Ytr;
                        Ytest=CV(ival).Ytest;
                        [B,Xloadings,MLKP] = PlsTrain(Xtr, Ytr, MLKP);
                        MLK.BCoeff = B;
                        YpredTrain = PlsPred(Xtr, Ytr, B, MLKP);
                        YpredTest = PlsPred(Xtest, Ytest, B, MLKP);
                        [MLK] = EvaluateOutput(Ytr,YpredTrain,Ytest,YpredTest,1,MLK,MLKP);
                        if (upper(MLKP.ProblemType) == 'CLS')
                            PercCorrTrain(ival)=MLK.PercCorrTrain;
                            PercCorrTest(ival)=MLK.PercCorrTest;
                        else
                            PercCorrTrain(ival)=mean(MLK.CorrTrain);
                            PercCorrTest(ival)=mean(MLK.CorrTest);
                            RmseTrain=RmseTrain+MLK.RmseTrain;
                            RmseTest=RmseTest+MLK.RmseTest;
                            CorrTrain=CorrTrain+MLK.CorrTrain;
                            CorrTest=CorrTest+MLK.CorrTest;
                        end
                    end
                    if (upper(MLKP.ProblemType) == 'CLS')
                        PCTrain=mean(PercCorrTrain);
                        PCTest=mean(PercCorrTest);
                    else
                        PCTrain=-mean(RmseTrain);
                        PCTest=-mean(RmseTest);
                    end
                    if(PCTest > OldPCTest)
                        OldPCTest=PCTest;
                        OPT.DistTypeX=DistTypeX;
                        OPT.Sigma=sigma(bulent);
                        OPT.Omega=omega(iomega);
                        OPT.Power=bulent;
                        OPT.B0=b0;
                        OPT.LV=MLKP.LatentVar;
                        NBestMLK.PercCorrTrain=PCTrain;
                        NBestMLK.PercCorrTest=PCTest;
                        if (upper(MLKP.ProblemType)=='REG')
                            NBestMLK.RmseTrain=RmseTrain/Nval;
                            NBestMLK.CorrTrain=CorrTrain/Nval;
                            NBestMLK.RmseTest=RmseTest/Nval;
                            NBestMLK.CorrTest=CorrTest/Nval;
                        end
                    end
                end
            end
        end
    end
end

Message=sprintf('\nOverall result with %s - %s with %d LVs',...
                 MapMethod, MLKP.TrainingMethod,OPT.LV);
disp(Message);
MeanPCTrain=NBestMLK.PercCorrTrain;
MeanPCTest=NBestMLK.PercCorrTest;
if (upper(MLKP.ProblemType) == 'CLS')
    Message=sprintf('Training: Percentage correct %g', MeanPCTrain);
    disp(Message);
    Message=sprintf('Test:     Percentage correct %g', MeanPCTest);
    disp(Message);
else
    for ivar=1:NvarY
        Message=sprintf('Training: RMSE %g   Corr %g',...
                        NBestMLK.RmseTrain(ivar),NBestMLK.CorrTrain(ivar));
        disp(Message);
        Message=sprintf('Test:     RMSE %g   Corr %g\n',...
                        NBestMLK.RmseTest(ivar),NBestMLK.CorrTest(ivar));
        disp(Message);
    end
end