## RNN基本↺

### 什么是RNN？↺

$\left({\mathbit{x}}_{1},{\mathbit{x}}_{2},\dots ,{\mathbit{x}}_{T}\right)$

，输出另外一个向量序列

$\left({\mathbit{y}}_{1},{\mathbit{y}}_{2},\dots ,{\mathbit{y}}_{T}\right)$

，并且满足如下递归关系

$\begin{array}{}\text{(1)}& {\mathbit{y}}_{t}=f\left({\mathbit{y}}_{t-1},{\mathbit{x}}_{t},t\right)\end{array}$

### 自己编写RNN↺

$\begin{array}{}\text{(2)}& {\mathbit{y}}_{t}=\mathrm{tanh}\left({\mathbit{W}}_{1}{\mathbit{y}}_{t-1}+{\mathbit{W}}_{2}{\mathbit{x}}_{t}+\mathbit{b}\right)\end{array}$

#! -*- coding: utf-8- -*-from keras.layers import Layerimport keras.backend as Kclass My_RNN(Layer):    def __init__(self, output_dim, **kwargs):        self.output_dim = output_dim # 输出维度        super(My_RNN, self).__init__(**kwargs)    def build(self, input_shape): # 定义可训练参数        self.kernel1 = self.add_weight(name='kernel1',                                      shape=(self.output_dim, self.output_dim),                                      initializer='glorot_normal',                                      trainable=True)        self.kernel2 = self.add_weight(name='kernel2',                                      shape=(input_shape[-1], self.output_dim),                                      initializer='glorot_normal',                                      trainable=True)        self.bias = self.add_weight(name='kernel',                                      shape=(self.output_dim,),                                      initializer='glorot_normal',                                      trainable=True)    def step_do(self, step_in, states): # 定义每一步的迭代        step_out = K.tanh(K.dot(states[0], self.kernel1) +                          K.dot(step_in, self.kernel2) +                          self.bias)        return step_out, [step_out]    def call(self, inputs): # 定义正式执行的函数        init_states = [K.zeros((K.shape(inputs)[0],                                self.output_dim)                              )] # 定义初始态(全零)        outputs = K.rnn(self.step_do, inputs, init_states) # 循环执行step_do函数        return outputs[0] # outputs是一个tuple，outputs[0]为最后时刻的输出，                          # outputs[1]为整个输出的时间序列，output[2]是一个list，                          # 是中间的隐藏状态。    def compute_output_shape(self, input_shape):        return (input_shape[0], self.output_dim)

${\mathbit{x}}_{t}$

，而states是一个list，代表

${\mathbit{y}}_{t-1}$

${\mathbit{y}}_{t-1}$

${\mathbit{y}}_{t}$

## ODE基本↺

### 什么是ODE？↺

ODE就是“常微分方程（Ordinary Differential Equation）”，这里指的是一般的常微分方程组：

$\begin{array}{}\text{(3)}& \stackrel{˙}{\mathbit{x}}\left(t\right)=\mathbit{f}\left(\mathbit{x}\left(t\right),t\right)\end{array}$

ODE可以产生非常丰富的函数。比如

${e}^{t}$

$\stackrel{˙}{x}=x$

$\mathrm{sin}t$

$\mathrm{cos}t$

$\stackrel{¨}{x}+x=0$

$\stackrel{˙}{x}=x$

${e}^{t}$

### 数值解ODE↺

ODE的数值解已经是一门非常成熟的学科了，这里我们也不多做介绍，仅引入最基本的由数学家欧拉提出来的迭代公式：

$\begin{array}{}\text{(4)}& \mathbit{x}\left(t+h\right)=\mathbit{x}\left(t\right)+h\mathbit{f}\left(\mathbit{x}\left(t\right),t\right)\end{array}$

$h$

$\begin{array}{}\text{(5)}& \frac{\mathbit{x}\left(t+h\right)-\mathbit{x}\left(t\right)}{h}\end{array}$

$\stackrel{˙}{\mathbit{x}}\left(t\right)$

。只要给定初始条件

$\mathbit{x}\left(0\right)$

，我们就可以根据

$\left(4\right)$

## ODE与RNN↺

### ODE也是RNN↺

$\left(4\right)$

$\left(1\right)$

，发现有什么联系了吗？

$\left(1\right)$

$t$

$\left(4\right)$

$t$

$\left(4\right)$

$\left(1\right)$

$\left(4\right)$

$h$

$t=nh$

，那么

$\left(4\right)$

$\begin{array}{}\text{(6)}& \mathbit{x}\left(\left(n+1\right)h\right)=\mathbit{x}\left(nh\right)+h\mathbit{f}\left(\mathbit{x}\left(nh\right),nh\right)\end{array}$

$\left(6\right)$

$n$

$\left(4\right)$

### 用RNN解ODE↺

$\begin{array}{}\text{(7)}& \left\{\begin{array}{r}\frac{d{x}_{1}}{dt}={r}_{1}{x}_{1}\left(1-\frac{{x}_{1}}{{N}_{1}}\right)-{a}_{1}{x}_{1}{x}_{2}\\ \frac{d{x}_{2}}{dt}={r}_{2}{x}_{2}\left(1-\frac{{x}_{2}}{{N}_{2}}\right)-{a}_{2}{x}_{1}{x}_{2}\end{array}\end{array}$

#! -*- coding: utf-8- -*-from keras.layers import Layerimport keras.backend as Kclass ODE_RNN(Layer):    def __init__(self, steps, h, **kwargs):        self.steps = steps        self.h = h        super(ODE_RNN, self).__init__(**kwargs)    def step_do(self, step_in, states): # 定义每一步的迭代        x = states[0]        r1,r2,a1,a2,iN1,iN2 = 0.1,0.3,0.0001,0.0002,0.002,0.003        _1 = r1 * x[:,0] * (1 - iN1 * x[:,0]) - a1 * x[:,0] * x[:,1]        _2 = r2 * x[:,1] * (1 - iN2 * x[:,1]) - a2 * x[:,0] * x[:,1]        _1 = K.expand_dims(_1, 1)        _2 = K.expand_dims(_2, 1)        _ = K.concatenate([_1, _2], 1)        step_out = x + self.h * _        return step_out, [step_out]    def call(self, inputs): # 这里的inputs就是初始条件        init_states = [inputs]        zeros = K.zeros((K.shape(inputs)[0],                         self.steps,                         K.shape(inputs)[1])) # 迭代过程用不着外部输入，所以                                              # 指定一个全零输入，只为形式上的传入        outputs = K.rnn(self.step_do, zeros, init_states) # 循环执行step_do函数        return outputs[1] # 这次我们输出整个结果序列    def compute_output_shape(self, input_shape):        return (input_shape[0], self.steps, input_shape[1])from keras.models import Sequentialimport numpy as npimport matplotlib.pyplot as pltsteps,h = 1000,0.1M = Sequential()M.add(ODE_RNN(steps, h, input_shape=(2,)))M.summary()# 直接前向传播就输出解了result = M.predict(np.array([[100, 150]]))[0] # 以[100, 150]为初始条件进行演算times = np.arange(1, steps+1) * h# 绘图plt.plot(times, result[:,0])plt.plot(times, result[:,1])plt.savefig('test.png')

$\left(7\right)$

RNN解两物种的竞争模型

### 反推ODE参数↺

$\overline{)\begin{array}{cccccccc}\text{时间}& 0& 10& 15& 30& 36& 40& 42\\ {x}_{1}& 100& 165& 197& 280& 305& 318& 324\\ {x}_{2}& 150& 283& 290& 276& 269& 266& 264\end{array}}$

$\left(7\right)$

#! -*- coding: utf-8- -*-from keras.layers import Layerimport keras.backend as Kdef my_init(shape, dtype=None): # 需要定义好初始化，这相当于需要实验估计参数的量级    return K.variable([0.1, 0.1, 0.001, 0.001, 0.001, 0.001])class ODE_RNN(Layer):        def __init__(self, steps, h, **kwargs):        self.steps = steps        self.h = h        super(ODE_RNN, self).__init__(**kwargs)        def build(self, input_shape): # 将原来的参数设为可训练的参数        self.kernel = self.add_weight(name='kernel',                                       shape=(6,),                                      initializer=my_init,                                      trainable=True)    def step_do(self, step_in, states): # 定义每一步的迭代        x = states[0]        r1,r2,a1,a2,iN1,iN2 = (self.kernel[0], self.kernel[1],                               self.kernel[2], self.kernel[3],                               self.kernel[4], self.kernel[5])        _1 = r1 * x[:,0] * (1 - iN1 * x[:,0]) - a1 * x[:,0] * x[:,1]        _2 = r2 * x[:,1] * (1 - iN2 * x[:,1]) - a2 * x[:,0] * x[:,1]        _1 = K.expand_dims(_1, 1)        _2 = K.expand_dims(_2, 1)        _ = K.concatenate([_1, _2], 1)        step_out = x + self.h * K.clip(_, -1e5, 1e5) # 防止梯度爆炸        return step_out, [step_out]        def call(self, inputs): # 这里的inputs就是初始条件        init_states = [inputs]        zeros = K.zeros((K.shape(inputs)[0],                         self.steps,
K.shape(inputs)[1])) # 迭代过程用不着外部输入，所以                                              # 指定一个全零输入，只为形式上的传入        outputs = K.rnn(self.step_do, zeros, init_states) # 循环执行step_do函数        return outputs[1] # 这次我们输出整个结果序列        def compute_output_shape(self, input_shape):        return (input_shape[0], self.steps, input_shape[1])from keras.models import Sequentialfrom keras.optimizers import Adamimport numpy as npimport matplotlib.pyplot as pltsteps,h = 50, 1 # 用大步长，减少步数，削弱长时依赖，也加快推断速度series = {0: [100, 150],          10: [165, 283],          15: [197, 290],          30: [280, 276],          36: [305, 269],          40: [318, 266],          42: [324, 264]}M = Sequential()M.add(ODE_RNN(steps, h, input_shape=(2,)))M.summary()# 构建训练样本# 其实就只有一个样本序列，X为初始条件，Y为后续时间序列X = np.array([series[0]])Y = np.zeros((1, steps, 2))for i,j in series.items():    if i != 0:        Y[0, int(i/h)-1] += series[i]# 自定义loss# 在训练的时候，只考虑有数据的几个时刻，没有数据的时刻被忽略def ode_loss(y_true, y_pred):    T = K.sum(K.abs(y_true), 2, keepdims=True)    T = K.cast(K.greater(T, 1e-3), 'float32')    return K.sum(T * K.square(y_true - y_pred), [1, 2])M.compile(loss=ode_loss,          optimizer=Adam(1e-4))M.fit(X, Y, epochs=10000) # 用低学习率训练足够多轮# 用训练出来的模型重新预测，绘图，比较结果result = M.predict(np.array([[100, 150]]))[0]times = np.arange(1, steps+1) * hplt.clf()plt.plot(times, result[:,0], color='blue')plt.plot(times, result[:,1], color='green')plt.plot(series.keys(), [i[0] for i in series.values()], 'o', color='blue')plt.plot(series.keys(), [i[1] for i in series.values()], 'o', color='green')plt.savefig('test.png')

RNN做ODE的参数估计效果（散点：有限的实验数据，曲线：估计出来的模型）