Ngày 18 tháng 11 năm 2019 Có nhiều phương pháp học máy nâng cao với độ chính xác dự đoán mạnh mẽ. Mặc dù các mô hình phức tạp có thể vượt trội so với các mô hình đơn giản trong việc dự đoán biến phản hồi, các mô hình đơn giản tốt hơn để hiểu được tác động và tầm quan trọng của từng tính năng trên một biến phản hồi. Khi nhiệm vụ trong tay có thể được mô tả bằng mô hình tuyến tính, hồi quy tuyến tính chiến thắng trên tất cả các phương pháp học máy khác trong giải thích tính năng do tính đơn giản của nó. Bài đăng này cố gắng giúp bạn hiểu về hồi quy tuyến tính trong không gian tính năng đa chiều, đánh giá độ chính xác của mô hình và cung cấp đoạn mã cho hồi quy tuyến tính nhiều trong Python. Chúng tôi sẽ sử dụng một dữ liệu mẫu trong suốt bài đăng này. Dữ liệu mẫu có liên quan đến ngành công nghiệp dầu khí. Nó có nguồn gốc từ Tiến sĩ Michael Pyrcz, Giáo sư Kỹ thuật Dầu khí tại Đại học Texas ở Austin. Dữ liệu gốc có thể được tìm thấy từ repo GitHub của anh ấy. Bạn cũng có thể sử dụng tải xuống trực tiếp hoặc truy cập trực tiếp bằng URL Pandas như bên dưới: In [1]:Danh mục> Học máy
0. Mô tả dữ liệu mẫu
import pandas as pd
file = '//aegis4048.github.io/downloads/notebooks/sample_data/unconv_MV_v5.csv'
df = pd.read_csv[file]
Out[2]:
| 0 | 1 | 12.08 | 2.92 | 2.80 | 81.40 | 1.16 | 2.31 | 4165.196191 |
| 1 | 2 | 12.38 | 3.53 | 3.22 | 46.17 | 0.89 | 1.88 | 3561.146205 |
| 2 | 3 | 14.02 | 2.59 | 4.01 | 72.80 | 0.89 | 2.72 | 4284.348574 |
| 3 | 4 | 17.67 | 6.75 | 2.63 | 39.81 | 1.08 | 1.88 | 5098.680869 |
| 4 | 5 | 17.52 | 4.57 | 3.18 | 10.94 | 1.51 | 1.90 | 3406.132832 |
| 5 | 6 | 14.53 | 4.81 | 2.69 | 53.60 | 0.94 | 1.67 | 4395.763259 |
| 6 | 7 | 13.49 | 3.60 | 2.93 | 63.71 | 0.80 | 1.85 | 4104.400989 |
| 7 | 8 | 11.58 | 3.03 | 3.25 | 53.00 | 0.69 | 1.93 | 3496.742701 |
| 8 | 9 | 12.52 | 2.72 | 2.43 | 65.77 | 0.95 | 1.98 | 4025.851153 |
| 9 | 10 | 13.25 | 3.94 | 3.71 | 66.20 | 1.14 | 2.65 | 4285.026122 |
Mô tả các tiêu đề
- Vâng: Chỉ số tốt
- POR: Độ xốp trung bình tốt [%]
- Perm: tính thấm [MD]
- AI: trở kháng accoustic [kg/m2s*10^6]
- Brittle: Tỷ lệ chống nổ [%]
- TOC: Tổng lượng carbon hữu cơ [%]
- VR: Độ phản xạ vitrinite [%]
- Prod: Sản xuất khí mỗi ngày [MCFD] - Biến phản hồiResponse Variable
Chúng tôi có sáu tính năng [POR, PERM, AI, giòn, TOC, VR] để dự đoán biến phản hồi [prod]. Dựa trên mức độ quan trọng của tính năng hoán vị được hiển thị trong Hình [1], POR là tính năng quan trọng nhất và giòn là tính năng quan trọng thứ hai.
Xếp hạng tính năng hoán vị nằm ngoài phạm vi của bài đăng này và sẽ không được thảo luận chi tiết. Tính năng đạt được đạt được với thư viện RFPIMP Python. Để biết thêm thông tin về Xếp hạng tính năng hoán vị, hãy tham khảo bài viết này: Hãy coi chừng Hiệu quả rừng ngẫu nhiên mặc định
Hình 1: Xếp hạng tính năng hoán vị
Mã nguồn cho Hình [1]
import rfpimp
import pandas as pd
import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
######################################## Data preparation #########################################
file = '//aegis4048.github.io/downloads/notebooks/sample_data/unconv_MV_v5.csv'
df = pd.read_csv[file]
features = ['Por', 'Perm', 'AI', 'Brittle', 'TOC', 'VR', 'Prod']
######################################## Train/test split #########################################
df_train, df_test = train_test_split[df, test_size=0.20]
df_train = df_train[features]
df_test = df_test[features]
X_train, y_train = df_train.drop['Prod',axis=1], df_train['Prod']
X_test, y_test = df_test.drop['Prod',axis=1], df_test['Prod']
################################################ Train #############################################
rf = RandomForestRegressor[n_estimators=100, n_jobs=-1]
rf.fit[X_train, y_train]
############################### Permutation feature importance #####################################
imp = rfpimp.importances[rf, X_test, y_test]
############################################## Plot ################################################
fig, ax = plt.subplots[figsize=[6, 3]]
ax.barh[imp.index, imp['Importance'], height=0.8, facecolor='grey', alpha=0.8, edgecolor='k']
ax.set_xlabel['Importance score']
ax.set_title['Permutation feature importance']
ax.text[0.8, 0.15, 'aegis4048.github.io', fontsize=12, ha='center', va='center',
transform=ax.transAxes, color='grey', alpha=0.5]
plt.gca[].invert_yaxis[]
fig.tight_layout[]
1. Hồi quy tuyến tính đa
Mô hình hồi quy tuyến tính có cấu trúc sau:
$$ y = \ beta_1 x_1 + \ beta_2 x_2 + \ cdots + \ beta_n x_n + \ beta_0 \ tag {1} $$
ở đâu
$ \ beta_n $
: hệ số hồi quy [trọng lượng] của tính năng $ n $ -th$n$-th feature
Mô hình hồi quy tuyến tính bivarate [có thể được hiển thị trong không gian 2D] là sự đơn giản hóa EQ [1]. Mô hình bivariate có cấu trúc sau:
$$ y = \ beta_1 x_1 + \ beta_0 \ tag {2} $$
Một bưc tranh đang gia ngan lơi noi. Chúng ta hãy cố gắng hiểu các thuộc tính của nhiều mô hình hồi quy tuyến tính với trực quan hóa. Đầu tiên, mô hình hồi quy tuyến tính 2D bivariate được hiển thị trong Hình [2], sử dụng POR làm một tính năng duy nhất. Mặc dù độ xốp là tính năng quan trọng nhất liên quan đến sản xuất khí, nhưng chỉ có độ xốp chỉ chiếm 74% phương sai của dữ liệu.
Hình 2: Mô hình hồi quy tuyến tính 2D
Mã nguồn cho Hình [2]
import pandas as pd
import matplotlib.pyplot as plt
from sklearn import linear_model
######################################## Data preparation #########################################
file = '//aegis4048.github.io/downloads/notebooks/sample_data/unconv_MV_v5.csv'
df = pd.read_csv[file]
X = df['Por'].values.reshape[-1,1]
y = df['Prod'].values
################################################ Train #############################################
ols = linear_model.LinearRegression[]
model = ols.fit[X, y]
response = model.predict[X]
############################################## Evaluate ############################################
r2 = model.score[X, y]
############################################## Plot ################################################
plt.style.use['default']
plt.style.use['ggplot']
fig, ax = plt.subplots[figsize=[8, 4]]
ax.plot[X, response, color='k', label='Regression model']
ax.scatter[X, y, edgecolor='k', facecolor='grey', alpha=0.7, label='Sample data']
ax.set_ylabel['Gas production [Mcf/day]', fontsize=14]
ax.set_xlabel['Porosity [%]', fontsize=14]
ax.text[0.8, 0.1, 'aegis4048.github.io', fontsize=13, ha='center', va='center',
transform=ax.transAxes, color='grey', alpha=0.5]
ax.legend[facecolor='white', fontsize=11]
ax.set_title['$R^2= %.2f$' % r2, fontsize=18]
fig.tight_layout[]
Làm thế nào mô hình sẽ trông như thế nào trong không gian 3D? Chúng ta hãy xem hình [3]. Do tính chất 3D của cốt truyện, nhiều lô được tạo ra từ các góc độ khác nhau. Hai tính năng [por và giòn] đã được sử dụng để dự đoán biến phản ứng prod. Với sự trợ giúp của tính năng bổ sung giòn, mô hình tuyến tính trải nghiệm độ chính xác đáng kể, hiện nắm bắt độ biến thiên của dữ liệu 93%.
Hình 3: Mô hình hồi quy tuyến tính 3D với các tính năng mạnh mẽ
Mã nguồn cho Hình [3]
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn import linear_model
from mpl_toolkits.mplot3d import Axes3D
######################################## Data preparation #########################################
file = '//aegis4048.github.io/downloads/notebooks/sample_data/unconv_MV_v5.csv'
df = pd.read_csv[file]
X = df[['Por', 'Brittle']].values.reshape[-1,2]
Y = df['Prod']
######################## Prepare model data point for visualization ###############################
x = X[:, 0]
y = X[:, 1]
z = Y
x_pred = np.linspace[6, 24, 30] # range of porosity values
y_pred = np.linspace[0, 100, 30] # range of brittleness values
xx_pred, yy_pred = np.meshgrid[x_pred, y_pred]
model_viz = np.array[[xx_pred.flatten[], yy_pred.flatten[]]].T
################################################ Train #############################################
ols = linear_model.LinearRegression[]
model = ols.fit[X, Y]
predicted = model.predict[model_viz]
############################################## Evaluate ############################################
r2 = model.score[X, Y]
############################################## Plot ################################################
plt.style.use['default']
fig = plt.figure[figsize=[12, 4]]
ax1 = fig.add_subplot[131, projection='3d']
ax2 = fig.add_subplot[132, projection='3d']
ax3 = fig.add_subplot[133, projection='3d']
axes = [ax1, ax2, ax3]
for ax in axes:
ax.plot[x, y, z, color='k', zorder=15, linestyle='none', marker='o', alpha=0.5]
ax.scatter[xx_pred.flatten[], yy_pred.flatten[], predicted, facecolor=[0,0,0,0], s=20, edgecolor='#70b3f0']
ax.set_xlabel['Porosity [%]', fontsize=12]
ax.set_ylabel['Brittleness', fontsize=12]
ax.set_zlabel['Gas Prod. [Mcf/day]', fontsize=12]
ax.locator_params[nbins=4, axis='x']
ax.locator_params[nbins=5, axis='x']
ax1.text2D[0.2, 0.32, 'aegis4048.github.io', fontsize=13, ha='center', va='center',
transform=ax1.transAxes, color='grey', alpha=0.5]
ax2.text2D[0.3, 0.42, 'aegis4048.github.io', fontsize=13, ha='center', va='center',
transform=ax2.transAxes, color='grey', alpha=0.5]
ax3.text2D[0.85, 0.85, 'aegis4048.github.io', fontsize=13, ha='center', va='center',
transform=ax3.transAxes, color='grey', alpha=0.5]
ax1.view_init[elev=28, azim=120]
ax2.view_init[elev=4, azim=114]
ax3.view_init[elev=60, azim=165]
fig.suptitle['$R^2 = %.2f$' % r2, fontsize=20]
fig.tight_layout[]
Mô hình tuyến tính 3D sẽ trông như thế nào nếu các tính năng ít mạnh hơn được chọn? Hãy chọn POR và VR làm tính năng mới của chúng tôi và phù hợp với mô hình tuyến tính. Trong hình [4] bên dưới, chúng ta thấy rằng R bình phương giảm so với hình [3] ở trên. Ảnh hưởng của hiệu suất mô hình giảm có thể được quan sát trực quan bằng cách so sánh các lô giữa của chúng; Các sơ đồ phân tán trong hình [3] có dày đặc hơn xung quanh mặt phẳng mô hình 2D so với các ô phân tán trong hình [4].
Hình 4: Mô hình hồi quy tuyến tính 3D với các tính năng yếu
Mã nguồn cho Hình [4]
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn import linear_model
from mpl_toolkits.mplot3d import Axes3D
######################################## Data preparation #########################################
file = '//aegis4048.github.io/downloads/notebooks/sample_data/unconv_MV_v5.csv'
df = pd.read_csv[file]
X = df[['Por', 'VR']].values.reshape[-1,2]
Y = df['Prod']
######################## Prepare model data point for visualization ###############################
x = X[:, 0]
y = X[:, 1]
z = Y
x_pred = np.linspace[6, 24, 30] # range of porosity values
y_pred = np.linspace[0.93, 2.9, 30] # range of VR values
xx_pred, yy_pred = np.meshgrid[x_pred, y_pred]
model_viz = np.array[[xx_pred.flatten[], yy_pred.flatten[]]].T
################################################ Train #############################################
ols = linear_model.LinearRegression[]
model = ols.fit[X, Y]
predicted = model.predict[model_viz]
############################################## Evaluate ############################################
r2 = model.score[X, Y]
############################################## Plot ################################################
plt.style.use['default']
fig = plt.figure[figsize=[12, 4]]
ax1 = fig.add_subplot[131, projection='3d']
ax2 = fig.add_subplot[132, projection='3d']
ax3 = fig.add_subplot[133, projection='3d']
axes = [ax1, ax2, ax3]
for ax in axes:
ax.plot[x, y, z, color='k', zorder=15, linestyle='none', marker='o', alpha=0.5]
ax.scatter[xx_pred.flatten[], yy_pred.flatten[], predicted, facecolor=[0,0,0,0], s=20, edgecolor='#70b3f0']
ax.set_xlabel['Porosity [%]', fontsize=12]
ax.set_ylabel['VR', fontsize=12]
ax.set_zlabel['Gas Prod. [Mcf/day]', fontsize=12]
ax.locator_params[nbins=4, axis='x']
ax.locator_params[nbins=5, axis='x']
ax1.text2D[0.2, 0.32, 'aegis4048.github.io', fontsize=13, ha='center', va='center',
transform=ax1.transAxes, color='grey', alpha=0.5]
ax2.text2D[0.3, 0.42, 'aegis4048.github.io', fontsize=13, ha='center', va='center',
transform=ax2.transAxes, color='grey', alpha=0.5]
ax3.text2D[0.85, 0.85, 'aegis4048.github.io', fontsize=13, ha='center', va='center',
transform=ax3.transAxes, color='grey', alpha=0.5]
ax1.view_init[elev=27, azim=112]
ax2.view_init[elev=16, azim=-51]
ax3.view_init[elev=60, azim=165]
fig.suptitle['$R^2 = %.2f$' % r2, fontsize=20]
fig.tight_layout[]
Chế độ xem đầy đủ của các mô hình tuyến tính được xây dựng dưới đây dưới dạng GIF. Lưu ý rằng mặt phẳng màu xanh luôn được chiếu tuyến tính, bất kể góc độ. Đây là lý do mà chúng tôi gọi đây là mô hình hồi quy "tuyến tính" nhiều. Mô hình sẽ luôn luôn là tuyến tính, bất kể tính chiều của các tính năng của bạn. Khi bạn có nhiều hơn 3 tính năng, mô hình sẽ rất khó được hiển thị, nhưng bạn có thể mong đợi rằng các mô hình tuyến tính có chiều cao cũng sẽ thể hiện xu hướng tuyến tính trong không gian tính năng của chúng.
Hình 5: Độ xốp và độ giòn mô hình tuyến tính GIF
Hình 6: Độ xốp và mô hình tuyến tính VR GIF
GIF được tạo bằng cách tạo 360 lô khác nhau được xem từ các góc khác nhau với đoạn mã sau và kết hợp thành một GIF từ IMGFLIP.
In [ ]:
for ii in np.arange[0, 360, 1]: ax.view_init[elev=32, azim=ii] fig.savefig['gif_image%d.png' % ii]
Ghi chú: Mã hóa dữ liệu - hồi quy với các biến phân loại
Hồi quy đòi hỏi các tính năng phải liên tục. Điều gì xảy ra nếu bạn có các tính năng phân loại quan trọng? Bạn có phải bỏ qua các biến phân loại và chỉ chạy hồi quy với các biến liên tục? Chúng ta có thể mã hóa các biến phân loại thành các biến số để tránh vấn đề này. Hãy xem hình dưới đây. Cấp độ tính năng ban đầu là một biến thể loại với ba loại thời gian. Điều này có nghĩa là có hệ thống phân cấp giữa các loại [ví dụ: thấp