EdmondFrank's 时光足迹

この先は暗い夜道だけかもしれない それでも信じて進むんだ。星がその道を少しでも照らしてくれるのを。
或许前路永夜,即便如此我也要前进,因为星光即使微弱也会我为照亮前途。
——《四月は君の嘘》

贝叶斯景象图



贝叶斯景象图

理论说明

均匀分布

对于一个含有N个未知元素的贝叶斯推断问题,我们隐式地为其先验分布创建了一个N维空间。先验分布上某一点的概率,都投射到某个高维的面或曲线上,其形状由先验分布决定。比如,假定有两个未知元素 ,其先验分布都是(0,5)上的均匀分布,那么先验分布存在于一个边长为5的正方形空间,而其概率面就是正方形上方的一个平面(由于假定了均匀分布,因此每一点概率相同)。

代码绘图演示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
%matplotlib inline
import scipy.stats as stats
from IPython.core.pylabtools import figsize
import numpy as np
figsize(12.5, 4)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

jet = plt.cm.jet
fig = plt.figure()
x = y = np.linspace(0, 5, 100)
X, Y = np.meshgrid(x, y)

plt.subplot(121)
uni_x = stats.uniform.pdf(x, loc=0, scale=5)
uni_y = stats.uniform.pdf(y, loc=0, scale=5)
M = np.dot(uni_x[:, None], uni_y[None, :])
im = plt.imshow(M, interpolation='none', origin='lower',
                cmap=jet, vmax=1, vmin=-.15, extent=(0, 5, 0, 5))

plt.xlim(0, 5)
plt.ylim(0, 5)
plt.title("Landscape formed by Uniform priors.")

ax = fig.add_subplot(122, projection='3d')
ax.plot_surface(X, Y, M, cmap=plt.cm.jet, vmax=1, vmin=-.15)
ax.view_init(azim=390)
plt.title("Uniform prior landscape; alternate view");

指数分布

再者,如果 的先验分布为Exp(3)和Exp(10),那么对应的空间便是二维平面上,各维都取正值确定的范围,而对应的概率面的形状就是一个从(0,0)点向正值方向流淌的瀑布。

以下的示例图就描绘了这样的情形,其中颜色越是趋向于暗红的位置,其先验概率就越高。反过来,颜色越是趋向于深蓝的位置,其先验概率就越低。

代码绘图演示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
figsize(12.5, 5)
fig = plt.figure()
plt.subplot(121)

exp_x = stats.expon.pdf(x, scale=3)
exp_y = stats.expon.pdf(x, scale=10)
M = np.dot(exp_x[:, None], exp_y[None, :])
CS = plt.contour(X, Y, M)
im = plt.imshow(M, interpolation='none', origin='lower',
                cmap=jet, extent=(0, 5, 0, 5))
#plt.xlabel("prior on $p_1$")
#plt.ylabel("prior on $p_2$")
plt.title("$Exp(3), Exp(10)$ prior landscape")

ax = fig.add_subplot(122, projection='3d')
ax.plot_surface(X, Y, M, cmap=jet)
ax.view_init(azim=390)
plt.title("$Exp(3), Exp(10)$ prior landscape; \nalternate view");

这些二维空间的例子很简单,我们的大脑可以轻易想象得到。但实际中,先验分布所在的空间和其概率面往往具有更高的维度。

观测值对先验分布的影响 

在实际中,观测样本对空间不会有影响,但它会改变概率面的形状,将其在某些局部区域拉伸或挤压,以表明参数的真实性在哪里。更多的数据意味着对概率面更多的拉伸和挤压,使得最初的概率面形状变得十分奇怪。反之数据越少,那么最初的形状就保留得越好。不管如何,最后得到的概率面描述了后验分布的形状。

假如我们现在想对两个参数为的泊松分布进行估计。那么我们将要分别比较用均匀分布和指数分布来对的先验分布进行假设的不同效果。

代码绘图演示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
# create the observed data

# sample size of data we observe, trying varying this (keep it less than 100 ;)
N = 1

# the true parameters, but of course we do not see these values...
lambda_1_true = 1
lambda_2_true = 3

#...we see the data generated, dependent on the above two values.
data = np.concatenate([
    stats.poisson.rvs(lambda_1_true, size=(N, 1)),
    stats.poisson.rvs(lambda_2_true, size=(N, 1))
], axis=1)
print("observed (2-dimensional,sample size = %d):" % N, data)

# plotting details.
x = y = np.linspace(.01, 5, 100)
likelihood_x = np.array([stats.poisson.pmf(data[:, 0], _x)
                        for _x in x]).prod(axis=1)
likelihood_y = np.array([stats.poisson.pmf(data[:, 1], _y)
                        for _y in y]).prod(axis=1)
L = np.dot(likelihood_x[:, None], likelihood_y[None, :])
observed (2-dimensional,sample size = 1): [[0 2]]
In [4]:
figsize(12.5, 12)
# matplotlib heavy lifting below, beware!
plt.subplot(221)
uni_x = stats.uniform.pdf(x, loc=0, scale=5)
uni_y = stats.uniform.pdf(x, loc=0, scale=5)
M = np.dot(uni_x[:, None], uni_y[None, :])
im = plt.imshow(M, interpolation='none', origin='lower',
                cmap=jet, vmax=1, vmin=-.15, extent=(0, 5, 0, 5))
plt.scatter(lambda_2_true, lambda_1_true, c="k", s=50, edgecolor="none")
plt.xlim(0, 5)
plt.ylim(0, 5)
plt.title("Landscape formed by Uniform priors on $p_1, p_2$.")

plt.subplot(223)
plt.contour(x, y, M * L)
im = plt.imshow(M * L, interpolation='none', origin='lower',
                cmap=jet, extent=(0, 5, 0, 5))
plt.title("Landscape warped by %d data observation;\n Uniform priors on $p_1, p_2$." % N)
plt.scatter(lambda_2_true, lambda_1_true, c="k", s=50, edgecolor="none")
plt.xlim(0, 5)
plt.ylim(0, 5)

plt.subplot(222)
exp_x = stats.expon.pdf(x, loc=0, scale=3)
exp_y = stats.expon.pdf(x, loc=0, scale=10)
M = np.dot(exp_x[:, None], exp_y[None, :])

plt.contour(x, y, M)
im = plt.imshow(M, interpolation='none', origin='lower',
                cmap=jet, extent=(0, 5, 0, 5))
plt.scatter(lambda_2_true, lambda_1_true, c="k", s=50, edgecolor="none")
plt.xlim(0, 5)
plt.ylim(0, 5)
plt.title("Landscape formed by Exponential priors on $p_1, p_2$.")

plt.subplot(224)
# This is the likelihood times prior, that results in the posterior.
plt.contour(x, y, M * L)
im = plt.imshow(M * L, interpolation='none', origin='lower',
                cmap=jet, extent=(0, 5, 0, 5))

plt.scatter(lambda_2_true, lambda_1_true, c="k", s=50, edgecolor="none")
plt.title("Landscape warped by %d data observation;\n Exponential priors on \
$p_1, p_2$." % N)
plt.xlim(0, 5)
plt.ylim(0, 5);

四张图里的黑点代表了参数的真实取值,左下图为均匀先验得到的后验分布图。虽然观测值相同,但是两种假设下的后验分布形状是不一样的。其主要原因是因为观测点的位置在两者的假设的前提先验概率是不一样的。这样,我们可以知道,即便只有一个观测值,形成的山峰也试图要包括参数值的真实位置。当然,在真正的推断中,仅用一个观测值显然也是十分不科学的,这里仅仅为了方便阐述而已。

本文参考自《Probabilistic-Programming-and-Bayesian-Methods-for-Hackers》

Python Keras + LSTM 进行单变量时间序列预测

Python Keras + LSTM 进行单变量时间序列预测

首先,时间序列预测问题是一个复杂的预测模型问题,它不像一般的回归预测模型。时间序列预测的输入变量是一组按时间顺序的数字序列。它既具有延续性又具有随机性,所以在建模难度上相对回归预测更大。

但同时,正好有一种强大的神经网络适合处理这种存在依赖关系的序列问题:RNN(Recurrent neural networks)。在过去几年中,应用 RNN 在语音识别,语言建模,翻译,图片描述等问题上已经取得一定成功,并且应用领域还在扩展。

LSTM网络

Long Short-Term Memory 网络亦称LSTM 网络,是一种在深度学习中应用的循环神经网络。可以学习长期依赖信息。LSTM 由Hochreiter & Schmidhuber (1997)提出,并在近期被Alex Graves进行了改良和推广。在很多问题,LSTM 都取得相当巨大的成功,并得到了广泛的使用。LSTM 通过刻意的设计来避免长期依赖问题。记住长期的信息在实践中是 LSTM 的默认行为,而非需要付出很大代价才能获得的能力。

具体应用

下面以一个洗发水销售的例子,来实现LSTM。 首先,你可以在这里下载到本文需要用的数据集。这是一个描述了3年内洗发水的月度销售数量的数据集。

数据读取

1
2
3
4
5
6
7
8
9
10
11
from pandas import read_csv
from pandas import datetime
from matplotlib import pyplot as plt
def parser(x):
    return datetime.strptime('190'+x,"%Y-%m")
series = read_csv('sales-of-shampoo-over-a-three-ye.csv',
                 header=0,parse_dates=[0],index_col=0,
                  squeeze=True, date_parser=parser)
print(series.head())
series.plot()
plt.show()

数据划分

首先我们把数据集划分成两个部分即:训练集和测试集。 那么我们该如何划分呢?因为我们今天研究的是时间序列分析,所以在数据集的划分上我们也应该按照时间来划分。我们可以将前两年的数据作为我们的训练集而将最后一年的数据作为测试集。

1
2
3
# split data into train and test
X = series.values
train, test = X[0:-12], X[-12:]

这里我们假设一个滚动预测的情景,又称前向模型验证(walk-forward model validation)。其原理很简单,举例来说就像当公司的预测期长达一年时,预测会将已过去的月份排除,而将预测期末的月份补上。好比一月份过去后,我们将其从预测中移除,同时次年的一月份就会作为收尾被添加到预测中以便预测总能保持12个月的完整性。

这样通过使用每月新的洗发水销售量来进行下个月的预测,我们就像模拟了一个更接近于真实世界的场景。

最后,我们将所有在测试集上的预测结果收集起来并计算出他们与真实值的均方根误差(RMSE)以此来作为评估我们模型的基准。

持续模型预测(Persistence Model Forecast)

持续性预测的基本思路就是从先前的(t-1)时间序列的结果用于预测当前时间(t)的取值。 那么根据以上的思路,我们可以通过滚动预测的原理从训练集的历史数据中获取最后一次观察值并使用它来预测当前时间的可能取值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
from pandas import read_csv
from pandas import datetime
from sklearn.metrics import mean_squared_error
from math import sqrt
from matplotlib import pyplot
# load dataset
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('sales-of-shampoo-over-a-three-ye.csv',
                  header=0, parse_dates=[0], index_col=0,
                  squeeze=True, date_parser=parser)
# split data into train and test
X = series.values
train, test = X[0:-12], X[-12:]
# walk-forward validation
history = [x for x in train]
predictions = list()
for i in range(len(test)):
  # make prediction
  predictions.append(history[-1])
  # observation
  history.append(test[i])
# report performance
rmse = sqrt(mean_squared_error(test, predictions))
print('RMSE: %.3f' % rmse)
# line plot of observed vs predicted
pyplot.plot(test)
pyplot.plot(predictions)
pyplot.show()

persistence_rmse.png

通过持续模型的预测,我们得到了一个最基础的预测模型以及RMSE(baseline)为了提升我们预测模型的效果,下面让我们进入正题来构建LSTM模型来对数据集进行时间序列预测。

数据处理

为了能够构建一个LSTM模型对训练集进行训练,我们首先要对数据进行一下处理:

  1. 将时间序列问题转化成监督学习问题
  2. 平稳时间序列
  3. 数据标准化

将时间序列转换成监督学习

对于一个时间序列问题,我们可以通过使用从最后一个(t-1)时刻的观测值作为输入的特征X和当前时刻(t)的观测值作为输出Y来实现转换。

因为,需要转换的是一组时间序列数据,所以无法组合成像真正的监督学习那样有明确一对一映射的输入输出关系。尤其是在数据集的最开始或最后时,两个位置总有一个位置无法在训练集中找到对应关系。为了解决这样的问题,我们通常的做法是,在最开始时将输入特征置为0,而它对应的输出就是时间序列的第一个元素。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
from pandas import read_csv
from pandas import datetime
from pandas import DataFrame
from pandas import concat

# frame a sequence as a supervised learning problem
def timeseries_to_supervised(data, lag=1):
  df = DataFrame(data)
  columns = [df.shift(i) for i in range(1, lag+1)]
  columns.append(df)
  df = concat(columns, axis=1)
  df.fillna(0, inplace=True)
  return df

# load dataset
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('sales-of-shampoo-over-a-three-ye.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
# transform to supervised learning
X = series.values
supervised = timeseries_to_supervised(X, 1)
print(supervised)

输出结果: 0 0 0 0.0 266.0 1 266.0 145.9 2 145.9 183.1 3 183.1 119.3 4 119.3 180.3 5 180.3 168.5 6 168.5 231.8 7 231.8 224.5 8 224.5 192.8 9 192.8 122.9 10 122.9 336.5 11 336.5 185.9 12 185.9 194.3 13 194.3 149.5 14 149.5 210.1 15 210.1 273.3 16 273.3 191.4 17 191.4 287.0 18 287.0 226.0 19 226.0 303.6 20 303.6 289.9 21 289.9 421.6 22 421.6 264.5 23 264.5 342.3 24 342.3 339.7 25 339.7 440.4 26 440.4 315.9 27 315.9 439.3 28 439.3 401.3 29 401.3 437.4 30 437.4 575.5 31 575.5 407.6 32 407.6 682.0 33 682.0 475.3 34 475.3 581.3 35 581.3 646.9

平稳时间序列

虽然不明显,但我们仍可以看出这个洗发水销售数据集在时间上呈上升趋势。因此我们说这个时间序列数据是非平稳的。那么,不平稳怎么办?

答案就是:差分。(有关差分的介绍点击此处

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
from pandas import read_csv
from pandas import datetime
from pandas import Series

# create a differenced series
def difference(dataset, interval=1):
  diff = list()
  for i in range(interval, len(dataset)):
      value = dataset[i] - dataset[i - interval]
      diff.append(value)
  return Series(diff)

# invert differenced value
def inverse_difference(history, yhat, interval=1):
  return yhat + history[-interval]

# load dataset
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('sales-of-shampoo-over-a-three-ye.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
print(series.head())
# transform to be stationary
differenced = difference(series, 1)
print(differenced.head())
# invert transform
inverted = list()
for i in range(len(differenced)):
  value = inverse_difference(series, differenced[i], len(series)-i)
  inverted.append(value)
inverted = Series(inverted)
print(inverted.head())
differenced.plot()
plt.show()

diff.png

经过一阶差分处理后,从图上看还是挺平稳的。

标准化数据

在数据输入前进行标准化可以非常有效的提升收敛速度和效果。尤其如果我们的激活函数是sigmoid或者tanh,其梯度最大的区间是0附近,当输入值很大或者很小的时候,sigmoid或者tanh的变化就基本平坦了(sigmoid的导数sig(1-sig)会趋于0),也就是进行梯度下降进行优化的时候,梯度会趋于0,而倒是优化速度很慢。

如果输入不进行归一化,由于我们初始化的时候一般都是0均值的的正太分布或者小范围的均匀分布(Xavier),如果输入中存在着尺度相差很大的特征,例如(10000,0.001)这样的,很容易导致激活函数的输入w1x1+w2x2+b变的很大或者很小,从而引起梯度趋于0。

而LSTM的默认激活函数就是tanh函数,它的输出范围在-1 到 1 之间,同时这是时间序列数据的首选范围。因此我们可以使用MinMaxScaler类将数据集转换到范围[-1,1]。像其他scikit用于转换数据的方法类一样,它需要以行和列的矩阵格式提供的数据。因此,在转换之前,我们必须重塑NumPy数组。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from pandas import read_csv
from pandas import datetime
from pandas import Series
from sklearn.preprocessing import MinMaxScaler
# load dataset
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('sales-of-shampoo-over-a-three-ye.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
print(series.head())
# transform scale
X = series.values
X = X.reshape(len(X), 1)
scaler = MinMaxScaler(feature_range=(-1, 1))
scaler = scaler.fit(X)
scaled_X = scaler.transform(X)
scaled_series = Series(scaled_X[:, 0])
print(scaled_series.head())
# invert transform
inverted_X = scaler.inverse_transform(scaled_X)
inverted_series = Series(inverted_X[:, 0])
print(inverted_series.head())

输出结果: Month 1901-01-01 266.0 1901-02-01 145.9 1901-03-01 183.1 1901-04-01 119.3 1901-05-01 180.3 Name: Sales of shampoo over a three year period, dtype: float64 0 -0.478585 1 -0.905456 2 -0.773236 3 -1.000000 4 -0.783188 dtype: float64 0 266.0 1 145.9 2 183.1 3 119.3 4 180.3 dtype: float64

构建LSTM模型

长短期记忆网络(LSTM)是一种递归神经网络(RNN)。 这类网络的的优点是它能学习并记住较长序列,并不依赖预先指定的窗口滞后观察值作为输入。 在Keras中,这被称为stateful,在定义LSTM网络层时将“stateful”语句设定为“True”。

LSTM层要求输入矩阵格式为:[样本,时间步长,特征]

鉴于训练数据集的形式定义为X输入和y输出,必须先将其转化为样本/时间步长/特征的形式。

完整代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
from pandas import DataFrame
from pandas import Series
from pandas import concat
from pandas import read_csv
from pandas import datetime
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import MinMaxScaler
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from math import sqrt
from matplotlib import pyplot
import numpy

# date-time parsing function for loading the dataset
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')

# frame a sequence as a supervised learning problem
def timeseries_to_supervised(data, lag=1):
  df = DataFrame(data)
  columns = [df.shift(i) for i in range(1, lag+1)]
  columns.append(df)
  df = concat(columns, axis=1)
  df.fillna(0, inplace=True)
  return df

# create a differenced series
def difference(dataset, interval=1):
  diff = list()
  for i in range(interval, len(dataset)):
      value = dataset[i] - dataset[i - interval]
      diff.append(value)
  return Series(diff)

# invert differenced value
def inverse_difference(history, yhat, interval=1):
  return yhat + history[-interval]

# scale train and test data to [-1, 1]
def scale(train, test):
  # fit scaler
  scaler = MinMaxScaler(feature_range=(-1, 1))
  scaler = scaler.fit(train)
  # transform train
  train = train.reshape(train.shape[0], train.shape[1])
  train_scaled = scaler.transform(train)
  # transform test
  test = test.reshape(test.shape[0], test.shape[1])
  test_scaled = scaler.transform(test)
  return scaler, train_scaled, test_scaled

# inverse scaling for a forecasted value
def invert_scale(scaler, X, value):
  new_row = [x for x in X] + [value]
  array = numpy.array(new_row)
  array = array.reshape(1, len(array))
  inverted = scaler.inverse_transform(array)
  return inverted[0, -1]

# fit an LSTM network to training data
def fit_lstm(train, batch_size, nb_epoch, neurons):
  X, y = train[:, 0:-1], train[:, -1]
  X = X.reshape(X.shape[0], 1, X.shape[1])
  model = Sequential()
  model.add(LSTM(neurons, batch_input_shape=(batch_size, X.shape[1], X.shape[2]), stateful=True))
  model.add(Dense(1))
  model.compile(loss='mean_squared_error', optimizer='adam')
  for i in range(nb_epoch):
      model.fit(X, y, epochs=1, batch_size=batch_size, verbose=0, shuffle=False)
      model.reset_states()
  return model

# make a one-step forecast
def forecast_lstm(model, batch_size, X):
  X = X.reshape(1, 1, len(X))
  yhat = model.predict(X, batch_size=batch_size)
  return yhat[0,0]

# load dataset
series = read_csv('sales-of-shampoo-over-a-three-ye.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)

# transform data to be stationary
raw_values = series.values
diff_values = difference(raw_values, 1)

# transform data to be supervised learning
supervised = timeseries_to_supervised(diff_values, 1)
supervised_values = supervised.values

# split data into train and test-sets
train, test = supervised_values[0:-12], supervised_values[-12:]

# transform the scale of the data
scaler, train_scaled, test_scaled = scale(train, test)

# fit the model
lstm_model = fit_lstm(train_scaled, 1, 3000, 4)
# forecast the entire training dataset to build up state for forecasting
train_reshaped = train_scaled[:, 0].reshape(len(train_scaled), 1, 1)
lstm_model.predict(train_reshaped, batch_size=1)

# walk-forward validation on the test data
predictions = list()
for i in range(len(test_scaled)):
  # make one-step forecast
  X, y = test_scaled[i, 0:-1], test_scaled[i, -1]
  yhat = forecast_lstm(lstm_model, 1, X)
  # invert scaling
  yhat = invert_scale(scaler, X, yhat)
  # invert differencing
  yhat = inverse_difference(raw_values, yhat, len(test_scaled)+1-i)
  # store forecast
  predictions.append(yhat)
  expected = raw_values[len(train) + i + 1]
  print('Month=%d, Predicted=%f, Expected=%f' % (i+1, yhat, expected))

#report performance
rmse = sqrt(mean_squared_error(raw_values[-12:], predictions))
print('Test RMSE: %.3f' % rmse)
# line plot of observed vs predicted
pyplot.plot(raw_values[-12:])
pyplot.plot(predictions)
pyplot.show()

lstm_pred.png

最后运行结果打印出测试数据集12个月份中每个月份的预期和预测销量。示例还打印了所有预测值得均方根误差。该模型显示洗发水月度销量的均方根误差为111.925,好于持续性模型得出的对应结果136.761。

另外,神经网络的一个难题是初始条件不同,它们给出结果就不同。一种解决办法是修改Keras使用的随机数种子值以确保结果可复制。另一种办法是使用不同的实验设置控制随机初始条件。

如何免费使用谷歌GPU训练神经网络



完全云端运行:免费使用谷歌GPU训练神经网络

背景

对,你没有听错,高大上的GPU,现在不花钱也能用上了。这是Google的一项免费云端机器学习服务,全名Colaboratory。

Colaboratory 是一个 Google 研究项目,旨在帮助传播机器学习培训和研究成果。它是一个 Jupyter 笔记本环境,不需要进行任何设置就可以使用,并且完全在云端运行。Colaboratory 笔记本存储在 Google 云端硬盘中,并且可以共享,就如同您使用 Google 文档或表格一样。Colaboratory 可免费使用,而且最重要的还提供免费的英伟达Tesla K80 GPU。还有这等好事?事不宜迟,本文马上介绍如何使用 Google CoLaboratory 训练神经网络。

准备工作

在Google Drive上创建文件夹

Colab用的数据都存储在Google Drive云端硬盘上,所以,我们需要先指定要在Google Drive上用的文件夹。

比如说,可以在Google Drive上创建一个“app”文件夹,或者其他什么名字,也可以选择Colab笔记本默认的文件夹。

新建Colab笔记本

在刚刚创建的app文件夹里点击右键,选择“更多”,然后从菜单里选择“Colaboratory”,这样就新建出了一个Colab笔记本。

若是更多选项中没有“Colaboratory”选项,可以点击“关联更多应用”选项,然后在打开的页面中,搜索“Colaboratory”,然后再点关联应用,再次点击右键就可以在“更多”选项中看到“Colaboratory”选项了。

设置免费GPU

新建Colaboratory成功后,在笔记本里点Edit>Notebook settings(编辑>笔记本设置),或者Runtime>Change runtime type(运行时>改变运行时类型),然后在Hardware accelerator(硬件加速器)一栏选择GPU。

然后,Google Colab就可以用了。

关联Google Drive

为了能让Colaboratory使用到你的Google Drive的文件,我们需要先运行下面这些代码,来安装必要的库、执行授权。

!apt-get install -y -qq software-properties-common python-software-properties module-init-tools
!add-apt-repository -y ppa:alessandro-strada/ppa 2>&1 > /dev/null
!apt-get update -qq 2>&1 > /dev/null
!apt-get -y install -qq google-drive-ocamlfuse fuse
from google.colab import auth
auth.authenticate_user()
from oauth2client.client import GoogleCredentials
creds = GoogleCredentials.get_application_default()
import getpass
!google-drive-ocamlfuse -headless -id={creds.client_id} -secret={creds.client_secret} < /dev/null 2>&1 | grep URL
vcode = getpass.getpass()
!echo {vcode} | google-drive-ocamlfuse -headless -id={creds.client_id} -secret={creds.client_secret}

运行的时候应该会看到下图所示的结果:

看见那个链接之后,点击它,复制验证码并粘贴到文本框里。(这里其实是调用了Google Drive的SDK来访问你的Google Drive,而这个验证码就相当于access_key了)

授权完成后,就可以挂载Google Drive了:

!mkdir -p drive
!google-drive-ocamlfuse drive

测试GPU

这时,我们在本地电脑上创建一个.py文件来测试一下,挂载是否成功以及GPU是否在工作吧。

echo "import tensorflow as tf\nprint(tf.test.gpu_device_name())" > test.py

然后将test.py上传到我们开始时创建的app的文件夹里。

然后在Colaboratory笔记本中运行一下代码:

!python3 drive/app/test.py

不出意外的话,就会输出类似以下的结果:

/usr/local/lib/python3.6/dist-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters
2018-02-18 12:37:05.172726: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:898] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2018-02-18 12:37:05.172988: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1208] Found device 0 with properties: 
name: Tesla K80 major: 3 minor: 7 memoryClockRate(GHz): 0.8235
pciBusID: 0000:00:04.0
totalMemory: 11.17GiB freeMemory: 503.62MiB
2018-02-18 12:37:05.173016: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1308] Adding visible gpu devices: 0
2018-02-18 12:37:05.457665: I tensorflow/core/common_runtime/gpu/gpu_device.cc:989] Creating TensorFlow device (/device:GPU:0 with 243 MB memory) -> physical GPU (device: 0, name: Tesla K80, pci bus id: 0000:00:04.0, compute capability: 3.7)
/device:GPU:0

到这里的话,那么恭喜你,你的GPU环境基本可以用了,只要把你的项目文件夹上传到你的app文件夹下,搭建好深度学习的库环境,就可以通过类似上面的操作进行神经网络训练了。

Tips

如何安装库?

安装Keras:

!pip install -q keras
import keras

安装PyTorch:

!pip install -q http://download.pytorch.org/whl/cu75/torch-0.2.0.post3-cp27-cp27mu-manylinux1_x86_64.whl torchvision
import torch

安装OpenCV:

!apt-get -qq install -y libsm6 libxext6 && pip install -q -U opencv-python
import cv2

安装XGBoost:

!pip install -q xgboost==0.4a30
import xgboost

安装GraphViz:

!apt-get -qq install -y graphviz && pip install -q pydot
import pydot

安装7zip Reader:

!apt-get -qq install -y libarchive-dev && pip install -q -U libarchive
import libarchive

安装其他库:

用!pip install或者!apt-get install命令。

元胞自动机与生命游戏(Game of Life)



元胞自动机与《生命游戏》(Game of Life)

背景

笔者最近读到一篇交通流仿真的论文里,提到了一个挺有意思的模型 -- 元胞自动机

在好奇心的驱使之下查询了不少资料,所以今天就来跟大家来分享一下“元胞自动机”这个模型以及它和“康威《生命游戏》”的关系。

为了让话题更加有趣,我们先从《生命游戏》开始谈起。

什么是《生命游戏》?

生命游戏由英国数学家约翰·何顿·康威提出,它其实是一个零玩家游戏,它包括一个二维矩形世界,这个世界中的每个方格居住着一个活着的或死了的细胞。

而整个《生命游戏》是贯彻着一条生命游戏定律的,即:如果一个生命,其周围的同类生命太少,会因为得不到帮助而死亡;如果太多,则会因为得不到足够的生命资源而死亡。 ——英国数学家约翰·康威

一个细胞在下一个时刻生死取决于相邻八个方格中活着的或死了的细胞的数量。如果相邻方格活着的细胞数量过多,这个细胞会因为资源匮乏而在下一个时刻死去;相反,如果周围活细胞过少,这个细胞会因太孤单而死去。实际中,你可以设定周围活细胞的数目怎样时才适宜该细胞的生存。如果这个数目设定过低,世界中的大部分细胞会因为找不到太多的活的邻居而死去,直到整个世界都没有生命;如果这个数目设定过高,世界中又会被生命充满而没有什么变化。

实际中,这个数目一般选取2或者3;这样整个生命世界才不至于太过荒凉或拥挤,而是一种动态的平衡。

这样的话,游戏的规则就是:当一个方格周围有2或3个活细胞时,方格中的活细胞在下一个时刻继续存活;即使这个时刻方格中没有活细胞,在下一个时刻也会“诞生”活细胞。在这个游戏中,还可以设定一些更加复杂的规则,例如当前方格的状况不仅由父一代决定,而且还考虑祖父一代的情况。你还可以作为这个世界的上帝,随意设定某个方格细胞的死活,以观察对世界的影响。

在游戏的进行中,杂乱无序的细胞会逐渐演化出各种精致、有形的结构;这些结构往往有很好的对称性,而且每一代都在变化形状。一些形状已经锁定,不会逐代变化。

有时,一些已经成形的结构会因为一些无序细胞的“入侵”而被破坏。但是形状和秩序经常能从杂乱中产生出来。

《生命游戏》和元胞自动机的关系以及什么是元胞自动机?

生命游戏的原理就是基于元胞自动机的,或者说《生命游戏》就是元胞自动机的一个展示。

元胞自动机(Cellular Automata,简称CA,也有人译为细胞自动机、点格自动机、分子自动机或单元自动机)。是一时间和空间都离散的动力系统。

散布在规则格网 (Lattice Grid)中的每一元胞(Cell)取有限的离散状态,遵循同样的作用规则,依据确定的局部规则作同步更新。大量元胞通过简单的相互作用而构成动态系统的演化。不同于一般的动力学模型,元胞自动机不是由严格定义的物理方程或函数确定,而是用一系列模型构造的规则构成。凡是满足这些规则的模型都可以算作是元胞自动机模型。

因此,元胞自动机是一类模型的总称,或者说是一个方法框架。其特点是时间、空间、状态都离散,每个变量只取有限多个状态,且其状态改变的规则在时间和空间上都是局部的。

初等元胞自动机( Elementary Cellular Automata, ECA)的基本要素如下

  • 空间:一维直线上等间距的点。可为某区间上的整数点的集合。
  • 状态集:S={s1,s2} 即只有两种不同的状态。这两种不同的状态可将其分别编码为0 与 1;若用图形表示,则可对应“黑”与“白” 或者其他两种不同的颜色。
  • 邻居:取邻居半径r=1,即每个元胞最多只有“左邻右舍”两个邻居。
  • 演化规则:任意设定, 最多2^8=256

元胞以相邻的8个元胞为邻居。即Moore邻居;一个元胞的生死由其在该时刻本身的生死状态和周围八个邻居的状态。

为了解释它,我们可以把计算机中的宇宙想象成是一堆方格子构成的封闭空间,尺寸为N的空间就有N*N个格子。而每一个格子都可以看成是一个生命体,每个生命都有生和死两种状态,如果该格子生就显示蓝色,死则显示白色。每一个格子旁边都有邻居格子存在,如果我们把3x3的9个格子构成的正方形看成一个基本单位的话,那么这个正方形中心的格子的邻居就是它旁边的8个格子。

每个格子的生死遵循下面的原则:

  1. 如果一个细胞周围有3个细胞为生(一个细胞周围共有8个细胞),则该细胞为生(即该细胞若原先为死,则转为生,若原先为生,则保持不变)。
  2. 如果一个细胞周围有2个细胞为生,则该细胞的生死状态保持不变;
  3. 在其它情况下,该细胞为死(即该细胞若原先为生,则转为死,若原先为死,则保持不变)

设定图像中每个像素的初始状态后依据上述的游戏规则演绎生命的变化,由于初始状态和迭代次数不同,将会得到令人叹服的优美图案。

这样就把这些若干个格子(生命体)构成了一个复杂的动态世界。运用简单的3条作用规则构成的群体会涌现出很多意想不到的复杂行为,这就是复杂性科学的研究焦点。

元胞自动机的应用

在实际应用过程中,有的元胞自动机模型对其中的某些特征进行了扩展,有的在规则设计中引入随机因素,如:森林火灾模型。 又如,在交通、通讯发达的今天, 研究流行病或计算机病毒的传播问题时, 我们还可以将空间背景换成复杂网络的结点,用网络邻接点作为邻居。

这样的调整显然比仍旧使用二维欧氏空间、采用欧氏距离的模型更加符合实际情况。 在大型场所人群紧急疏散问题模拟研究中,可以考虑年龄、性别等因素,即元胞不是同质的,更加有利于使模拟系统接近真实系统。

《生命游戏实现》Python版

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
import numpy as np
import matplotlib.pyplot as plt


class GameOfLife(object):

  def __init__(self, cells_shape):
    """
    Parameters
    ----------
    cells_shape : 一个元组,表示画布的大小。

    Examples
    --------
    建立一个高20,宽30的画布
    game = GameOfLife((20, 30))
    
    """

    # 矩阵的四周不参与运算
    self.cells = np.zeros(cells_shape)

    real_width = cells_shape[0] - 2
    real_height = cells_shape[1] - 2

    self.cells[1:-1, 1:-1] = np.random.randint(2, size=(real_width, real_height))
    self.timer = 0
    self.mask = np.ones(9)
    self.mask[4] = 0

  def update_state(self):
    """更新一次状态"""
    buf = np.zeros(self.cells.shape)
    cells = self.cells
    for i in range(1, cells.shape[0] - 1):
      for j in range(1, cells.shape[0] - 1):
        # 计算该细胞周围的存活细胞数
        neighbor = cells[i-1:i+2, j-1:j+2].reshape((-1, ))
        neighbor_num = np.convolve(self.mask, neighbor, 'valid')[0]
        if neighbor_num == 3:
          buf[i, j] = 1
        elif neighbor_num == 2:
          buf[i, j] = cells[i, j]
        else:
          buf[i, j] = 0
    self.cells = buf
    self.timer += 1

  def plot_state(self):
    """画出当前的状态"""
    plt.title('Iter :{}'.format(self.timer))
    plt.imshow(self.cells)
    plt.show()

  def update_and_plot(self, n_iter):
    """更新状态并画图
    Parameters
    ----------
    n_iter : 更新的轮数
    """
    plt.ion()
    for _ in range(n_iter):
      plt.title('Iter :{}'.format(self.timer))
      plt.imshow(self.cells)
      self.update_state()
      plt.pause(0.2)
    plt.ioff()


if __name__ == '__main__':
  game = GameOfLife(cells_shape=(60, 60))
  game.update_and_plot(200)

效果图