Add files via upload

RBF.py

@@ -0,0 +1,207 @@
+import numpy as np
+from numpy.linalg import norm
+import matplotlib.pyplot as plt
+import math
+import pandas as pd
+import random as rd
+plt.close('all')
+
+'''IBM stock data'''
+IBM = pd.read_csv('IBM_stock_close_price.csv')
+ibm = np.array(IBM)
+plt.figure(7)
+plt.plot(ibm[:,1])
+plt.title('IBM close price 1961~1962')
+plt.xlabel('days')
+plt.ylabel('close price')
+
+Max = max(IBM.iloc[:,1])
+Min = min(IBM.iloc[:,1])
+
+for i in range(ibm.shape[0]):
+    ibm[i,1] = (ibm[i,1] - Min) / (Max - Min)
+
+ibm_train = ibm[0:250,1]
+ibm_test = ibm[250:,1]
+
+
+'''initialized center'''
+Xmax = 1
+Xmin = 0
+hidden_neuron = 15; '''hidden neuron for encoding to a single input variable'''
+center = []
+for neuron in range(hidden_neuron):
+    center.append( Xmin + (Xmax - Xmin)/(2*hidden_neuron) + (neuron - 1)*((Xmax - Xmin)/hidden_neuron) )
+    center[neuron] = float('%.4f' % center[neuron]) # get float number by 4 decimal point
+    '''10 RBF centers'''
+
+print('Initial center =', center)
+#dmax = (2*(hidden_neuron-1) - 1)/(2*hidden_neuron)
+dmax = max(center) - min(center)
+width = dmax/(2*hidden_neuron)**0.5
+Width = width*np.ones((1,hidden_neuron)) # width array for optimization
+beta = - 1 / (2 * (width ** 2))
+delta = width*np.ones((hidden_neuron,1))
+
+class RBFNN:
+
+    def __init__(self, n_input, n_hidden, w_in_hidden = None, w_out = None):
+        self.n_input = n_input  # practically, input neurons depend on sample features. Set up first is unnecessary
+        self.n_hidden = n_hidden
+        self.w_in_hidden = w_in_hidden  # usually set = 1
+        '''3 parameters to be optimized during iteration, thus give them self attribute'''
+        self.w_out = w_out  # randomly initialized weight
+        self.center = center
+        self.delta = delta
+        self.output = []
+        self.hidden_out = []
+
+    def gaussian(self, input, Center):  # calculate 1 * 5 array input mapping to 10 * 5 output matrix
+        '''input x dimension = number of features'''
+        encode_neuron = self.n_hidden
+        '''I don't know how to set the width theoretically'''
+        hidden_out = np.zeros((encode_neuron, 1))  # 10*1 matrix
+        for j in range(encode_neuron):
+            hidden_out[j, 0] = np.exp(-self.Euclidean(input, Center[j]) ** 2 / 2 * (self.delta[j] ** 2))
+
+        return hidden_out  # 10 outputs of hidden neuron
+
+    def Euclidean(self, a, b):
+        return norm(a - b)  # return euclidean distance of 2 vector
+
+    def HiddenLayer(self, input):
+        centroid = self.center
+        hidden_output = self.gaussian(input,centroid)  # hidden neurons output of 1 input sample
+        self.hidden_out = hidden_output
+        return hidden_output # 10*1 matrix
+
+    def feed_forward(self, input):  # feed forward calculation from input to output layer, calculate one sample at a time
+        hidden_out = self.HiddenLayer(input) # 10*1 array
+        weight = self.w_out # 1*10 array
+        output = np.dot(weight, hidden_out)
+        return output # output value of single sample
+
+    def train(self, train_sample, sample_target):
+
+        # Gradient descent learning: weight, center, width
+        # update self attribute, keep the update to next iteration
+        # Weight update:
+        lr = 0.2  # learning rate
+        num_hidden = self.n_hidden
+
+        samples = train_sample.shape[0]
+        w_gradient = np.zeros((num_hidden, 1))
+        c_gradient = np.zeros((num_hidden, 1))
+        delta_gradient = np.zeros((num_hidden, 1))
+        beta = np.zeros((num_hidden,1))
+        Error = 0 # total error
+
+        for sample in range(samples):
+            output = self.feed_forward(train_sample[sample,:]) # calculate one output at a time
+            Error += output - sample_target[sample]
+            # Batch learning
+
+            for j in range(num_hidden):
+                beta[j] = ( -1 / (2*(self.delta[j]**2)) )
+
+                w_gradient[j] = w_gradient[j] + (output - sample_target[sample])*\
+                                                np.exp( beta[j] * (self.Euclidean(train_sample[sample,:], self.center[j]) ** 2) )
+                c_gradient[j] = c_gradient[j] + (output - sample_target[sample])*self.w_out[0,j]*beta[j]*( self.Euclidean(train_sample[sample,:],self.center[j] )) \
+                                  *np.exp( beta[j] * (self.Euclidean(train_sample[sample,:], self.center[j]) ** 2) )
+                delta_gradient[j] = delta_gradient[j] + (output - sample_target[sample])*(self.Euclidean(train_sample[sample,:],self.center[j])**2)\
+                                                        *self.w_out[0,j]*( (self.delta[j])**(-3) )* \
+                                                        np.exp(beta[j] * (
+                                                        self.Euclidean(train_sample[sample, :], self.center[j]) ** 2))
+                if self.delta[j] > 0:
+                   if delta_gradient[j] < 0:
+                      print('error =',output - sample_target[sample])
+
+        for i in range(num_hidden):
+           self.w_out[0,i] += -lr*w_gradient[i]/samples
+           self.center[i] += 2*lr*c_gradient[i]/samples
+           self.delta[i] += -lr*delta_gradient[i]/samples
+
+
+        #print('w_out =', self.w_out)
+        #print('center =', self.center)
+        print('delta =', self.delta)
+
+
+        return Error # scale in [0,1]
+
+    def predict(self, test_data):
+        # predict for one test data
+        estimate = self.feed_forward(test_data)
+        # decode
+        estimate = Min + estimate * (Max - Min)
+        return estimate
+
+if __name__ == '__main__':
+
+   Max = max(IBM.iloc[:,1])
+   Min = min(IBM.iloc[:,1])
+   # network parameters
+   # Use close price as prediction
+   slide_window = 3  # days of data to predict one future value, which is seen as feature
+   duration = ibm_train.shape[0] - slide_window # days duration of whole predict process
+   #duration = 1
+   input_neurons = slide_window # number of features = input neurons number
+   hidden_neurons = 15
+   weight = np.zeros((1,hidden_neurons)) # weight = 1*10 matrix
+   for j in range(hidden_neurons):
+       weight[0,j] = rd.uniform(0,1)
+
+   rbf = RBFNN(input_neurons, hidden_neurons, None, weight) # build an RBF network archictecture
+
+   # create training set time series
+   samples = duration
+   train_set = np.zeros((duration,slide_window))
+   target = np.zeros((samples,1))
+
+   for sample in range(samples):
+       train_set[sample,:] = ibm_train[sample:(sample + slide_window)]
+       target[sample] = ibm[sample + slide_window,1]
+
+
+   Iteration = 100
+   Error = 0
+   Performance = np.zeros((Iteration,1))
+   Predict = np.zeros((Iteration,1))
+   for iterate in range(Iteration):
+       error = 0
+       # Batch learning
+       error = rbf.train(train_set, target)  # output = total error
+       E = error * (Max - Min) / samples + Min
+       print('Error in epoch: ', E)
+       Performance[iterate] = E
+       Predict[iterate] = rbf.predict(train_set[0, :])
+
+   #print(Performance)
+   plt.figure()
+   plt.title('Error Performance')
+   plt.plot(Performance)
+
+   # training set performance
+   check = np.zeros((train_set.shape[0], 1))
+   for i in range(train_set.shape[0]):
+       check[i] = rbf.predict(train_set[i,:])
+       target[i] = Min + (Max - Min)*target[i]
+
+   # test set performance
+   test_duration = ibm_test.shape[0] - slide_window
+   test_set = np.zeros((test_duration,slide_window))
+   for sample in range(test_set.shape[0]):
+       test_set[sample,:] = ibm_test[sample:sample + slide_window]
+   test_target = ibm_test[slide_window:]
+   test_out = np.zeros((train_set.shape[0], 1))
+   for i in range(test_set.shape[0]):
+       test_out[i] = rbf.predict(test_set[i, :])
+       test_target[i] = Min + (Max - Min) * test_target[i]
+
+   plt.figure()
+   plt.title('test performance')
+   plt.plot(test_target, 'r-', label='Real output')
+   plt.plot(test_out, 'b-', label='Predict output')
+   plt.legend(loc='upper right')
+
+   plt.show()