In [1]:
import pandas
data = pandas.read_csv("iris.csv", header=None,
names=["lod.dl.", "lod.sz.", "pl.dl.", "pl.sz.", "Gatunek"])
data[:8]
Out[1]:
In [2]:
import numpy as np
m = len(data)
X = np.matrix(data[["lod.dl.", "lod.sz.", "pl.dl.", "pl.sz."]])
X = np.hstack((np.ones(m).reshape(m,1), X))
y = np.matrix(data[["Gatunek"]]).reshape(m,1)
def mapY(y, c):
n = len(y)
yBi = np.matrix(np.zeros(n)).reshape(n, 1)
yBi[y == c] = 1.
return yBi
def indicatorMatrix(y):
classes = np.unique(y.tolist())
Y = mapY(y, classes[0])
for c in classes[1:]:
Y = np.hstack((Y, mapY(y, c)))
Y = np.matrix(Y).reshape(len(y), len(classes))
return Y
#print(X[:10])
#print(y[:10])
Y = indicatorMatrix(y)
#print(Y[:10])
# Dzielimy dane na trenujące i testowe
XTrain, XTest = X[:100], X[100:]
YTrain, YTest = Y[:100], Y[100:]
In [3]:
#%matplotlib inline
#import matplotlib.pyplot as plt
#import seaborn as sns
#g = sns.pairplot(data, hue="Gatunek")
from IPython.display import Image
Image(filename='iris1.png')
Out[3]:
In [4]:
def logistic(x):
return 1.0/(1.0 + np.exp(-x))
x = np.linspace(-5,5,200)
y = logistic(x)
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)
plt.ylim(-.1,1.1)
ax.plot(x, y, linewidth='2');
Funkcja regresji logistycznej:
$$h_\theta(x) = g(\theta^Tx) = \dfrac{1}{1+e^{-\theta^Tx}}$$Funkcja regresji logistycznej (wersja macierzowa, można zastosować do wielu wierszy):
$$h_\theta(X) = g(X\theta) = \dfrac{1}{1+e^{-X\theta}}$$
In [6]:
def h(theta, X):
return 1.0/(1.0 + np.exp(-X*theta))
thetaTemp = np.zeros(5).reshape(5,1)
h(thetaTemp, XTrain[:10])
Out[6]:
Funkcja kosztu dla regresji logistycznej: $$\small$$ $$ J(\theta) = -\dfrac{1}{m} [\sum_{i=1}^{m} y^{(i)} \log h_\theta(x^{(i)})+ (1-y^{(i)}) \log (1-h_\theta(x^{(i)}))]$$
Zapis macierzowy gradientu (z wyjątkiem $h_\theta$ taki sam jak dla regresji liniowej):
$$\nabla J(\theta) = \frac{1}{|\vec y|} X^T\left(h_\theta(X)-\vec y\right)$$
In [7]:
def J(h,theta,X,y):
m = len(y)
s1 = np.multiply(y, np.log(h(theta,X)))
s2 = np.multiply((1 - y), np.log(1 - h(theta,X)))
return -np.sum(s1 + s2, axis=0)/m
YTrain0 = YTrain[:,0]
print(J(h, thetaTemp, XTrain, YTrain0))
def dJ(h,theta,X,y):
return 1.0/len(y)*(X.T*(h(theta,X)-y))
print(dJ(h, thetaTemp, XTrain, YTrain0))
In [8]:
# Zmienić tutaj test class na 0,1,2
testClass = 0
YTrainBi = YTrain[:,testClass]
#data0 = np.concatenate((XTrain, YTrainBi), axis=1)
#df = pandas.DataFrame(data0[:,1:])
#df.columns = data.columns
#g = sns.pairplot(df, vars=[c for c in data.columns if c != "Gatunek"],
# hue="Gatunek")
from IPython.display import Image
Image(filename='iris2.png')
Out[8]:
In [16]:
def mbSGD(h, fJ, fdJ, theta, X, y,
alpha=0.01, maxSteps=10000, batchSize = 10):
#errorCurr = fJ(h, theta, X, y)
i = 0
b = batchSize
m = X.shape[0]
while i < maxSteps:
start = (i*b) % m
end = ((i+1)*b) % m
if(end <= start):
end = m
Xbatch = X[start:end]
Ybatch = y[start:end]
theta = theta - alpha * fdJ(h, theta, Xbatch, Ybatch)
i += 1
return theta, []
# Liczy dla klasy ustawionej wcześniej w testClass
thetaBest, errors = mbSGD(h, J, dJ, thetaTemp, XTrain, YTrainBi,
alpha=0.001, maxSteps=10000, batchSize = 10)
print(thetaBest)
In [13]:
YTestBi = YTest[:,testClass]
def classifyBi(X):
prob = h(thetaBest, X).item()
return (1, prob) if prob > 0.5 else (0, prob)
acc = 0.0
for i, rest in enumerate(YTestBi):
cls, prob = classifyBi(XTest[i])
print int(YTestBi[i].item()), "<=>", cls, "-- prob:", round(prob, 4)
acc += cls == YTestBi[i].item()
print "Accuracy:", acc/len(XTest)
In [ ]:
#g = sns.pairplot(data, hue="Gatunek")
from IPython.display import Image
Image(filename='iris1.png')
In [ ]:
# Zapis macierzowy, dlaczego to jest poprawne?
def softmax(X):
return np.exp(X)/np.sum(np.exp(X))
X = np.matrix([2.1, 0.5, 0.8, 0.9, 3.2]).reshape(5,1)
P = softmax(X)
print(X)
print "Suma X =", np.sum(X), "\n"
print(P)
print "Suma P =", np.sum(P)
In [ ]:
def trainMaxEnt(X, Y):
n = X.shape[1]
thetas = []
for c in range(Y.shape[1]):
YBi = Y[:,c]
theta = np.matrix(np.random.random(n)).reshape(n,1)
thetaBest, errors = GD(h, J, dJ, theta,
X, YBi, alpha=0.1, eps=10**-4)
thetas.append(thetaBest)
return thetas
thetas = trainMaxEnt(XTrain, YTrain);
for theta in thetas:
print(theta)
In [ ]:
def classify(thetas, X):
regs = np.array([(X*theta).item() for theta in thetas])
probs = softmax(regs)
return np.argmax(probs), probs
print(YTest[:10])
YTestCls = YTest * np.matrix((0,1,2)).T
print(YTestCls[:10])
acc = 0.0
for i in range(len(YTestCls)):
cls, probs = classify(thetas, XTest[i])
correct = int(YTestCls[i].item())
print correct, "<=>", cls, " - ", correct == cls, np.round(probs, 4).tolist()
acc += correct == cls
print "Accuracy =", acc/float(len(XTest))
Mając $X$ oraz $\vec{y}$, liczymy $\theta$ np. za pomocą macierzy normalnej (nie ma żadnej modyfikacji algorytmu uczącego czy funkji kosztu):
$$ \theta = (X^TX)^{-1}X^T\vec{y} $$lub możemy obliczyć wszystkie zestawy parametrów $\theta$ dla wszystkich klas za pomocą jednej operacji, gdzie Y to macierz indykatorów wszystkich klas:
$$ \Theta = (X^TX)^{-1}X^TY $$
In [ ]:
def norm(X, y):
return np.linalg.pinv(X.T*X)*X.T*y
# Macierz wszystkich parametrów
thetas = norm(XTrain, YTrain)
print(thetas)
# W porównianiu pojedyncze zestawy parametrów
theta0 = norm(XTrain, YTrain[:,0])
print(theta0)
theta1 = norm(XTrain, YTrain[:,1])
print(theta1)
theta2 = norm(XTrain, YTrain[:,2])
print(theta2)
Mając macierz $\Theta$ potrafimy teraz obliczyć maciesz kosztów $K$ dla wszystjich wierszy testowych i dla wszystkich klas za pomocą jednej operacji, czyli
$$ K = X\Theta $$Aby zaklasyfikować wszystkie wiersze ze zbioru testowego $X$, wystarczy teraz wykonać $\mathop{\textrm{arg}\,\textrm{max}}$ po wierszach macierzy $K=X\Theta$.
$$ c_i = \mathop{\textrm{arg}\,\textrm{max}}_{j \in \{1, \ldots ,k\}} \{ (X\Theta)_{ij} \} , \textrm{dla } i = 1,\ldots,m $$
In [ ]:
def h(theta,X):
return X*theta
def classify(thetas, X):
costs = h(thetas, X)
best = np.argmax(costs, axis=1)
return best, costs
Res, Costs = classify(thetas, XTest)
acc = sum(Res == YTestCls).item()/float(YTestCls.shape[0])
print "Accuracy:", acc