模拟：fMRI中的GLM模型到底要不要加入所有刺激的onset和duration

%% GLM分析
%SPM glm的分析原理
%输入: onset duration condition(/stimulus)
%输出: 估计的beta(分条件)
%
%过程: 
%    > 通过onset和duration估计出每个condition可能的bold响应——也就是X,
%    > 通过onset和duration提取出观测到的bold响应——也就是Y (非必须)
%    > 通过GLM得出估计出的beta——也就是每个condition对Y的影响的强度
% 
% SPM的GLM的特别之处:
%    不知道真实的X, X是由设计矩阵卷积HRF得到的.
%
% 探索的结论:
%     > 尽可能在设计矩阵中加入所有condition

clear
%% data 总共5个条件
% 无法观测到的数据
B = [1;0;0.5;-2;0] % 真正的beta(理论上存在,实际上无法观测)
X = randn(100,5);  % X的每一列相当于一个condition/stimulus所引起的bold响应
                   % (此数据理论上存在,实际上无法观测)
                   % SPM建模中的X是由设计矩阵卷积得到的,与真实的X可能存在差距
% 可观测数据
BOLD = X*B + randn(100,1)*.5  % 所有stimulus的bold响应合并起来, 得到观测到的bold响应

% 
conditions = [1 3 5]
X1= X(:,conditions)     % 只在模型中加入一个condition的情况                   

%% 建模过程
% 最完善模型
[b,~,~] = glmfit(X,BOLD,'normal') % 
% mean(stats.resid)
% 实践中可能出现的模型
[b,~,~] = glmfit(X1,BOLD,'normal')
% mean(stats.resid)
B(conditions)


%% 迭代平均
% 用迭代算beta平均值来比较一下
for ii = 1:100
[b,~,~] = glmfit(X,BOLD,'normal') ;% 

B0(:,ii) = b(conditions+1);

[b,~,~] = glmfit(X1,BOLD,'normal');

B1(:,ii) = b(2:end);

end

B0m = mean(B0,2)
B1m = mean(B1,2)