40
def gradient(X_norm,y,theta,alpha,m,n,num_it): 
    temp=np.array(np.zeros_like(theta,float)) 
    for i in range(0,num_it): 
     h=np.dot(X_norm,theta) 
     #temp[j]=theta[j]-(alpha/m)*( np.sum((h-y)*X_norm[:,j][np.newaxis,:]) ) 
     temp[0]=theta[0]-(alpha/m)*(np.sum(h-y)) 
     temp[1]=theta[1]-(alpha/m)*(np.sum((h-y)*X_norm[:,1])) 
     theta=temp 
    return theta 



X_norm,mean,std=featureScale(X) 
#length of X (number of rows) 
m=len(X) 
X_norm=np.array([np.ones(m),X_norm]) 
n,m=np.shape(X_norm) 
num_it=1500 
alpha=0.01 
theta=np.zeros(n,float)[:,np.newaxis] 
X_norm=X_norm.transpose() 
theta=gradient(X_norm,y,theta,alpha,m,n,num_it) 
print theta 

Moja theta z powyższego kodu jest 100.2 100.2, ale powinno być 100.2 61.09 w Matlab, która jest poprawna.metoda gradientu prostego za pomocą Python i numpy

+5

średnik są ignorowane w python i wcięcia, jeśli są zasadnicze. –

Odpowiedz

101

Myślę, że twój kod jest nieco zbyt skomplikowany i potrzebuje więcej struktury, ponieważ w przeciwnym razie stracisz wszystkie równania i operacje. W końcu ten regresji sprowadza się cztery operacje:

  1. Oblicz hipotezy h = X * teta
  2. Obliczenie strat = h - Y i może kwadratu kosztów (strata^2)/2M
  3. obliczyć gradient = X”* utrata/m
  4. Aktualizacja parametry theta = theta - alfa * Gradient

W twoim przypadku, myślę, że masz mylić m z n. Tutaj m oznacza liczbę przykładów w twoim zestawie treningowym, a nie liczbę funkcji.

Rzućmy okiem na moje zmianę kodu:

import numpy as np 
import random 

# m denotes the number of examples here, not the number of features 
def gradientDescent(x, y, theta, alpha, m, numIterations): 
    xTrans = x.transpose() 
    for i in range(0, numIterations): 
     hypothesis = np.dot(x, theta) 
     loss = hypothesis - y 
     # avg cost per example (the 2 in 2*m doesn't really matter here. 
     # But to be consistent with the gradient, I include it) 
     cost = np.sum(loss ** 2)/(2 * m) 
     print("Iteration %d | Cost: %f" % (i, cost)) 
     # avg gradient per example 
     gradient = np.dot(xTrans, loss)/m 
     # update 
     theta = theta - alpha * gradient 
    return theta 


def genData(numPoints, bias, variance): 
    x = np.zeros(shape=(numPoints, 2)) 
    y = np.zeros(shape=numPoints) 
    # basically a straight line 
    for i in range(0, numPoints): 
     # bias feature 
     x[i][0] = 1 
     x[i][1] = i 
     # our target variable 
     y[i] = (i + bias) + random.uniform(0, 1) * variance 
    return x, y 

# gen 100 points with a bias of 25 and 10 variance as a bit of noise 
x, y = genData(100, 25, 10) 
m, n = np.shape(x) 
numIterations= 100000 
alpha = 0.0005 
theta = np.ones(n) 
theta = gradientDescent(x, y, theta, alpha, m, numIterations) 
print(theta) 

w pierwszej chwili stworzyć mały losowy zbiór danych, który powinien wyglądać tak:

Linear Regression

Jak widać I dodano również wygenerowaną linię regresji i formułę obliczoną przez Excel.

Musisz dbać o intuicję regresji za pomocą gradientowego zniżania. Gdy wykonasz pełne przetwarzanie wsadowe danych X, musisz zmniejszyć m-straty każdego przykładu do pojedynczej aktualizacji wagi. W tym przypadku jest to średnia z sumy na gradientach, a więc podział przez m.

Następną rzeczą, na którą należy zwrócić uwagę, jest śledzenie zbieżności i dostosowanie szybkości uczenia się. Jeśli o to chodzi, zawsze powinieneś śledzić swój koszt w każdej iteracji, a może nawet spiskuj.

Jeśli uruchomić mój przykład, theta zwrócony będzie wyglądać następująco:

Iteration 99997 | Cost: 47883.706462 
Iteration 99998 | Cost: 47883.706462 
Iteration 99999 | Cost: 47883.706462 
[ 29.25567368 1.01108458] 

Która jest całkiem blisko do równania, które zostały obliczone przez program Excel (y = x + 30). Zwróć uwagę, że kiedy przekazywaliśmy wartość odchylenia do pierwszej kolumny, pierwsza wartość theta oznacza wagę obciążenia.

+3

W gradientDescent, '/ 2 * m' ma być'/(2 * m) '? – physicsmichael

+1

@ vgm64 ouch, tak masz rację! Poprawiono to. –

+2

Użycie słowa 'loss' dla absolutnej różnicy nie jest dobrym pomysłem, ponieważ" strata "jest zwykle synonimem" kosztu ". Nie musisz też w ogóle podawać 'm', tablice NumPy znają swój własny kształt. –

8

Poniżej znajduje się moja implementacja gradientu zejścia dla problemu regresji liniowej.

Najpierw należy obliczyć gradient, taki jak X.T * (X * w - y)/N, i zaktualizować bieżący teta równoczesnym gradientem.

  • X: cecha matrycy
  • r: wartość docelową
  • w: wag/wartości
  • N. wielkość zbiorze treningowym

Oto kod pyton:

import pandas as pd 
import numpy as np 
from matplotlib import pyplot as plt 
import random 

def generateSample(N, variance=100): 
    X = np.matrix(range(N)).T + 1 
    Y = np.matrix([random.random() * variance + i * 10 + 900 for i in range(len(X))]).T 
    return X, Y 

def fitModel_gradient(x, y): 
    N = len(x) 
    w = np.zeros((x.shape[1], 1)) 
    eta = 0.0001 

    maxIteration = 100000 
    for i in range(maxIteration): 
     error = x * w - y 
     gradient = x.T * error/N 
     w = w - eta * gradient 
    return w 

def plotModel(x, y, w): 
    plt.plot(x[:,1], y, "x") 
    plt.plot(x[:,1], x * w, "r-") 
    plt.show() 

def test(N, variance, modelFunction): 
    X, Y = generateSample(N, variance) 
    X = np.hstack([np.matrix(np.ones(len(X))).T, X]) 
    w = modelFunction(X, Y) 
    plotModel(X, Y, w) 


test(50, 600, fitModel_gradient) 
test(50, 1000, fitModel_gradient) 
test(100, 200, fitModel_gradient) 

test1 test2 test2

+3

Niepotrzebne polecenie importu: importowanie pand jako pd – fviktor

+0

@Muatik Nie rozumiem, jak można uzyskać gradient z wewnętrznym produktem błędu i zestawu szkoleń: 'gradient = x.T * błąd/N' Jaka jest logika tego? – eggie5

1

wiem, to pytanie już była odpowiedź, ale zrobiłem jakąś aktualizację funkcji GD:

### COST FUNCTION 

def cost(theta,X,y): 
    ### Evaluate half MSE (Mean square error) 
    error = np.dot(X,theta) - y 
    J = np.sum(error ** 2)/(2*m) 
    return J 

cost(theta,X,y) 



def GD(X,y,theta,alpha): 

    cost_histo = [0] 
    theta_histo = [0] 

    # an arbitrary gradient, to pass the initial while() check 
    delta = [np.repeat(1,len(X))] 
    # Initial theta 
    old_cost = cost(theta,X,y) 

    while (np.max(np.abs(delta)) > 1e-6): 
     error = np.dot(X,theta) - y 
     delta = np.dot(np.transpose(X),error)/len(y) 
     trial_theta = theta - alpha * delta 
     trial_cost = cost(trial_theta,X,y) 
     while (trial_cost >= old_cost): 
      trial_theta = (theta +trial_theta)/2 
      trial_cost = cost(trial_theta,X,y) 
      cost_histo = cost_histo + trial_cost 
      theta_histo = theta_histo + trial_theta 
     old_cost = trial_cost 
     theta = trial_theta 
    Intercept = theta[0] 
    Slope = theta[1] 
    return [Intercept,Slope] 

res = GD(X,y,theta,alpha) 

Ta funkcja zmniejszenia alfa nad iteracji dokonywania funkcję zbyt zbiegają się szybciej zobaczyć Estimating linear regression with Gradient Descent (Steepest Descent) dla przykład w R. Stosuję tę samą logikę, ale w Pythonie.