粒度端元分析的matlab程序m文件

% % End-member modelling script 
% % by E. Dietze A. Borchers & K. Hartmann （2009/2010）
% 
% % Preconditions of input data: 
% % grain size classes as variables in columns withheader row in 1st row 
% % samples in rows with header column as 1st column
% % optional water depth （or e.g. depth in core） in 2nd column
% 
% % Locations of necessary parameters
% % 1. Constant sum value: c line 34
% % 2. Percentile weight for scaling and transformation: l line 43
% % 3. Loop （9A） to find optimal EMmodel: start line 90 end line 202
% % Alternative: number of end-members for manual choice （9B）: q line 94. 
% % If 9A is active put loop lines 90 and 202 in comment（%）（also lines 204
% % to 231） and remove comment signs from line 234 on 
% % 4. Visualisation steps （no. 20 and 21 to 22） only reasonable for loop
% % routine （9A） and manually chosen end-member number （9B） respectively%% 
clear all
data = load（‘WLG Surface sample GS 001.txt‘）; % data import including labels
X = data（2:end 3:end）;
[mn] = size（X）; % size of the matrix
% X: （mxn） matrix with original data
% m rows: observations （e.g. sample ID/ sites/ depths）
% n columns: variables （e.g. phi classes of grain sizes）
 
% for labelling
phi = data（13:end）‘; % variables
dcs = data（2:end1）; % observations
depth = data（2:end2）; % additional information from sampling
 
% Scaling of data to constant sum
c = 100;
% c: scalar value of constant sum         
 
for j = 1:m
    X（j:） = X（j:） / sum（X（j:）‘） * c; 
end
% X: matrix with data scaled to constant sum
 
%% 2. Weight transformation
l = 10;
% l: percentile weight after Manson & Imbrie （1964） Klovan & Imbrie （1971） 
 for i = 1:m
    for j = 1:n
        h = prctile（X（:j） l）;
        g = prctile（X（:j） 100-l）;
        W（ij） = （X（ij） - h） / （g - h）;
    end
end
% W: （mxn） matrix with weight transformed data
 
%% 3. Definition of the similarity matrix
A = W‘ * W; 
% A: （mxm） matrix Gamma i.e. minor product matrix after Weltje （1997） 
 
%% 4. Extraction of the eigenspace （Reyment & J?reskog 1993）
[VD] = eig（A）; 
% V: （nxn） matrix with eigenvectors 
% D: （nxn） matrix with eigenvalues in the diagonal 
 
%% 5. Eigenvalues in column vector
L = diag（D）; 
% L: vector Lambda of size n which contains all eigenvalues from D 
 
%% 6. Calculation of eigenvalue proportions
Ln = L / sum（L）; 
% Ln: vector Lambda normalised which gives the fraction of variance 
% of the principal components
 
%% 7. Proportion of variance for decision on minimum number of eigenvectors
Lv = cumsum（flipud（Ln））; 
% Lv: vector Lambda cumulative with cumulative fractions of the 
% eigenvector?s variance concerning the absolute variance of W
 
%% 8. A plot for visualisation of Lv
% serves for evaluating the proportions and how many eigenvectors are 
% necessary for explaining a certain level of the data’s variance
 
figure（1）
tickstep = 1:1:n;
plot（tickstep Lv） xlabel（‘Number of eigenvectors‘）