该算法需利用高斯随机分布来制造预测,因此另封装一个randn()函数制造伪高斯随机数,然后在Forecast()函数中对样本进行网络训练,并利用训练结果进行网络预测,
#gauss.h/////////////////////////////////////////////////////////////////////
#define pi 3.1415926535897
#define rd (rand()/(RAND_MAX+1.0))
double randn(int type)
{
srand((unsigned)time(NULL));
//按照12个均匀分布之和减去6得到正态分布函数的x值
if (type==1)
return rd+rd+rd+rd+rd+rd+rd+rd+rd+rd+rd+rd-6.0;
//按照计算公式y=sqrt(-2*ln(U))*cos(2*PI*V)计算得到x
else if(type==2)
return sqrt(-2*log(rand()/(RAND_MAX+1.0)))*cos(2*pi*rand()/(RAND_MAX+1.0));
else
return randn(0.0,1.0,-10.0,10.0);
}
///////////////////////////////////////////////////////////////////////
//////////预测算法///////////////////////////////////////////////////////
vector< vector< double > > Input;////////////训练输入样本
vector< vector< double > > Output;//////////训练输出样本
vector< vector< double > > PreInput;////////预测输入样本
vector< vector< double > > PreOutput; ///////预测输出结果
int input_col;//定义训练输入样本列,数据计算以列为单位,与MATLAB中列向量同意义
int input_row;//定义训练输入样本行
int inputtest_row;//定义训练输出样本行
int output_col;//定义预测输入样本列,即需要预测的数据,1列
int output_row;//定义预测输入样本行
int outputtest_row;//定义预测输出样本行,1列
int maxgen; //迭代次数
int n; //隐形节点个数
double Input[20][1];
///////////////以下为个训练向量与预测向量结构初始化,预测时需进行数据有效化填充////////////////////////////////
int i=0;
Input.resize(input_row);
for(i=0;i<input_row;i++)
{
Input[i].resize(input_col);
}
Output.resize(inputtest_row);
for(i=0;i<inputtest_row;i++)
{
Output[i].resize(input_col);
}
PreInput.resize(output_row);
for(i=0;i<output_row;i++)
{
PreInput[i].resize(output_col);
}
PreOutput.resize(outputtest_row);
for(i=0;i<outputtest_row;i++)
{
PreOutput[i].resize(output_col);
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////本人为动态改变预测体,采用vector向量结构计算,可随时改变数组结构///////////////////////////
/////////vector的使用需添加标准库头文件<vertor>,并定义命名空间,using namespace std;///
///////////////规避double数组只能静态定义的特性///////////////////////////////////////////////////////////////////////////////////
///////////调试中发现double二维及以上常规数组无法进行指针与数组名之间的转化/////////////////////
/////////////具体错误cannot convert double** to double[][]////////////////////////////////////////
void Forecast()
{
int M=input_col;%输入节点个数
int N=output_col;%输出节点个数
double lr1=0.01; //%学习概率
double lr2=0.001; //%学习概率
int i,j,k,kk,kkk;
srand((unsigned)time( NULL ));
double value=0.0;
CString str = _T("");
vector< vector< double > > Wjk(n);
for(i=0;i<n;i++)
{
Wjk[i].resize(input_col);
}
vector< vector< double > > Wjk_1(n);
for(i=0;i<n;i++)
{
Wjk_1[i].resize(input_col);
}
vector< vector< double > > Wjk_2(n);
for(i=0;i<n;i++)
{
Wjk_2[i].resize(input_col);
}
for(i = 0;i <n;i++)
{
for(j = 0;j < input_col;j++)
{
Wjk[i][j] = randn(2);
Wjk_1[i][j] = Wjk[i][j];
Wjk_2[i][j] = Wjk_1[i][j];
}
}
vector< vector< double > > Wij(output_col);
for(i=0;i<output_col;i++)
{
Wij[i].resize(n);
}
vector< vector< double > > Wij_1(output_col);
for(i=0;i<output_col;i++)
{
Wij_1[i].resize(n);
}
vector< vector< double > > Wij_2(output_col);
for(i=0;i<output_col;i++)
{
Wij_2[i].resize(n);
}
for(i = 0;i <output_col;i++)
{
for(j = 0;j < n;j++)
{
Wij[i][j] = randn(2);
Wij_1[i][j] = Wij[i][j];
Wij_2[i][j] = Wij_1[i][j];
}
}
vector< vector< double > > a(1);
for(i=0;i<1;i++)
{
a[i].resize(n);
}
vector< vector< double > > a_1(1);
for(i=0;i<1;i++)
{
a_1[i].resize(n);
}
vector< vector< double > > a_2(1);
for(i=0;i<1;i++)
{
a_2[i].resize(n);
}
for(i = 0;i <1;i++)
{
for(j = 0;j < n;j++)
{
a[i][j] = randn(2);
a_1[i][j] = a[i][j];
a_2[i][j] = a_1[i][j];
}
}
vector< vector< double > > b(1);
for(i=0;i<1;i++)
{
b[i].resize(n);
}
vector< vector< double > > b_1(1);
for(i=0;i<1;i++)
{
b_1[i].resize(n);
}
vector< vector< double > > b_2(1);
for(i=0;i<1;i++)
{
b_2[i].resize(n);
}
for(i = 0;i <1;i++)
{
for(j = 0;j < n;j++)
{
b[i][j] = randn(2);
b_1[i][j] = b[i][j];
b_2[i][j] = b_1[i][j];
}
}
vector< vector< double > > y(1);
for(i=0;i<1;i++)
{
y[i].resize(output_col);
}
vector< vector< double > > net(1);
for(i=0;i<1;i++)
{
net[i].resize(n);
}
vector< vector< double > > net_ab(1);
for(i=0;i<1;i++)
{
net_ab[i].resize(n);
}
vector< vector< double > > d_Wjk(n);
for(i=0;i<n;i++)
{
d_Wjk[i].resize(input_col);
}
vector< vector< double > > d_Wij(output_col);
for(i=0;i<output_col;i++)
{
d_Wij[i].resize(n);
}
vector< vector< double > > d_a(1);
for(i=0;i<1;i++)
{
d_a[i].resize(n);
}
vector< vector< double > > d_b(1);
for(i=0;i<1;i++)
{
d_b[i].resize(n);
}
//%% 输入输出数据归一化
double tempmax = 0.0,tempmin = 0.0;
vector< double > input_max(input_col);
vector< double > input_min(input_col);
vector< vector< double > > Input1(input_row);
for(i=0;i<input_row;i++)
{
Input1[i].resize(input_col);
}
for(i = 0;i <input_col;i++)
{
tempmax=Input[0][i];
tempmin=Input[0][i];
for(j = 0;j < input_row;j++)
{
if(Input[j][i] > tempmax)
tempmax = Input[j][i];
if(Input[j][i] < tempmin)
tempmin = Input[j][i];
}
input_max[i] = tempmax;
input_min[i] = tempmin;
for(j = 0;j < input_row;j++)
{
if(input_max[i] == input_min[i])
Input1[j][i]=0;
else
Input1[j][i] = 2*(Input[j][i] - input_min[i])/(input_max[i] - input_min[i])-1;
}
}
vector< double > output_max(output_col);
vector< double > output_min(output_col);
vector< vector< double > > PreInput1(output_row);
for(i=0;i<output_row;i++)
{
PreInput1[i].resize(output_col);
}
for(i = 0;i <output_col;i++)
{
tempmax=PreInput[0][i];
tempmin=PreInput[0][i];
for(j = 0;j < output_row;j++)
{
if(PreInput[j][i] > tempmax)
tempmax = PreInput[j][i];
if(PreInput[j][i] < tempmin)
tempmin = PreInput[j][i];
}
output_max[i] = tempmax;
output_min[i] = tempmin;
for(j = 0;j < output_row;j++)
{
if(output_max[i] == output_min[i])
PreInput1[j][i] = 0;
else
PreInput1[j][i] = 2*(PreInput[j][i] - output_min[i])/(output_max[i] - output_min[i])-1;
}
}
//%% 网络训练
vector< double > x(input_col);
vector< double > yqw(output_col);
vector< double > error(maxgen);
double temp = 0.0;
for(i = 0; i < maxgen; i++)
{
//%误差累计
error[i]=0.0;
for(kk = 0;kk <input_row;kk++)
{
for(kkk = 0;kkk <input_col ; kkk++)
{
x[kkk] = Input1[kk][kkk];
}
for(kkk = 0;kkk <output_col ; kkk++)
{
yqw[kkk] = PreInput1[kk][kkk];
}
for(j = 0;j <n; j++ )
{
for(k = 0;k<input_col ; k++)
{
net[0][j]=net[0][j]+Wjk[j][k]*x[k];
net_ab[0][j]=(net[0][j]-b[0][j])/a[0][j];
}
temp = mymorlet(net_ab[0][j]);
for(k = 0;k<output_col ; k++)
{
y[0][k] = y[0][k] + Wij[k][j]*temp;
}
}
for(j = 0;j <output_col ;j++)
{
temp = temp+abs(yqw[j] - y[0][j]);
}
error[i] = error[i]+temp;
for(j = 0;j < n;j++)
{
temp = mymorlet(net_ab[0][j]);
for(k = 0;k < output_col;k++)
{
d_Wij[k][j] = d_Wij[k][j] - (yqw[k] - y[0][k])*temp;
}
temp = d_mymorlet(net_ab[0][j]);
for(k = 0;k< input_col ;k++)
{
for(kkk = 0; kkk< output_col ;kkk++)
{
d_Wjk[j][k] = d_Wjk[j][k] + (yqw[kkk] - y[0][kkk])*Wij[kkk][j];
}
d_Wjk[j][k] = -d_Wjk[j][k]*temp*x[k]/a[0][j];
}
for(k = 0;k <output_col ;k++)
{
d_b[0][j] = d_b[0][j] + (yqw[k] - y[0][k])*Wij[k][j];
}
d_b[0][j] = d_b[0][j]*temp/a[0][j];
for(k = 0; k <output_col; k++)
{
d_a[0][j] = d_a[0][j] + (yqw[k] - y[0][k])*Wij[k][j];
}
d_a[0][j] = d_a[0][j]*temp*((net[0][j] - b[0][j])/b[0][j])/a[0][j];
}
///权值参数更新
for(j = 0;j <n;j++)
{
for(k = 0;k < input_col;k++)
{
Wjk[j][k] = Wjk[j][k] - lr1*d_Wjk[j][k];
d_Wjk[j][k] = 0.0;
Wjk_1[j][k] =Wjk[j][k];
Wjk_2[j][k] =Wjk_1[j][k];
}
b[0][j] = b[0][j] - lr2*d_b[0][j];
a[0][j] = a[0][j] - lr2*d_a[0][j];
d_a[0][j] = 0.0;
d_b[0][j] = 0.0;
net[0][j] = 0.0;
net_ab[0][j] = 0.0;
a_1[0][j] = a[0][j];
a_2[0][j] = a_1[0][j];
b_1[0][j] = b[0][j];
b_1[0][j] = b_1[0][j];
}
for(j = 0;j <output_col;j++)
{
for(k = 0;k < n;k++)
{
Wij[j][k] = Wij[j][k] - lr1*d_Wij[j][k];
d_Wij[j][k] = 0.0;
Wij_1[j][k] =Wij[j][k];
Wij_2[j][k] =Wij_1[j][k];
}
y[0][j] = 0.0;
}
}
}
////// 网络预测
//////预测输入归一化
vector< double > inputtest_max(input_col);
vector< double > inputtest_min(input_col);
vector< vector< double > > OutputAsInput1(inputtest_row);
for(i=0;i<inputtest_row;i++)
{
OutputAsInput1[i].resize(input_col);
}
//double inputtest_min[input_col] = {{0.0}};
for(i = 0;i <input_col;i++)
{
tempmax=Output[0][i];
tempmin=Output[0][i];
for(j = 0;j < inputtest_row;j++)
{
if(Output[j][i] > tempmax)
tempmax = Output[j][i];
if(Output[j][i] < tempmin)
tempmin = Output[j][i];
}
inputtest_max[i] = tempmax;
inputtest_min[i] = tempmin;
for(j = 0;j < inputtest_row;j++)
{
if(inputtest_max[i] == inputtest_min[i])
OutputAsInput1[j][i] = 0.5*(input_max[i] + input_min[i]);
else
OutputAsInput1[j][i] = (Output[j][i] - inputtest_min[i])/(inputtest_max[i] - inputtest_min[i])*(input_max[i] - input_min[i]) + input_min[i];
}
}
/////网络预测
vector< double > x_test(input_col);
vector< vector< double > > yuce(1);
for(i=0;i<1;i++)
{
yuce[i].resize(inputtest_row);
}
for(i = 0; i <inputtest_row; i++)
{
for(j = 0; j <input_col; j++)
{
x_test[j] = OutputAsInput1[i][j];
}
for(j = 0; j <n; j++)
{
for(k = 0; k < input_col; k++)
{
net[0][j] = net[0][j] + Wjk[j][k]*x_test[k];
net_ab[0][j] = (net[0][j] - b[0][j])/a[0][j];
}
temp = mymorlet(net_ab[0][j]);
for(k = 0; k < output_col; k++)
{
y[0][k] = y[0][k] + Wij[k][j]*temp;
}
}
yuce[0][i] = y[0][k-1];
for(j = 0;j <output_col;j++)
{
y[0][j] = 0.0;
}
for(j = 0;j < n;j++)
{
net[0][j] = 0.0;
net_ab[0][j] = 0.0;
}
}
vector< double > yun(outputtest_row);
for(i = 0; i < outputtest_row; i ++)
{
yun[i] = (yuce[0][i]+1)/2*(output_max[0] - output_min[0]) + output_min[0];
PreOutput[i][output_col-1]=yun[i];
}
return;
}