HW1 Script 4

Script 4

Analytical_model_generation.m

%% Analytical moodel generation

clear all
load('Al_LHS_v1.mat')
fID_Progress = fopen('Progress_Script', 'a');
fprintf(fID_Progress, 'Analytical_Model_Generation_Begin\n');
tic
% 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
fprintf(fID_Progress, '\t%f\n', toc);

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
fprintf(fID_Progress, '\t%f\n', toc);
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

fprintf(fID_Progress, '\t%f\n', toc);
fprintf(fID_Progress, 'Analytical_Model_Generation_End\n');

fclose('all');