function run_te_pipeline_validation()
%% ========= 配置 =========
rng(2025);
% 被试列表
subject_folders = {'sub1','sub2','sub5','sub6','sub7','sub9','sub11', ...
'sub12','sub13','sub15','sub16','sub17','sub18','sub19', ...
'sub21','sub23','sub24'};
sessions = {'Session1','Session2','Session3'};
lang_by_session = {'Color','Japanese','Russian'};
conds = {'RScon','RSinc','Scon','Sinc'};
% 使用跨平台路径
base_outdir_WEAK = fullfile('Data','Liucc','GW','TEdata_SIM_WEAK');
base_outdir_STRONG = fullfile('Data','Liucc','GW','TEdata_SIM_STRONG');
% 计算参数（调试时可减小N_ROI）
N_ROI = 20; % 快速测试用20，正式运行改回68
p_order = 2;
TR = 0.72;
fs_neural = 2;
T_neural = 900;
N_sur = 60;
maxlag_gc = 2;
%% ========= 生成 Weak 数据 =========
weak_cfg = struct();
weak_cfg.base_outdir = base_outdir_WEAK;
weak_cfg.N_ROI = N_ROI;
weak_cfg.p_order = p_order;
weak_cfg.TR = TR;
weak_cfg.fs_neural = fs_neural;
weak_cfg.T_neural = T_neural;
weak_cfg.N_sur = N_sur;
weak_cfg.maxlag_gc = maxlag_gc;
weak_cfg.density = 0.05;
weak_cfg.w_mean = 0.25;
weak_cfg.w_std = 0.05;
weak_cfg.sigma_noise = 1.0;
weak_cfg.add_global_signal= true;
weak_cfg.global_SNR = 2.5;
weak_cfg.session_scaler = [1.00,1.03,1.03];
weak_cfg.delta_RScon = 0.03;
weak_cfg.delta_Scon = 0.03;
weak_cfg.delta_RSinc = 0.04;
weak_cfg.delta_Sinc = 0.04;
weak_cfg.jitter_A_eps = 0.08;
weak_cfg.flip_rate = 0.40;
weak_cfg.HRF_jitter = 0.12;
simulate_dataset(subject_folders, sessions, lang_by_session, conds, weak_cfg);
%% ========= 生成 Strong 数据 =========
strong_cfg = weak_cfg;
strong_cfg.base_outdir = base_outdir_STRONG;
strong_cfg.delta_RScon = 0.10;
strong_cfg.delta_Scon = 0.10;
strong_cfg.delta_RSinc = 0.12;
strong_cfg.delta_Sinc = 0.12;
strong_cfg.jitter_A_eps = 0.03;
strong_cfg.flip_rate = 0.00;
strong_cfg.global_SNR = 0.0;
strong_cfg.session_scaler = [1.00,1.05,1.05];
simulate_dataset(subject_folders, sessions, lang_by_session, conds, strong_cfg);
%% ========= 分析 =========
target_session = 'Session2';
target_lang = 'Japanese';
cond_A = 'RSinc';
cond_B = 'RScon';
alpha = 0.05;
fprintf('\n================= 分析 WEAK 数据 =================\n');
run_group_analysis(base_outdir_WEAK, subject_folders, target_session, target_lang, cond_A, cond_B, alpha);
fprintf('\n================= 分析 STRONG 数据 =================\n');
run_group_analysis(base_outdir_STRONG, subject_folders, target_session, target_lang, cond_A, cond_B, alpha);
fprintf('\n✅ 完成\n');
end
%% =======================================================================
function simulate_dataset(subject_folders, sessions, lang_by_session, conds, C)
% 确保输出目录存在
if ~exist(C.base_outdir, 'dir')
mkdir(C.base_outdir);
end
% 生成基础MVAR
A_base = gen_sparse_stable_mvar(C.N_ROI, C.p_order, C.density, C.w_mean, C.w_std);
% 定义条件差异
frontals = 1:10; parietals = 11:20; visuals = 21:30;
mod = struct();
mod.RScon = sparse(C.N_ROI,C.N_ROI);
mod.Scon = sparse(C.N_ROI,C.N_ROI);
mod.RSinc = sparse(C.N_ROI,C.N_ROI);
mod.Sinc = sparse(C.N_ROI,C.N_ROI);
mod.RScon(parietals, visuals) = C.delta_RScon;
mod.Scon( parietals, visuals) = C.delta_Scon;
mod.RSinc(parietals, frontals)= C.delta_RSinc;
mod.Sinc( parietals, frontals)= C.delta_Sinc;
% 主循环
for s = 1:numel(subject_folders)
outdir = fullfile(C.base_outdir, subject_folders{s}, 'output_with_HRF');
if ~exist(outdir, 'dir')
mkdir(outdir);
end
A_indiv = jitter_A(A_base, C.jitter_A_eps);
for ss = 1:numel(sessions)
lang = lang_by_session{ss};
for c = 1:numel(conds)
cond = conds{c};
% 应用条件差异
A_eff = apply_condition(A_indiv, mod.(cond));
if C.flip_rate > 0
A_eff = maybe_flip_task_direction(A_eff, mod.(cond), C.flip_rate);
end
A_eff = scale_spectral_radius(A_eff, 0.95);
A_eff = cellfun(@(Ak) Ak * C.session_scaler(ss), A_eff, 'uni', 0);
% 生成时间序列
Xn = simulate_mvar(A_eff, C.T_neural, C.sigma_noise);
% 添加全局信号
if C.add_global_signal && C.global_SNR > 0
gsig = bandlimited_noise(C.T_neural, 0.01, 0.1);
Xn = Xn + C.global_SNR * gsig;
end
% HRF卷积
Xb = neural_to_bold_jitterHRF(Xn, C.fs_neural, C.TR, C.HRF_jitter);
% 估计GC
[GCval, GCr] = estimate_gc_surrogate(Xb, C.maxlag_gc, C.N_sur);
% 保存结果
g = struct();
g.GCval = GCval;
g.GCr = GCr;
g.info = struct('subject', subject_folders{s}, 'session', sessions{ss}, ...
'language', lang, 'condition', cond, 'TR', C.TR, ...
'p_order', C.p_order, 'N_sur', C.N_sur);
gt = struct('A', A_eff, 'A_base', A_base, 'delta', full(mod.(cond)));
fname = sprintf('GCall_%s_%s_%s_with_HRF_Onset_ITER100.mat', sessions{ss}, lang, cond);
save(fullfile(outdir, fname), 'g', 'gt', '-v7.3');
end
end
end
fprintf('✅ 模拟数据生成完成\n');
end
%% =======================================================================
function run_group_analysis(base_outdir, subject_folders, session, lang, condA, condB, alpha)
fprintf('\n--- 单条件分析: %s ---\n', condA);
[GC_stack, GT_edges] = load_stack(base_outdir, subject_folders, session, lang, condA);
[SigA, tA, pA] = onesample_t_FDR(GC_stack, alpha, 'right');
report_stats(SigA, tA, pA, GT_edges, '单条件');
fprintf('\n--- 条件对比: %s - %s ---\n', condA, condB);
[GC_A, GT_A] = load_stack(base_outdir, subject_folders, session, lang, condA);
[GC_B, GT_B] = load_stack(base_outdir, subject_folders, session, lang, condB);
Delta = GC_A - GC_B;
[SigD, tD, pD] = onesample_t_FDR(Delta, alpha, 'right');
% 修正真值差异计算
GT_delta = abs(mean(GT_A, 3) - mean(GT_B, 3)) > 0.01;
report_stats(SigD, tD, pD, GT_delta, '条件对比');
end
%% =================== 核心功能函数 ===================
function [GC_stack, GT_stack] = load_stack(base_outdir, subject_folders, session, lang, cond)
GCs = []; GTs = [];
for s = 1:numel(subject_folders)
f = fullfile(base_outdir, subject_folders{s}, 'output_with_HRF', ...
sprintf('GCall_%s_%s_%s_with_HRF_Onset_ITER100.mat', session, lang, cond));
try
S = load(f, 'g', 'gt');
GCs(s,:,:) = S.g.GCval;
GTs(:,:,s) = A_to_single(S.gt.A);
catch
fprintf('警告: 文件加载失败 %s\n', f);
end
end
GC_stack = GCs;
GT_stack = GTs;
end
function [Sig, tmat, pmat] = onesample_t_FDR(GC_stack, alpha, tail)
[Nsub, R, ~] = size(GC_stack);
tmat = zeros(R,R);
pmat = ones(R,R);
for i = 1:R
for j = 1:R
if i == j, continue; end
x = squeeze(GC_stack(:,i,j));
[~, p, ~, stats] = ttest(x, 0, 'Tail', tail);
tmat(i,j) = stats.tstat;
pmat(i,j) = p;
end
end
mask = ~eye(R);
pvec = pmat(mask);
qvec = fdr_bh(pvec, alpha);
Sig = false(R,R);
Sig(mask) = qvec < alpha;
end
function q = fdr_bh(p, alpha)
p = p(:);
[ps, ix] = sort(p);
m = numel(p);
rank = (1:m)';
q_s = ps .* m ./ rank;
for k = m-1:-1:1
q_s(k) = min(q_s(k), q_s(k+1));
end
q = zeros(size(p));
q(ix) = q_s;
end
%% =================== 模拟工具函数 ===================
function A = gen_sparse_stable_mvar(N, p, density, w_mean, w_std)
A = cell(1,p);
for k = 1:p
mask = rand(N,N) < density;
vals = w_mean + w_std*randn(N,N);
A{k} = mask .* vals;
A{k}(eye(N)==1) = 0;
end
A = scale_spectral_radius(A, 0.95);
end
function A2 = jitter_A(A, eps_scale)
A2 = cell(size(A));
for k = 1:numel(A)
A2{k} = A{k} + eps_scale*randn(size(A{k}));
A2{k}(eye(size(A2{k},1))==1) = 0;
end
A2 = scale_spectral_radius(A2, 0.95);
end
function A2 = apply_condition(A, delta)
A2 = A;
A2{1} = A2{1} + full(delta);
A2 = scale_spectral_radius(A2, 0.95);
end
function A = scale_spectral_radius(A, rho_target)
N = size(A{1},1);
M = zeros(N);
for k = 1:numel(A)
M = M + abs(A{k});
end
e = max(abs(eig(M)));
if e > 0
sc = rho_target / e;
for k = 1:numel(A)
A{k} = A{k} * sc;
end
end
end
function A2 = maybe_flip_task_direction(A以下是针对MATLAB 2018a兼容性优化的完整修改代码，解决了函数定义、路径兼容性和关键函数缺失问题：
```matlab
function run_te_pipeline_validation()
%% ========= 配置 =========
rng(2025);
% 被试列表
subject_folders = {'sub1','sub2','sub5','sub6','sub7','sub9','sub11', ...
'sub12','sub13','sub15','sub16','sub17','sub18','sub19', ...
'sub21','sub23','sub24'};
sessions = {'Session1','Session2','Session3'};
lang_by_session = {'Color','Japanese','Russian'};
conds = {'RScon','RSinc','Scon','Sinc'};
% 兼容Windows/Mac的路径
base_outdir_WEAK = fullfile('Data', 'Liucc', 'GW', 'TEdata_SIM_WEAK');
base_outdir_STRONG = fullfile('Data', 'Liucc', 'GW', 'TEdata_SIM_STRONG');
% 计算参数
N_ROI = 20; % 快速测试模式 (正式运行改为68)
p_order = 2;
TR = 0.72;
fs_neural = 2;
T_neural = 900;
N_sur = 60;
maxlag_gc = 2;
%% ========= 生成 Weak 数据 =========
weak_cfg = struct();
weak_cfg.base_outdir = base_outdir_WEAK;
weak_cfg.N_ROI = N_ROI;
weak_cfg.p_order = p_order;
weak_cfg.TR = TR;
weak_cfg.fs_neural = fs_neural;
weak_cfg.T_neural = T_neural;
weak_cfg.N_sur = N_sur;
weak_cfg.maxlag_gc = maxlag_gc;
weak_cfg.density = 0.05;
weak_cfg.w_mean = 0.25;
weak_cfg.w_std = 0.05;
weak_cfg.sigma_noise = 1.0;
weak_cfg.add_global_signal= true;
weak_cfg.global_SNR = 2.5;
weak_cfg.session_scaler = [1.00,1.03,1.03];
weak_cfg.delta_RScon = 0.03;
weak_cfg.delta_Scon = 0.03;
weak_cfg.delta_RSinc = 0.04;
weak_cfg.delta_Sinc = 0.04;
weak_cfg.jitter_A_eps = 0.08;
weak_cfg.flip_rate = 0.40;
weak_cfg.HRF_jitter = 0.12;
simulate_dataset(subject_folders, sessions, lang_by_session, conds, weak_cfg);
%% ========= 生成 Strong 数据 =========
strong_cfg = weak_cfg;
strong_cfg.base_outdir = base_outdir_STRONG;
strong_cfg.delta_RScon = 0.10;
strong_cfg.delta_Scon = 0.10;
strong_cfg.delta_RSinc = 0.12;
strong_cfg.delta_Sinc = 0.12;
strong_cfg.jitter_A_eps = 0.03;
strong_cfg.flip_rate = 0.00;
strong_cfg.global_SNR = 0.0;
strong_cfg.session_scaler = [1.00,1.05,1.05];
simulate_dataset(subject_folders, sessions, lang_by_session, conds, strong_cfg);
%% ========= 分析 =========
target_session = 'Session2';
target_lang = 'Japanese';
cond_A = 'RSinc';
cond_B = 'RScon';
alpha = 0.05;
fprintf('\n================= 分析 WEAK 数据 =================\n');
run_group_analysis(base_outdir_WEAK, subject_folders, target_session, target_lang, cond_A, cond_B, alpha);
fprintf('\n================= 分析 STRONG 数据 =================\n');
run_group_analysis(base_outdir_STRONG, subject_folders, target_session, target_lang, cond_A, cond_B, alpha);
fprintf('\n✅ 完成\n');
end
%% =================== 核心函数 ===================
function simulate_dataset(subject_folders, sessions, lang_by_session, conds, C)
for s = 1:numel(subject_folders)
outdir = fullfile(C.base_outdir, subject_folders{s}, 'output_with_HRF');
if ~exist(outdir, 'dir')
mkdir(outdir);
end
% 生成基础稳定MVAR
A_base = gen_sparse_stable_mvar(C.N_ROI, C.p_order, C.density, C.w_mean, C.w_std);
A_indiv = jitter_A(A_base, C.jitter_A_eps);
% 定义条件差异
frontals = 1:floor(C.N_ROI/3);
parietals = (floor(C.N_ROI/3)+1):floor(2*C.N_ROI/3);
visuals = (floor(2*C.N_ROI/3)+1):C.N_ROI;
mod = struct();
mod.RScon = sparse(C.N_ROI, C.N_ROI);
mod.Scon = sparse(C.N_ROI, C.N_ROI);
mod.RSinc = sparse(C.N_ROI, C.N_ROI);
mod.Sinc = sparse(C.N_ROI, C.N_ROI);
mod.RScon(parietals, visuals) = C.delta_RScon;
mod.Scon(parietals, visuals) = C.delta_Scon;
mod.RSinc(parietals, frontals) = C.delta_RSinc;
mod.Sinc(parietals, frontals) = C.delta_Sinc;
for ss = 1:numel(sessions)
lang = lang_by_session{ss};
for c = 1:numel(conds)
cond = conds{c};
% 应用条件差异
A_eff = apply_condition(A_indiv, mod.(cond));
if C.flip_rate > 0
A_eff = maybe_flip_task_direction(A_eff, mod.(cond), C.flip_rate);
end
A_eff = scale_spectral_radius(A_eff, 0.95);
for k = 1:numel(A_eff)
A_eff{k} = A_eff{k} * C.session_scaler(ss);
end
% 生成神经信号
Xn = simulate_mvar(A_eff, C.T_neural, C.sigma_noise);
% 添加全局信号
if C.add_global_signal && C.global_SNR > 0
gsig = bandlimited_noise(C.T_neural, 0.01, 0.1);
Xn = Xn + C.global_SNR * gsig(:);
end
% HRF卷积
Xb = neural_to_bold_jitterHRF(Xn, C.fs_neural, C.TR, C.HRF_jitter);
% 估计GC
[GCval, GCr] = estimate_gc_surrogate(Xb, C.maxlag_gc, C.N_sur);
% 保存结果
g = struct();
g.GCval = GCval;
g.GCr = GCr;
g.info = struct('subject', subject_folders{s}, 'session', sessions{ss}, ...
'language', lang, 'condition', cond, 'TR', C.TR, ...
'p_order', C.p_order, 'N_sur', C.N_sur);
gt = struct('A', {A_eff}, 'A_base', {A_base}, 'delta', full(mod.(cond)));
fname = sprintf('GCall_%s_%s_%s_with_HRF_Onset_ITER100.mat', sessions{ss}, lang, cond);
save(fullfile(outdir, fname), 'g', 'gt', '-v7.3');
end
end
end
fprintf('✅ 模拟数据已生成: %s\n', C.base_outdir);
end
%% =================== 辅助函数 ===================
function A = gen_sparse_stable_mvar(N, p, density, w_mean, w_std)
A = cell(1, p);
for k = 1:p
mask = rand(N) < density;
mask = triu(mask, 1) + tril(mask, -1); % 避免对角线
M = mask .* (w_mean + w_std * randn(N));
A{k} = M;
end
A = scale_spectral_radius(A, 0.95);
end
function A_out = jitter_A(A, eps_scale)
A_out = cell(size(A));
for k = 1:numel(A)
J = eps_scale * randn(size(A{k}));
J(logical(eye(size(J,1)))) = 0; % 保持对角线为0
A_out{k} = A{k} + J;
end
A_out = scale_spectral_radius(A_out, 0.95);
end
function A_out = apply_condition(A, delta)
A_out = A;
A_out{1} = A_out{1} + full(delta);
A_out = scale_spectral_radius(A_out, 0.95);
end
function A = scale_spectral_radius(A, rho_target)
N = size(A{1}, 1);
M = zeros(N);
for k = 1:numel(A)
M = M + abs(A{k});
end
rho = max(abs(eig(M)));
if rho > 0
scale_factor = rho_target / rho;
for k = 1:numel(A)
A{k} = A{k} * scale_factor;
end
end
end
function A_out = maybe_flip_task_direction(A, delta, flip_rate)
A_out = A;
if rand > flip_rate
return;
end
[rows, cols] = find(delta);
for i = 1:numel(rows)
r = rows(i); c = cols(i);
w = delta(r, c);
A_out{1}(r, c) = A_out{1}(r, c) - w;
A_out{1}(c, r) = A_out{1}(c, r) + w;
end
end
function X = simulate_mvar(A, T, sigma)
N = size(A{1}, 1);
p = numel(A);
X = zeros(T, N);
X(1:p, :) = randn(p, N);
for t = (p+1):T
Xt = zeros(1, N);
for k = 1:p
Xt = Xt + X(t-k, :) * A{k}';
end
X(t, :) = Xt + sigma * randn(1, N);
end
end
function y = bandlimited_noise(T, fl, fh)
t = (0:T-1)';
noise = randn(T, 1);
f = fft(noise);
freq = (0:T-1)'/T;
freq(freq > 0.5) = 1 - freq(freq > 0.5);
bandpass = (freq >= fl) & (freq <= fh);
f(~bandpass) = 0;
y = real(ifft(f));
y = y - mean(y);
end
function Xb = neural_to_bold_jitterHRF(Xn, fs_neural, TR, jitter)
[T, N] = size(Xn);
dt = 1/fs_neural;
hrf = spm_hrf(dt);
t_hrf = (0:length(hrf)-1)' * dt;
% 应用抖动
sc = max(0.5, 1 + jitter * randn());
t_scaled = t_hrf * sc;
hrf_scaled = interp1(t_hrf, hrf, t_scaled, 'linear', 0);
hrf_scaled = hrf_scaled / sum(hrf_scaled);
% 卷积
Xb_full = zeros(T + length(hrf_scaled) - 1, N);
for i = 1:N
Xb_full(:, i) = conv(Xn(:, i), hrf_scaled, 'full');
end
% 降采样
ds_factor = round(TR * fs_neural);
Xb = Xb_full(1:ds_factor:end, :);
Xb = detrend(Xb);
Xb = Xb ./ std(Xb, 0, 1);
end
function hrf = spm_hrf(dt)
u = 0:dt:32;
hrf = gampdf(u, 6, 1) - gampdf(u, 16, 1)/6;
hrf = hrf / sum(hrf);
end
function y = gampdf(x, a, b)
y = (x.^(a-1) .* exp(-x/b)) / (gamma(a) * b^a);
y(x <= 0) = 0;
end
function [GCval, GCr] = estimate_gc_surrogate(X, maxlag, Nsur)
[T, N] = size(X);
GCval = zeros(N);
GCr = zeros(N);
shifts = randi([round(0.1*T), round(0.9*T)], Nsur, 1);
for i = 1:N
for j = 1:N
if i == j
continue;
end
% 真实GC
gc_real = pairwise_gc(X(:, i), X(:, j), maxlag);
GCval(i, j) = gc_real;
% 置换检验
sur_vals = zeros(Nsur, 1);
for s = 1:Nsur
Xj_sur = circshift(X(:, j), shifts(s));
sur_vals(s) = pairwise_gc(X(:, i), Xj_sur, maxlag);
end
mu_sur = mean(sur_vals);
sd_sur = std(sur_vals);
GCr(i, j) = (gc_real - mu_sur) / (sd_sur + eps);
end
end
end
function F = pairwise_gc(y, x, p)
y = y(:);
x = x(:);
T = length(y);
Y = y(p+1:end);
% 受限模型（仅y的历史）
Ry = zeros(T-p, p);
for k = 1:p
Ry(:, k) = y(p+1-k:T-k);
end
b_r = Ry \ Y;
res_r = Y - Ry * b_r;
var_r = var(res_r);
% 非受限模型（y和x的历史）
R_all = [Ry, zeros(T-p, p)];
for k = 1:p
R_all(:, p+k) = x(p+1-k:T-k);
end
b_u = R_all \ Y;
res_u = Y - R_all * b_u;
var_u = var(res_u);
F = log(var_r / var_u);
end
%% =================== 分析函数 ===================
function run_group_analysis(base_dir, subjects, session, lang, condA, condB, alpha)
fprintf('\n--- [%s | %s | %s] 单条件分析 ---\n', session, lang, condA);
[GC_stack_A, GT_edges_A] = load_stack(base_dir, subjects, session, lang, condA);
[SigA, tA, pA] = onesample_t_FDR(GC_stack_A, alpha, 'right');
report_stats(SigA, tA, pA, GT_edges_A, '单条件');
fprintf('\n--- [%s | %s | %s - %s] 条件对比 ---\n', session, lang, condA, condB);
[GC_stack_B, ~] = load_stack(base_dir, subjects, session, lang, condB);
Delta = GC_stack_A - GC_stack_B;
[SigD, tD, pD] = onesample_t_FDR(Delta, alpha, 'right');
GT_delta = abs(mean(GT_edges_A, 3) - mean(load_gt(base_dir, subjects, session, lang, condB), 3)) > 0.01;
report_stats(SigD, tD, pD, GT_delta, '条件对比');
end
function [GC_stack, GT_stack] = load_stack(base_dir, subjects, session, lang, cond)
nSubs = numel(subjects);
first = true;
for s = 1:nSubs
fpath = fullfile(base_dir, subjects{s}, 'output_with_HRF', ...
sprintf('GCall_%s_%s_%s_with_HRF_Onset_ITER100.mat', session, lang, cond));
data = load(fpath, 'g', 'gt');
if first
N = size(data.g.GCval, 1);
GC_stack = zeros(nSubs, N, N);
GT_stack = zeros(N, N, nSubs);
first = false;
end
GC_stack(s, :, :) = data.g.GCval;
GT_stack(:, :, s) = A_to_single(data.gt.A);
end
end
function GT_mean = load_gt(base_dir, subjects, session, lang, cond)
nSubs = numel(subjects);
first = true;
for s = 1:nSubs
fpath = fullfile(base_dir, subjects{s}, 'output_with_HRF', ...
sprintf('GCall_%s_%s_%s_with_HRF_Onset_ITER100.mat', session, lang, cond));
data = load(fpath, 'gt');
if first
N = size(data.gt.A{1}, 1);
GT_mean = zeros(N, N, nSubs);
first = false;
end
GT_mean(:, :, s) = A_to_single(data.gt.A);
end
end
function A1 = A_to_single(Acell)
A1 = zeros(size(Acell{1}));
for k = 1:numel(Acell)
A1 = A1 + Acell{k};
end
end
function [Sig, tmat, pmat] = onesample_t_FDR(GC_stack, alpha, tail)
[nSubs, N, ~] = size(GC_stack);
tmat = zeros(N);
pmat = ones(N);
mask = ~eye(N);
for i = 1:N
for j = 1:N
if i == j, continue; end
vals = squeeze(GC_stack(:, i, j));
[~, p, ~, stats] = ttest(vals, 0, 'Tail', tail);
tmat(i, j) = stats.tstat;
pmat(i, j) = p;
end
end
% FDR校正
pvec = pmat(mask);
[~, q] = mafdr(pvec, 'BHFDR', true); % MATLAB内置FDR
Sig = false(N);
Sig(mask) = q < alpha;
end
function report_stats(Sig, tmat, pmat, GT, tag)
N = size(Sig, 1);
n_sig = nnz(Sig);
fprintf('显著边数量: %d\n', n_sig);
if ~isempty(GT) && isequal(size(GT), [N, N])
gt_mask = GT & ~eye(N);
overlap = nnz(Sig & gt_mask);
gt_edges = nnz(gt_mask);
fprintf('真值重叠率: %.2f%% (%d/%d)\n', 100*overlap/gt_edges, overlap, gt_edges);
end
if n_sig > 0
[~, idx] = sort(tmat(:), 'descend');
topk = min(5, n_sig);
fprintf('Top %d 显著边:\n', topk);
for k = 1:topk
[i, j] = ind2sub([N, N], idx(k));
if ~Sig(i, j), continue; end
fprintf(' %2d -> %2d: t=%.2f, p=%.3g\n', i, j, tmat(i,j), pmat(i,j));
end
end
fprintf('\n');
end
这是完整的代码吗？我跑了，错误: 文件:valition1.m 行:271 列:42
文本字符无效。请检查不受支持的符号、不可见的字符或非 ASCII 字符的粘贴。
最新发布
关于符号~的探讨
博客内容提及了符号~，但信息较少，未明确关于该符号的具体信息技术相关内容。
扫一扫
1705
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
