回声状态网络ESN（含python、MATLAB代码数据）

时间序列预测：把获得的数据按时间顺序排序，分析其变化方向和程度，用来预测未来若干时期可能达到的水平。

非线性自回归时间序列：特性是连续的观测数据比较长时间间隔的观测数据更具有相关性

e.g ：Mackey-Glass(MG)时间序列，其表达式如下：

回声状态网络ESN：Echo State Network

• 可用于时间序列预测
• 由输入层、储备池、输出层构成
• 用一个储备池模式来代替传统神经网络全接连的隐层
• 储备池中神经元的连接状态是随机的，即储备池中哪些神经元会相连，连接权重都是随机初始化的，且不会再改变
• 因此只需训练输出层来拟合数据，无需像传统的多层感知机一样通过梯度下降不断修改权重。
• 输入权重在开始时随机生成不变

实例：

• 数据集：MackeyGlass_t17.txt
• 步骤：
• 将序列前1000用于训练，后1000测试

1. 随机初始化输入权值与储备池权值
2. 依次输入数据，该时刻储备池接收两个方向输入：t时刻的输入及上一时刻储备池输出状态。加上动量因子（可看做储备池的更新速率，不用就设为1）   
3. 由目标输出Yt和计算好的储备池状态矩阵计算输出权值
4. ’用训练得的输出矩阵直接用于测试（前一时刻状态训练时已得）

• 效果图

matlab:

% A minimalistic sparse Echo State Networks demo with Mackey-Glass (delay 17) data % in "plain" Octave/Matlab. % by Mantas Lukosevicius 2012-2018 % http://mantas.info % load the data trainLen = 2000; testLen = 2000; initLen = 100; data = load('MackeyGlass_t17.txt'); % plot some of it figure(10); plot(data(1:1000)); title('A sample of data'); % generate the ESN reservoir inSize = 1; outSize = 1; resSize = 1000; a = 0.3; % leaking rate rand( 'seed', 42 ); Win = (rand(resSize,1+inSize)-0.5) .* 1; % dense W: W = rand(resSize,resSize)-0.5; % sparse W: % W = sprand(resSize,resSize,0.01); % W_mask = (W~=0); % W(W_mask) = (W(W_mask)-0.5); % normalizing and setting spectral radius disp 'Computing spectral radius...'; opt.disp = 0; rhoW = abs(eigs(W,1,'LM',opt)); disp 'done.' W = W .* ( 1.25 /rhoW); % allocated memory for the design (collected states) matrix X = zeros(1+inSize+resSize,trainLen-initLen); % set the corresponding target matrix directly Yt = data(initLen+2:trainLen+1)'; % run the reservoir with the data and collect X x = zeros(resSize,1); for t = 1:trainLen u = data(t); x = (1-a)*x + a*tanh( Win*[1;u] + W*x ); if t > initLen X(:,t-initLen) = [1;u;x]; end end % train the output by ridge regression reg = 1e-8; % regularization coefficient % direct equations from texts: %X_T = X'; %Wout = Yt*X_T * inv(X*X_T + reg*eye(1+inSize+resSize)); % using Matlab solver: Wout = ((X*X' + reg*eye(1+inSize+resSize)) \ (X*Yt'))'; % run the trained ESN in a generative mode. no need to initialize here, % because x is initialized with training data and we continue from there. Y = zeros(outSize,testLen); u = data(trainLen+1); for t = 1:testLen x = (1-a)*x + a*tanh( Win*[1;u] + W*x ); y = Wout*[1;u;x]; Y(:,t) = y; % generative mode: u = y; % this would be a predictive mode: %u = data(trainLen+t+1); end errorLen = 500; mse = sum((data(trainLen+2:trainLen+errorLen+1)'-Y(1,1:errorLen)).^2)./errorLen; disp( ['MSE = ', num2str( mse )] ); % plot some signals figure(1); plot( data(trainLen+2:trainLen+testLen+1), 'color', [0,0.75,0] ); hold on; plot( Y', 'b' ); hold off; axis tight; title('Target and generated signals y(n) starting at n=0'); legend('Target signal', 'Free-running predicted signal'); figure(2); plot( X(1:20,1:200)' ); title('Some reservoir activations x(n)'); figure(3); bar( Wout' ) title('Output weights W^{out}');

python: 

from numpy import * from matplotlib.pyplot import * import scipy.linalg import matplotlib.pyplot as plt # 加载数据 # 前2000个数据用来训练，2001-4000的数据用来测试。训练数据中，前100项用来初始化储备池，以让储备池中形成良好的回声之后再开始训练。 trainLen = 2000 testLen = 2000 initLen = 100#前100项用来初始化储备池 data = loadtxt(r'data\MackeyGlass_t17.txt') # 绘制前1000条数据 figure(0).clear() plot(data[0:1000]) title('A sample of data') # 生成ESN储层 inSize = outSize = 1 # inSize 输入维数 K resSize = 1000 # 储备池规模 N a = 0.3 # 可以看作储备池更新的速度，可不加，即设为1. random.seed(42) # 随机初始化 Win 和 W Win = (random.rand(resSize, 1 + inSize) - 0.5) * 1 # 输入矩阵 N * 1+K W = random.rand(resSize, resSize) - 0.5 # 储备池连接矩阵 N * N # 对W进行防缩，以满足稀疏的要求。 # 方案 1 - 直接缩放 (快且有脏数据, 特定储层): W *= 0.135 # 方案 2 - 归一化并设置谱半径 (正确, 慢): print('计算谱半径...') rhoW = max(abs(linalg.eig(W)[0])) # linalg.eig(W)[0]:特征值 linalg.eig(W)[1]:特征向量 W *= 1.25 / rhoW # 为设计（收集状态）矩阵分配内存 X = zeros((1 + inSize + resSize, trainLen - initLen)) # 储备池的状态矩阵x(t)：每一列是每个时刻的储备池状态。后面会转置 # 直接设置相应的目标矩阵 Yt = data[None, initLen + 1:trainLen + 1] # 输出矩阵:每一行是一个时刻的输出 # 输入所有的训练数据，然后得到每一时刻的输入值和储备池状态。 x = zeros((resSize, 1)) for t in range(trainLen): u = data[t] x = (1 - a) * x + a * tanh(dot(Win, vstack((1, u))) + dot(W, x)) # vstack((1, u)):将偏置量1加入输入序列 if t >= initLen: # 空转100次后，开始记录储备池状态 X[:, t - initLen] = vstack((1, u, x))[:, 0] # 使用Wout根据输入值和储备池状态去拟合目标值，这是一个简单的线性回归问题，这里使用的是岭回归(Ridge Regression)。 reg = 1e-8 # 正则化系数 X_T = X.T # Wout: 1 * 1+K+N Wout = dot(dot(Yt, X_T), linalg.inv(dot(X, X_T) + \ reg * eye( 1 + inSize + resSize))) # linalg.inv矩阵求逆；numpy.eye()生成对角矩阵，规模:1+inSize+resSize，默认对角线全1，其余全0 # Wout = dot( Yt, linalg.pinv(X) ) # 使用训练数据进行前向处理得到结果 # run the trained ESN in a generative mode. no need to initialize here, # because x is initialized with training data and we continue from there. Y = zeros((outSize, testLen)) u = data[trainLen] for t in range(testLen): x = (1 - a) * x + a * tanh(dot(Win, vstack((1, u))) + dot(W, x)) y = dot(Wout, vstack((1, u, x))) # 输出矩阵(1 * 1+K+N)*此刻状态矩阵(1+K+N * 1)=此刻预测值 Y[:, t] = y # t时刻的预测值 Y: 1 * testLen # 生成模型 u = y # 预测模型 # u = data[trainLen+t+1] # 计算第一个errorLen时间步长的MSE errorLen = 500 mse = sum(square(data[trainLen + 1:trainLen + errorLen + 1] - Y[0, 0: errorLen])) / errorLen print('MSE = {0}'.format(str(mse))) # 绘制测试集的真实数据和预测数据 figure(1).clear() plot(data[trainLen + 1:trainLen + testLen + 1], 'g') plot(Y.T, 'b') title('Target and generated signals $y(n)$ starting at $n=0$') legend(['Target signal', 'Free-running predicted signal']) # 绘制储备池中前200个时刻状态(x(t))的前20个储层结点值 figure(2).clear() plot(X[0:20, 0:200].T) title('Some reservoir activations $\mathbf{x}(n)$') # 绘制输出矩阵 figure(3).clear() #bar(np.arange(1 + inSize + resSize), Wout.T, 8) plot(np.arange(1 + inSize + resSize),Wout.T) title('Output weights $\mathbf{W}^{out}$') show() 

 数据：MackeyGlass_t17.txt

https://blog.csdn.net/wangqianqianya/article/details/