230517 / BSA11. python에서 통계 모델링

BSA08-Regress-Whitewine.ipynb

 

패키지 호출

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Lasso, Ridge
from sklearn.metrics import mean_squared_error
from sklearn.metrics import confusion_matrix, accuracy_score
from sklearn.linear_model import LogisticRegression
from sklearn import datasets

import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.discrete.discrete_model import MNLogit

 

데이터 전처리

whitedata = pd.read_csv("white.csv")

# 결측치 확인
whitedata.isna().sum()

# 데이터 구조 파악
whitedata.head()

 

1. 단순선형모형

# 단순선형모형 : 설명변수 1개
y = whitedata.y  # whitedata["y"]와 같은 의미
x = whitedata.x1
plt.scatter(x,y,color="blue"); plt.show()  # 산점도

# whitedata["x1"] : 판다스의 시리즈를 데이터프레임 형태로 만듦
x = whitedata["x1"].to_frame()
y = whitedata.y

# simple linear regression
slm = LinearRegression().fit(x,y) 
w, b = slm.coef_[0], slm.intercept_  # w=베타1, b=베타0
print(w,b)

x0, x1 = x["x1"].min(), x["x1"].max()
print(x0, x1)

plt.scatter(x,y,color="blue")
plt.plot([x0,x1],[w*x0+b,w*x1+b],'r')  # (x0, w*x0+b), (x1, w*x1+b) 두 점을 잇는 빨간 선분 추가

# statsmodels : ordinal least square
model = sm.OLS(y, sm.add_constant(x))   # sm.OLS(y,x)와 비교
results = model.fit()   # sm.OLS(y, sm.add_constant(x)).fit()
print(results.summary())
plt.show()

 

2. 다중선형회귀분석

1) sklearn

## Multiple linear regression
X = whitedata.drop(["y"],axis=1)
y = whitedata["y"]

# 학습, 검정데이터 분리
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=100)

# 모델 적용
# n_jobs=-1 : 컴퓨터의 코어의 수를 최대로 활용 (모두 사용해서 분석)
# # fit_intercept=True (default) : 절편이 있는 모형
lrm = LinearRegression(n_jobs=-1)
result = lrm.fit(X_train, y_train)
print(result.coef_)
print(result.feature_names_in_)

# 예측
forecast = lrm.predict(X_test)
print(forecast[0:10])

# 모델 성능 평가
# y hat과 y 비교
accuracy = lrm.score(X_test, y_test)    # 결정계수
print(accuracy)

 

2) statsmodels

model = sm.OLS(y_train, sm.add_constant(X_train))
results = model.fit() 
print(results.summary())
# 유의하지 않은 변수는 제외하고 분석을 진행해야 함

 

3. LASSO, Ridge

1) sklearn

# alpha = 람다값
lasso = Lasso(alpha=0.1) 
ridge = Ridge(alpha=0.1)  # alpha=0: OLS (0보다 크면 다 보여줌)


lasso.fit(X_train, y_train)
print(lasso.intercept_ , lasso.coef_ )
# 추정한 베타값들이 0.1보다 작아서 0으로 만듦 (즉, 모형에 포함시킬 필요가 없다고 판단)
# but. 표준화해서 계산한 값이므로 실제 베타값들과는 차이가 있음


forecast = lasso.predict(X_test)
print(forecast[0:10])


# 가장 설명력이 뛰어난 람다값(alpha)을 찾기 위한 작업
# MSE가 작은 모델(설명력이 좋은)을 찾음
## training-validation procedure(반복 with train-test-split)
from sklearn.metrics import mean_squared_error


result = np.zeros((10,5))
a = [0.001, 0.002, 0.003, 0.004, 0.005]
for step in range(10):
    X0, X1, y0, y1 = train_test_split(X, y, test_size=0.3)
    for choice in range(len(a)):
        lasso = Lasso(alpha=a[choice])      # sklearn.preprocessing.StandardScaler
        lasso.fit(X0, y0); forecast = lasso.predict(X1)
        result[step,choice] = mean_squared_error(y1,forecast)
print(result)


result = pd.DataFrame(result)
result.mean()
# alpha=0.001을 선택하는 것이 적합함 (MSE가 가장 작으므로)

 

2) statsmodels

mlogit_sm = MNLogit(Y,sm.add_constant(X))
mlogit_smL = mlogit_sm.fit_regularized(method="l1", alpha=0.5)  # lasso 방법을 적용하고 싶다면
mlogit_smL.summary()

yhatL = np.argmax(mlogit_smL.predict(sm.add_constant(X)),1)
print(yhatL)

# sklearn에서는 다범주 데이터에도 로지스틱 회귀모형 적용 가능
# multinomial : 세 범주의 확률의 합이 1이 되도록
# ovr : class1, class2, class3일 확률 각각 계산 -> 각 확률의 합이 1이 되도록 보정
mlogitR = LogisticRegression(C=1e5)  # multi_class="auto" {'ovr','multinomial','auto'}  
mlogitR.fit(X,Y)

muR = mlogitR.predict_proba(X)
yhatR = np.argmax(muR,1)
print(muR[0:10])
print(yhatR)


## 그림으로 표시
h = .02  # step size in the mesh
x_min, x_max = X[:, 0].min() - .5, X[:, 0].max() + .5
y_min, y_max = X[:, 1].min() - .5, X[:, 1].max() + .5
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
Z = mlogitR.predict(np.c_[xx.ravel(), yy.ravel()])
# put the result into a color plot
Z = Z.reshape(xx.shape)
plt.figure(1, figsize=(8, 7))
plt.pcolormesh(xx, yy, Z, cmap=plt.cm.Paired)
# Plot also the training points
plt.scatter(X[:, 0], X[:, 1], c=Y, edgecolors='k', cmap=plt.cm.Paired)
plt.xlabel('Sepal length'); plt.ylabel('Sepal width')
plt.xlim(xx.min(), xx.max());  plt.ylim(yy.min(), yy.max())
plt.xticks(()); plt.yticks(())
plt.show()

 

4. 로지스틱 회귀분석

1) sklearn

logitS = LogisticRegression(penalty=None)  ## C=1/alpha=1, penality="l2" {"l1","l2","elasticnet",None}
logitS.fit(X_train, y_train)  # 별도로 베타0을 넣지 않아도 있는 것처럼 분석
# 특별히 언급하지 않으면 Ridge regression 진행
# penalty=None 옵션 지정해야 기본적인 로지스틱 모형 분석 가능

print(logitS.coef_, logitS.intercept_)


# Lasso의 로지스틱 모형
logitL = LogisticRegression(penalty="l1",solver="saga")  # LASSO , C = 1/alpha    
logitL.fit(X_train, y_train)  

muL = logitL.predict_proba(X_test)  # predict_probability
yhat = logitL.predict(X_test)
print(muL[0:10])
print(yhat[0:10])

logitL.score(X_test, y_test)

 

2) statsmodels

X = whitedata.drop(["y"],axis=1)
y = whitedata["y"]

whitedata["good"] = (whitedata.y > 5).astype(float)
y_binary = whitedata["good"]
X_train, X_test, y_train, y_test = train_test_split(X, y_binary, test_size=0.3, random_state=100)

# statmodels
logit_sm = sm.Logit(y_train, sm.add_constant(X_train)).fit()
print(logit_sm.summary())

print("Parameters:",logit_sm.params)  # 파라미터값 출력
print("Eta:", logit_sm.fittedvalues[0:10])  # 예측값 
# Eta값을 시그모이드 함수에 넣어서 0.5보다 크면 좋은 것

muhat = logit_sm.predict(sm.add_constant(X_test))
print("prob:", muhat[0:10])
# Eta값을 시그모이드 함수에 넣은 값을 출력함
# 0.5보다 크면 good, 0.5보다 작으면 bad

yhat = (muhat > 0.5).astype(int)  # True=1, False=0
print(yhat[0:10])

# 혼동행렬 
cm = confusion_matrix(y_test, yhat)  
print ("Confusion Matrix : \n", cm) 

# accuracy score of the model
print('Test accuracy = ', accuracy_score(y_test, yhat))

# 변수 지정
model = "good ~ x2+x4+x6+x8+x9+x10+x11"  # 유의한 설명변수를 가지고 모형 설정
logit_smf = smf.logit(formula=str(model),data=whitedata).fit()
logit_smf.summary()
# 특별히 지정하지 않아도 베타0(Intercept) 포함해서 분석

 

5. 다범주 로짓 모형

# iris data
iris = datasets.load_iris()
X = iris.data[:, :2]  # 2개의 설명변수(0, 1)만 사용
Y = iris.target

mlogit_sm = MNLogit(Y,sm.add_constant(X)).fit(maxiter=100,method="bfgs")
mlogit_sm.summary()

prob = mlogit_sm.predict(sm.add_constant(X))
print(prob[0:10])

yhat = np.argmax(prob,1)
print(yhat)