ICME 2012 POT-GEN-4

< Homework 1

Author(s): Laalitha Liyanage

License

By using the codes provided here you accept the the Mississippi State University's license agreement. Please read the agreement carefully before usage.

Overview

The following script takes the sampled space generated from script-3 and builds analytical models which represents the parameter space. No modification of the script is needed except the name of the database generated in script-3.

Running the script

This should run fast on a workstation. Save the file as any_name.m. Do

swsetup matlab

and then to execute script do

matlab -nodisplay -nosplash < any_name.m

Ouput of the script

The script will output R² values. These should be close to 1.0. If not note down those target properties (same order as your goal array) and weight them zero in the next step (script-5)

Analytical model generation script

%% Analytical moodel generation

clear all
load('Al_LHS_v1.mat')

% Generate cubic model
clear cub
cub(1,1:size(X,2)) = 0;
count = 1;
for i = 1:size(X,2)
    count = count + 1;
    cub(count,i) = 1;
end
models{1} = cub; % Linear model

for i = 1:size(X,2)
    count = count + 1;
    cub(count,i) = 2;
end
models{2} = cub; % Pure quadratic model

cub = models{1};
count = size(cub,1);
for i = 1:size(X,2)-1
    for j = i+1:size(X,2)
        count = count + 1;
        cub(count,i) = 1;
        cub(count,j) = 1;
    end
end
models{3} = cub; % Quadratic interactions model

for i = 1:size(X,2)
    count = count + 1;
    cub(count,i) = 2;
end
models{4} = cub; % Full quadratic model

for i = 1:size(X,2)-1
    for j = i+1:size(X,2)
        count = count + 1;
        cub(count,i) = 2;
        cub(count,j) = 1;
        count = count + 1;
        cub(count,i) = 1;
        cub(count,j) = 2;
    end
end
for i = 1:size(X,2)-2
    for j = i+1:size(X,2)-1
        for k = j+1:size(X,2)
            count = count + 1;
            cub(count,i) = 1;
            cub(count,j) = 1;
            cub(count,k) = 1;
        end
    end
end
models{5} = cub; % Quadratic with cubic interactions

for i = 1:size(X,2)
    count = count + 1;
    cub(count,i) = 3;
end
models{6} = cub; % Full cubic model

for i = 1:size(Y,2)

    % Trying another way
    y = Y(:,i);
    x = X; % don't add column of ones for this

    rsquare = 0;

    for j = 1:length(models)
        if size(models{j},1) < size(y,1)
            stats = regstats(y,x,models{j});
            % Calculate PRESS
            r = stats.r;
            hii = diag(stats.hatmat);
            PRESS = sum((r./(1-hii)).^2);

            % Calculate the R^2 of the prediction (Montgomery, p. 419)
            SS_total = stats.fstat.sse + stats.fstat.ssr;
            Rsquare = 1-PRESS/SS_total;
            if Rsquare >= rsquare
                rsquare = Rsquare;
                model = models{j};
                beta = stats.beta;
            end
        end
    end
    disp(['Response: ' num2str(i)  '   r2: ' num2str(rsquare)])

    save(['objfun' num2str(i) '.mat'],'models','beta')

end