matlab实习实验报告.docx

实验一 时间序列ARMA学习
一、基本形式

ARMA模型分为以下三种：
1.自回归模型（AR）
如果时间序列 满足

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

以及 E(  ) = 0，则称时间序列  为服从p阶的自回归模型。
AR模型是用过去的观测值和现在的干扰值的线性组合来预测未来值，式子中的p为阶数，β为待定参数，y为一个平稳序列。
自回归模型的平稳条件：
滞后算子多项式

的根均在单位圆外，即φ(B) = 0的根大于1。
2.移动平均模型（MA）
如果时间序列  满足

则称时间序列  为服从q阶移动平均模型；
MA模型是用过去的干扰值和现在的干扰值的线性组合预测未来的值，式中q为阶数，α为待定参数，y为一个平稳序列。
移动平均模型平稳条件：任何条件下都平稳。
3.自回归滑动平均模型（ARMA）
如果时间序列  满足：

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

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

96.7
126.5
150
116.7
39.2
99.5
73.3
28.3
41.7
81.7
71.7
104.3
125.8
166.7
72.5
38.7
66
53.3
15.7
85.5
120.5
71.7
116.7
264.3
142.3
75.5
47.5
130.7
76.2
23.5
66.2
77.3
80.5
92.8
142
171.7
94
73.3
48.8
63.3
35.3
54.2
75
73.3
141.7
122.2
152
101.2
58.3
45.2
60

43.7
107.8
73.3
78
139.2
126.5
109.5
84.5
83.3
77.7
52.8
50
55.8
64.5
78.3
158
148.7
105.5
110.5
118.3
62.7
36.7
63.5
62.7
104.2
81.7
110.5
147.2
125.7
99.7
98.8
66.7
65
21.3
86.7
62.8
83.3
这里只给出部分数据，原表格附在文档最后，此案例中选取最近50多年的数据局进行分析。
1.AR模型求解
首先整理数据，进行平稳性判断，使用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)%²
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('³õʼֵ','Ô¤²âÖµ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('³õʼֵ','Ô¤²âÖµ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)%²自相关
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('³õʼֵ','Ô¤²âÖµMA')
bzc=sum((yp-c).^2)
结果如下：
比较图像可看出，ARMA数据的整体误差小于AR和MA数据的误差，因而ARMA方法的提出是有用的。
再用误差平方和的比较定量的显示模型的差距。
bzc=sum((yp-c).^2)
用matlab计算得出数据：
AR
MA
ARMA
3.1636e+05
7.5764e+05
2.4427e+05
比较数据可得ARMA模型更加精确。
三、改进方法
通过讨论，我们决定对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 keras.layers.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()
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()
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://deeplearning4j.org/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 定义，但数值计算软件如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，…，Mn的线性方程组.
2.将Mi(i=0，1，…，n)代入方程组(1)就得到s(x)关于各子区间的表达式.
特别指出，若第2种端点条件取为
M0=Mn=0(s″(x0)=s″(xn)=0)，
据此得到的样条插值函数称为自然样条，它在理论上，计算实践上都是很重要的.上面求解三次样条插值的方法称为三弯矩法，是三次样条插值解算方法中最常用的一种。 
三、勒让德高斯积分
高斯积分法是精度最高的插值型数值积分，具有2n+12n+1阶精度，并且高斯积分总是稳定。而高斯求积系数，可以由Lagrange多项式插值系数进行积分得到.高斯-勒让德求积公式是构造高精度差值积分的最好方法之一。他是通过让节点和积分系数待定让函数f(x)f(x)以此取i=0,1,2....ni=0,1,2....n次多项式使其尽可能多的能够精确成立来求出积分节点和积分系数。高斯积分的代数精度是2n−12n−1，而且是最高的。通常运用的是(−1,1)(−1,1)的积分节点和积分系数，其他积分域是通过一定的变换变换到-1到1之间积分。

四、问题解决方案
使用三次样条插值的方法，逐次取点进行曲线拟合，构成比较符合的曲线，再用勒让德高斯求积方法求出曲线的长度，下面是三次样条插值的流程图：

C语言求解三次样条的代码如下，默认边界条件为二阶导数=0，只给出主要部分函数（源代码test222.cpp）。
#include
#include
using namespace std;

const int MAX=50;
float x[MAX],y[MAX],h[MAX];//x为横坐标，y为纵坐标，h为步长

float c[MAX],a[MAX],fxym[MAX];
float f(int x1, int x2,int x3)//求差分函数
{
float a=(y[x3]-y[x2])/(x[x3]-x[x2]);
float b=(y[x2]-y[x1])/(x[x2]-x[x1]);
return (a-b)/(x[x3]-x[x1]);
}
void cal_m(int n)//追赶法求解弯矩向量
{
float B[MAX];
B[0]=c[0]/2;
for (int i=1;i=0;i--)
fxym[i]=fxym[i]-B[i]*fxym[i+1];
}
void printout(int n);
用一组测试数据得出的结果，与原结果一致。

混合编程输出四个系数参数：
#include "8spl.h"
#include "mex.h"
void chazhi(double *p)
{
int k, n;
double t;
static double x[11] = { -1.0,-0.95,-0.75,-0.55,-0.3,0.0,0.2,0.45,0.6,0.8,1.0 };
static double y[11] = { 0.0384615,0.0424403,0.06639,0.116788,0.307692,1.0,0.5,0.164948,0.1,0.0588236,0.0384615 };
static double s[4];
k = -1; n = 11;
t = 0.85;
spl(x, y, n, k, t, s);
for (int i = 0; i <4; i++)
{
*(p + i) = *(s+i);
}
}
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
double *p;
double q[2][11];
plhs[0] = mxCreateDoubleMatrix(1,4, mxREAL);
p = mxGetPr(plhs[0]);

chazhi(p);
}

高斯勒让德求积法的C语言实现：
#include "math.h"
#include "stdio.h"
#include "9lrgs.c"//函数头文件
main()
{ double a,b,eps,g,lrgsf(double);
a=0.75; b=0.85; eps=0.000001;
g=lrgs(a,b,eps,lrgsf);
printf("\n");
printf("g=%e\n",g);
printf("\n");
}

double lrgsf(x)
double x;
{ double y;
y=0.390152*pow(30-x,3)-2.9*pow(x-29,2)+2.60985*(30-x)+5.9;//
return(y);
}
#include "math.h"
double lrgs(a,b,eps,f)
double a,b,eps,(*f)();
{ int m,i,j;
double s,p,ep,h,aa,bb,w,x,g;
static double t[5]={-0.9061798459,-0.5384693101,0.0,
0.5384693101,0.9061798459};
static double c[5]={0.2369268851,0.4786286705,0.5688888889,
0.4786286705,0.2369268851};
m=1;
h=b-a; s=fabs(0.001*h);
p=1.0e+35; ep=eps+1.0
while ((ep>=eps)&&(fabs(h)>s))
{ g=0.0;
for (i=1;i<=m;i++)
{ aa=a+(i-1.0)*h; bb=a+i*h;
w=0.0;
for (j=0;j<=4;j++)
{ x=((bb-aa)*t[j]+(bb+aa))/2.0;
w=w+(*f)(x)*c[j];
}
g=g+w;
}
g=g*h/2.0;
ep=fabs(g-p)/(1.0+fabs(g));
p=g; m=m+1; h=(b-a)/m;
}
return(g);
}
用y=0.390152*pow(30-x,3)-2.9*pow(x-29,2)+2.60985*(30-x)+5.9函数，代入C代码中，计算结果与mathematic计算做比较，精度达到小数点后六位，可见算法的准确性。
计算结果为：1755公里
与实际长度1566公里有较大误差，主要是计算过程比较麻烦，用