Prob 3-1¶
In [207]:
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
Import Data,建立training set and testing set¶
In [208]:
# x1:T, x2:V, x3:AP, x4:RH, x0:1
data3 = pd.read_excel("data.xlsx")
#data = data3.head()
data = data3
data['x0'] = 1 #增加一排常數項當作x0
data['x1'] = data3['T']
data['x2'] = data3['V']
data['x3'] = data3['AP']
data['x4'] = data3['RH']
data['y'] = data3['EP']
data = data.drop('T',1).drop('V',1).drop('AP',1).drop('RH',1).drop('EP',1)
data.head()
Out[208]:
In [210]:
train_x = np.matrix( data.drop('y', 1)[0:400] )
train_y = np.matrix( data.y[0:400] ).reshape(-1,1)
test_x = np.matrix( data.drop('y', 1)[400:500] )
test_y = np.matrix( data.y[400:500] ).reshape(-1,1)
train_df = data[0:400]
test_df = data[400:500]
D = 5 # dim (T,V,AP,RH) 維度是4, 加上x0維度是5
#train_y
建立model¶
In [211]:
def basis_func_M2(x): #x=[x0,x1,x2,x3,x4] np.array #一筆資料
D = x.size
X = []
feature_names = []
for i in range(D):
for j in range(i,D):
X.append(x[i]*x[j])
feature_names.append('x'+str(i)+'*x'+str(j))
return X , feature_names
def basis_func_M3(x): #x=[x0,x1,x2,x3,x4] np.array #一筆資料
D = x.size
X = []
feature_names = []
for i in range(D):
for j in range(i,D):
for k in range(j,D):
X.append(x[i]*x[j]*x[k])
feature_names.append('x'+str(i)+'*x'+str(j)+'*x'+str(k))
return X , feature_names
def basis_func(x,deg): #x=[x0,x1,x2,x3,x4] np.array #回傳基底函數作用後的list
if deg==2:
return basis_func_M2(x)
elif deg==3:
return basis_func_M3(x)
else:
return []
In [212]:
def buildup_basis_func_matrix(train_x,deg): # deg:degree(order)
N = train_x.size / train_x[0].size
#N=2
MX = [] #matrix of x
F = [] #feature_names
for i in range(N):
r = np.array(train_x[i])[0] #取出某一row, 因是matrix形式要取list
#print r
X = []
X,F = basis_func(r,deg) #丟進基底函數
MX.append(X)
return np.matrix(MX),F # matrix of fi(x), feature_names
#buildup_basis_func_matrix(train_x,3)
In [213]:
def buildup_model_no_regularization(train_x,train_y,deg):
MX,F = buildup_basis_func_matrix(train_x,deg)
#return MX,F
psX = np.linalg.pinv(MX) #算 pseudo-inverse
Y = train_y #Y=[[y0],[y1]...] target/label,column matrix
#print psX
W = np.dot(psX,Y) #weight W = [ [w0] ,[w1]...] , column matrix
return W,F
#W,F = buildup_model_no_regularization(train_x,train_y,3)
#print W
In [214]:
def predict_model_no_regularization(test_x,test_y,W,deg): # weight W = [ [w0] ,[w1]...] 是train好的 Weight, column matrix
MX,F = buildup_basis_func_matrix(test_x,deg)
pred_y = np.dot(MX,W)
return pred_y # pred_y is column matrix
#predict_model_no_regularization(test_x,test_y,W,3)
$ L1-norm: 1/2 * lambda * |w| $
$ L2-norm: 1/2 * lambda * |w|^2 $
In [215]:
def colMatrix2Array(colM): # column matrix (np.marix) --> (np.array)
L = np.array(colM.reshape(1,colM.size).tolist()[0])
return L, L.size
def L1_norm(W): # W: column matrix (np.marix) -> |W|
return np.sum(np.abs(W))
print L1_norm(np.matrix([[-12],[11]]))
def L2_norm(W): # W: column matrix (np.marix) -> |W|^2
return np.sum(np.square(W))
print L2_norm(np.matrix([[-12],[11]]))
def calError_model(pred_y,test_y,norm,lenda,W): #pred_y,test_y, W 都是column matrix
py, N = colMatrix2Array(pred_y)
ty, N = colMatrix2Array(test_y)
sum_square_diff = np.sum( np.square( py-ty ) )
E_D = 0.5* sum_square_diff
RMS = np.sqrt(sum_square_diff/N)
#print N,sum_square_diff
if norm == 0: # norm = 0 or lenda=0 means no_regularization
return E_D, E_D, 0, RMS
elif norm == 1:
E_W = 0.5*lenda*L1_norm(W)
return (E_D+E_W), E_D, E_W, RMS
elif norm == 2:
E_W = 0.5*lenda*L2_norm(W)
return (E_D+E_W), E_D, E_W, RMS
else:
return -1, -1, -1, -1
calError_model(np.matrix([[12],[11]]),np.matrix([[10],[11]]),1,1,W)
Out[215]:
In [216]:
def Reg(train_x,train_y,test_x,test_y,degree): #main
#---------build up model & find coeff---------------------------
W,feature_names = buildup_model_no_regularization(train_x,train_y,degree)
coef,fn = colMatrix2Array(W)
coef_pair = zip(np.abs(coef),range(fn))
#---------plot coeff---------------------------
print 'plot:'
plt.plot(range(fn),np.abs(coef),'b-')
plt.xlabel('coef order')
plt.ylabel('abs(w)')
plt.title('Coef Weight')
plt.show()
#plt.close()
#---------show coeff---------------------------
coef_pair.sort()
rank = [ coef_pair[-i] for i in range(1,4) ]
print 'coef=',rank
print 'feature_names=', [ feature_names[pair[1]] for pair in rank ]
print '*'*64
#---------predict y from test_x---------------------
#s_name = 'pred_y_deg'+str(degree)+'_byTest'
pred_y = predict_model_no_regularization(test_x,test_y,W,degree)
E, E_D, E_W, RMS_test = calError_model(pred_y,test_y,0,0,W)# norm = 0 , lenda=0 means no_regularization
#---------predict y from train_x---------------------
#s_name = 'pred_y_deg'+str(degree)+'_byTrain'
pred_y = predict_model_no_regularization(train_x,train_y,W,degree)
E, E_D, E_W, RMS_train = calError_model(pred_y,train_y,0,0,W)# norm = 0 , lenda=0 means no_regularization
return [RMS_train,RMS_test]
In [217]:
deg =2
RMS = Reg(train_x,train_y,test_x,test_y,deg) #算出來=4.17413184555, 跟原本w/o sklearn的答案4.17413184559極為近似
print '*'*64
print 'RMS=', RMS
In [218]:
deg =3
RMS = Reg(train_x,train_y,test_x,test_y,deg) #算出來 4.08913247008, 跟原本4.0766289726稍有落差
print '*'*64
print 'RMS=', RMS
Prob 3-2 ¶
For M = 3 (order=3)¶
In [219]:
cols = ['x1','x2','x3','x4']
RMS_matrix = [0,0,0,0]
for i in range(4):
col_name = cols[i]
train_x0 = np.array( data['x0'])
train_xn = np.array( data[col_name])
train_x = np.matrix( zip(train_x0,train_xn))
train_y = np.matrix( data.y ).reshape(-1,1)
print train_x[0].size
RMS_matrix[i]= Reg(train_x,train_y,train_x,train_y,3)[0]
RMS_df = pd.DataFrame(np.matrix( zip(cols,RMS_matrix)) )
RMS_df #印出每一個feature的RMS
Out[219]:
In [ ]:
In [ ]:
沒有留言:
張貼留言