回歸分析--多元回歸
介紹一下多元回歸分析中的統計量
案例分析及python實踐
# 導入相關包import pandas as pdimport numpy as npimport mathimport scipyimport matplotlib.pyplot as pltfrom scipy.stats import t
# 構建數據columns = {'A':"分行編號", 'B':"不良貸款(億元)", 'C':"貸款餘額(億元)", 'D':"累計應收貸款(億元)", 'E':"貸款項目個數", 'F':"固定資產投資額(億元)"}data={"A":[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25], "B":[0.9,1.1,4.8,3.2,7.8,2.7,1.6,12.5,1.0,2.6,0.3,4.0,0.8,3.5,10.2,3.0,0.2,0.4,1.0,6.8,11.6,1.6,1.2,7.2,3.2], "C":[67.3,111.3,173.0,80.8,199.7,16.2,107.4,185.4,96.1,72.8,64.2,132.2,58.6,174.6,263.5,79.3,14.8,73.5,24.7,139.4,368.2,95.7,109.6,196.2,102.2], "D":[6.8,19.8,7.7,7.2,16.5,2.2,10.7,27.1,1.7,9.1,2.1,11.2,6.0,12.7,15.6,8.9,0,5.9,5.0,7.2,16.8,3.8,10.3,15.8,12.0], "E":[5,16,17,10,19,1,17,18,10,14,11,23,14,26,34,15,2,11,4,28,32,10,14,16,10], "F":[51.9,90.9,73.7,14.5,63.2,2.2,20.2,43.8,55.9,64.3,42.7,76.7,22.8,117.1,146.7,29.9,42.1,25.3,13.4,64.3,163.9,44.5,67.9,39.7,97.1] }df = pd.DataFrame(data)X = df[["C", "D", "E", "F"]]Y = df[["B"]]
# 構建多元線性回歸模型from sklearn.linear_model import LinearRegressionlreg = LinearRegression()lreg.fit(X, Y)x = Xy_pred = lreg.predict(X)y_true = np.array(Y).reshape(-1,1)coef = lreg.coef_[0]intercept = lreg.intercept_[0]
# 自定義函數def log_like(y_true, y_pred): """ y_true: 真實值 y_pred:預測值 """ sig = np.sqrt(sum((y_true - y_pred)**2)[0] / len(y_pred)) # 殘差標準差δ y_sig = np.exp(-(y_true - y_pred) ** 2 / (2 * sig ** 2)) / (math.sqrt(2 * math.pi) * sig) loglik = sum(np.log(y_sig)) return loglikdef param_var(x): """ x:只含自變量寬表 """ n = len(x) beta0 = np.ones((n,1)) df_to_matrix = x.as_matrix() concat_matrix = np.hstack((beta0, df_to_matrix)) # 矩陣合併 transpose_matrix = np.transpose(concat_matrix) # 矩陣轉置 dot_matrix = np.dot(transpose_matrix, concat_matrix) # (X.T X)^(-1) inv_matrix = np.linalg.inv(dot_matrix) # 求(X.T X)^(-1) 逆矩陣 diag = np.diag(inv_matrix) # 獲取矩陣對角線,即每個參數的方差 return diagdef param_test_stat(x, Se, intercept, coef, alpha=0.05): n = len(x) k = len(x.columns) beta_array = param_var(x) beta_k = beta_array.shape[0] coef = [intercept] + list(coef) std_err = [] t_Stat = [] P_value = [] t_intv = [] coefLower = [] coefupper = [] for i in range(beta_k): se_belta = np.sqrt(Se**2 * beta_array[i]) # 回歸係數的抽樣標準誤差 t = coef[i] / se_belta # 用於檢驗回歸係數的t統計量, 即檢驗統計量t p_value = scipy.stats.t.sf(np.abs(t), n-k-1)*2 # 用於檢驗回歸係數的P值(P_value) t_score = scipy.stats.t.isf(alpha/2, df = n-k-1) # t臨界值 coef_lower = coef[i] - t_score * se_belta # 回歸係數(斜率)的置信區間下限 coef_upper = coef[i] + t_score * se_belta # 回歸係數(斜率)的置信區間上限 std_err.append(round(se_belta, 3)) t_Stat.append(round(t,3)) P_value.append(round(p_value,3)) t_intv.append(round(t_score,3)) coefLower.append(round(coef_lower,3)) coefupper.append(round(coef_upper,3)) dict_ = {"coefficients":list(map(lambda x:round(x, 4), coef)), 'std_err':std_err, 't_Stat':t_Stat, 'P_value':P_value, 't臨界值':t_intv, 'Lower_95%':coefLower, 'Upper_95%':coefupper} index = ["intercept"] + list(x.columns) stat = pd.DataFrame(dict_, index=index) return stat
# 自定義函數(計算輸出各回歸分析統計量)def get_lr_stats(x, y_true, y_pred, coef, intercept, alpha=0.05): n = len(x) k = len(x.columns) ssr = sum((y_pred - np.mean(y_true))**2)[0] # 回歸平方和 SSR sse = sum((y_true - y_pred)**2)[0] # 殘差平方和 SSE sst = ssr + sse # 總平方和 SST msr = ssr / k # 均方回歸 MSR mse = sse / (n-k-1) # 均方殘差 MSE R_square = ssr / sst # 判定係數R^2 Adjusted_R_square = 1-(1-R_square)*((n-1) / (n-k-1)) # 調整的判定係數 Multiple_R = np.sqrt(R_square) # 復相關係數 Se = np.sqrt(sse/(n - k - 1)) # 估計標準誤差 loglike = log_like(y_true, y_pred)[0] AIC = 2*(k+1) - 2 * loglike # (k+1) 代表k個回歸參數或係數和1個截距參數 BIC = -2*loglike + (k+1)*np.log(n) # 線性關係的顯著性檢驗 F = (ssr / k) / (sse / ( n - k - 1 )) # 檢驗統計量F (線性關係的檢驗) pf = scipy.stats.f.sf(F, k, n-k-1) # 用於檢驗的顯著性F,即Significance F Fa = scipy.stats.f.isf(alpha, dfn=k, dfd=n-k-1) # F臨界值 # 回歸係數的顯著性檢驗 stat = param_test_stat(x, Se, intercept, coef, alpha=alpha) # 輸出各回歸分析統計量 print('='*80) print('df_Model:{} df_Residuals:{}'.format(k, n-k-1), '\n') print('loglike:{} AIC:{} BIC:{}'.format(round(loglike,3), round(AIC,1), round(BIC,1)), '\n') print('SST:{} SSR:{} SSE:{} MSR:{} MSE:{} Se:{}'.format(round(sst,4), round(ssr,4), round(sse,4), round(msr,4), round(mse,4), round(Se,4)), '\n') print('Multiple_R:{} R_square:{} Adjusted_R_square:{}'.format(round(Multiple_R,4), round(R_square,4), round(Adjusted_R_square,4)), '\n') print('F:{} pf:{} Fa:{}'.format(round(F,4), pf, round(Fa,4))) print('='*80) print(stat) print('='*80) return 0
輸出結果如下:
對比statsmodels下ols結果:
參考資料:
【1】https://www.zhihu.com/question/328568463
【2】https://blog.csdn.net/qq_38998213/article/details/83480147