import math import numpy as np import matplotlib.pyplot as plt import matplotlib.colors as pltcolors from sklearn import metrics as metrics from sklearn.linear_model import LogisticRegression as SkLinReg import scipy as sy import seaborn as sns import pandas import tensorflow as tf 

# Single feature, Gaussian noise nFeatures = 1 def generateBatch(N): # xMin = 0 xMax = 1 b = 0.2 std = 0.2 # Threshold from 0 to 1 threshold = 1 # x = np.random.uniform(xMin, xMax, N) # 4th degree relation between y and x yClean = 2*(x**4 + (x-0.3)**3 + b) labels = yClean + np.random.normal(0, std, N) > threshold return (x, yClean, labels) 

N = 2000 # x and y have 1 dim in R, label has 1 dim in B xTrain, yCleanTrain, labelTrain = generateBatch(N) colors = ['blue','red'] fig = plt.figure(figsize=(8,4)) plt.subplot(1,2,1) plt.scatter(xTrain, yCleanTrain, c=labelTrain, cmap=pltcolors.ListedColormap(colors), marker=',', alpha=0.01) plt.xlabel('x') plt.ylabel('y') plt.grid() plt.subplot(1,2,2) plt.scatter(xTrain, labelTrain, marker=',', alpha=0.01) plt.xlabel('x') plt.ylabel('label') plt.grid() 

count, bins, ignored = plt.hist(labelTrain*1.0, 10, density=True, alpha=0.5) p = np.mean(labelTrain) print('Bernouilli parameter of the distribution:', p) 

xTest, yTest, labelTest = generateBatch(N) 

def plotHeatMap(X, classes, title=None, fmt='.2g', ax=None, xlabel=None, ylabel=None): """ Fix heatmap plot from Seaborn with pyplot 3.1.0, 3.1.1  https://stackoverflow.com/questions/56942670/matplotlib-seaborn-first-and-last-row-cut-in-half-of-heatmap-plot  """ ax = sns.heatmap(X, xticklabels=classes, yticklabels=classes, annot=True, fmt=fmt, cmap=plt.cm.Blues, ax=ax) #notation: "annot" not "annote" bottom, top = ax.get_ylim() ax.set_ylim(bottom + 0.5, top - 0.5) if title: ax.set_title(title) if xlabel: ax.set_xlabel(xlabel) if ylabel: ax.set_ylabel(ylabel) def plotConfusionMatrix(yTrue, yEst, classes, title=None, fmt='.2g', ax=None): plotHeatMap(metrics.confusion_matrix(yTrue, yEst), classes, title, fmt, ax, xlabel='Estimations', \ ylabel='True values'); 

def logistic(X): return (1+(np.exp(-(X))))**-1 

xx = np.linspace(-10, 10) xlogistic = logistic(xx) plt.figure(figsize=(10,5)) plt.subplot(1, 2, 1) plt.plot(xx, xlogistic) plt.grid() plt.subplot(1, 2, 2) plt.plot(xx, np.log(xlogistic)) plt.grid() 

b = np.linspace(-5, 5) w = 1 px = np.zeros(len(b)) for i in range(len(b)): fx = logistic(b[i] + w*xTrain) px[i] = 1/N * np.sum(np.log(labelTrain*fx + (1-labelTrain)*(1-fx))) plt.plot(b, px); plt.xlabel('$b$') plt.ylabel('l(b, X)') plt.grid() bHat = b[np.argmax(px)] print('Estimate b =', bHat) 

w = np.linspace(-20, 20) px = np.zeros(len(w)) for i in range(len(w)): fx = logistic(bHat + w[i]*xTrain) px[i] = 1/N * np.sum(np.log(labelTrain*fx + (1-labelTrain)*(1-fx))) plt.plot(w, px); plt.xlabel('w') plt.ylabel('l(w, X)') plt.grid() wHat = w[np.argmax(px)] print('Estimate w =', wHat) 

pXTest0 = logistic(bHat + wHat * xTest) labelEst0 = pXTest0 > 0.5 

plt.scatter(xTest, pXTest0, c=labelEst0, cmap=pltcolors.ListedColormap(colors), marker=',', alpha=0.01); plt.scatter(xTest, yTest/np.max(yTest), c = labelTest, cmap=pltcolors.ListedColormap(colors), marker='x', alpha=0.01); plt.xlabel('x') plt.legend(('Estimated probability', 'Normalized model')); 

plt.hist(labelEst0*1.0, 10, density=True) print('Bernouilli parameter =', np.mean(labelEst0)) 

accuracy0 = np.sum(labelTest == labelEst0)/N print('Accuracy =', accuracy0) 

print('Precision =', np.sum(labelTest[labelEst0 == 1])/np.sum(labelEst0)) 

print('Recall =', np.sum(labelTest[labelEst0 == 1])/np.sum(labelTest)) 

plotConfusionMatrix(labelTest, labelEst0, np.array(['Blue', 'Red'])); 

print(metrics.classification_report(labelTest, labelEst0)) 

model1 = SkLinReg(solver='lbfgs') model1.fit(xTrain.reshape(-1,1), labelTrain) model1.coef_ 

labelEst1 = model1.predict(xTest.reshape(-1,1)) print('Accuracy =',model1.score(xTest.reshape(-1,1), labelTest)) 

plt.hist(labelEst1*1.0, 10, density=True) print('Bernouilli parameter =', np.mean(labelEst1)) 

plotConfusionMatrix(labelTest, labelEst1, np.array(['Blue', 'Red'])) 

print(metrics.classification_report(labelTest, labelEst1)) 

logit_roc_auc = metrics.roc_auc_score(labelTest, labelEst1) fpr, tpr, thresholds = metrics.roc_curve(labelTest, model1.predict_proba(xTest.reshape(-1,1))[:,1]) plt.plot(fpr, tpr, label='Logistic Regression (area = %0.2f)' % logit_roc_auc) plt.plot([0, 1], [0, 1],'r--') plt.xlim([0.0, 1.0]) plt.ylim([0.0, 1.05]) plt.xlabel('False Positive Rate') plt.ylabel('True Positive Rate') plt.title('Receiver operating characteristic') plt.legend(loc="lower right"); 

# Labels as floats {0., 1.} labelTrainF = np.multiply(labelTrain, 1.0) labelTrainF.dtype, labelTrainF.shape 

# (Mini) Batch size nBatch = 100 # Number of batches per Epoch nBatchPerEpoch =20 # Number of epochs nEpochMax = 1000 

# Simple custom layer exposing the linear regression model class MyLogisticRegressionLayer(tf.keras.layers.Layer): def __init__(self, *args, **kwargs): super(MyLogisticRegressionLayer, self).__init__(*args, **kwargs) def build(self, input_shape): self.w = self.add_weight( shape=input_shape[0], dtype=self.dtype, initializer=tf.keras.initializers.ones(), #regularizer=tf.keras.regularizers.l2(0.02), trainable=True) self.b = self.add_weight( shape=1, dtype=self.dtype, initializer=tf.keras.initializers.ones(), #regularizer=tf.keras.regularizers.l2(0.02), trainable=True) @tf.function def call(self, x, training=None): return tf.math.sigmoid(tf.math.add(tf.math.multiply(x, self.w), self.b)) 

# Using TensorFlow 2.0 style of metrics to implement accuracy class MyBinaryAccuracy(tf.keras.metrics.Metric): def __init__(self, name='my_accuracy', **kwargs): super(MyBinaryAccuracy, self).__init__(name=name, **kwargs) self.accuracySum = self.add_weight(name='accuracySum', initializer='zeros') self.accuracyCount = self.add_weight(name='accuracyCount', initializer='zeros') def update_state(self, labels, yEst): labels = tf.cast(labels, tf.bool) labelEst = tf.greater(yEst, 0.5) values = tf.cast(tf.equal(labels, labelEst), self.dtype) self.accuracySum.assign_add(tf.reduce_sum(values)) self.accuracyCount.assign_add(values.get_shape()[0]) def result(self): return self.accuracySum / self.accuracyCount 

# Model 1, instantiate the custom layer model1 = tf.keras.Sequential([MyLogisticRegressionLayer(input_shape=[nFeatures], dtype="float64")]) # Stochastic Gradient Descent Optimizer optim1 = tf.keras.optimizers.SGD(0.01) # Perform a train step on a mini-batch # This function's code is rewritten by TensorFlow 2.0 and shall be compiled at every execution of the optimizer @tf.function def trainStep1(x, labels): with tf.GradientTape() as tape: predictions = model1(x, training=True) loss = -tf.reduce_sum(tf.math.log((labels * predictions) + ((1 - labels) * (1 - predictions)))) #loss = tf.keras.losses.categorical_crossentropy(labels, predictions) gradients = tape.gradient(loss, model1.trainable_variables) optim1.apply_gradients(zip(gradients, model1.trainable_variables)) return loss, predictions # Initialize values and loop on epochs and mini batch epoch = 0 cost_epoch = 1 histo = [] accuracy = MyBinaryAccuracy() for epoch in range(nEpochMax): cost_cumul = 0 accuracy.reset_states() for b in range(0, nBatchPerEpoch*nBatch, nBatch): cost, predictions = trainStep1(xTrain[b : b + nBatch], labelTrainF[b : b + nBatch]) cost_cumul += cost accuracy.update_state(labelTrainF[b : b + nBatch], predictions) cost_epoch = cost_cumul / nBatchPerEpoch W = model1.get_weights() histo.append((cost_epoch.numpy(), accuracy.result().numpy(), W[1][0], W[0])) print("Predicted model: {b:.3f} + {w:.3f} x, num epochs={c}".format(w=W[0], b=W[1][0], c=len(histo))) # Save history as a Panda Data Frame df = pandas.DataFrame(histo, columns = ('cost', 'accuracy', 'b', 'w0')) 

plt.scatter(df['b'], df['w0'], marker='.', alpha=0.2); plt.xlabel('intercept') plt.ylabel('weight'); 

fig, ax = plt.subplots(1,2, figsize=(16, 4)) ax[0].plot(df['cost']) ax[0].grid() ax[1].plot(df['accuracy']) ax[1].grid()