📄 bayes.cpp
字号:
#include <iostream.h>
#include <math.h>
#include <fstream.h>
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <time.h>
#define NN 27//
#define N1 11//
#define N2 16//
#define predict_sample 24//
#define dd 11//
#define N 6//
//global variables
//int N= //样本数
//int d= //维数
char shuru[100];//样本数据存放的文件名及其路径
char parameter[100];//
char result[100];//存放对预测样本的函数输出值，以及预测精度结果等的文件名及其存放路径
double t0=14.5,grad=3.0,D1=1954.6,srm18=3.5,h1=4.2,ri1=35.0,ssp=-62.0;
double rw75=0.911,gmin=40.0,gmax=140.0,gcur=3.7,xm=1.61,xn=1.73;
double rsh=8.0,xx=7.779,fep0=2.43,A=1.04;
int G=2,swop=3,sg=1,mn[2]={N1,N2};
double oil[N1][dd],water[N2][dd],unk[predict_sample][dd];
void spcorr(double,double,double *,double*);
void can(double *,double *);
void sum1(double *,int,double *);
void sum2(double *,int,double *,double *);
void guige(double *,int,double *,double *,double *);
void ave(int,double *,double *,double *);
void def(int,double *,double *,double *);
void ivsnc(double *);
void pm(double *,double *,double *,double *,double *,int *);
void setX();
/*G表示样本类别数*/
/*N表示规格化后得样本属性个数*/
/*swop＝3表示计算sw时所用公式为第三个*/
/*mn*/
/***************************主函数********************************/
void main()
{
double m1,u,dep,fep,h,d;
double v[N],zz[N],x11[N1][N],x22[N2][N],v1[N],v2[N],u1[N],u2[N];
double uu[N],x1[N1][N],x2[N2][N],t[N][N],aa[N],b1[N1][N];
double w[6][6],s[6][6],b2[N2][N],co[2],c[2][6],ff[2],pmax,z[6];
double tave[6],a[2][6],sum;
double z1[100][6],z2[100][6];
int i,j,k,m,l,i1,ii,jj,ie,ip,IG,ih,ni;
/************第一部分****对样本数据进行处理求出判别公式***********/
setX();//从文件中将数据读入数组
cout<<endl<<endl;
cout<<" ********************对样本数据进行处理***********************"<<endl;
cout<<endl;
cout<<"原始样本数据如下:"<<endl;
cout<<endl;
cout<<"G NO1 VSH POR SO FRB FRC FRE"<<endl<<endl;
for(k=0;k<N;k++)
{
v1[k]=0.0;
u1[k]=0.0;
v2[k]=0.0;
u2[k]=0.0;
}
m=0;
for(i=0;i<G;i++)
{
ni=mn[i];
m=m+ni;
if(i==0)
{
for(j=0;j<ni;j++)
{
printf("%2d",i+1);
printf("%4d",j+1);
can(oil[j],zz);
for(k=0;k<N;k++)
{ x11[j][k]=zz[k];
printf("%10.3f",x11[j][k]);
}
printf("\n");
}
sum1(*x11,ni,v1);
}
printf("\n");
if(i==1)
{
for(j=0;j<ni;j++)
{
printf("%2d",i+1);
printf("%4d",j+1);
can(water[j],zz);
for(k=0;k<N;k++)
{
x22[j][k]=zz[k];
printf("%10.3f",x22[j][k]);
}
printf("\n");
}
sum1(*x22,ni,v2);
}
}
for(k=0;k<N;k++)
{
v[k]=(v1[k]+v2[k])/m;
}
for(i=0;i<G;i++)
{
ni=mn[i];
if(i==0)
{
sum2(*x11,ni,v,u1);
}
if(i==1)
{
sum2(*x22,ni,v,u2);
}
}
m1=1;
for(k=0;k<N;k++)
{
uu[k]=(u1[k]+u2[k])/(m-m1);
}
cout<<endl<<endl;
cout<<"规格化数据数据如下："<<endl<<endl;
cout<<"G NO1 VSH POR SO FRB FRC FRE"<<endl<<endl;
for(i=0;i<G;i++)
{
ni=mn[i];
if(i==0)
{
guige(*x11,ni,v,uu,*x1);
for(j=0;j<ni;j++)
{
printf("%2d",i+1);
printf("%4d",j+1);
for(k=0;k<N;k++)
printf("%10.3f",x1[j][k]);
printf("\n");
}
printf("\n");
}
if(i==1)
{
guige(*x22,ni,v,uu,*x2);
for(j=0;j<ni;j++)
{
printf("%2d",i+1);
printf("%4d",j+1);
for(k=0;k<N;k++)
printf("%10.3f",x2[j][k]);
printf("\n");
}
}
}
for(i=0;i<N;i++)
for(j=0;j<N;j++)
t[i][j]=0.0;
for(i=0;i<G;i++)
{
ni=mn[i];
if(i==0)
{
ave(ni,*x1,aa,*b1);
def(ni,*b1,*w,*s);
}
if(i==1)
{
ave(ni,*x2,aa,*b2);
def(ni,*b2,*w,*s);
}
for(j=0;j<N;j++)
a[i][j]=aa[j];
for(j=0;j<N;j++)
for(k=0;k<N;k++)
t[j][k]=t[j][k]+s[j][k];
}
for(j=0;j<N;j++)
{
sum=0.0;
m=0;
m=m+mn[i];
for(i=0;i<G;i++)
{
sum=sum+mn[i]*a[i][j];
m=m+mn[i];
}
tave[j]=sum/m;
}
cout<<endl;
cout<<"均值为:"<<endl<<endl;
for(i=0;i<G;i++)
{
for(j=0;j<N;j++)
printf("%12.4f",a[i][j]);
printf("\n");
}
cout<<endl;
cout<<"总均值为:"<<endl<<endl;
for(j=0;j<N;j++)
printf("%13.4e",tave[j]);
cout<<endl<<endl;
cout<<"协方差矩阵为:"<<endl<<endl;
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
printf("%12.5f",t[i][j]);
printf("\n");
}
cout<<endl;
cout<<"协方差逆矩阵为:"<<endl<<endl;
ivsnc(*t);
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
printf("%12.5f",t[i][j]);
printf("\n");
}
for(i=0;i<G;i++)
{
sum=0.0;
for(k=0;k<N;k++)
{
c[i][k]=0.0;
for(l=0;l<N;l++)
{
c[i][k]=c[i][k]+t[k][l]*a[i][l];
sum=sum+t[l][k]*a[i][k]*a[i][l];
}
}
co[i]=(-0.5)*sum;
}
cout<<endl;
cout<<"判别系数为:"<<endl<<endl;
for(i=0;i<G;i++)
{
for(k=0;k<N;k++)
{
printf("%11.5f",c[i][k]);
}
printf("%11.5f",co[i]);
printf("\n");
}
printf("\n");
cout<<"判别结果:"<<endl<<endl;
cout<<" G NO1 VSH POR SO FRB FRC FRE PMAX NO2"<<endl;
for(i=0;i<G;i++)
{
ni=mn[i];
for(ii=0;ii<ni;ii++)
{
for(jj=0;jj<N;jj++)
{
if(i==0)z[jj]=x1[ii][jj];
if(i==1)z[jj]=x2[ii][jj];
}
pm(*c,co,z,ff,&pmax,&i1);
printf("%3d%4d",i+1,ii+1);
for(jj=0;jj<N;jj++)
printf("%8.3f ",z[jj]);
printf("%8.4f%3d\n",pmax,i1+1);
}
printf("\n");
}
u=0.0;
for(k=0;k<N;k++)
{
h=0.0;
for(l=0;l<N;l++)
{
d=0.0;
for(i=0;i<G;i++)
d=d+ni*a[i][k]-tave[k]*(a[i][l]-tave[l]);
h=h+t[k][l]*d;
}
u=u+h;
}
cout<<" 显著性检验结果:"<<endl<<endl;
if(u<=xx)cout<<" 组内差异小"<<endl;
else cout<<" 组间差异不显著"<<endl;
IG=G-1;
for(ie=0;ie<IG;ie++)
{
ih=ie+1;
for(ip=ih-1;ip<G;ip++)
{
for(k=0;k<N;k++)
{
dep=(c[k][ie]-c[k][ip])*(a[ie][k]-a[ip][k]);
}
fep=((m-G-N+1)*mn[ie]*mn[ip]*dep*dep)/((m-G)*N*mn[ie]+mn[ip]);
if(fep<=fep0)
{
cout<<" 组间差异不显著"<<endl;
goto end1;
}
}
}
cout<<" 组间差异显著"<<endl;
end1:cout<<endl;
cout<<" *********************对未知样本进行预测**************************"<<endl<<endl;
cout<<" 提示:输入样本以输入深度(D)为零表示输入结束"<<endl<<endl;
cout<<" 按顺序依次输入 D H RM18 DM RLLD RLLS GR AC RXO RI SP"<<endl<<endl;
cout<<" 以空格分开各个样本属性值,回车然后输入下一个未知样本"<<endl<<endl;
cout<<" 请输入未知样本:\n"<<endl;
for(i=0;i<predict_sample;i++)
{
can(unk[i],zz);
for(j=0;j<N;j++)
{
z1[i][j]=zz[j];
}
}
guige(*z1,predict_sample,v,uu,*z2);
for(i=0;i<predict_sample;i++)
{
printf("%3d",i+1);
for(j=0;j<N;j++)
{
z[j]=z2[i][j];
printf("%10.4f",z[j]);
}
pm(*c,co,z,ff,&pmax,&i1);
printf("%10.4f%4d",pmax,i1+1);
printf("\n");
}
printf("预测结束\n");
}
void spcorr(double h,double ri,double *rm,double *yy)
{
double y[5],x[5]={5.0,20.0,50.0,100.0,200.0},p[4],q[3],r[2];
double x1,s;
int i=0;
x1=ri/(*rm);
h=h/0.3048;
h=log10(h);
y[0]=0.7287-1.6228*h+1.2741*h*h-0.3564*pow(h,3);
y[1]=1.2735-2.8651*h+2.2902*h*h-0.6383*pow(h,3);
y[2]=1.4156-2.6646*h+1.7730*h*h-0.4146*pow(h,3);
y[3]=1.9524-3.7751*h+2.6247*h*h-0.6449*pow(h,3);
y[4]=3.0503-6.5901*h+5.1026*h*h-1.3703*pow(h,3);
for(i=0;i<5;i++)
y[i]=pow(10.0,y[i]);
for(i=0;i<4;i++)
p[i]=(y[i+1]-y[i])/(x[i+1]-x[i]);
for(i=0;i<3;i++)
q[i]=(p[i+1]-p[i])/(x[i+2]-x[i]);
for(i=0;i<2;i++)
r[i]=(q[i+1]-q[i])/(x[i+3]-x[i]);
s=(r[1]-r[0])/(x[4]-x[0]);
*yy=y[0]+p[0]*(x1-x[0])+q[0]*(x1-x[0])*(x1-x[1])+r[0]*(x1-x[0])*(x1-x[1])*(x1-x[2])+s*(x1-x[0])*(x1-x[1])*(x1-x[2])*(x1-x[3]);
return;
}
void can(double *shu,double *zz)
{
double d,h,rm18,dm,rlld,rlls,gr,ac,rxo,ri,sp,t,psp,k,c,rw,rwn,rmn,rmfn;
double rmfen,rwen,rmf,rt,sh1,sh,sw,so,sxo,fre,frc,frb,por,alfa,rm,yy;
int ki;
d=shu[0];h=shu[1];rm18=shu[2];dm=shu[3];rlld=shu[4];rlls=shu[5];
gr=shu[6];ac=shu[7];rxo=shu[8];ri=shu[9];sp=shu[10];
t=t0+grad*d/100;
rm=71.4*rm18/(1.8*t+39.0);
spcorr(h,ri,&rm,&yy);
psp=yy*sp;
/*计算rw*/
t=1.8*t+32.0;
k=60.0+0.133*t;
c=1.78363+4.94273*0.1*dm-1.76949*pow(dm,2)+5.89192*0.1*pow(dm,3);
if(rw75!=0.0)rwn=rw75;
else
{
rmn=71.4*rm18/82.0;
rmfn=c*pow(rmn,1.07);
if(rmfn>0.1)rmfen=0.85*rmfn;
else rmfen=(146.0*rmfn-5.0)/(337.0*rmfn+77.0);
rwen=rmfen/pow(10,(-psp/k));
if(rwen>0.12)rwn=-0.58+pow(10.0,(0.69*rwen-0.24));
else rwn=(77.0*rwen+5.0)/(146.0-337.0*rwen);
}
rw=82.0*rwn/(t+7.0);
rmf=c*pow(rm,1.07);
/*计算rt*/
rt=(rlld-0.29*rxo)*pow((rlld/rlls),0.1)/0.71;
if((rlld<=rlls)||(rlld<=rxo))rt=rlld;
/*计算vsh*/
sh1=(gr-gmin)/(gmax-gmin);
if(sh1>1.0)sh1=1.0;
if(sh1<0.0)sh1=0.0;
sh=(pow(2.0,(gcur*sh1))-1.0)/(pow(2.0,gcur)-1.0);
/*计算por*/
ac=ac/0.3048;
if(sg!=1)por=1.64*pow(10,-3.0)*ac-0.232;
else por=1.64*pow(10,-3.0)*ac-0.35*sh-0.202;
if(por<0.001)por=0.001;
/*计算sw,sxo*/
ki=swop;
if(ki==1)
{
sw=(sqrt(0.81*rw/rt)-sh*(rw/(0.4*rsh)))/por;
sxo=(1.0/por)*(sqrt(0.81*rmf/rxo)-sh*(rmf/(0.4*rsh)));
}
if(ki==2)
{
xm=1.87+0.019/por;
if(por>0.1)xm=2.1;
if(xm>4.0)xm=4.0;
A=1.0;
sw=pow(A*rw/(pow(por,xm)*rt),(1.0/xn));
sxo=pow(A*rmf/(pow(por,xm)*rxo),1.0/xn);
}
if(ki==3)
{
sw=pow((A*rw/(pow(por,xm)*rt)),(1.0/xn));
sxo=pow((A*rmf)/(pow(por,xm)*rxo),(1.0/xn));
}
if(sw>0.999)sw=0.999;
if(sw<0.001)sw=0.001;
if(sxo>0.8)sxo=sw;
if(sxo<sw)sxo=sw;
fre=sw/sxo;
so=1.0-sw;
frc=rxo/rt;
alfa=psp/ssp;
frb=alfa*rt/rw;
zz[0]=sh;
zz[1]=por;
zz[2]=so;
zz[3]=frb;
zz[4]=frc;
zz[5]=fre;
return;
}
void sum1(double *x,int ni,double *v)
{
double sum;
int k,j;
for(k=0;k<N;k++)
{
sum=0.0;
for(j=0;j<ni;j++)
sum=sum+*(x+j*N+k);
v[k]=sum;
}
return;
}
void sum2(double *x,int ni,double *v,double *u)
{
int j,k;
double sum;
for(k=0;k<N;k++)
{
sum=0.0;
for(j=0;j<ni;j++)
{
sum=sum+pow((*(x+j*6+k)-v[k]),2.0);
}
u[k]=sum;
if(u[k]<0.000001)u[k]=0.000001;
u[k]=sqrt(u[k]);
}
return;
}
void guige(double *x,int ni,double *v,double *uu,double *y)
{
int j,k;
for(j=0;j<ni;j++)
for(k=0;k<N;k++)
*(y+k+j*6)=(*(x+k+j*6)-v[k])/uu[k];
return;
}
void ave(int ni,double *x,double *aa,double *b)
{
int j,k;
double sum;
for(j=0;j<N;j++)
{
sum=0.0;
for(k=0;k<ni;k++)
sum=sum+*(x+6*k+j);
aa[j]=sum/ni;
for(k=0;k<ni;k++)
*(b+6*k+j)=*(x+6*k+j)-aa[j];
}
return;
}
void def(int ni,double *b,double *w,double *s)
{
int k,l,j;
double sum;
for(k=0;k<N;k++)
for(l=0;l<N;l++)
{
sum=0.0;
for(j=0;j<ni;j++)
sum=sum+(*(b+k+6*j))*(*(b+l+6*j));
*(w+l+6*k)=sum;
*(s+l+6*k)=*(w+l+6*k)/(ni-1);
}
return;
}
void ivsnc(double *t)
{
int i,j,k;
double tkk,tit;
for(k=0;k<N;k++)
{
tkk=*(t+k+N*k);
for(j=0;j<N;j++)
{
*(t+j+N*k)=(*(t+j+N*k))/tkk;
}
*(t+k+N*k)=1.0/tkk;
for(i=0;i<N;i++)
{
if(i==k)continue;
tit=*(t+N*i+k);
for(j=0;j<N;j++)
{
*(t+j+N*i)=*(t+j+N*i)-tit*(*(t+j+N*k));
}
*(t+k+N*i)=(-1.0)*tit*(*(t+k+N*k));
}
}
return;
}
void pm(double *c,double *co,double *z,double *ff,double *pmax,int *i1)
{
double fma,sum,h,p,q;
int i,k,j;
fma=0.0;
for(i=0;i<G;i++)
{
sum=0.0;
for(k=0;k<N;k++)
{
sum=sum+(*(c+k+N*i))*z[k];
}
ff[i]=sum+co[i];
if(fma>=ff[i])continue;
fma=ff[i];
}
*pmax=0.0;
h=0.0;
for(j=0;j<G;j++)
{
h=h+exp(ff[j]-fma);
}
for(i=0;i<G;i++)
{
q=exp(ff[i]-fma);
p=q/h;
if((*pmax)>=p)continue;
*pmax=p;
*i1=i;
}
return;
}
//计算样本X[]
void setX()
{
char shuru[100]={"e:\\data\\total.txt"};
cout<<"请输入样本数据的文件名及其路径"<<endl;
FILE *fp=fopen(shuru,"r");
if(fp == NULL)
{
cout<<"can't open the file "<<shuru<<endl;
exit(1);
}
int i1=0,i2=0,j,d1;
float value;
for(int i=0;i<NN;i++)
{
j=0;
int c=fgetc(fp);
if(c==-1)break;
if(char(c)=='+')
{
c=fgetc(fp);
while(1)
{
do
{
c=fgetc(fp);
if(c=='\n')
{ i1++;
goto out2;
}
}while(isspace(c));
ungetc(c,fp);
fscanf(fp,"%d:%f",&d1,&value);
oil[i1][j]=value;
j++;
}
}
else if(char(c)=='-')
{
c=fgetc(fp);
while(1)
{
do
{
c=fgetc(fp);
if(c=='\n')
{
i2++;
goto out2;
}
}while(isspace(c));
ungetc(c,fp);
fscanf(fp,"%d:%f",&d1,&value);
water[i2][j]=value;
j++;
}
}
out2:continue;
}
for(i=0;i<predict_sample;i++)
{
j=0;
int c=fgetc(fp);
if(c==-1)break;
c=fgetc(fp);
while(1)
{
do
{
c=fgetc(fp);
if(c=='\n')
{ i1++;
goto out3;
}
}while(isspace(c));
ungetc(c,fp);
fscanf(fp,"%d:%f",&d1,&value);
unk[i][j]=value;
j++;
}
out3:continue;
}
fclose(fp);
for(i=0;i<predict_sample;i++)
{
for(j=0;j<dd;j++)
{
cout<<unk[i][j]<<" ";
}
cout<<endl;
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -