深度学习算法实践6---逻辑回归算法应用

在上篇博文中，我们介绍了深度学习算法的实现，并且以MNIST手写数字识别为例，验证了该算法的有效性。

但是我们学习逻辑回归算法的目的是解决我们的实际问题，而不是学习算法本身。逻辑回归算法在实际中的应用还是很广泛的，例如在医学领域的疾病预测中，我们就可以列出一系疾病相关因素，然后根据某位患者的具体情况，应用逻辑回归算法，判断该患者是否患有某种疾病。当然，逻辑回归算法还是有局限性的，其比较适合于处理线性可分的分类问题，但是对于线性不可分的分类问题，这种算法的价值就会大打折扣了。但是我们可以将逻辑回归算法，视为没有隐藏层的前馈网络，通过增加隐藏层，就可以处理各种线性不可分问题了。借助于Theano的框架，在后面博文中我们会介绍BP网络、多层卷积网络（LeNet），大家可以看到，在Theano中，实现这些模型是一件非常简单的事情。

言归正传，如果我们要用逻辑回归算法解决实际问题，我们主要需要改变的就是load_data函数，使其从我们规定的数据源中读取数据。在此，我们先设计一个训练数据读入的工具类SegLoader，文件名为seg_loader.py，代码如下所示：

[python] view plain copy

 

1. from __future__ import print_function  
2.   
3. __docformat__ = 'restructedtext en'  
4.   
5. import six.moves.cPickle as pickle  
6. import gzip  
7. import os  
8. import sys  
9. import timeit  
10.   
11. import numpy  
12.   
13. import theano  
14. import theano.tensor as T  
15.   
16. class SegLoader(object):  
17.     def load_data(self, dataset):  
18.         samplesNumber = 6  
19.         features = 2  
20.         train_set = (numpy.ndarray(shape=(samplesNumber, features), dtype=numpy.float32), numpy.ndarray(shape=(samplesNumber), dtype=int))  
21.         self.prepare_dataset(train_set)  
22.         valid_set = (train_set[0].copy(), train_set[1].copy())  
23.         test_set = (train_set[0].copy(), train_set[1].copy())  
24.         test_set_x, test_set_y = self.shared_dataset(test_set)  
25.         valid_set_x, valid_set_y = self.shared_dataset(valid_set)  
26.         train_set_x, train_set_y = self.shared_dataset(train_set)  
27.         rval = [(train_set_x, train_set_y), (valid_set_x, valid_set_y),  
28.             (test_set_x, test_set_y)]  
29.         return rval  
30.   
31.     def shared_dataset(self, data_xy, borrow=True):  
32.         data_x, data_y = data_xy  
33.         shared_x = theano.shared(numpy.asarray(data_x,  
34.                                                dtype=theano.config.floatX),  
35.                                  borrow=borrow)  
36.         shared_y = theano.shared(numpy.asarray(data_y,  
37.                                                dtype=theano.config.floatX),  
38.                                  borrow=borrow)  
39.         return shared_x, T.cast(shared_y, 'int32')  
40.   
41.     def prepare_dataset(self, dataset):  
42.         dataset[0][0][0] = 1.0  
43.         dataset[0][0][1] = 1.0  
44.         dataset[1][0] = 1  
45.   
46.         dataset[0][1][0] = 2.0  
47.         dataset[0][1][1] = 2.0  
48.         dataset[1][1] = 1  
49.   
50.         dataset[0][2][0] = 3.0  
51.         dataset[0][2][1] = 3.0  
52.         dataset[1][2] = 1  
53.   
54.         dataset[0][3][0] = 1.5  
55.         dataset[0][3][1] = 2.0  
56.         dataset[1][3] = 0  
57.   
58.         dataset[0][4][0] = 2.5  
59.         dataset[0][4][1] = 4.0  
60.         dataset[1][4] = 0  
61.   
62.         dataset[0][5][0] = 3.5  
63.         dataset[0][5][1] = 7.0  
64.         dataset[1][5] = 0  

上面的代码非常简单，生成一个元组train_set，包含两个元素，第一个元素是一个类型为float32的二维数组，每行代表一个样本，第一列代表X坐标，第二列代表Y坐标，train_set元组的第二个元素为一维整数数组，每个元素代表一个样本的分类结果，这里有两个大类，1代表在Y=X的直线上，0代表不在该直线上，prepare_dataset准备了6个训练样。因为这个问题非常简单，所以6个样本基本就够用了，但是对实际问题而言，显然需要相当大的样本量。

接着我们定义这个线性分割的执行引擎LrSegEngine，源码文件为lr_seg_engine.py，代码如下所示：

[python] view plain copy

 

1. from __future__ import print_function  
2.   
3. __docformat__ = 'restructedtext en'  
4.   
5. import six.moves.cPickle as pickle  
6. import gzip  
7. import os  
8. import sys  
9. import timeit  
10.   
11. import numpy  
12.   
13. import theano  
14. import theano.tensor as T  
15. from logistic_regression import LogisticRegression  
16. from seg_loader import SegLoader  
17.   
18. class LrSegEngine(object):  
19.     def __init__(self):  
20.         print("Logistic Regression MNIST Engine")  
21.         self.learning_rate = 0.13  
22.         self.n_epochs = 1000  
23.         self.batch_size = 1 # 600  
24.         self.dataset = 'mnist.pkl.gz'  
25.   
26.     def train(self):  
27.         print("Yantao:train the model")  
28.         loader = SegLoader()  
29.         datasets = loader.load_data(self.dataset)  
30.         train_set_x, train_set_y = datasets[0]  
31.         valid_set_x, valid_set_y = datasets[1]  
32.         test_set_x, test_set_y = datasets[2]  
33.         n_train_batches = train_set_x.get_value(borrow=True).shape[0] // self.batch_size  
34.         n_valid_batches = valid_set_x.get_value(borrow=True).shape[0] // self.batch_size  
35.         n_test_batches = test_set_x.get_value(borrow=True).shape[0] // self.batch_size  
36.         index = T.lscalar()  
37.         x = T.matrix('x')  
38.         y = T.ivector('y')  
39.         # in:x,y out: 1 in y=x otherwise 0  
40.         classifier = LogisticRegression(input=x, n_in=2, n_out=2)  
41.         cost = classifier.negative_log_likelihood(y)  
42.         test_model = theano.function(  
43.             inputs=[index],  
44.             outputs=classifier.errors(y),  
45.             givens={  
46.                 x: test_set_x[index * self.batch_size: (index + 1) * self.batch_size],  
47.                 y: test_set_y[index * self.batch_size: (index + 1) * self.batch_size]  
48.             }  
49.         )  
50.         validate_model = theano.function(  
51.             inputs=[index],  
52.             outputs=classifier.errors(y),  
53.             givens={  
54.                 x: valid_set_x[index * self.batch_size: (index + 1) * self.batch_size],  
55.                 y: valid_set_y[index * self.batch_size: (index + 1) * self.batch_size]  
56.             }  
57.         )  
58.         g_W = T.grad(cost=cost, wrt=classifier.W)  
59.         g_b = T.grad(cost=cost, wrt=classifier.b)  
60.         updates = [(classifier.W, classifier.W - self.learning_rate * g_W),  
61.                (classifier.b, classifier.b - self.learning_rate * g_b)]  
62.         train_model = theano.function(  
63.             inputs=[index],  
64.             outputs=cost,  
65.             updates=updates,  
66.             givens={  
67.                 x: train_set_x[index * self.batch_size: (index + 1) * self.batch_size],  
68.                 y: train_set_y[index * self.batch_size: (index + 1) * self.batch_size]  
69.             }  
70.         )  
71.         patience = 5000    
72.         patience_increase = 2    
73.         improvement_threshold = 0.995    
74.         validation_frequency = min(n_train_batches, patience // 2)  
75.         best_validation_loss = numpy.inf  
76.         test_score = 0.  
77.         start_time = timeit.default_timer()  
78.         done_looping = False  
79.         epoch = 0  
80.         while (epoch < self.n_epochs) and (not done_looping):  
81.             epoch = epoch + 1  
82.             for minibatch_index in range(n_train_batches):  
83.                 minibatch_avg_cost = train_model(minibatch_index)  
84.                 # iteration number  
85.                 iter = (epoch - 1) * n_train_batches + minibatch_index  
86.                 if (iter + 1) % validation_frequency == 0:  
87.                     # compute zero-one loss on validation set  
88.                     validation_losses = [validate_model(i)  
89.                                      for i in range(n_valid_batches)]  
90.                     this_validation_loss = numpy.mean(validation_losses)  
91.                     print(  
92.                         'epoch %i, minibatch %i/%i, validation error %f %%' %  
93.                         (  
94.                             epoch,  
95.                             minibatch_index + 1,  
96.                             n_train_batches,  
97.                             this_validation_loss * 100.  
98.                         )  
99.                     )  
100.                     if this_validation_loss < best_validation_loss:  
101.                         #improve patience if loss improvement is good enough  
102.                         if this_validation_loss < best_validation_loss * improvement_threshold:  
103.                             patience = max(patience, iter * patience_increase)  
104.                         best_validation_loss = this_validation_loss  
105.                         # test it on the test set  
106.                         test_losses = [test_model(i)  
107.                                    for i in range(n_test_batches)]  
108.                         test_score = numpy.mean(test_losses)  
109.                         print(  
110.                             (  
111.                                 '     epoch %i, minibatch %i/%i, test error of'  
112.                                 ' best model %f %%'  
113.                             ) %  
114.                             (  
115.                                 epoch,  
116.                                 minibatch_index + 1,  
117.                                 n_train_batches,  
118.                                 test_score * 100.  
119.                             )  
120.                         )  
121.                         # save the best model  
122.                         with open('best_model.pkl', 'wb') as f:  
123.                             pickle.dump(classifier, f)  
124.                 if patience <= iter:  
125.                     done_looping = True  
126.                     break  
127.         end_time = timeit.default_timer()  
128.         print(  
129.             (  
130.                 'Optimization complete with best validation score of %f %%,'  
131.                 'with test performance %f %%'  
132.             )  
133.             % (best_validation_loss * 100., test_score * 100.)  
134.         )  
135.         print('The code run for %d epochs, with %f epochs/sec' % (  
136.             epoch, 1. * epoch / (end_time - start_time)))  
137.         print(('The code for file ' +  
138.                os.path.split(__file__)[1] +  
139.                ' ran for %.1fs' % ((end_time - start_time))), file=sys.stderr)  
140.   
141.     def run(self, data):  
142.         print("run the model")  
143.         classifier = pickle.load(open('best_model.pkl', 'rb'))  
144.         predict_model = theano.function(  
145.             inputs=[classifier.input],  
146.             outputs=classifier.y_pred  
147.         )  
148.         rst = predict_model(data)  
149.         print(rst)  

在这里的train方法，与上篇博文处理MNIST手写数字识别的代码基本一致，只需要注意以下几点：首先，由于我们只有6个样本，因此将样本批次的大小设置为1（在MNIST手写数字识别中，由于有6万个训练样本，所以批次大小为600）；其次，在初始化逻辑回归模型时，输入维度n_in，设置为2，表示样本只有两个特征即x,y坐标，输出维度也为2，表示有两个类别，1是在y=x线上，0代表不在线上。

接着我们定义逻辑回归模型类LogisticRegression，源码文件为logistic_regression.py，代码如下所示：

[python] view plain copy

 

1. from __future__ import print_function  
2.   
3. __docformat__ = 'restructedtext en'  
4.   
5. import six.moves.cPickle as pickle  
6. import gzip  
7. import os  
8. import sys  
9. import timeit  
10.   
11. import numpy  
12.   
13. import theano  
14. import theano.tensor as T  
15.   
16. class LogisticRegression(object):    
17.     def __init__(self, input, n_in, n_out):    
18.         self.W = theano.shared(    
19.             value=numpy.zeros(    
20.                 (n_in, n_out),    
21.                 dtype=theano.config.floatX    
22.             ),    
23.             name='W',    
24.             borrow=True    
25.         )    
26.         self.b = theano.shared(    
27.             value=numpy.zeros(    
28.                 (n_out,),    
29.                 dtype=theano.config.floatX    
30.             ),    
31.             name='b',    
32.             borrow=True    
33.         )    
34.         self.p_y_given_x = T.nnet.softmax(T.dot(input, self.W) + self.b)    
35.         self.y_pred = T.argmax(self.p_y_given_x, axis=1)    
36.         self.params = [self.W, self.b]    
37.         self.input = input    
38.         print("Yantao: ***********************************")  
39.     
40.     def negative_log_likelihood(self, y):    
41.         return -T.mean(T.log(self.p_y_given_x)[T.arange(y.shape[0]), y])    
42.     
43.     def errors(self, y):    
44.         if y.ndim != self.y_pred.ndim:    
45.             raise TypeError(    
46.                 'y should have the same shape as self.y_pred',    
47.                 ('y', y.type, 'y_pred', self.y_pred.type)    
48.             )    
49.         if y.dtype.startswith('int'):    
50.             return T.mean(T.neq(self.y_pred, y))    
51.         else:    
52.             raise NotImplementedError()    

上面的代码与上篇博文几乎没有变化，只是将其单独保存到一个文件中而已。

接下来是模型训练lr_train.py，代码如下所示：

[python] view plain copy

 

1. from __future__ import print_function  
2.   
3. __docformat__ = 'restructedtext en'  
4.   
5. import six.moves.cPickle as pickle  
6. import gzip  
7. import os  
8. import sys  
9. import timeit  
10.   
11. import numpy  
12.   
13. import theano  
14. import theano.tensor as T  
15.   
16. from logistic_regression import LogisticRegression  
17. from seg_loader import SegLoader  
18. from lr_seg_engine import LrSegEngine  
19.   
20. if __name__ == '__main__':  
21.     engine = LrSegEngine()  
22.     engine.train()  

上面代码只是简单调用逻辑回归分割的引擎类的train方法，完成对模型的训练，其会将最佳的结果保存到best_model.pkl文件中。

当模型训练好之后，我们就可以拿模型来进行分类了，lr_run.py的代码如下所示：

[python] view plain copy

 

1. from seg_loader import SegLoader  
2. from lr_seg_engine import LrSegEngine  
3.   
4. if __name__ == '__main__':  
5.     print("test program v1.0")  
6.     engine = LrSegEngine()  
7.     data = [[2.0, 2.0]]  
8.     print(data)  
9.     engine.run(data)  

上面代码首先初始化一个二维数组，其中只有一个样本元素，坐标为(2.0, 2.0)，然后调用逻辑回归分割引擎的run方法，其将给出分类结果，运行这个程序，会得到类似如下所示的结果：

test program v1.0

Logistic Regression MNIST Engine

[[2.0, 2.0]]

run the model

[1]