-
diff --git a/doc/pub/LogReg/html/._LogReg-bs017.html b/doc/pub/LogReg/html/._LogReg-bs017.html new file mode 100644 index 000000000..1698f4b23 --- /dev/null +++ b/doc/pub/LogReg/html/._LogReg-bs017.html @@ -0,0 +1,247 @@ + + + + + + + +
+ + + + +
+The Hamiltonian of the two-dimensional Ising model without an external field for a constant coupling constant \( J \) is given by +$$ +\begin{align} + H = -J \sum_{\langle ij\rangle} S_i S_j, +\tag{2} +\end{align} +$$ + +where \( S_i \in \{-1, 1\} \) and \( \langle ij \rangle \) signifies that we only iterate over the nearest neighbors in the lattice. We will be looking at a system of \( L = 40 \) spins in each dimension, i.e., \( L^2 = 1600 \) spins in total. Opposed to the one-dimensional Ising model we will get a phase transition from an ordered phase to a disordered phase at the critical temperature + +$$ +\begin{align} + \frac{T_c}{J} = \frac{2}{\log\left(1 + \sqrt{2}\right)} \approx 2.26, +\tag{3} +\end{align} +$$ + +as shown by Lars Onsager. + +
+Here we use logistic regression to predict when a phase transition
+occurs. The data we will look at is a set of spin configurations,
+i.e., individual lattices with spins, labeled ordered 1 or
+disordered 0. Our job is to build a model which will take in a
+spin configuration and predict whether or not the spin configuration
+constitutes an ordered or a disordered phase. To achieve this we will
+represent the lattices as flattened arrays with \( 1600 \) elements
+instead of a matrix of \( 40 \times 40 \) elements. As an extra test of
+the performance of the algorithms we will divide the dataset into
+three pieces. We will do a conventional train-test-split on a
+combination of totally ordered and totally disordered phases. The
+remaining "critical-like" states will be used as test data which we
+hope the model will be able to make good extrapolated predictions on.
+
+
+ + +
import pickle
+import os
+import glob
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import seaborn as sns
+import sklearn.model_selection as skms
+import sklearn.linear_model as skl
+import sklearn.metrics as skm
+import tqdm
+import copy
+import time
+from IPython.display import display
+
+%matplotlib inline
+
+sns.set(color_codes=True)
++
+ +
+ + +
+ + + + +
+Using the data from Mehta et al. (specifically the two datasets named Ising2DFM_reSample_L40_T=All.pkl and Ising2DFM_reSample_L40_T=All_labels.pkl) we have to unpack the data into numpy arrays.
+
+
+ + +
filenames = glob.glob(os.path.join("..", "dat", "*"))
+label_filename = list(filter(lambda x: "label" in x, filenames))[0]
+dat_filename = list(filter(lambda x: "label" not in x, filenames))[0]
+
+# Read in the labels
+with open(label_filename, "rb") as f:
+ labels = pickle.load(f)
+
+# Read in the corresponding configurations
+with open(dat_filename, "rb") as f:
+ data = np.unpackbits(pickle.load(f)).reshape(-1, 1600).astype("int")
+
+# Set spin-down to -1
+data[data == 0] = -1
++This dataset consists of \( 10000 \) samples, i.e., \( 10000 \) spin +configurations with \( 40 \times 40 \) spins each, for \( 16 \) temperatures +between \( 0.25 \) to \( 4.0 \). Next we create a train/test-split and keep +the data in the critical phase as a separate dataset for +extrapolation-testing. + +
+ + +
# Set up slices of the dataset
+ordered = slice(0, 70000)
+critical = slice(70000, 100000)
+disordered = slice(100000, 160000)
+
+X_train, X_test, y_train, y_test = skms.train_test_split(
+ np.concatenate((data[ordered], data[disordered])),
+ np.concatenate((labels[ordered], labels[disordered])),
+ test_size=0.95
+)
++Using a small training set yields a better accuracy. This will be discussed in the end. + +
+
+ +
+ + +
+ + + + +
+Logistic regression is a linear model for classification. Recalling +the cost function for ordinary least squares with both L2 (ridge) and +L1 (LASSO) penalties we will see that the logistic cost function is +very similar. In OLS we wish to predict a continuous variable +\( \hat{y} \) using +$$ +\begin{align} + \hat{y} = X\omega, +\tag{4} +\end{align} +$$ + +
+where \( X \in \mathbb{R}^{n \times p} \) is the input data and \( \omega^{p +\times d} \) are the weights of the regression. In a classification +setting (binary classification in our situation) we are interested in +a positive or negative answer. We can thus define either answer to be +above or below some threshold. But, in order to limit the size of the +answer and also to get a probability interpretation on how sure we are +for either answer we can compute the sigmoid function of OLS. That is, + +$$ +\begin{align} + f(X\omega) = \frac{1}{1 + \exp(-X\omega)}. +\tag{5} +\end{align} +$$ + +We are thus interested in minizming the following cost function +$$ +\begin{align} + C(X, \omega) = \sum_{i = 1}^n \left\{ + - y_i\log\left( f(x_i^T\omega) \right) + - (1 - y_i)\log\left[1 - f(x_i^T\omega)\right] + \right\}, +\tag{6} +\end{align} +$$ + +
+where we will restrict ourselves to a value for \( f(z) \) as the sigmoid +described above. We can also tack on a L2 (Ridge) or L1 (LASSO) +penalization to this cost function in the same manner we did for +linear regression. + +
+
+ +
+ + +
+ + + + +
+The penalization factor \( \lambda \) is inverted in the case of the +logistic regression model we use. We will explore several values of +\( \lambda \) using both L1 and L2 penalization. We do this using a grid +search over different parameters and run a 3-fold cross validation for +each configuration. In other words, we fit a model 3 times for each +configuration of the hyper parameters. + +
+ + +
lambdas = np.logspace(-7, -1, 7)
+
+param_grid = {
+ "C": list(1.0/lambdas),
+ "penalty": ["l1", "l2"]
+}
+clf = skms.GridSearchCV(
+ skl.LogisticRegression(),
+ param_grid=param_grid,
+ n_jobs=-1,
+ return_train_score=True
+)
+t0 = time.time()
+clf.fit(X_train, y_train)
+t1 = time.time()
+
+print (
+ "Time spent fitting GridSearchCV(LogisticRegression): {0:.3f} sec".format(
+ t1 - t0
+ )
+)
++We can see that logistic regression is quite slow and using the grid +search and cross validation results in quite a heavy +computation. Below we show the results of the different +configurations. + +
+ + +
logreg_df = pd.DataFrame(clf.cv_results_)
+
+display(logreg_df)
++
+ +
+ + +
+ + + + +
+To determine how well a classification model is performing we count +the number of correctly labeled classes and divide by the number of +classes in total. The accuracy is thus given by + +$$ +\begin{align} + a(y, \hat{y}) = \frac{1}{n}\sum_{i = 1}^{n} I(y_i = \hat{y}_i), +\tag{7} +\end{align} +$$ + +
+where \( I(y_i = \hat{y}_i) \) is the indicator function given by + +$$ +\begin{align} + I(x = y) = \begin{cases} + 1 & x = y, +\tag{8}\\ + 0 & x \neq y. + \end{cases} +\tag{9} +\end{align} +$$ + +
+This is the accuracy provided by Scikit-learn when using sklearn.metrics.accuracyscore. + +
+Below we compute the accuracy of the best fit model on the training data (which should give a good accuracy), the test data (which has not been shown to the model) and the critical data (completely new data that needs to be extrapolated). + +
+ + +
train_accuracy = skm.accuracy_score(y_train, clf.predict(X_train))
+test_accuracy = skm.accuracy_score(y_test, clf.predict(X_test))
+critical_accuracy = skm.accuracy_score(labels[critical], clf.predict(data[critical]))
+
+print ("Accuracy on train data: {0}".format(train_accuracy))
+print ("Accuracy on test data: {0}".format(test_accuracy))
+print ("Accuracy on critical data: {0}".format(critical_accuracy))
++We can see that we get quite good accuracy on the training data, but gradually worsening accuracy on the test and critical data. + +
+
+ +
+ + +
+ + + + +
+Below we show a different metric for determining the quality of our +model, namely the reciever operating characteristic (ROC). The ROC +curve tells us how well the model correctly classifies the different +labels. We plot the true positive rate (the rate of predicted +positive classes that are positive) versus the false positive rate +(the rate of predicted positive classes that are negative). The ROC +curve is built by computing the true positive rate and the false +positive rate for varying thresholds, i.e, which probability we +should acredit a certain class. + +
+By computing the area under the curve (AUC) of the ROC curve we get an estimate of how well our model is performing. Pure guessing will get an AUC of \( 0.5 \). A perfect score will get an AUC of \( 1.0 \). + +
+ + +
fig = plt.figure(figsize=(20, 14))
+
+for (_X, _y), label in zip(
+ [
+ (X_train, y_train),
+ (X_test, y_test),
+ (data[critical], labels[critical])
+ ],
+ ["Train", "Test", "Critical"]
+):
+ proba = clf.predict_proba(_X)
+ fpr, tpr, _ = skm.roc_curve(_y, proba[:, 1])
+ roc_auc = skm.auc(fpr, tpr)
+
+ print ("LogisticRegression AUC ({0}): {1}".format(label, roc_auc))
+
+ plt.plot(fpr, tpr, label="{0} (AUC = {1})".format(label, roc_auc), linewidth=4.0)
+
+plt.plot([0, 1], [0, 1], "--", label="Guessing (AUC = 0.5)", linewidth=4.0)
+
+plt.title(r"The ROC curve for LogisticRegression", fontsize=18)
+plt.xlabel(r"False positive rate", fontsize=18)
+plt.ylabel(r"True positive rate", fontsize=18)
+plt.axis([-0.01, 1.01, -0.01, 1.01])
+plt.xticks(fontsize=18)
+plt.yticks(fontsize=18)
+plt.legend(loc="best", fontsize=18)
+plt.show()
++We can see that this plot of the ROC looks very strange. This tells us +that logistic regression is quite inept at predicting the Ising model +transition and is therefore highly non-linear. The ROC curve for the +training data looks quite good, but as the testing data is so far off +we see that we are dealing with an overfit model. + +
+A previous run with \( 50\% \) of the data used for training yielded a +worse performance than using a smaller training set. This again gives +confidence to the fact that logistic regression is not able to +correctly fit the Ising model as it is not a linear model. + +
+ +
+ + ++The Hamiltonian of the two-dimensional Ising model without an external field for a constant coupling constant \( J \) is given by +
+$$
+\begin{align}
+ H = -J \sum_{\langle ij\rangle} S_i S_j,
+\tag{2}
+\end{align}
+$$
+
+
+where \( S_i \in \{-1, 1\} \) and \( \langle ij \rangle \) signifies that we only iterate over the nearest neighbors in the lattice. We will be looking at a system of \( L = 40 \) spins in each dimension, i.e., \( L^2 = 1600 \) spins in total. Opposed to the one-dimensional Ising model we will get a phase transition from an ordered phase to a disordered phase at the critical temperature
+
+
+$$
+\begin{align}
+ \frac{T_c}{J} = \frac{2}{\log\left(1 + \sqrt{2}\right)} \approx 2.26,
+\tag{3}
+\end{align}
+$$
+
+
+as shown by Lars Onsager.
+
+
+Here we use logistic regression to predict when a phase transition
+occurs. The data we will look at is a set of spin configurations,
+i.e., individual lattices with spins, labeled ordered 1 or
+disordered 0. Our job is to build a model which will take in a
+spin configuration and predict whether or not the spin configuration
+constitutes an ordered or a disordered phase. To achieve this we will
+represent the lattices as flattened arrays with \( 1600 \) elements
+instead of a matrix of \( 40 \times 40 \) elements. As an extra test of
+the performance of the algorithms we will divide the dataset into
+three pieces. We will do a conventional train-test-split on a
+combination of totally ordered and totally disordered phases. The
+remaining "critical-like" states will be used as test data which we
+hope the model will be able to make good extrapolated predictions on.
+
+
+ + +
import pickle
+import os
+import glob
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import seaborn as sns
+import sklearn.model_selection as skms
+import sklearn.linear_model as skl
+import sklearn.metrics as skm
+import tqdm
+import copy
+import time
+from IPython.display import display
+
+%matplotlib inline
+
+sns.set(color_codes=True)
+
+Using the data from Mehta et al. (specifically the two datasets named Ising2DFM_reSample_L40_T=All.pkl and Ising2DFM_reSample_L40_T=All_labels.pkl) we have to unpack the data into numpy arrays.
+
+
+ + +
filenames = glob.glob(os.path.join("..", "dat", "*"))
+label_filename = list(filter(lambda x: "label" in x, filenames))[0]
+dat_filename = list(filter(lambda x: "label" not in x, filenames))[0]
+
+# Read in the labels
+with open(label_filename, "rb") as f:
+ labels = pickle.load(f)
+
+# Read in the corresponding configurations
+with open(dat_filename, "rb") as f:
+ data = np.unpackbits(pickle.load(f)).reshape(-1, 1600).astype("int")
+
+# Set spin-down to -1
+data[data == 0] = -1
++This dataset consists of \( 10000 \) samples, i.e., \( 10000 \) spin +configurations with \( 40 \times 40 \) spins each, for \( 16 \) temperatures +between \( 0.25 \) to \( 4.0 \). Next we create a train/test-split and keep +the data in the critical phase as a separate dataset for +extrapolation-testing. + +
+ + +
# Set up slices of the dataset
+ordered = slice(0, 70000)
+critical = slice(70000, 100000)
+disordered = slice(100000, 160000)
+
+X_train, X_test, y_train, y_test = skms.train_test_split(
+ np.concatenate((data[ordered], data[disordered])),
+ np.concatenate((labels[ordered], labels[disordered])),
+ test_size=0.95
+)
++Using a small training set yields a better accuracy. This will be discussed in the end. +
+Logistic regression is a linear model for classification. Recalling +the cost function for ordinary least squares with both L2 (ridge) and +L1 (LASSO) penalties we will see that the logistic cost function is +very similar. In OLS we wish to predict a continuous variable +\( \hat{y} \) using +
+$$
+\begin{align}
+ \hat{y} = X\omega,
+\tag{4}
+\end{align}
+$$
+
+
+
+where \( X \in \mathbb{R}^{n \times p} \) is the input data and \( \omega^{p +\times d} \) are the weights of the regression. In a classification +setting (binary classification in our situation) we are interested in +a positive or negative answer. We can thus define either answer to be +above or below some threshold. But, in order to limit the size of the +answer and also to get a probability interpretation on how sure we are +for either answer we can compute the sigmoid function of OLS. That is, + +
+$$
+\begin{align}
+ f(X\omega) = \frac{1}{1 + \exp(-X\omega)}.
+\tag{5}
+\end{align}
+$$
+
+
+We are thus interested in minizming the following cost function
+
+$$
+\begin{align}
+ C(X, \omega) = \sum_{i = 1}^n \left\{
+ - y_i\log\left( f(x_i^T\omega) \right)
+ - (1 - y_i)\log\left[1 - f(x_i^T\omega)\right]
+ \right\},
+\tag{6}
+\end{align}
+$$
+
+
+
+where we will restrict ourselves to a value for \( f(z) \) as the sigmoid +described above. We can also tack on a L2 (Ridge) or L1 (LASSO) +penalization to this cost function in the same manner we did for +linear regression. +
+The penalization factor \( \lambda \) is inverted in the case of the +logistic regression model we use. We will explore several values of +\( \lambda \) using both L1 and L2 penalization. We do this using a grid +search over different parameters and run a 3-fold cross validation for +each configuration. In other words, we fit a model 3 times for each +configuration of the hyper parameters. + +
+ + +
lambdas = np.logspace(-7, -1, 7)
+
+param_grid = {
+ "C": list(1.0/lambdas),
+ "penalty": ["l1", "l2"]
+}
+clf = skms.GridSearchCV(
+ skl.LogisticRegression(),
+ param_grid=param_grid,
+ n_jobs=-1,
+ return_train_score=True
+)
+t0 = time.time()
+clf.fit(X_train, y_train)
+t1 = time.time()
+
+print (
+ "Time spent fitting GridSearchCV(LogisticRegression): {0:.3f} sec".format(
+ t1 - t0
+ )
+)
++We can see that logistic regression is quite slow and using the grid +search and cross validation results in quite a heavy +computation. Below we show the results of the different +configurations. + +
+ + +
logreg_df = pd.DataFrame(clf.cv_results_)
+
+display(logreg_df)
++To determine how well a classification model is performing we count +the number of correctly labeled classes and divide by the number of +classes in total. The accuracy is thus given by + +
+$$
+\begin{align}
+ a(y, \hat{y}) = \frac{1}{n}\sum_{i = 1}^{n} I(y_i = \hat{y}_i),
+\tag{7}
+\end{align}
+$$
+
+
+
+where \( I(y_i = \hat{y}_i) \) is the indicator function given by + +
+$$
+\begin{align}
+ I(x = y) = \begin{cases}
+ 1 & x = y,
+\tag{8}\\
+ 0 & x \neq y.
+ \end{cases}
+\tag{9}
+\end{align}
+$$
+
+
+
+This is the accuracy provided by Scikit-learn when using sklearn.metrics.accuracyscore. + +
+Below we compute the accuracy of the best fit model on the training data (which should give a good accuracy), the test data (which has not been shown to the model) and the critical data (completely new data that needs to be extrapolated). + +
+ + +
train_accuracy = skm.accuracy_score(y_train, clf.predict(X_train))
+test_accuracy = skm.accuracy_score(y_test, clf.predict(X_test))
+critical_accuracy = skm.accuracy_score(labels[critical], clf.predict(data[critical]))
+
+print ("Accuracy on train data: {0}".format(train_accuracy))
+print ("Accuracy on test data: {0}".format(test_accuracy))
+print ("Accuracy on critical data: {0}".format(critical_accuracy))
++We can see that we get quite good accuracy on the training data, but gradually worsening accuracy on the test and critical data. +
+Below we show a different metric for determining the quality of our +model, namely the reciever operating characteristic (ROC). The ROC +curve tells us how well the model correctly classifies the different +labels. We plot the true positive rate (the rate of predicted +positive classes that are positive) versus the false positive rate +(the rate of predicted positive classes that are negative). The ROC +curve is built by computing the true positive rate and the false +positive rate for varying thresholds, i.e, which probability we +should acredit a certain class. + +
+By computing the area under the curve (AUC) of the ROC curve we get an estimate of how well our model is performing. Pure guessing will get an AUC of \( 0.5 \). A perfect score will get an AUC of \( 1.0 \). + +
+ + +
fig = plt.figure(figsize=(20, 14))
+
+for (_X, _y), label in zip(
+ [
+ (X_train, y_train),
+ (X_test, y_test),
+ (data[critical], labels[critical])
+ ],
+ ["Train", "Test", "Critical"]
+):
+ proba = clf.predict_proba(_X)
+ fpr, tpr, _ = skm.roc_curve(_y, proba[:, 1])
+ roc_auc = skm.auc(fpr, tpr)
+
+ print ("LogisticRegression AUC ({0}): {1}".format(label, roc_auc))
+
+ plt.plot(fpr, tpr, label="{0} (AUC = {1})".format(label, roc_auc), linewidth=4.0)
+
+plt.plot([0, 1], [0, 1], "--", label="Guessing (AUC = 0.5)", linewidth=4.0)
+
+plt.title(r"The ROC curve for LogisticRegression", fontsize=18)
+plt.xlabel(r"False positive rate", fontsize=18)
+plt.ylabel(r"True positive rate", fontsize=18)
+plt.axis([-0.01, 1.01, -0.01, 1.01])
+plt.xticks(fontsize=18)
+plt.yticks(fontsize=18)
+plt.legend(loc="best", fontsize=18)
+plt.show()
++We can see that this plot of the ROC looks very strange. This tells us +that logistic regression is quite inept at predicting the Ising model +transition and is therefore highly non-linear. The ROC curve for the +training data looks quite good, but as the testing data is so far off +we see that we are dealing with an overfit model. + +
+A previous run with \( 50\% \) of the data used for training yielded a +worse performance than using a smaller training set. This again gives +confidence to the fact that logistic regression is not able to +correctly fit the Ising model as it is not a linear model. +
+ + +
+The Hamiltonian of the two-dimensional Ising model without an external field for a constant coupling constant \( J \) is given by +$$ +\begin{align} + H = -J \sum_{\langle ij\rangle} S_i S_j, +\label{_auto2} +\end{align} +$$ + +where \( S_i \in \{-1, 1\} \) and \( \langle ij \rangle \) signifies that we only iterate over the nearest neighbors in the lattice. We will be looking at a system of \( L = 40 \) spins in each dimension, i.e., \( L^2 = 1600 \) spins in total. Opposed to the one-dimensional Ising model we will get a phase transition from an ordered phase to a disordered phase at the critical temperature + +$$ +\begin{align} + \frac{T_c}{J} = \frac{2}{\log\left(1 + \sqrt{2}\right)} \approx 2.26, +\label{_auto3} +\end{align} +$$ + +as shown by Lars Onsager. + +
+Here we use logistic regression to predict when a phase transition
+occurs. The data we will look at is a set of spin configurations,
+i.e., individual lattices with spins, labeled ordered 1 or
+disordered 0. Our job is to build a model which will take in a
+spin configuration and predict whether or not the spin configuration
+constitutes an ordered or a disordered phase. To achieve this we will
+represent the lattices as flattened arrays with \( 1600 \) elements
+instead of a matrix of \( 40 \times 40 \) elements. As an extra test of
+the performance of the algorithms we will divide the dataset into
+three pieces. We will do a conventional train-test-split on a
+combination of totally ordered and totally disordered phases. The
+remaining "critical-like" states will be used as test data which we
+hope the model will be able to make good extrapolated predictions on.
+
+
+ + +
import pickle
+import os
+import glob
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import seaborn as sns
+import sklearn.model_selection as skms
+import sklearn.linear_model as skl
+import sklearn.metrics as skm
+import tqdm
+import copy
+import time
+from IPython.display import display
+
+%matplotlib inline
+
+sns.set(color_codes=True)
+
+
+
+
+Using the data from Mehta et al. (specifically the two datasets named Ising2DFM_reSample_L40_T=All.pkl and Ising2DFM_reSample_L40_T=All_labels.pkl) we have to unpack the data into numpy arrays.
+
+
+ + +
filenames = glob.glob(os.path.join("..", "dat", "*"))
+label_filename = list(filter(lambda x: "label" in x, filenames))[0]
+dat_filename = list(filter(lambda x: "label" not in x, filenames))[0]
+
+# Read in the labels
+with open(label_filename, "rb") as f:
+ labels = pickle.load(f)
+
+# Read in the corresponding configurations
+with open(dat_filename, "rb") as f:
+ data = np.unpackbits(pickle.load(f)).reshape(-1, 1600).astype("int")
+
+# Set spin-down to -1
+data[data == 0] = -1
++This dataset consists of \( 10000 \) samples, i.e., \( 10000 \) spin +configurations with \( 40 \times 40 \) spins each, for \( 16 \) temperatures +between \( 0.25 \) to \( 4.0 \). Next we create a train/test-split and keep +the data in the critical phase as a separate dataset for +extrapolation-testing. + +
+ + +
# Set up slices of the dataset
+ordered = slice(0, 70000)
+critical = slice(70000, 100000)
+disordered = slice(100000, 160000)
+
+X_train, X_test, y_train, y_test = skms.train_test_split(
+ np.concatenate((data[ordered], data[disordered])),
+ np.concatenate((labels[ordered], labels[disordered])),
+ test_size=0.95
+)
++Using a small training set yields a better accuracy. This will be discussed in the end. + +
+
+
+
+Logistic regression is a linear model for classification. Recalling +the cost function for ordinary least squares with both L2 (ridge) and +L1 (LASSO) penalties we will see that the logistic cost function is +very similar. In OLS we wish to predict a continuous variable +\( \hat{y} \) using +$$ +\begin{align} + \hat{y} = X\omega, +\label{_auto4} +\end{align} +$$ + +
+where \( X \in \mathbb{R}^{n \times p} \) is the input data and \( \omega^{p +\times d} \) are the weights of the regression. In a classification +setting (binary classification in our situation) we are interested in +a positive or negative answer. We can thus define either answer to be +above or below some threshold. But, in order to limit the size of the +answer and also to get a probability interpretation on how sure we are +for either answer we can compute the sigmoid function of OLS. That is, + +$$ +\begin{align} + f(X\omega) = \frac{1}{1 + \exp(-X\omega)}. +\label{_auto5} +\end{align} +$$ + +We are thus interested in minizming the following cost function +$$ +\begin{align} + C(X, \omega) = \sum_{i = 1}^n \left\{ + - y_i\log\left( f(x_i^T\omega) \right) + - (1 - y_i)\log\left[1 - f(x_i^T\omega)\right] + \right\}, +\label{_auto6} +\end{align} +$$ + +
+where we will restrict ourselves to a value for \( f(z) \) as the sigmoid +described above. We can also tack on a L2 (Ridge) or L1 (LASSO) +penalization to this cost function in the same manner we did for +linear regression. + +
+
+
+
+The penalization factor \( \lambda \) is inverted in the case of the +logistic regression model we use. We will explore several values of +\( \lambda \) using both L1 and L2 penalization. We do this using a grid +search over different parameters and run a 3-fold cross validation for +each configuration. In other words, we fit a model 3 times for each +configuration of the hyper parameters. + +
+ + +
lambdas = np.logspace(-7, -1, 7)
+
+param_grid = {
+ "C": list(1.0/lambdas),
+ "penalty": ["l1", "l2"]
+}
+clf = skms.GridSearchCV(
+ skl.LogisticRegression(),
+ param_grid=param_grid,
+ n_jobs=-1,
+ return_train_score=True
+)
+t0 = time.time()
+clf.fit(X_train, y_train)
+t1 = time.time()
+
+print (
+ "Time spent fitting GridSearchCV(LogisticRegression): {0:.3f} sec".format(
+ t1 - t0
+ )
+)
++We can see that logistic regression is quite slow and using the grid +search and cross validation results in quite a heavy +computation. Below we show the results of the different +configurations. + +
+ + +
logreg_df = pd.DataFrame(clf.cv_results_)
+
+display(logreg_df)
+
+
+
+
+To determine how well a classification model is performing we count +the number of correctly labeled classes and divide by the number of +classes in total. The accuracy is thus given by + +$$ +\begin{align} + a(y, \hat{y}) = \frac{1}{n}\sum_{i = 1}^{n} I(y_i = \hat{y}_i), +\label{_auto7} +\end{align} +$$ + +
+where \( I(y_i = \hat{y}_i) \) is the indicator function given by + +$$ +\begin{align} + I(x = y) = \begin{cases} + 1 & x = y, +\label{_auto8}\\ + 0 & x \neq y. + \end{cases} +\label{_auto9} +\end{align} +$$ + +
+This is the accuracy provided by Scikit-learn when using sklearn.metrics.accuracyscore. + +
+Below we compute the accuracy of the best fit model on the training data (which should give a good accuracy), the test data (which has not been shown to the model) and the critical data (completely new data that needs to be extrapolated). + +
+ + +
train_accuracy = skm.accuracy_score(y_train, clf.predict(X_train))
+test_accuracy = skm.accuracy_score(y_test, clf.predict(X_test))
+critical_accuracy = skm.accuracy_score(labels[critical], clf.predict(data[critical]))
+
+print ("Accuracy on train data: {0}".format(train_accuracy))
+print ("Accuracy on test data: {0}".format(test_accuracy))
+print ("Accuracy on critical data: {0}".format(critical_accuracy))
++We can see that we get quite good accuracy on the training data, but gradually worsening accuracy on the test and critical data. + +
+
+
+
+Below we show a different metric for determining the quality of our +model, namely the reciever operating characteristic (ROC). The ROC +curve tells us how well the model correctly classifies the different +labels. We plot the true positive rate (the rate of predicted +positive classes that are positive) versus the false positive rate +(the rate of predicted positive classes that are negative). The ROC +curve is built by computing the true positive rate and the false +positive rate for varying thresholds, i.e, which probability we +should acredit a certain class. + +
+By computing the area under the curve (AUC) of the ROC curve we get an estimate of how well our model is performing. Pure guessing will get an AUC of \( 0.5 \). A perfect score will get an AUC of \( 1.0 \). + +
+ + +
fig = plt.figure(figsize=(20, 14))
+
+for (_X, _y), label in zip(
+ [
+ (X_train, y_train),
+ (X_test, y_test),
+ (data[critical], labels[critical])
+ ],
+ ["Train", "Test", "Critical"]
+):
+ proba = clf.predict_proba(_X)
+ fpr, tpr, _ = skm.roc_curve(_y, proba[:, 1])
+ roc_auc = skm.auc(fpr, tpr)
+
+ print ("LogisticRegression AUC ({0}): {1}".format(label, roc_auc))
+
+ plt.plot(fpr, tpr, label="{0} (AUC = {1})".format(label, roc_auc), linewidth=4.0)
+
+plt.plot([0, 1], [0, 1], "--", label="Guessing (AUC = 0.5)", linewidth=4.0)
+
+plt.title(r"The ROC curve for LogisticRegression", fontsize=18)
+plt.xlabel(r"False positive rate", fontsize=18)
+plt.ylabel(r"True positive rate", fontsize=18)
+plt.axis([-0.01, 1.01, -0.01, 1.01])
+plt.xticks(fontsize=18)
+plt.yticks(fontsize=18)
+plt.legend(loc="best", fontsize=18)
+plt.show()
++We can see that this plot of the ROC looks very strange. This tells us +that logistic regression is quite inept at predicting the Ising model +transition and is therefore highly non-linear. The ROC curve for the +training data looks quite good, but as the testing data is so far off +we see that we are dealing with an overfit model. + +
+A previous run with \( 50\% \) of the data used for training yielded a +worse performance than using a smaller training set. This again gives +confidence to the fact that logistic regression is not able to +correctly fit the Ising model as it is not a linear model. diff --git a/doc/pub/LogReg/html/LogReg.html b/doc/pub/LogReg/html/LogReg.html index 3a86df613..17d00292d 100644 --- a/doc/pub/LogReg/html/LogReg.html +++ b/doc/pub/LogReg/html/LogReg.html @@ -54,7 +54,17 @@ div { text-align: justify; text-justify: inter-word; } ('Including more classes', 2, None, '___sec12'), ('The Softmax function', 2, None, '___sec13'), ('A _scikit-learn_ example', 2, None, '___sec14'), - ('A simple classification problem', 2, None, '___sec15')]} + ('A simple classification problem', 2, None, '___sec15'), + ('The two-dimensional Ising model, Predicting phase transition ' + 'of the two-dimensional Ising model', + 2, + None, + '___sec16'), + ('Reading in the data', 2, None, '___sec17'), + ('Logistic regression', 2, None, '___sec18'), + ('Exploring the logistic regression', 2, None, '___sec19'), + ('Accuracy of a classification model', 2, None, '___sec20'), + ('Analyzing the results', 2, None, '___sec21')]} end of tocinfo -->
@@ -515,6 +525,333 @@ plt.show() main()+ + +
+The Hamiltonian of the two-dimensional Ising model without an external field for a constant coupling constant \( J \) is given by +$$ +\begin{align} + H = -J \sum_{\langle ij\rangle} S_i S_j, +\label{_auto2} +\end{align} +$$ + +where \( S_i \in \{-1, 1\} \) and \( \langle ij \rangle \) signifies that we only iterate over the nearest neighbors in the lattice. We will be looking at a system of \( L = 40 \) spins in each dimension, i.e., \( L^2 = 1600 \) spins in total. Opposed to the one-dimensional Ising model we will get a phase transition from an ordered phase to a disordered phase at the critical temperature + +$$ +\begin{align} + \frac{T_c}{J} = \frac{2}{\log\left(1 + \sqrt{2}\right)} \approx 2.26, +\label{_auto3} +\end{align} +$$ + +as shown by Lars Onsager. + +
+Here we use logistic regression to predict when a phase transition
+occurs. The data we will look at is a set of spin configurations,
+i.e., individual lattices with spins, labeled ordered 1 or
+disordered 0. Our job is to build a model which will take in a
+spin configuration and predict whether or not the spin configuration
+constitutes an ordered or a disordered phase. To achieve this we will
+represent the lattices as flattened arrays with \( 1600 \) elements
+instead of a matrix of \( 40 \times 40 \) elements. As an extra test of
+the performance of the algorithms we will divide the dataset into
+three pieces. We will do a conventional train-test-split on a
+combination of totally ordered and totally disordered phases. The
+remaining "critical-like" states will be used as test data which we
+hope the model will be able to make good extrapolated predictions on.
+
+
+ + +
import pickle
+import os
+import glob
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import seaborn as sns
+import sklearn.model_selection as skms
+import sklearn.linear_model as skl
+import sklearn.metrics as skm
+import tqdm
+import copy
+import time
+from IPython.display import display
+
+%matplotlib inline
+
+sns.set(color_codes=True)
+
+
+
+
+Using the data from Mehta et al. (specifically the two datasets named Ising2DFM_reSample_L40_T=All.pkl and Ising2DFM_reSample_L40_T=All_labels.pkl) we have to unpack the data into numpy arrays.
+
+
+ + +
filenames = glob.glob(os.path.join("..", "dat", "*"))
+label_filename = list(filter(lambda x: "label" in x, filenames))[0]
+dat_filename = list(filter(lambda x: "label" not in x, filenames))[0]
+
+# Read in the labels
+with open(label_filename, "rb") as f:
+ labels = pickle.load(f)
+
+# Read in the corresponding configurations
+with open(dat_filename, "rb") as f:
+ data = np.unpackbits(pickle.load(f)).reshape(-1, 1600).astype("int")
+
+# Set spin-down to -1
+data[data == 0] = -1
++This dataset consists of \( 10000 \) samples, i.e., \( 10000 \) spin +configurations with \( 40 \times 40 \) spins each, for \( 16 \) temperatures +between \( 0.25 \) to \( 4.0 \). Next we create a train/test-split and keep +the data in the critical phase as a separate dataset for +extrapolation-testing. + +
+ + +
# Set up slices of the dataset
+ordered = slice(0, 70000)
+critical = slice(70000, 100000)
+disordered = slice(100000, 160000)
+
+X_train, X_test, y_train, y_test = skms.train_test_split(
+ np.concatenate((data[ordered], data[disordered])),
+ np.concatenate((labels[ordered], labels[disordered])),
+ test_size=0.95
+)
++Using a small training set yields a better accuracy. This will be discussed in the end. + +
+
+
+
+Logistic regression is a linear model for classification. Recalling +the cost function for ordinary least squares with both L2 (ridge) and +L1 (LASSO) penalties we will see that the logistic cost function is +very similar. In OLS we wish to predict a continuous variable +\( \hat{y} \) using +$$ +\begin{align} + \hat{y} = X\omega, +\label{_auto4} +\end{align} +$$ + +
+where \( X \in \mathbb{R}^{n \times p} \) is the input data and \( \omega^{p +\times d} \) are the weights of the regression. In a classification +setting (binary classification in our situation) we are interested in +a positive or negative answer. We can thus define either answer to be +above or below some threshold. But, in order to limit the size of the +answer and also to get a probability interpretation on how sure we are +for either answer we can compute the sigmoid function of OLS. That is, + +$$ +\begin{align} + f(X\omega) = \frac{1}{1 + \exp(-X\omega)}. +\label{_auto5} +\end{align} +$$ + +We are thus interested in minizming the following cost function +$$ +\begin{align} + C(X, \omega) = \sum_{i = 1}^n \left\{ + - y_i\log\left( f(x_i^T\omega) \right) + - (1 - y_i)\log\left[1 - f(x_i^T\omega)\right] + \right\}, +\label{_auto6} +\end{align} +$$ + +
+where we will restrict ourselves to a value for \( f(z) \) as the sigmoid +described above. We can also tack on a L2 (Ridge) or L1 (LASSO) +penalization to this cost function in the same manner we did for +linear regression. + +
+
+
+
+The penalization factor \( \lambda \) is inverted in the case of the +logistic regression model we use. We will explore several values of +\( \lambda \) using both L1 and L2 penalization. We do this using a grid +search over different parameters and run a 3-fold cross validation for +each configuration. In other words, we fit a model 3 times for each +configuration of the hyper parameters. + +
+ + +
lambdas = np.logspace(-7, -1, 7)
+
+param_grid = {
+ "C": list(1.0/lambdas),
+ "penalty": ["l1", "l2"]
+}
+clf = skms.GridSearchCV(
+ skl.LogisticRegression(),
+ param_grid=param_grid,
+ n_jobs=-1,
+ return_train_score=True
+)
+t0 = time.time()
+clf.fit(X_train, y_train)
+t1 = time.time()
+
+print (
+ "Time spent fitting GridSearchCV(LogisticRegression): {0:.3f} sec".format(
+ t1 - t0
+ )
+)
++We can see that logistic regression is quite slow and using the grid +search and cross validation results in quite a heavy +computation. Below we show the results of the different +configurations. + +
+ + +
logreg_df = pd.DataFrame(clf.cv_results_)
+
+display(logreg_df)
+
+
+
+
+To determine how well a classification model is performing we count +the number of correctly labeled classes and divide by the number of +classes in total. The accuracy is thus given by + +$$ +\begin{align} + a(y, \hat{y}) = \frac{1}{n}\sum_{i = 1}^{n} I(y_i = \hat{y}_i), +\label{_auto7} +\end{align} +$$ + +
+where \( I(y_i = \hat{y}_i) \) is the indicator function given by + +$$ +\begin{align} + I(x = y) = \begin{cases} + 1 & x = y, +\label{_auto8}\\ + 0 & x \neq y. + \end{cases} +\label{_auto9} +\end{align} +$$ + +
+This is the accuracy provided by Scikit-learn when using sklearn.metrics.accuracyscore. + +
+Below we compute the accuracy of the best fit model on the training data (which should give a good accuracy), the test data (which has not been shown to the model) and the critical data (completely new data that needs to be extrapolated). + +
+ + +
train_accuracy = skm.accuracy_score(y_train, clf.predict(X_train))
+test_accuracy = skm.accuracy_score(y_test, clf.predict(X_test))
+critical_accuracy = skm.accuracy_score(labels[critical], clf.predict(data[critical]))
+
+print ("Accuracy on train data: {0}".format(train_accuracy))
+print ("Accuracy on test data: {0}".format(test_accuracy))
+print ("Accuracy on critical data: {0}".format(critical_accuracy))
++We can see that we get quite good accuracy on the training data, but gradually worsening accuracy on the test and critical data. + +
+
+
+
+Below we show a different metric for determining the quality of our +model, namely the reciever operating characteristic (ROC). The ROC +curve tells us how well the model correctly classifies the different +labels. We plot the true positive rate (the rate of predicted +positive classes that are positive) versus the false positive rate +(the rate of predicted positive classes that are negative). The ROC +curve is built by computing the true positive rate and the false +positive rate for varying thresholds, i.e, which probability we +should acredit a certain class. + +
+By computing the area under the curve (AUC) of the ROC curve we get an estimate of how well our model is performing. Pure guessing will get an AUC of \( 0.5 \). A perfect score will get an AUC of \( 1.0 \). + +
+ + +
fig = plt.figure(figsize=(20, 14))
+
+for (_X, _y), label in zip(
+ [
+ (X_train, y_train),
+ (X_test, y_test),
+ (data[critical], labels[critical])
+ ],
+ ["Train", "Test", "Critical"]
+):
+ proba = clf.predict_proba(_X)
+ fpr, tpr, _ = skm.roc_curve(_y, proba[:, 1])
+ roc_auc = skm.auc(fpr, tpr)
+
+ print ("LogisticRegression AUC ({0}): {1}".format(label, roc_auc))
+
+ plt.plot(fpr, tpr, label="{0} (AUC = {1})".format(label, roc_auc), linewidth=4.0)
+
+plt.plot([0, 1], [0, 1], "--", label="Guessing (AUC = 0.5)", linewidth=4.0)
+
+plt.title(r"The ROC curve for LogisticRegression", fontsize=18)
+plt.xlabel(r"False positive rate", fontsize=18)
+plt.ylabel(r"True positive rate", fontsize=18)
+plt.axis([-0.01, 1.01, -0.01, 1.01])
+plt.xticks(fontsize=18)
+plt.yticks(fontsize=18)
+plt.legend(loc="best", fontsize=18)
+plt.show()
++We can see that this plot of the ROC looks very strange. This tells us +that logistic regression is quite inept at predicting the Ising model +transition and is therefore highly non-linear. The ROC curve for the +training data looks quite good, but as the testing data is so far off +we see that we are dealing with an overfit model. + +
+A previous run with \( 50\% \) of the data used for training yielded a +worse performance than using a smaller training set. This again gives +confidence to the fact that logistic regression is not able to +correctly fit the Ising model as it is not a linear model. diff --git a/doc/pub/LogReg/ipynb/LogReg.ipynb b/doc/pub/LogReg/ipynb/LogReg.ipynb index cf3532b4d..835cad9c5 100644 --- a/doc/pub/LogReg/ipynb/LogReg.ipynb +++ b/doc/pub/LogReg/ipynb/LogReg.ipynb @@ -616,6 +616,489 @@ "if __name__ == \"__main__\":\n", " main()" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "## The two-dimensional Ising model, Predicting phase transition of the two-dimensional Ising model\n", + "\n", + "The Hamiltonian of the two-dimensional Ising model without an external field for a constant coupling constant $J$ is given by" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + " H = -J \\sum_{\\langle ij\\rangle} S_i S_j,\n", + "\\label{_auto2} \\tag{2}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $S_i \\in \\{-1, 1\\}$ and $\\langle ij \\rangle$ signifies that we only iterate over the nearest neighbors in the lattice. We will be looking at a system of $L = 40$ spins in each dimension, i.e., $L^2 = 1600$ spins in total. Opposed to the one-dimensional Ising model we will get a phase transition from an **ordered** phase to a **disordered** phase at the critical temperature" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "$$\n", + "\\begin{equation}\n", + " \\frac{T_c}{J} = \\frac{2}{\\log\\left(1 + \\sqrt{2}\\right)} \\approx 2.26,\n", + "\\label{_auto3} \\tag{3}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "as shown by Lars Onsager.\n", + "\n", + "\n", + "Here we use **logistic regression** to predict when a phase transition\n", + "occurs. The data we will look at is a set of spin configurations,\n", + "i.e., individual lattices with spins, labeled **ordered** `1` or\n", + "**disordered** `0`. Our job is to build a model which will take in a\n", + "spin configuration and predict whether or not the spin configuration\n", + "constitutes an ordered or a disordered phase. To achieve this we will\n", + "represent the lattices as flattened arrays with $1600$ elements\n", + "instead of a matrix of $40 \\times 40$ elements. As an extra test of\n", + "the performance of the algorithms we will divide the dataset into\n", + "three pieces. We will do a conventional train-test-split on a\n", + "combination of totally ordered and totally disordered phases. The\n", + "remaining \"critical-like\" states will be used as test data which we\n", + "hope the model will be able to make good extrapolated predictions on." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import pickle\n", + "import os\n", + "import glob\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import seaborn as sns\n", + "import sklearn.model_selection as skms\n", + "import sklearn.linear_model as skl\n", + "import sklearn.metrics as skm\n", + "import tqdm\n", + "import copy\n", + "import time\n", + "from IPython.display import display\n", + "\n", + "%matplotlib inline\n", + "\n", + "sns.set(color_codes=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Reading in the data\n", + "\n", + "Using the data from [Mehta et al.](https://physics.bu.edu/~pankajm/ML-Review-Datasets/isingMC/) (specifically the two datasets named `Ising2DFM_reSample_L40_T=All.pkl` and `Ising2DFM_reSample_L40_T=All_labels.pkl`) we have to unpack the data into numpy arrays." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "filenames = glob.glob(os.path.join(\"..\", \"dat\", \"*\"))\n", + "label_filename = list(filter(lambda x: \"label\" in x, filenames))[0]\n", + "dat_filename = list(filter(lambda x: \"label\" not in x, filenames))[0]\n", + "\n", + "# Read in the labels\n", + "with open(label_filename, \"rb\") as f:\n", + " labels = pickle.load(f)\n", + "\n", + "# Read in the corresponding configurations\n", + "with open(dat_filename, \"rb\") as f:\n", + " data = np.unpackbits(pickle.load(f)).reshape(-1, 1600).astype(\"int\")\n", + "\n", + "# Set spin-down to -1\n", + "data[data == 0] = -1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This dataset consists of $10000$ samples, i.e., $10000$ spin\n", + "configurations with $40 \\times 40$ spins each, for $16$ temperatures\n", + "between $0.25$ to $4.0$. Next we create a train/test-split and keep\n", + "the data in the critical phase as a separate dataset for\n", + "extrapolation-testing." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Set up slices of the dataset\n", + "ordered = slice(0, 70000)\n", + "critical = slice(70000, 100000)\n", + "disordered = slice(100000, 160000)\n", + "\n", + "X_train, X_test, y_train, y_test = skms.train_test_split(\n", + " np.concatenate((data[ordered], data[disordered])),\n", + " np.concatenate((labels[ordered], labels[disordered])),\n", + " test_size=0.95\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Using a small training set yields a better accuracy. This will be discussed in the end.\n", + "\n", + "## Logistic regression\n", + "\n", + "Logistic regression is a linear model for classification. Recalling\n", + "the cost function for ordinary least squares with both L2 (ridge) and\n", + "L1 (LASSO) penalties we will see that the logistic cost function is\n", + "very similar. In OLS we wish to predict a continuous variable\n", + "$\\hat{y}$ using" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "$$\n", + "\\begin{equation}\n", + " \\hat{y} = X\\omega,\n", + "\\label{_auto4} \\tag{4}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $X \\in \\mathbb{R}^{n \\times p}$ is the input data and $\\omega^{p\n", + "\\times d}$ are the weights of the regression. In a classification\n", + "setting (binary classification in our situation) we are interested in\n", + "a positive or negative answer. We can thus define either answer to be\n", + "above or below some threshold. But, in order to limit the size of the\n", + "answer and also to get a probability interpretation on how sure we are\n", + "for either answer we can compute the sigmoid function of OLS. That is," + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "$$\n", + "\\begin{equation}\n", + " f(X\\omega) = \\frac{1}{1 + \\exp(-X\\omega)}.\n", + "\\label{_auto5} \\tag{5}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We are thus interested in minizming the following cost function" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "$$\n", + "\\begin{equation}\n", + " C(X, \\omega) = \\sum_{i = 1}^n \\left\\{\n", + " - y_i\\log\\left( f(x_i^T\\omega) \\right)\n", + " - (1 - y_i)\\log\\left[1 - f(x_i^T\\omega)\\right]\n", + " \\right\\},\n", + "\\label{_auto6} \\tag{6}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where we will restrict ourselves to a value for $f(z)$ as the sigmoid\n", + "described above. We can also tack on a L2 (Ridge) or L1 (LASSO)\n", + "penalization to this cost function in the same manner we did for\n", + "linear regression.\n", + "\n", + "## Exploring the logistic regression\n", + "\n", + "The penalization factor $\\lambda$ is inverted in the case of the\n", + "logistic regression model we use. We will explore several values of\n", + "$\\lambda$ using both L1 and L2 penalization. We do this using a grid\n", + "search over different parameters and run a 3-fold cross validation for\n", + "each configuration. In other words, we fit a model 3 times for each\n", + "configuration of the hyper parameters." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "lambdas = np.logspace(-7, -1, 7)\n", + "\n", + "param_grid = {\n", + " \"C\": list(1.0/lambdas),\n", + " \"penalty\": [\"l1\", \"l2\"]\n", + "}\n", + "clf = skms.GridSearchCV(\n", + " skl.LogisticRegression(),\n", + " param_grid=param_grid,\n", + " n_jobs=-1,\n", + " return_train_score=True\n", + ")\n", + "t0 = time.time()\n", + "clf.fit(X_train, y_train)\n", + "t1 = time.time()\n", + "\n", + "print (\n", + " \"Time spent fitting GridSearchCV(LogisticRegression): {0:.3f} sec\".format(\n", + " t1 - t0\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that logistic regression is quite slow and using the grid\n", + "search and cross validation results in quite a heavy\n", + "computation. Below we show the results of the different\n", + "configurations." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "logreg_df = pd.DataFrame(clf.cv_results_)\n", + "\n", + "display(logreg_df)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Accuracy of a classification model\n", + "\n", + "To determine how well a classification model is performing we count\n", + "the number of correctly labeled classes and divide by the number of\n", + "classes in total. The accuracy is thus given by" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "$$\n", + "\\begin{equation}\n", + " a(y, \\hat{y}) = \\frac{1}{n}\\sum_{i = 1}^{n} I(y_i = \\hat{y}_i),\n", + "\\label{_auto7} \\tag{7}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $I(y_i = \\hat{y}_i)$ is the indicator function given by" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "$$\n", + "\\begin{equation}\n", + " I(x = y) = \\begin{cases}\n", + " 1 x = y, \n", + "\\label{_auto8} \\tag{8}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "$$\n", + "\\begin{equation} \n", + " 0 x \\neq y.\n", + " \\end{cases}\n", + "\\label{_auto9} \\tag{9}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is the accuracy provided by Scikit-learn when using **sklearn.metrics.accuracyscore**.\n", + "\n", + "Below we compute the accuracy of the best fit model on the training data (which should give a good accuracy), the test data (which has not been shown to the model) and the critical data (completely new data that needs to be extrapolated)." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "train_accuracy = skm.accuracy_score(y_train, clf.predict(X_train))\n", + "test_accuracy = skm.accuracy_score(y_test, clf.predict(X_test))\n", + "critical_accuracy = skm.accuracy_score(labels[critical], clf.predict(data[critical]))\n", + "\n", + "print (\"Accuracy on train data: {0}\".format(train_accuracy))\n", + "print (\"Accuracy on test data: {0}\".format(test_accuracy))\n", + "print (\"Accuracy on critical data: {0}\".format(critical_accuracy))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that we get quite good accuracy on the training data, but gradually worsening accuracy on the test and critical data.\n", + "\n", + "## Analyzing the results\n", + "\n", + "Below we show a different metric for determining the quality of our\n", + "model, namely the **reciever operating characteristic** (ROC). The ROC\n", + "curve tells us how well the model correctly classifies the different\n", + "labels. We plot the **true positive rate** (the rate of predicted\n", + "positive classes that are positive) versus the **false positive rate**\n", + "(the rate of predicted positive classes that are negative). The ROC\n", + "curve is built by computing the true positive rate and the false\n", + "positive rate for varying **thresholds**, i.e, which probability we\n", + "should acredit a certain class.\n", + "\n", + "By computing the **area under the curve** (AUC) of the ROC curve we get an estimate of how well our model is performing. Pure guessing will get an AUC of $0.5$. A perfect score will get an AUC of $1.0$." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(20, 14))\n", + "\n", + "for (_X, _y), label in zip(\n", + " [\n", + " (X_train, y_train),\n", + " (X_test, y_test),\n", + " (data[critical], labels[critical])\n", + " ],\n", + " [\"Train\", \"Test\", \"Critical\"]\n", + "):\n", + " proba = clf.predict_proba(_X)\n", + " fpr, tpr, _ = skm.roc_curve(_y, proba[:, 1])\n", + " roc_auc = skm.auc(fpr, tpr)\n", + "\n", + " print (\"LogisticRegression AUC ({0}): {1}\".format(label, roc_auc))\n", + "\n", + " plt.plot(fpr, tpr, label=\"{0} (AUC = {1})\".format(label, roc_auc), linewidth=4.0)\n", + "\n", + "plt.plot([0, 1], [0, 1], \"--\", label=\"Guessing (AUC = 0.5)\", linewidth=4.0)\n", + "\n", + "plt.title(r\"The ROC curve for LogisticRegression\", fontsize=18)\n", + "plt.xlabel(r\"False positive rate\", fontsize=18)\n", + "plt.ylabel(r\"True positive rate\", fontsize=18)\n", + "plt.axis([-0.01, 1.01, -0.01, 1.01])\n", + "plt.xticks(fontsize=18)\n", + "plt.yticks(fontsize=18)\n", + "plt.legend(loc=\"best\", fontsize=18)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that this plot of the ROC looks very strange. This tells us\n", + "that logistic regression is quite inept at predicting the Ising model\n", + "transition and is therefore highly non-linear. The ROC curve for the\n", + "training data looks quite good, but as the testing data is so far off\n", + "we see that we are dealing with an overfit model.\n", + "\n", + "A previous run with $50\\%$ of the data used for training yielded a\n", + "worse performance than using a smaller training set. This again gives\n", + "confidence to the fact that logistic regression is not able to\n", + "correctly fit the Ising model as it is not a linear model." + ] } ], "metadata": {}, diff --git a/doc/pub/LogReg/ipynb/ipynb-LogReg-src.tar.gz b/doc/pub/LogReg/ipynb/ipynb-LogReg-src.tar.gz index bd49f72aa..81585084e 100644 Binary files a/doc/pub/LogReg/ipynb/ipynb-LogReg-src.tar.gz and b/doc/pub/LogReg/ipynb/ipynb-LogReg-src.tar.gz differ diff --git a/doc/pub/LogReg/pdf/LogReg-beamer-handouts2x3.pdf b/doc/pub/LogReg/pdf/LogReg-beamer-handouts2x3.pdf index e7e36e3de..2aa78af95 100644 Binary files a/doc/pub/LogReg/pdf/LogReg-beamer-handouts2x3.pdf and b/doc/pub/LogReg/pdf/LogReg-beamer-handouts2x3.pdf differ diff --git a/doc/pub/LogReg/pdf/LogReg-beamer.pdf b/doc/pub/LogReg/pdf/LogReg-beamer.pdf index ce0339144..d5c376136 100644 Binary files a/doc/pub/LogReg/pdf/LogReg-beamer.pdf and b/doc/pub/LogReg/pdf/LogReg-beamer.pdf differ diff --git a/doc/pub/LogReg/pdf/LogReg-minted.pdf b/doc/pub/LogReg/pdf/LogReg-minted.pdf index 66f878004..2ffffaae2 100644 Binary files a/doc/pub/LogReg/pdf/LogReg-minted.pdf and b/doc/pub/LogReg/pdf/LogReg-minted.pdf differ diff --git a/doc/src/LogisticRegression/LogReg.do.txt b/doc/src/LogisticRegression/LogReg.do.txt index cb5c7ef3f..a22273019 100644 --- a/doc/src/LogisticRegression/LogReg.do.txt +++ b/doc/src/LogisticRegression/LogReg.do.txt @@ -388,3 +388,292 @@ def main(): if __name__ == "__main__": main() !ec + +!split +===== The two-dimensional Ising model, Predicting phase transition of the two-dimensional Ising model ===== + +The Hamiltonian of the two-dimensional Ising model without an external field for a constant coupling constant $J$ is given by +!bt +\begin{align} + H = -J \sum_{\langle ij\rangle} S_i S_j, +\end{align} +!et +where $S_i \in \{-1, 1\}$ and $\langle ij \rangle$ signifies that we only iterate over the nearest neighbors in the lattice. We will be looking at a system of $L = 40$ spins in each dimension, i.e., $L^2 = 1600$ spins in total. Opposed to the one-dimensional Ising model we will get a phase transition from an _ordered_ phase to a _disordered_ phase at the critical temperature + +!bt +\begin{align} + \frac{T_c}{J} = \frac{2}{\log\left(1 + \sqrt{2}\right)} \approx 2.26, +\end{align} +!et +as shown by Lars Onsager. + + +Here we use _logistic regression_ to predict when a phase transition +occurs. The data we will look at is a set of spin configurations, +i.e., individual lattices with spins, labeled _ordered_ `1` or +_disordered_ `0`. Our job is to build a model which will take in a +spin configuration and predict whether or not the spin configuration +constitutes an ordered or a disordered phase. To achieve this we will +represent the lattices as flattened arrays with $1600$ elements +instead of a matrix of $40 \times 40$ elements. As an extra test of +the performance of the algorithms we will divide the dataset into +three pieces. We will do a conventional train-test-split on a +combination of totally ordered and totally disordered phases. The +remaining "critical-like" states will be used as test data which we +hope the model will be able to make good extrapolated predictions on. + + +!bc pycod +import pickle +import os +import glob +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import seaborn as sns +import sklearn.model_selection as skms +import sklearn.linear_model as skl +import sklearn.metrics as skm +import tqdm +import copy +import time +from IPython.display import display + +%matplotlib inline + +sns.set(color_codes=True) +!ec + +!split +===== Reading in the data ===== + +Using the data from "Mehta et al.":"https://physics.bu.edu/~pankajm/ML-Review-Datasets/isingMC/" (specifically the two datasets named `Ising2DFM_reSample_L40_T=All.pkl` and `Ising2DFM_reSample_L40_T=All_labels.pkl`) we have to unpack the data into numpy arrays. + + +!bc pycod +filenames = glob.glob(os.path.join("..", "dat", "*")) +label_filename = list(filter(lambda x: "label" in x, filenames))[0] +dat_filename = list(filter(lambda x: "label" not in x, filenames))[0] + +# Read in the labels +with open(label_filename, "rb") as f: + labels = pickle.load(f) + +# Read in the corresponding configurations +with open(dat_filename, "rb") as f: + data = np.unpackbits(pickle.load(f)).reshape(-1, 1600).astype("int") + +# Set spin-down to -1 +data[data == 0] = -1 +!ec + +This dataset consists of $10000$ samples, i.e., $10000$ spin +configurations with $40 \times 40$ spins each, for $16$ temperatures +between $0.25$ to $4.0$. Next we create a train/test-split and keep +the data in the critical phase as a separate dataset for +extrapolation-testing. + + +!bc pycod +# Set up slices of the dataset +ordered = slice(0, 70000) +critical = slice(70000, 100000) +disordered = slice(100000, 160000) + +X_train, X_test, y_train, y_test = skms.train_test_split( + np.concatenate((data[ordered], data[disordered])), + np.concatenate((labels[ordered], labels[disordered])), + test_size=0.95 +) +!ec + +Using a small training set yields a better accuracy. This will be discussed in the end. + +!split +===== Logistic regression ===== + +Logistic regression is a linear model for classification. Recalling +the cost function for ordinary least squares with both L2 (ridge) and +L1 (LASSO) penalties we will see that the logistic cost function is +very similar. In OLS we wish to predict a continuous variable +$\hat{y}$ using +!bt +\begin{align} + \hat{y} = X\omega, +\end{align} +!et + +where $X \in \mathbb{R}^{n \times p}$ is the input data and $\omega^{p +\times d}$ are the weights of the regression. In a classification +setting (binary classification in our situation) we are interested in +a positive or negative answer. We can thus define either answer to be +above or below some threshold. But, in order to limit the size of the +answer and also to get a probability interpretation on how sure we are +for either answer we can compute the sigmoid function of OLS. That is, + +!bt +\begin{align} + f(X\omega) = \frac{1}{1 + \exp(-X\omega)}. +\end{align} +!et +We are thus interested in minizming the following cost function +!bt +\begin{align} + C(X, \omega) = \sum_{i = 1}^n \left\{ + - y_i\log\left( f(x_i^T\omega) \right) + - (1 - y_i)\log\left[1 - f(x_i^T\omega)\right] + \right\}, +\end{align} +!et + +where we will restrict ourselves to a value for $f(z)$ as the sigmoid +described above. We can also tack on a L2 (Ridge) or L1 (LASSO) +penalization to this cost function in the same manner we did for +linear regression. + +!split +===== Exploring the logistic regression ===== + +The penalization factor $\lambda$ is inverted in the case of the +logistic regression model we use. We will explore several values of +$\lambda$ using both L1 and L2 penalization. We do this using a grid +search over different parameters and run a 3-fold cross validation for +each configuration. In other words, we fit a model 3 times for each +configuration of the hyper parameters. + + +!bc pycod +lambdas = np.logspace(-7, -1, 7) + +param_grid = { + "C": list(1.0/lambdas), + "penalty": ["l1", "l2"] +} +clf = skms.GridSearchCV( + skl.LogisticRegression(), + param_grid=param_grid, + n_jobs=-1, + return_train_score=True +) +t0 = time.time() +clf.fit(X_train, y_train) +t1 = time.time() + +print ( + "Time spent fitting GridSearchCV(LogisticRegression): {0:.3f} sec".format( + t1 - t0 + ) +) +!ec + +We can see that logistic regression is quite slow and using the grid +search and cross validation results in quite a heavy +computation. Below we show the results of the different +configurations. + + +!bc pycod +logreg_df = pd.DataFrame(clf.cv_results_) + +display(logreg_df) +!ec + +!split +===== Accuracy of a classification model ===== + +To determine how well a classification model is performing we count +the number of correctly labeled classes and divide by the number of +classes in total. The accuracy is thus given by + +!bt +\begin{align} + a(y, \hat{y}) = \frac{1}{n}\sum_{i = 1}^{n} I(y_i = \hat{y}_i), +\end{align} +!et + +where $I(y_i = \hat{y}_i)$ is the indicator function given by + +!bt +\begin{align} + I(x = y) = \begin{cases} + 1 & x = y, \\ + 0 & x \neq y. + \end{cases} +\end{align} +!et + +This is the accuracy provided by Scikit-learn when using _sklearn.metrics.accuracyscore_. + +Below we compute the accuracy of the best fit model on the training data (which should give a good accuracy), the test data (which has not been shown to the model) and the critical data (completely new data that needs to be extrapolated). + + +!bc pycod +train_accuracy = skm.accuracy_score(y_train, clf.predict(X_train)) +test_accuracy = skm.accuracy_score(y_test, clf.predict(X_test)) +critical_accuracy = skm.accuracy_score(labels[critical], clf.predict(data[critical])) + +print ("Accuracy on train data: {0}".format(train_accuracy)) +print ("Accuracy on test data: {0}".format(test_accuracy)) +print ("Accuracy on critical data: {0}".format(critical_accuracy)) +!ec + +We can see that we get quite good accuracy on the training data, but gradually worsening accuracy on the test and critical data. + +!split +===== Analyzing the results ===== + +Below we show a different metric for determining the quality of our +model, namely the _reciever operating characteristic_ (ROC). The ROC +curve tells us how well the model correctly classifies the different +labels. We plot the _true positive rate_ (the rate of predicted +positive classes that are positive) versus the _false positive rate_ +(the rate of predicted positive classes that are negative). The ROC +curve is built by computing the true positive rate and the false +positive rate for varying _thresholds_, i.e, which probability we +should acredit a certain class. + +By computing the _area under the curve_ (AUC) of the ROC curve we get an estimate of how well our model is performing. Pure guessing will get an AUC of $0.5$. A perfect score will get an AUC of $1.0$. + + +!bc pycod +fig = plt.figure(figsize=(20, 14)) + +for (_X, _y), label in zip( + [ + (X_train, y_train), + (X_test, y_test), + (data[critical], labels[critical]) + ], + ["Train", "Test", "Critical"] +): + proba = clf.predict_proba(_X) + fpr, tpr, _ = skm.roc_curve(_y, proba[:, 1]) + roc_auc = skm.auc(fpr, tpr) + + print ("LogisticRegression AUC ({0}): {1}".format(label, roc_auc)) + + plt.plot(fpr, tpr, label="{0} (AUC = {1})".format(label, roc_auc), linewidth=4.0) + +plt.plot([0, 1], [0, 1], "--", label="Guessing (AUC = 0.5)", linewidth=4.0) + +plt.title(r"The ROC curve for LogisticRegression", fontsize=18) +plt.xlabel(r"False positive rate", fontsize=18) +plt.ylabel(r"True positive rate", fontsize=18) +plt.axis([-0.01, 1.01, -0.01, 1.01]) +plt.xticks(fontsize=18) +plt.yticks(fontsize=18) +plt.legend(loc="best", fontsize=18) +plt.show() +!ec + +We can see that this plot of the ROC looks very strange. This tells us +that logistic regression is quite inept at predicting the Ising model +transition and is therefore highly non-linear. The ROC curve for the +training data looks quite good, but as the testing data is so far off +we see that we are dealing with an overfit model. + +A previous run with $50\%$ of the data used for training yielded a +worse performance than using a smaller training set. This again gives +confidence to the fact that logistic regression is not able to +correctly fit the Ising model as it is not a linear model. +