matlab实习实验报告

实验一时间序列ARMA学习

一、基本形式

ARMA模型分为以下三种：

1.自回归模型（AR）

如果时间序列满足

其中是独立同分布的随机变量序列，且满足：

以及E( ) = 0，则称时间序列为服从p阶的自回归模型。

AR模型是用过去的观测值和现在的干扰值的线性组合来预测未来值，式子中的p为阶数，β为待定参数，y为一个平稳序列。

自回归模型的平稳条件：

滞后算子多项式

的根均在单位圆外，即φ(B) = 0的根大于1。

2.移动平均模型（MA）

如果时间序列满足

则称时间序列为服从q阶移动平均模型；

MA模型是用过去的干扰值和现在的干扰值的线性组合预测未来的值，式中q为阶数，α为待定参数，y为一个平稳序列。

移动平均模型平稳条件：任何条件下都平稳。

3.自回归滑动平均模型（ARMA）

如果时间序列满足：

则称时间序列为服从(p,q)阶自回归滑动平均混合模型。

模型求解流程图：

二、实例分析

太阳的光球表面有时会出现一些暗的区域，它是磁场聚集的地方，这就是太阳黑子。黑子是太阳表面可以看到的最突出的现象。一个中等大小的黑子大概和地球的大小差不多。长期的观测发现，黑子多的时候，其他太阳活动现象也会比较频繁。从外文网站查询到太阳黑子数量的变化，从而预测接下来数量的走势。

首先整理数据，进行平稳性判断，使用matlab语言中的adftest()函数即可。

数据为平稳数据，因而无需做平稳化处理，使用AR模型对原数据作出预测。

Matlab代码为：

clc

clear

a=csvread('heizi.csv');

b=a(2700:3234,[1,4]);

c=b(:,2);

diff(c);

adftest(c)%

figure(2);

subplot(2,1,1)

autocorr(c,500)%2

set(gca,'Xlim',[0 500]);

subplot(2,1,2)

parcorr(c,50);

set(gca,'Xlim',[0 50]);

z=iddata(c);%×

m=armax(z(1:end),'na',20,'nc',0);%

yp = predict(m,c,1);

figure(3)

plot(c,'-.');

hold on

plot(yp,'r')

grid

legend('3?ê??μ','?¤2a?μAR')

bzc=sum((yp-c).^2)

结果如下：

2.MA模型求解

Matlab代码：

clc

clear

a=csvread('heizi.csv');

b=a(2700:3234,[1,4]);

c=b(:,2);

diff(c);

adftest(c)

figure(2);

subplot(2,1,1)

autocorr(c,500)

set(gca,'Xlim',[0 500]);

subplot(2,1,2)

parcorr(c,50);

set(gca,'Xlim',[0 50]);

z=iddata(c);

m=armax(z(1:end),'na',0,'nc',4); yp = predict(m,c,1);

figure(3)

plot(c,'-.');

hold on

plot(yp,'r')

grid

legend('3?ê??μ','?¤2a?μMA')

bzc=sum((yp-c).^2)

结果如下：

3.ARMA模型求解

对于ARMA模型，需要做出阶数判断，用ACF和PACF确定p、q阶数,在使用ARMA模型求解。

Matlab代码：

clc

clear

a=csvread('heizi.csv');

b=a(2700:3234,[1,4]);

c=b(:,2);

diff(c);

adftest(c)

figure(2);

subplot(2,1,1)

autocorr(c,500)%2自相关

set(gca,'Xlim',[0 500]);

subplot(2,1,2)

parcorr(c,50); %偏自相关

set(gca,'Xlim',[0 50]);

结果如下：

由图可以看出ACF和PACF都是拖尾，故可用ARMA模型求解，模拟数据matlab代码如下：

a=csvread('heizi.csv');

b=a(2700:3234,[1,4]);

c=b(:,2);

diff(c);

adftest(c)%

figure(2);

subplot(2,1,1)

autocorr(c,500)%

set(gca,'Xlim',[0 500]);

subplot(2,1,2)

parcorr(c,50); %

set(gca,'Xlim',[0 50]);

z=iddata(c);

m=armax(z(1:end),'na',19,'nc',19);

yp = predict(m,c,1);

figure(3)

plot(c,'-.');

hold on

plot(yp,'r')

grid

legend('3?ê??μ','?¤2a?μMA')

bzc=sum((yp-c).^2)

结果如下：

比较图像可看出，ARMA数据的整体误差小于AR和MA数据的误差，因而

三、改进方法

通过讨论，我们决定对ARMA模型做改进，主要是在定阶和模型的误差方面作出改进。传统的ARMA模型定阶方法是比较繁琐的，对于有些不太好处理的数据，很难判断出符合的阶数，误差也比较大。对于上述问题，我们做了合理的改进，主要改进方法如下：

1.定阶的方法（确定模型p,q）

在传统的ARMA模型中，定阶方法有观察自相关偏相关图像确定阶数，aic 最小原则确定阶数，等等。但是在实际操作中，我们发现这些步骤不好做，对于上述例子来说，自相关图像就不好观察，无法确定确切的拖尾截尾的点，因而无法判断。使用aic最小原则通常比较繁琐，数据量大的时候计算又耗费时间，在例题中为了得到aic最小的pq,我们做了实际计算，pq均在20以内时，计算完需要一分钟以上，效率不高，如果有更多的数据，将无法处理。对于这个问题，我们想到用神经网络训练让其判断数据的p，q。需要大量的原始数据来训练网络，用R语言生成符合时间序列阶数的大量数据，操作如下：

R语言生成数据：arima.sim()

分别取p,q为0,1,2,3,4,5生成符合的平稳数据，各生成1100个数据，1000个用来训练网络，100用来检测。

Python代码如下：

from __future__ import print_function

import numpy as np

import random as r

from https://www.360docs.net/doc/9514700465.html,yers.core import Dense, Dropout, Activation

from keras.models import Sequential

from keras.models import load_model

def getdata(p, q):

'''

生成数据的函数

输入阶数p，q，输出对应的平稳序列

'''

if p == 0 and q == 0:

a = r.random()

return np.array([a for i in range(20)])

elif p == 0 and q != 0:

elps, e = np.random.random(size=20), np.random.normal(size=20)

e = e / np.max(e)

cita = np.random.normal(0, 1, q)

cita = cita / np.max(cita)

y = np.zeros(20)

for i in range(q, 20):

y[i] = elps[i] - np.dot(cita, e[i - q:i])

y = y / np.max(y)

return y

elif p != 0 and q == 0:

elps, e = np.random.random(size=20), np.random.normal(size=20)

e = e / np.max(e)

fai = np.random.random(size=p)

if p >= q:

y = np.append(fai[0:p], np.zeros(20 - p))

for i in range(p, 20):

y[i] = -np.dot(fai, y[i - p:i]) + elps[i]

y = y / np.max(y)

return y

else:

elps, e = np.random.random(size=20), np.random.normal(size=20)

e = e / np.max(e)

fai, cita = np.random.random(size=p), np.random.normal(0, 1, q)

cita = cita / np.max(cita)

if p >= q:

y = np.append(fai[0:p], np.zeros(20 - p))

for i in range(p, 20):

y[i] = -np.dot(fai, y[i - p:i]) + elps[i] - \

np.dot(cita, e[i - q:i])

else:

y = np.append(fai[0:p], np.zeros(20 - p))

for i in range(q, 20):

y[i] = -np.dot(fai, y[i - p:i]) + elps[i] - \

np.dot(cita, e[i - q:i])

y = y / np.max(y)

return y

def getlabel(p, q):

'''

对标签进行独热编码

'''

y = np.zeros(36)

y[p * 6 + q] = 1

return y

# 生成数据并归一化

x_train, y_train, x_test, y_test = [], [], [], []

x, y, x1 = [], [], []

with open('data1.txt', 'r') as f:

lines = f.readlines()

for line in lines:

x1.append(float(line))

for i in range(500):

x_train.append(np.array(x1[i * 20:i * 20 + 20])) y_train.append(getlabel(0, 2))

x, y, x1 = [], [], []

with open('data2.txt', 'r') as f:

lines = f.readlines()

for line in lines:

x1.append(float(line))

for i in range(500):

x_train.append(np.array(x1[i * 20:i * 20 + 20])) y_train.append(getlabel(2, 0))

ps, qs = [0, 1, 2, 3, 4, 5], [0, 1, 2, 3, 4, 5]

for p in ps:

for q in qs:

for i in range(1000):

x_train.append(getdata(p, q))

y_train.append(getlabel(p, q))

x_train = np.array(x_train)

y_train = np.array(y_train)

for p in ps:

for q in qs:

for i in range(100):

x_test.append(getdata(p, q))

y_test.append(getlabel(p, q))

x_test = np.array(x_test)

y_test = np.array(y_test)

# 设置初始随机数种子

np.random.seed(1337)

# 设置batch的大小

batch_size = 100

# 设置类别的个数

nb_classes = 36

# 设置迭代的次数

nb_epoch = 5000

train_data, train_labels, test_data, test_labels = x_train, y_train, x_test, y_test

'''

模型一

'''

# 建立顺序型模型

model1 = Sequential()

# 输入层有20个神经元

# 第一个隐层有15个神经元，激活函数为ReLu，Dropout比例为0.2

model1.add(Dense(15, input_shape=(20,)))

model1.add(Activation('relu'))

model1.add(Dropout(0.2))

# 第二个隐层有18个神经元，激活函数为ReLu，Dropout比例为0.2

model11.add(Dense(18))

model11.add(Activation('relu'))

model11.add(Dropout(0.2))

# 第三个隐层有15个神经元，激活函数为ReLu，Dropout比例为0.2

model1.add(Dense(15))

model1.add(Activation('relu'))

model1.add(Dropout(0.2))

# 输出层有36个神经元，激活函数为SoftMax，得到分类结果

model1.add(Dense(36))

model1.add(Activation('softmax'))

# 输出模型的整体信息

model1.summary()

https://www.360docs.net/doc/9514700465.html,pile(loss='categorical_crossentropy',

optimizer='adam',

metrics=['accuracy'])

history = model1.fit(train_data, train_labels,

batch_size=batch_size,

epochs=nb_epoch,

verbose=1,

validation_data=(test_data, test_labels))

score = model1.evaluate(test_data, test_labels, verbose=0)

# model1.save('model1.h5')

# 输出训练好的模型在测试集上的表现

print('Test score:', score[0])

print('Test accuracy:', score[1])

'''

模型二

'''

# 建立顺序型模型

model2 = Sequential()

# 输入层有20个神经元

# 隐层有100个神经元，激活函数为ReLu，Dropout比例为0.2

model2.add(Dense(100, input_shape=(20,)))

model2.add(Activation('relu'))

model2.add(Dropout(0.2))

# 输出层有36个神经元，激活函数为SoftMax，得到分类结果

model2.add(Dense(36))

model2.add(Activation('softmax'))

# 输出模型的整体信息

model2.summary()

https://www.360docs.net/doc/9514700465.html,pile(loss='categorical_crossentropy',

optimizer='adam',

metrics=['accuracy'])

history = model2.fit(train_data, train_labels,

batch_size=batch_size,

epochs=nb_epoch,

verbose=1,

validation_data=(test_data, test_labels))

score = model2.evaluate(test_data, test_labels, verbose=0)

# model2.save('model2.h5')

# 输出训练好的模型在测试集上的表现

print('Test score:', score[0])

print('Test accuracy:', score[1])

训练完成之后，用余下的100组数据进行检验，模型一正确率为0.65，模型二正确率为0.82，可见深网络比宽网络更准确。

由此，p,q的阶数确定将变得简单，缺点是准确率为0.82，有出错的可能。四、参考文献

[1]司守奎，孙兆亮.数学建模算法与应用.北京：国防工业出版社，2017.

[2]李维国，同登科.数值计算方法.山东：中国石油大学出版社，2016.

[3]郑阿奇，曹戈.MATLAB实用教程.北京：电子工业出版社，2016.

[4]高俊芳,吴清.时间序列ARMA模型及其应用[J].上海工程技术大学学

报,1996(04):68-73.

[5]胡代平,刘豹.利用神经网络确定ARMA模型的结构[J].1998(04):34-36+30.

[6]colah's blog.Understanding LSTM Networks.

http://colah.github.io/posts/2015-08-Understanding-LSTMs/

[7]DL4J.LSTM和循环网络基础教程.https://https://www.360docs.net/doc/9514700465.html,/cn/lstm

实验二MATLAB混合编程

一、Mex混合编程技术学习

为matlab配置混合编程的软件，选用C++.

mex文件的编写可以分为三个步骤，数据传入，数据处理，数据导出

1、传入矩阵用mxGetPr，传入数值用mxGetScalar

2、获取2维矩阵维度用mxGetM，mxGetN

3、数据在内存上都是连续的，访问的时候按列的顺序访问

4、数据传出用mxCreateDoubleMatrix，先创建一个double类型的指针用来存放要导出的数据，然后再拷贝到plhs[i]对应的指针上

二、Mex混合编程实现

通过一个案例来说明一中的具体操作，

#include"mex.h"

//计算过程

void juzhenplus(double *y ,double *p,int n,int m)

{

int i, j;

for (i = 0; i < n; i++)

for (j = 0; j < m; j++)

*(y+j+n*i)=*(p+i*n+j) * 5;

}

//接口过程

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])

{

double x, *y, *p;

int n,t;

x = mxGetM(prhs[0]);

t = mxGetN(prhs[0]);

plhs[0] = mxCreateDoubleMatrix(x, t, mxREAL);

y = mxGetPr(plhs[0]);

p = mxGetPr(prhs[0]);

juzhenplus(y,p,x,t);

}

在VS2017中编写上面的代码，实现的是从matlab传入一个矩阵，在c++中进行乘以常数的操作，再输出到matlab中。

结果显示，混用是成功的。

实验三三次样条插值法与勒让德高斯

积分法的应用

一、问题提出

在实际生活中我们会遇到这种问题，如何求得一段非常不规则的曲线长度，曲曲折折我们无法用测量工具来准确地测得长度。但是有时我们需要测量这些数据，比如海岸线的长度。本方案将对此提出解决方法。

问题来源：台湾省边界长度

用matlab识别出图片轮廓，提取边界上的点。

matlab边缘检测即可得出所有需要的点，通过插值的方法使计算误差变小. 二、三次样条插值

三次样条插值（Cubic Spline Interpolation）简称Spline插值，是通过一系

列形值点的一条光滑曲线，数学上通过求解三弯矩方程组得出曲线函数组的过程。

实际计算时还需要引入边界条件才能完成计算。一般的计算方法书上都没有说明非扭结边界的定义，但数值计算软件如Matlab都把非扭结边界条件作为默认的边界条件。

三次样条函数:

定义:函数S(x)∈C2[a,b] ，且在每个小区间[ xj,xj+1 ]上是三次多项式，其中

a =x0

b 是给定节点，则称S(x)是节点x0,x1,...xn上的三次样条函数。

若在节点x j 上给定函数值Yj= f (Xj).( j =0, 1, , n) ，并成立

S(xj ) =yj .( j= 0, 1, , n) ，则称S(x)为三次样条插值函数。

实际计算时还需要引入边界条件才能完成计算。边界通常有自然边界（边界点的二阶导为0），夹持边界（边界点导数给定），非扭结边界（使两端点的三阶导与这两端点的邻近点的三阶导相等）。一般的计算方法书上都没有说明非扭结边界的定义，但数值计算软件如Matlab都把非扭结边界条件作为默认的边界条件。

已知:

a. n+1个数据点[xi, yi], i = 0, 1, …, n

b. 每一分段都是三次多项式函数曲线

c. 节点达到二阶连续

d. 左右两端点处特性（自然边界，固定边界，非节点边界）

根据定点，求出每段样条曲线方程中的系数，即可得到每段曲线的具体表达式。

插值和连续性：

, 其中i = 0, 1, …, n-1

微分连续性：

, 其中i = 0, 1, …, n-2

样条曲线的微分式：

将步长带入样条曲线的条件：a. 由(i = 0, 1, …, n-1)推出

b. 由(i = 0, 1, …, n-1)推出

c. 由(i = 0, 1, …, n-2)推出

由此可得：

d. 由(i = 0, 1, …, n-2)推出

设，则

a. 可写为：

，推出

b. 将ci, di带入可得：

c. 将bi, ci, di带入(i = 0, 1, …, n-2)可得：

这样我们可以构造一个以m为未知数的线性方程组，而且在端点条件已知的情况下我们是知道其中几个未知数的值的

端点条件

由i的取值范围可知，共有n-1个公式，但却有n+1个未知量m 。要想求解该方程组，还需另外两个式子。所以需要对两端点x0和xn的微分加些限制。选择不是唯一的，3种比较常用的限制如下。

a. 自由边界(Natural)

首尾两端没有受到任何让它们弯曲的力，即。具体表示为

和

则要求解的方程组可写为：

b. 固定边界(Clamped)

首尾两端点的微分值是被指定的，这里分别定为A和B。则可以推出

将上述两个公式带入方程组，新的方程组左侧为

c. 非节点边界(Not-A-Knot)

指定样条曲线的三次微分匹配，即

根据和，则上述条件变为

新的方程组系数矩阵可写为：

求解过程：

1.将方程组(2)与三种端点条件的任何一种联合，解关于M0，M1，…，M n 的线性方程组.

2.将M i(i=0，1，…，n)代入方程组(1)就得到s(x)关于各子区间的表达式.

特别指出，若第2种端点条件取为

M0=M n=0(s″(x0)=s″(x n)=0)，

据此得到的样条插值函数称为自然样条，它在理论上，计算实践上都是很重要的.上面求解三次样条插值的方法称为三弯矩法，是三次样条插值解算方法中最常用的一种。

三、勒让德高斯积分

高斯积分法是精度最高的插值型数值积分，具有2n+12n+1阶精度，并且高斯积分总是稳定。而高斯求积系数，可以由Lagrange多项式插值系数进行积分得到.高斯-勒让德求积公式是构造高精度差值积分的最好方法之一。他是通过让节点和积分系数待定让函数f(x)f(x)以此取i=0,1,2....n i=0,1,2....n次多项式使其尽可能多的能够精确成立来求出积分节点和积分系数。高斯积分的代数精度是

2n?12n?1，而且是最高的。通常运用的是(?1,1)(?1,1)的积分节点和积分系数，其他积分域是通过一定的变换变换到-1到1之间积分。