function [index_hat,cc,rece,year_com,expl,pc1_expl_percent,expl_ite,rece_ite] = pcreg4(sampWithDate,indexWithDate,calYear,valYear,delCol,...
    selPCs,selectPC2,regressMth,stepCriterion,lmmodType,r_weight0,r_p)
% selPCs: 'expl80','first1','eig1'
% selectPC2: [0.2, 0.1] (pc与index相关绝对值达0.2，显著度小于0.1）
% stepCriterion: 'SSE','aic''rsquared'
% lmmodType: 'linear','purequadratic'
% delMode（'fast','low'）
%        删除样本模式，fast为delChronSeqWithStep函数，step指定步长，结果的每个cell包含前一个cell的所有值，
%        low为delChronSeq函数结果，逐个删除样本，结果的每个cell不包含前一个cell的所有值，
% calYear = [1921,1970];
% valYear = [1871,1920];

%% 1. 分离年序列 与 样本 
samp_year = sampWithDate(:,1);   % 507-2013共1507年
index_year = indexWithDate(:,1);
samp_0 = sampWithDate(:,2:end);
index = indexWithDate(:,2:end);

%%  2. 初始化
% delcol nest过程删除列的序号，由delChronSeqWithStep()或delChronSeq()计算
% 对于非nest过程，delCol为空
if ~isempty(r_weight0), r_weight0 = abs(r_weight0(:)); end
if isempty(stepCriterion), stepCriterion = 'aic'; end
if isempty(lmmodType), lmmodType = 'linear'; end

cc = []; 
index_hat = nan(size(samp_0,1),length(delCol)+1);   % 储存指数预测结果
delCol_vec = [];
npc = 10; %只计算前几个pc

%% 3. PCA
for ii = 1:length(delCol)+1    
% if mod(ii,50) == 0, disp(ii); end   % 输出计算进程

% 3.1 进入一次nest, 删除一部分样本（列），删除后样本储存到samp中用于计算，
    samp = samp_0;
    if ~isempty(r_weight0), r_weight = r_weight0; end
    % 3.1.1 第一次不删样本，第2次以后开始删除样本
    if ii>=2   
        samp(:,delCol{ii-1}) = [];
        if ~isempty(r_weight0), r_weight(delCol{ii-1}) = []; end
    end

% 3.2 选择 年表 共有时段做pca
    % 3.2.1 选择 年表 共有时段。 共同段的行号idx、年序列year_com、样品矩阵samp_com
    idx = commonSeg(samp);   % 年表 共同时段在 samp 中的行索引idx
    year_com = samp_year(idx);  % 年表 共同时段的year序列    
    samp_com = samp(idx,:);
    % 3.2.2 共同时段的samp标准化,并计算PCA
    % samp_com = zscore(samp_com);
    samp_com = mapminmax(samp_com',0,1);
    samp_com = samp_com';
    %
    if ~isempty(r_weight0)
        % flag = samp_com<0;
        samp_com = bsxfun(@power,samp_com,(r_weight.^r_p)');
        % samp_com(flag) = samp_com(flag).*(-1);
    end
    [eofs,pcs,expl] = pca(samp_com,'VariableWeights','variance');  % 'NumComponents',8  ,
    t = cumsum(expl)./sum(expl);
    pc1_expl_percent(ii,1) = t(1);
% % **************** 筛选PCs条件**************
% 对于260个以上样本，'expl'选择的更少，随着nest继续，样本量减少，减小到小于260时,'eig1'选择的更少。
% 样本量大时，'eig1'选择得过多，如997个样本，选择129个，（expl选择66个），造成所选的后面主成分波动小，共线高，回归时误差大。
% 4. 选择主成分pc
selPCs = char(selPCs);
selPCs_ab = selPCs(1:3);
switch selPCs_ab
    case 'fir' % 标准二：第一主成分
        n_pc = str2num(selPCs(6:end));
    case 'exp'  % 标准一：解释度>80%
        expl = cumsum(expl)./sum(expl);
        n_pc = find(expl> str2num(selPCs(5:end))./100,1);
    case 'eig' % 标准三：特征值大于1
        n_pc = find(expl>=1,1,'last');   % pcs(:,n_pc) = pcs(:,n_pc)./[sqrt(expl(n_pc))]';
end
n_pc = 1:min(size(pcs,2),n_pc);
% *******************************************

% 确定储存 预测结果/第一主成分，在矩阵index_hat中的起始行
    idx_store_beg = find(samp_year == year_com(1));  
%     pc1(idx_store_beg: idx_store_beg+length(year_com)-1, ii) = pcs(:,1);  %  储存第一主成分到pc1中
    
% 主成分 与 气候指数 的共同时段， 求所选的主成分与气候指数的相关系数cc
    year_end = min([calYear(end),year_com(end),index_year(end)]);
    lm_x = pcs(year_com>=calYear(1) & year_com<=year_end,  n_pc);   %pcs每行对应year_com
    lm_y = index( index_year>=calYear(1) & index_year<=year_end );   % 线性回归y
    [pcs_cc,pcs_cc_sig] = corr(lm_x,lm_y);  
    % 这两行输出 第1主成分与指数的相关系数及pvalue
    cc(ii,1) = pcs_cc(1);
    cc(ii,2) = pcs_cc_sig(1);
% （Optional）根据pcs与index的相关系数，继续筛选PCs,选择与index相关显著的 (如  'expl' ，'eig1'时）
if ~isempty(selectPC2)
    n_pc = n_pc(abs(pcs_cc)>= selectPC2(1) &  pcs_cc_sig <= selectPC2(2) ) ; %n_pc(pcs_cc_sig<0.05 & abs(pcs_cc)>0.2)
    lm_x = pcs(find(year_com==calYear(1)) : find(year_com==year_end),  n_pc);
end


%%
%    
% %主成分 与 气候指数 的共同时段 做回归，求出回归参数lm_par
% % 方案一：原方案 线性回归
% %     
%     lm_mod = fitlm(lm_x,lm_y);  % [lm_par,~,~,~,stats] = regress(lm_y,regress_x(lm_x));   % 线性回归
%     lm_par = table2array(lm_mod.Coefficients(:,1));
%     stats = anova(lm_mod,'summary');  % pvalue
%     cc(ii,7:13) = [lm_mod.Rsquared.Ordinary,lm_mod.Rsquared.Adjusted,stats{2,5},...
%         [lm_mod.ModelCriterion.AIC,lm_mod.ModelCriterion.AICc,...
%         lm_mod.ModelCriterion.BIC,lm_mod.ModelCriterion.CAIC]];  %r2,r2adjust,p,  AIC,AICc,BIC,CAIC
%     cc(ii,14) = length(n_pc); % 选择的主成分个数
% % 根据回归参数lm_par，求出指数预测值 lm_yhat
%     lm_yhat = [ones(size(pcs,1),1),pcs(:,n_pc)]*lm_par;   % 回归预测结果
%     index_hat(idx_store_beg: idx_store_beg+length(year_com)-1, ii) = lm_yhat;  % 储存指数的预测值到index_hat  

% debug绘图
%     plot(samp_year,index_hat(:,1)); hold on
%     plot(index_year,index);
if strcmp(regressMth,'linear')
    [lm_par,~,~,~,stats] = regress(lm_y, [ones(size(lm_x,1),1),lm_x] );
    lm_yhat = [ones(size(pcs,1),1),pcs(:,n_pc)]*lm_par;   % 回归预测结果
    index_hat(idx_store_beg: idx_store_beg+length(year_com)-1, ii) = lm_yhat;  % 预测结果储存指数的预测值到index_hat 
    cc(ii,7:8) = stats([1,3]);   % 回归预测calYear阶段的相关系数和pvalue
    
% % % 方案二：逐步回归 （速度慢）
elseif strcmp(regressMth,'step')
    lm_mod = stepwiselm(lm_x,lm_y,'Criterion',stepCriterion, 'Upper',lmmodType );  % ,'Verbose',0
    %'Criterion','aic''rsquared' 
    % 'Upper','linear','purequadratic'    
    % 'Verbose',0 (0-2) 不显示信息到显示详细信息
    lm_yhat = predict(lm_mod,pcs(:,n_pc));
    index_hat(idx_store_beg: idx_store_beg+length(year_com)-1, ii) = lm_yhat;
    cc(ii,7) = lm_mod.Rsquared.Ordinary.^0.5;
    cc(ii,8) = lm_mod.NumPredictors;  % 使用了多少pc
    if ~isempty(selectPC2)
        if selectPC2(2)<1 || selectPC2(1)>0
            disp('           ↑ 存在根据r或p对pcs进行二次筛选，stepwise给出的"Adding X_"不是真正的序号！')
        end
    end
end   
    
% 其他统计量      
    cc(ii,3:4) = [year_com(1),year_com(end)];  % 开始结束年
    cc(ii,5) = length(year_com);   %  总年数
    cc(ii,6) = size(samp,2); % 样本量 

    % 调均值方差，暂时一列
    t = index_hat(samp_year>=calYear(1) & samp_year<=year_end,ii);
    index_hat(:,ii) = (index_hat(:,ii) - mean(t))./std(t).*std(lm_y)+mean(lm_y);
      
    % 检验
    if ~isempty(valYear)
        val_x = index_hat( samp_year>=valYear(1) & samp_year<=valYear(end), ii);
        val_y = index( index_year>=valYear(1) & index_year<=valYear(end) );
        [cc(ii,8),cc(ii,9)] = corr(val_x,val_y); 
    end
    

end

% 每次迭代的预测值的expl
r = corr(index_hat(samp_year>=calYear(1) & samp_year<=calYear(end),:), ...
    index(index_year>=calYear(1) & index_year<=calYear(end)));
expl_ite = r.^2;

% 每次迭代的rece
if ~isempty(valYear)
    rece_ite = nan(length(delCol)+1,2);
    for ii = 1:length(delCol)+1
        yr1 = samp_year(find(~isnan(index_hat(:,ii)),1,'last'));
        yr2 = indexWithDate(find(~isnan(indexWithDate(:,2)),1,'last'),1);
        rece_ite(ii,1:2) = RECE([samp_year,index_hat(:,ii)],indexWithDate,[calYear(1),min(yr1,yr2)],valYear) ;
    end
else
    rece_ite = [];
end

% 合并重建
% index_hat = nanmean(index_hat,2);
index_hat = spliceRecon(index_hat);

if ~isempty(valYear)
    yr1 = samp_year(find(~isnan(index_hat),1,'last'));
    yr2 = indexWithDate(find(~isnan(indexWithDate(:,2)),1,'last'),1);
    rece = RECE([samp_year,index_hat],indexWithDate,[calYear(1),min(yr1,yr2)],valYear) ;
else
    rece = [];
end

r = corr(index_hat(samp_year>=calYear(1) & samp_year<=calYear(end)), ...
    index(index_year>=calYear(1) & index_year<=calYear(end)));
expl = r.^2;

% % 预测值加1列年份，并移除nan值
% index_hat = [samp_year,index_hat];
% index_hat = index_hat(~isnan(index_hat(:,2)),:);


% 将预测值 与 器测值相减 求 残差resid
%     com_year = intersect(samp_year,index_year);  % 年表与指数的共同的year序列
%     posi_samp = find(ismember(samp_year,com_year));   % 共同时段 在年表中的位置 
%     posi_index = find(ismember(index_year,com_year));  % 共同时段 在指数中的位置
%     resid = index_hat( posi_samp, :) - index(posi_index);
%     % 去除残差全为nan的年，如后面年份并没有预测值
%     logi_onePred = ~all(isnan(resid),2);   % year With At Least One Predictor
%     com_year = com_year(logi_onePred);
%     resid = resid(logi_onePred,:);
% %
% RMSE_eachYear = sqrt(nanmean(resid.^2,2));
% RMSE_medi = median(RMSE_eachYear);