% read in the data

fname=('simobs');

ysim=textread(strcat(fname,'_norm.txt'));

ysim=ysim';

% normalize

ysimmean=mean(ysim,2);

ysimStd=ysim-repmat(ysimmean,1,m);

ysimsd=std(ysimStd,0,2);

ysimStd=ysimStd./repmat(ysimsd,1,m);

% SVD and reduction

[U,S,V]=svd(ysimStd,0);

lam=diag(S).^2/sum(diag(S).^2);

lam=cumsum(lam);

pu=[];

pu=sum(lam<0.99999)+1;

Ksim=U( :,1: pu)*S(1: pu,1: pu)./sqrt(m);

There is an accompanying paper for this code (DOI 10.1198/016214507000000888 ) that discusses the computational expense of having a large number of output values being the reason for using PCA to reduce the outputs but I don't understand why the reduction would be on the observations and not the number of features. Does this make sense to anyone or can anybody point me in the right direction?