HW1 Script 4

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

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');
```