📄 rbf网络.cpp
字号:
// RBF网络.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#define EP 0.01 //中心变化量预设值
#define X_Num 6 //输入层神经元个数
//#define I 51 //隐含层神经元个数
#define pp 0.85 //分割因子
#define Y_Num 1 //输出层神经元个数
#define T_Num 34 //样本组数
#define STUDY_SPEED 0.01 //学习速率
#define Max 500 //计算中心时，最大迭代次数
int Ix,N_Max;
int No[T_Num];
int NC[T_Num];
double LL[T_Num];
double *w; //输出层权值,包含阀值
double *t; //聚类中心
double *t1; //训练中前一次的聚类中心
double x[X_Num][T_Num]; //训练样本，即输入
double *G; //隐含层输出矩阵，包含阀值系数1
double Y[T_Num][1]; //训练时的网络输出值
double g; //方差——高斯因子
double xx[X_Num][1]; //预测时的输入向量
double *uu; //预测时隐含层输出
double yy[1][1]; //预测输出值
double *GT,*GM,*GD; //GT存放G的转置，GM存放(GT*G)及其逆阵，GD存放(GT*G)的逆阵与GT的乘积
//////////**********用到的矩阵运算子程序**********//////////
void multiply(double *A,double *B,double *M,int m,int k,int n) //求两矩阵A(m行k列)和B(k行n列)相乘,结果为M(m行n列)
{
int i,j,t;
for(i=0;i<m;i++) //M矩阵初始化为全0
{
for(j=0;j<n;j++)
// *(M+i*n+j) = 0.0; //*******也可以这样表示*******
M[i*n+j] = 0.0; //*******这样表示较直观*******
}
for(i=0;i<m;i++)
{
for(j=0;j<n;j++)
{
for(t=0;t<k;t++)
M[i*n+j] = M[i*n+j]+A[i*k+t]*B[t*n+j]; //获得M中一个元素的算法
}
}
}
void transpose(double *A,double *B,int m,int n) //求矩阵A(n行m列)的转秩，结果为B(m行n列)
{
int i,j;
for(i=0;i<m;i++)
{
for(j=0;j<n;j++)
B[i*n+j] = A[j*m+i]; //获得B中一个元素的算法
}
}
int reverse(double *a,int n) //全选主元高斯—约当(Gauss—Jordan)消去法求矩阵a的逆阵
{
int *is,*js,i,j,k,l,u,v;
double d,p;
is=(int *)malloc(n*sizeof(int));
js=(int *)malloc(n*sizeof(int));
for (k=0; k<=n-1; k++)
{ d=0.0;
for (i=k; i<=n-1; i++)
for (j=k; j<=n-1; j++)
{ l=i*n+j; p=fabs(a[l]);
if (p>d) { d=p; is[k]=i; js[k]=j;}
}
if (d+1.0==1.0)
{ free(is); free(js); printf("err**not inv\n");
return(0);
}
if (is[k]!=k)
for (j=0; j<=n-1; j++)
{ u=k*n+j; v=is[k]*n+j;
p=a[u]; a[u]=a[v]; a[v]=p;
}
if (js[k]!=k)
for (i=0; i<=n-1; i++)
{ u=i*n+k; v=i*n+js[k];
p=a[u]; a[u]=a[v]; a[v]=p;
}
l=k*n+k;
a[l]=1.0/a[l];
for (j=0; j<=n-1; j++)
if (j!=k)
{ u=k*n+j; a[u]=a[u]*a[l];}
for (i=0; i<=n-1; i++)
if (i!=k)
for (j=0; j<=n-1; j++)
if (j!=k)
{ u=i*n+j;
a[u]=a[u]-a[i*n+k]*a[k*n+j];
}
for (i=0; i<=n-1; i++)
if (i!=k)
{ u=i*n+k; a[u]=-a[u]*a[l];}
}
for (k=n-1; k>=0; k--)
{ if (js[k]!=k)
for (j=0; j<=n-1; j++)
{ u=k*n+j; v=js[k]*n+j;
p=a[u]; a[u]=a[v]; a[v]=p;
}
if (is[k]!=k)
for (i=0; i<=n-1; i++)
{ u=i*n+k; v=i*n+is[k];
p=a[u]; a[u]=a[v]; a[v]=p;
}
}
free(is); free(js);
return(1);
}
void compose(double *data,int *no,int n) //快速排序法
{
int i,j,x,t=0;
double xx;
double *data1;
int *no1;
data1=(double *)malloc(n*sizeof(double));
no1=(int *)malloc(n*sizeof(int));
for(i=0;i<n;i++)
data1[i] = data[i];
for(i=0;i<n;i++)
no1[i] = no[i];
for(j=0;j<n-1;j++)
{
for(i=t+1;i<n;i++)
{
if(data1[i]>data1[t])
{
xx = data1[t];
data1[t] = data1[i];
data1[i] = xx;
x = no1[t];
no1[t] = no1[i];
no1[i] = x;
}
}
t++;
}
for(i=0;i<n;i++)
{
data[i] = data1[n-i-1];
no[i] = no1[n-i-1];
}
free(data1);
free(no1);
}
//////////**********学习中心——K均值聚类算法**********//////////
double interval(double *A,double *B,int n)
{
int i;
double sum=0.0;
for(i=0;i<n;i++)
sum+=pow(fabs(A[i]-B[i]),2);
return (sqrt(sum));
}
////径向距离/////
int min_juli(int m,int *set_number) //m为样本序号，*set_number返回输入向量所在的聚类
{
int k,jj=0;
double sum_max=1000.0;
double sum=0.0;
for(jj=0;jj<Ix;jj++)
{
for (k=0;k<X_Num;k++)
sum+=pow(fabs(x[k][m]-t[k*Ix+jj]),2);
if(sum==0) //若样本即为所选的初始中心，则直接跳出
{
*set_number = jj;
break;
}
else
if(sum<sum_max)
{
*set_number = jj;
sum_max = sum;
}
sum=0;
}
return (*set_number);
}
//////////////k均值算法//////////////////////
void k_means()
{
int i,j;
int set_number=0; //聚类序号
int set_member_number = 0; //某一个聚类中输入向量序号
double min_distance;
double sum;
double max_change=10.0; //对应中心最大的变化量
N_Max = 0; //迭代次数初始化为0
////////******主循环******////////
while (max_change>EP&&N_Max<Max) //循环条件:中心误差未满足要求，并且未达到最大迭代次数(设为500)
{
for(i=0;i<X_Num;i++) //聚类前先记录原来的中心
for(j=0;j<Ix;j++)
t1[i*Ix+j] = t[i*Ix+j];
set_member_number = 0; //每次都是从第 0 组样本开始聚类
while(set_member_number<T_Num)
{
min_juli(set_member_number,&set_number); //返回一个距离最小的聚类序号
for (i=0;i<X_Num;i++) //调整中心
t[i*Ix+set_number]+=STUDY_SPEED*(x[i][set_member_number]-t[i*Ix+set_number]);
set_member_number++;
}
sum = 0.0;
min_distance=0.0;
for(j=0;j<Ix;j++) //计算中心的变化
{
for(i=0;i<X_Num;i++)
sum+=pow(fabs(t[i*Ix+j]-t1[i*Ix+j]),2);
if(sqrt(sum)>=min_distance)
min_distance = sqrt(sum);
sum = 0.0;
}
max_change = min_distance; //得到此次迭代中心的最大变化量
N_Max++;
}
//////////主循环结束////////////
}
//////////**********学习中心结束**********//////////
//////////**********确定方差**********//////////
void getError()
{
int i,j,n;
double dsum=0,dmax=0.0;
for(n=0;n<Ix-1;n++) //求所选中心之间的最大距离
for(j=n+1;j<Ix;j++)
for(i=0;i<X_Num;i++)
dsum+=pow(fabs(t[i*Ix+n]-t[i*Ix+j]),2);
if(dsum>dmax)
dmax=dsum;
dmax = sqrt(dmax);
g = dmax/sqrt(2*Ix); //求方差——高斯因子
}
//////////**********学习权值——用伪逆的方法求解**********//////////
void wStudy()
{
int i,j,n;
double q,sum=0;
for (i=0;i<T_Num;i++) //隐含层输出矩阵G初始化（全置1）
for (j=0;j<Ix+1;j++)
G[i*(Ix+1)+j] = 1.0;
for (j=0;j<Ix;j++) //计算隐含层输出矩阵G
for (i=0;i<T_Num;i++)
{
for (n=0;n<X_Num;n++)
sum+=pow(fabs(x[n][i]-t[n*Ix+j]),2);
q = sum/((-2.0)*g*g);
G[i*(Ix+1)+j] = exp(q);
sum = 0.0;
}
////////********以下用伪逆的方法求解w[][]********////////
transpose(G,GT,Ix+1,T_Num);
multiply(GT,G,GM,Ix+1,T_Num,Ix+1);
i = reverse(GM,Ix+1);
multiply(GM,GT,GD,Ix+1,Ix+1,T_Num);
multiply(GD,*Y,w,Ix+1,T_Num,1);
}
////////********预测输出值********////////
////****计算隐藏层输出——需要输入数据xx[X_Num][1]****////
void compute_uu()
{
int i,j;
double q,sum=0.0;
for(i=0;i<=Ix;i++)
uu[0*Ix+i] = 1.0;
for(j=0;j<Ix;j++)
{
for(i=0;i<X_Num;i++)
sum+=pow(fabs(xx[i][0]-t[i*Ix+j]),2);
q = sum/((-2.0)*g*g);
uu[0*Ix+j] = exp(q);
sum = 0.0;
}
}
void ones(double *A,double *range,int m,int n) //A为待归一化的矩阵（m行n列），range为矩阵每列上下限（2行n列）
{
int i,j;
for(j=0;j<n;j++)
{
for(i=0;i<m;i++)
A[i*n+j] = (A[i*n+j]-range[n+j])/(range[j]-range[n+j]);
}
}
void main()
{
int i,j,n,mm,jj=0,p=2;
FILE *fp;
char file[] = "0000.dat";
if((fp=fopen(file,"a+"))==NULL)
{
printf("can not open this file!\n");
}
double data[X_Num+T_Num];
double XX[308][2]={{30.9, 28.5},
{31.3, 28.6},
{31.1, 27.8},
{31.3, 29.1},
{31.3, 28.1},
{32.0, 29.5},
{30.2, 28.3},
{30.2, 29.0},
{30.8, 28.6},
{30.1, 28.6},
{30.2, 28.6},
{30.1, 28.6},
{30.9, 28.6},
{30.6, 29.1},
{31.0, 29.3},
{30.1, 29.1},
{30.5, 28.9},
{30.3, 28.7},
{30.5, 28.8},
{30.5, 28.8},
{32.3, 27.0},
{30.3, 25.2},
{28.9, 25.1},
{27.6, 25.2},
{26.9, 25.2},
{26.9, 24.8},
{34.7, 35.7},
{37.6, 38.0},
{41.4, 40.9},
{45.4, 46.9},
{49.2, 49.5},
{52.2, 52.6},
{55.9, 55.9},
{58.5, 59.6},
{61.3, 61.9},
{64.8, 65.8},
{68.4, 67.4},
{70.4, 70.1},
{73.4, 74.6},
{75.4, 76.2},
{77.4, 79.5},
{80.0, 79.8},
{81.7, 82.5},
{84.5, 84.9},
{85.8, 86.8},
{87.0, 90.0},
{87.4, 86.6},
{86.0, 84.5},
{84.1, 83.7},
{84.1, 82.7},
{82.4, 81.8},
{81.9, 79.1},
{81.0, 78.3},
{78.3, 79.3},
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -