westoshy
3/14/2017 - 12:50 PM

線形回帰分析

線形回帰分析

% http://home.hiroshima-u.ac.jp/ete131/index.cgi?page=Matlab_Excel_Regression

% 被説明変数ベクトルの作成
y = [1 3 4 2 1 0 4]'; 

% 標本数の取得
T = size(y, 1);

% 説明変数行列の作成
iota = ones(T, 1);
x = [2 4 3 0 0 1 2]';
X = [iota x];

% 説明変数の数の取得
k = size(X, 2);

% 回帰係数ベクトルのOLS推定値ベクトルの計算
b = (X'*X) \ (X'*y);
b1 = b(1);
b2 = b(2);

% 残差ベクトルの計算
e = y - X*b;

% TSS, RSS, ESSの計算
TSS = (y - mean(y)*iota)'*(y - mean(y)*iota);
RSS = e'*e;
ESS = TSS - RSS;

% 分散の不偏推定値の計算
s2 = RSS/(T - k);

% 標準誤差の計算 
s = sqrt(s2);

% OLS推定量の分散共分散行列の計算
cov = s2*inv(X'*X);

% 標準誤差の計算
s_b1 = sqrt(cov(1, 1));
s_b2 = sqrt(cov(2, 2));

% t値の計算
t_b1 = b(1) / s_b1;
t_b2 = b(2) / s_b2;

% P値の計算
p_b1 = tcdf(-abs(t_b1), T - k) + 1 - tcdf(abs(t_b1), T - k);
p_b2 = tcdf(-abs(t_b2), T - k) + 1 - tcdf(abs(t_b2), T - k);

% 信頼係数0.95の信頼区間の計算
tast = tinv(0.975, T-k);
cilb_b1 = b1 - tast*s_b1;
ciub_b1 = b1 + tast*s_b1;
cilb_b2 = b2 - tast*s_b2;
ciub_b2 = b2 + tast*s_b2;

% 決定係数,重相関係数,自由度調整済決定係数の計算 
R2 = 1 - RSS / TSS;
R = sqrt(R2);
adjR2 = 1 - (RSS/(T-k)) / (TSS/(T-1));

% F値の計算
F = ( (TSS - RSS)/(k-1) ) / ( RSS/(T-k) );
p_F = 1 - fcdf(F, k-1, T-k);

% 結果出力の準備
out_reg = [k-1, ESS, ESS/(k-1), F, p_F];
out_rss = [T-k, RSS, RSS/(T-k)];
out_tol = [T-1, TSS]; 
out_b1 = [b1, s_b1, t_b1, p_b1, cilb_b1, ciub_b1];
out_b2 = [b2, s_b2, t_b2, p_b2, cilb_b2, ciub_b2];

% -------------------- EXCEL 分析ツール風の出力 ----------------------   %

fprintf('\n')
fprintf('----------------------------------\n')
fprintf('     回帰統計\n')
fprintf('----------------------------------\n')
fprintf('重相関R\t %8.6f \n', R)
fprintf('重決定R2\t %8.6f \n', R2)
fprintf('補正R2\t %8.6f \n', adjR2)
fprintf('標準偏差\t %8.6f \n', s)
fprintf('観測数\t\t %d \n', T)
fprintf('----------------------------------\n')
fprintf('\n')
fprintf('分散分析表\n')
fprintf('-------------------------------------------------\n')
fprintf('     自由度   変動   分散  観測された分散比  有意F\n')
fprintf('-------------------------------------------------\n')
fprintf('回帰\t %d %9.6f %9.6f %9.6f %9.6f\n', out_reg)
fprintf('残差\t %d %9.6f %9.6f\n', out_rss)
fprintf('合計\t %d %9.5f \n', out_tol)
fprintf('-------------------------------------------------\n')
fprintf('\n')
fprintf( '-----------------------------------------------------------------\n')
fprintf('          係数      標準誤差      t       P-値     下限95%   上限95%\n')
fprintf( '-----------------------------------------------------------------\n')
fprintf('切片\t %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f \n', out_b1)
fprintf('X値1\t %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f \n', out_b2)
fprintf( '-----------------------------------------------------------------\n')
fprintf('\n')