📄 assignt.cpp
字号:
/****************************************************************************
* Computer Excersice Section 2.5 Problem 2 *
* Student: Lieni Xie ID: 4052420 *
* source code file : lieni.cpp *
* output file: lienioutput.txt *
* purpose: 1)Design adichotomizer for two categories *
* 2)Determine the empirical training error *
* 3)Use the Bhattacharyya bound to bound the error *
*****************************************************************************/
#include "stdio.h"
#include "math.h"
#include <stdlib.h>
void main()
{
struct Sample1
{ long double x[3];
};
struct Sigema1
{ long double s[3][3];
};
Sample1 samp_w1[10],samp_w2[10];
Sigema1 sigma_w1[10],sigma_w2[10];
static long double x[3];
int i,j,k,iCorrect=0,iWrong=0;
static long double u1[3],u2[3],x_u1[3],x_u2[3],seta1[3][3],seta2[3][3];
long double seta_w1_hl,seta_w2_hl;
long double inver_seta1[3][3],inver_seta2[3][3];
long double a11,a12,a13,a21,a22,a23,a31,a32,a33,b11,b12,b13,b21,b22,b23,b31,b32,b33;
long double x_uw1[3],x_uw2[3],dmatrixw1[3],dmatrixw2[3];
long double mma1,mma2,gxw1,gxw2;
long double dLogsigma1, dLogsigma2;
long double setaplus[3][3],setaplus_hl,inver_setaplus[3][3],JB1;
long double JB,JBln,JBUT,u2_u1[3],Pe;
FILE *fp;
long double hanglieshi(long double a[3][3]);
// Data of Samples
samp_w1[0].x[0]=-5.01; samp_w1[0].x[1]=-8.12; samp_w1[0].x[2]=-3.68;
samp_w1[1].x[0]=-5.43; samp_w1[1].x[1]=-3.48; samp_w1[1].x[2]=-3.54;
samp_w1[2].x[0]=1.08; samp_w1[2].x[1]=-5.52; samp_w1[2].x[2]=1.66;
samp_w1[3].x[0]=0.86; samp_w1[3].x[1]=-3.78; samp_w1[3].x[2]=-4.11;
samp_w1[4].x[0]=-2.67; samp_w1[4].x[1]=0.63; samp_w1[4].x[2]=7.39;
samp_w1[5].x[0]=4.94; samp_w1[5].x[1]=3.29; samp_w1[5].x[2]=2.08;
samp_w1[6].x[0]=-2.51; samp_w1[6].x[1]=2.09; samp_w1[6].x[2]=-2.59;
samp_w1[7].x[0]=-2.25; samp_w1[7].x[1]=-2.13; samp_w1[7].x[2]=-6.94;
samp_w1[8].x[0]=5.56; samp_w1[8].x[1]=2.86; samp_w1[8].x[2]=-2.26;
samp_w1[9].x[0]=1.03; samp_w1[9].x[1]=-3.33; samp_w1[9].x[2]=4.33;
samp_w2[0].x[0]=-0.91; samp_w2[0].x[1]=-0.18; samp_w2[0].x[2]=-0.05;
samp_w2[1].x[0]=1.30; samp_w2[1].x[1]=-2.06; samp_w2[1].x[2]=-3.53;
samp_w2[2].x[0]=-7.75; samp_w2[2].x[1]=-4.54; samp_w2[2].x[2]=-0.95;
samp_w2[3].x[0]=-5.47; samp_w2[3].x[1]=0.50; samp_w2[3].x[2]=3.92;
samp_w2[4].x[0]=6.14; samp_w2[4].x[1]=5.72; samp_w2[4].x[2]=-4.85;
samp_w2[5].x[0]=3.60; samp_w2[5].x[1]=1.26; samp_w2[5].x[2]=4.36;
samp_w2[6].x[0]=5.37; samp_w2[6].x[1]=-4.63; samp_w2[6].x[2]=-3.65;
samp_w2[7].x[0]=7.18; samp_w2[7].x[1]=1.46; samp_w2[7].x[2]=-6.66;
samp_w2[8].x[0]=-7.39; samp_w2[8].x[1]=1.17; samp_w2[8].x[2]=6.30;
samp_w2[9].x[0]=-7.50; samp_w2[9].x[1]=-6.32; samp_w2[9].x[2]=-0.31;
u1[0]=0;u1[1]=0;u1[2]=0;u2[0]=0;u2[1]=0;u2[2]=0;
for(i=0;i<10;i++) //caculate Mu
{
u1[0]=samp_w1[i].x[0]+u1[0];
u1[1]=samp_w1[i].x[1]+u1[1];
u1[2]=samp_w1[i].x[2]+u1[2];
u2[0]=samp_w2[i].x[0]+u2[0];
u2[1]=samp_w2[i].x[1]+u2[1];
u2[2]=samp_w2[i].x[2]+u2[2];
}
for(i=0;i<3;i++)
{
u1[i]=u1[i]/10;
u2[i]=u2[i]/10;
}
if((fp=fopen("lienioutput.txt","wr"))==NULL)
{ printf("can not open this file\n");
exit(0);
}
////////////////////////////////////////////////////////////
printf("mu1=[%6.2f,%6.2f,%6.2f]T ",u1[0],u1[1],u1[2]);
printf("mu2=[%6.2f,%6.2f,%6.2f]T\n",u2[0],u2[1],u2[2]);
fprintf(fp,"mu1=[%6.2f,%6.2f,%6.2f]T ",u1[0],u1[1],u1[2]);
fprintf(fp,"mu2=[%6.2f,%6.2f,%6.2f]T\n\n",u2[0],u2[1],u2[2]);
/////////////////////////////////////////////////////////
seta1[3][3]=((0,0,0),(0,0,0),(0,0,0));
seta2[3][3]=((0,0,0),(0,0,0),(0,0,0));
for(i=0;i<10;i++) //caculate sigma
{
x_u1[0]=samp_w1[i].x[0]-u1[0];
x_u1[1]=samp_w1[i].x[1]-u1[1];
x_u1[2]=samp_w1[i].x[2]-u1[2];
x_u2[0]=samp_w2[i].x[0]-u2[0];
x_u2[1]=samp_w2[i].x[1]-u2[1];
x_u2[2]=samp_w2[i].x[2]-u2[2];
sigma_w1[i].s[0][0]=x_u1[0]*x_u1[0];
sigma_w1[i].s[0][1]=x_u1[0]*x_u1[1];
sigma_w1[i].s[0][2]=x_u1[0]*x_u1[2];
sigma_w1[i].s[1][0]=x_u1[1]*x_u1[0];
sigma_w1[i].s[1][1]=x_u1[1]*x_u1[1];
sigma_w1[i].s[1][2]=x_u1[1]*x_u1[2];
sigma_w1[i].s[2][0]=x_u1[2]*x_u1[0];
sigma_w1[i].s[2][1]=x_u1[2]*x_u1[1];
sigma_w1[i].s[2][2]=x_u1[2]*x_u1[2];
sigma_w2[i].s[0][0]=x_u2[0]*x_u2[0];
sigma_w2[i].s[0][1]=x_u2[0]*x_u2[1];
sigma_w2[i].s[0][2]=x_u2[0]*x_u2[2];
sigma_w2[i].s[1][0]=x_u2[1]*x_u2[0];
sigma_w2[i].s[1][1]=x_u2[1]*x_u2[1];
sigma_w2[i].s[1][2]=x_u2[1]*x_u2[2];
sigma_w2[i].s[2][0]=x_u2[2]*x_u2[0];
sigma_w2[i].s[2][1]=x_u2[2]*x_u2[1];
sigma_w2[i].s[2][2]=x_u2[2]*x_u2[2];
seta1[0][0]=sigma_w1[i].s[0][0]+seta1[0][0];
seta1[0][1]=sigma_w1[i].s[0][1]+seta1[0][1];
seta1[0][2]=sigma_w1[i].s[0][2]+seta1[0][2];
seta1[1][0]=sigma_w1[i].s[1][0]+seta1[1][0];
seta1[1][1]=sigma_w1[i].s[1][1]+seta1[1][1];
seta1[1][2]=sigma_w1[i].s[1][2]+seta1[1][2];
seta1[2][0]=sigma_w1[i].s[2][0]+seta1[2][0];
seta1[2][1]=sigma_w1[i].s[2][1]+seta1[2][1];
seta1[2][2]=sigma_w1[i].s[2][2]+seta1[2][2];
seta2[0][0]=sigma_w2[i].s[0][0]+seta2[0][0];
seta2[0][1]=sigma_w2[i].s[0][1]+seta2[0][1];
seta2[0][2]=sigma_w2[i].s[0][2]+seta2[0][2];
seta2[1][0]=sigma_w2[i].s[1][0]+seta2[1][0];
seta2[1][1]=sigma_w2[i].s[1][1]+seta2[1][1];
seta2[1][2]=sigma_w2[i].s[1][2]+seta2[1][2];
seta2[2][0]=sigma_w2[i].s[2][0]+seta2[2][0];
seta2[2][1]=sigma_w2[i].s[2][1]+seta2[2][1];
seta2[2][2]=sigma_w2[i].s[2][2]+seta2[2][2];
}//end sigma matrix
printf("SIGMA1= SIGMA2=\n");
fprintf(fp,"SIGMA1= SIGMA2=\n");
for(j=0;j<3;j++)
{
for(k=0;k<3;k++)
{
seta1[j][k]=seta1[j][k]/10;
seta2[j][k]=seta2[j][k]/10;
}
printf(" |%6.2f,%6.2f,%6.2f| |%6.2f,%6.2f,%6.2f|\n",seta1[j][0],seta1[j][1],seta1[j][2],seta2[j][0],seta2[j][1],seta2[j][2]);
fprintf(fp," |%6.2f,%6.2f,%6.2f| |%6.2f,%6.2f,%6.2f|\n",seta1[j][0],seta1[j][1],seta1[j][2],seta2[j][0],seta2[j][1],seta2[j][2]);
}//end seta/10
//|sigma|
seta_w1_hl=hanglieshi(seta1);
seta_w2_hl=hanglieshi(seta2);
printf("|SIGMA1|=%8.2f, |SIGMA2|=%8.2f\n",seta_w1_hl,seta_w2_hl);
fprintf(fp,"\n|SIGMA1|=%8.2f, |SIGMA2|=%8.2f\n\n",seta_w1_hl,seta_w2_hl);
//inverse matrix
a11=seta1[0][0]; a12=seta1[0][1]; a13=seta1[0][2];
a21=seta1[1][0]; a22=seta1[1][1]; a23=seta1[1][2];
a31=seta1[2][0]; a32=seta1[2][1]; a33=seta1[2][2];
b11=a22*a33-a23*a32; inver_seta1[0][0]=b11;
b12=-(a21*a33-a23*a31); inver_seta1[1][0]=b12;
b13=a21*a32-a22*a31; inver_seta1[2][0]=b13;
b21=-(a12*a33-a13*a32); inver_seta1[0][1]=b21;
b22=a11*a33-a13*a31; inver_seta1[1][1]=b22;
b23=-(a11*a32-a12*a31); inver_seta1[2][1]=b23;
b31=a12*a23-a13*a22; inver_seta1[0][2]=b31;
b32=-(a11*a23-a13*a22); inver_seta1[1][2]=b32;
b33=a11*a22-a12*a21; inver_seta1[2][2]=b33;
a11=seta2[0][0]; a12=seta2[0][1]; a13=seta2[0][2];
a21=seta2[1][0]; a22=seta2[1][1]; a23=seta2[1][2];
a31=seta2[2][0]; a32=seta2[2][1]; a33=seta2[2][2];
b11=a22*a33-a23*a32; inver_seta2[0][0]=b11;
b12=-(a21*a33-a23*a31); inver_seta2[1][0]=b12;
b13=a21*a32-a22*a31; inver_seta2[2][0]=b13;
b21=-(a12*a33-a13*a32); inver_seta2[0][1]=b21;
b22=a11*a33-a13*a31; inver_seta2[1][1]=b22;
b23=-(a11*a32-a12*a31); inver_seta2[2][1]=b23;
b31=a12*a23-a13*a22; inver_seta2[0][2]=b31;
b32=-(a11*a23-a13*a22); inver_seta2[1][2]=b32;
b33=a11*a22-a12*a21; inver_seta2[2][2]=b33;
printf("inverse SIGMA1= inverse SIGMA2=\n");
fprintf(fp,"inverse SIGMA1= inverse SIGMA2=\n");
for(j=0;j<3;j++)
{ printf(" |%6.2f,%6.2f,%6.2f| |%6.2f,%6.2f,%6.2f|\n",inver_seta1[j][0]/seta_w1_hl,inver_seta1[j][1]/seta_w1_hl,inver_seta1[j][2]/seta_w1_hl,inver_seta2[j][0]/seta_w2_hl,inver_seta2[j][1]/seta_w2_hl,inver_seta2[j][2]/seta_w2_hl);
fprintf(fp," |%6.2f,%6.2f,%6.2f| |%6.2f,%6.2f,%6.2f|\n",inver_seta1[j][0]/seta_w1_hl,inver_seta1[j][1]/seta_w1_hl,inver_seta1[j][2]/seta_w1_hl,inver_seta2[j][0]/seta_w2_hl,inver_seta2[j][1]/seta_w2_hl,inver_seta2[j][2]/seta_w2_hl);
}
/// Aadichotomizer for two categories
printf("--------------------------------------------------------------------\n");
fprintf(fp,"\n--------------------------------------------------------------------\n");
printf("Sample g1(x) g2(x) Decision Correct or Not\n");
fprintf(fp,"Sample g1(x) g2(x) Decision Correct or Not\n");
printf("--------------------------------------------------------------------\n");
fprintf(fp,"--------------------------------------------------------------------\n");
int h;
for(h=0;h<20;h++)
{
if(h<10)
{
x_uw1[0]=samp_w1[h].x[0];
x_uw1[1]=samp_w1[h].x[1];
x_uw1[2]=samp_w1[h].x[2];
}else{
x_uw1[0]=samp_w2[h-10].x[0];
x_uw1[1]=samp_w2[h-10].x[1];
x_uw1[2]=samp_w2[h-10].x[2];
}
printf("%2d[%5.2f,%5.2f,%5.2f]",h,x_uw1[0],x_uw1[1],x_uw1[2]);
fprintf(fp,"%2d[%5.2f,%5.2f,%5.2f]",h,x_uw1[0],x_uw1[1],x_uw1[2]);
for(i=0;i<3;i++)
{
x_uw1[i]=x_uw1[i]-u1[i];
x_uw2[i]=x_uw1[i]-u2[i];
dmatrixw1[i]=x_uw1[i];
dmatrixw2[i]=x_uw2[i];
}
for(i=0;i<3;i++)
{
dmatrixw1[i]=x_uw1[0]*inver_seta1[0][i]+x_uw1[1]*inver_seta1[1][i]+x_uw1[2]*inver_seta1[2][i];
dmatrixw2[i]=x_uw2[0]*inver_seta2[0][i]+x_uw2[1]*inver_seta2[1][i]+x_uw2[2]*inver_seta2[2][i];
}
mma1=(dmatrixw1[0]*x_uw1[0]+dmatrixw1[1]*x_uw1[1]+dmatrixw1[2]*x_uw1[2])/seta_w1_hl;
mma2=(dmatrixw2[0]*x_uw2[0]+dmatrixw2[1]*x_uw2[1]+dmatrixw2[2]*x_uw2[2])/seta_w2_hl;
dLogsigma1=log(seta_w1_hl)/2;
dLogsigma2=log(seta_w2_hl)/2;
gxw1=-mma1*0.5-dLogsigma1;
gxw2=-mma2*0.5-dLogsigma2;
printf(" %6.2f %6.2f",gxw1,gxw2);
fprintf(fp," %6.2f %6.2f",gxw1,gxw2);
if(gxw1>gxw2)
{
printf(" w1");fprintf(fp," w1");
if(h<10){
iCorrect++;
printf(" Correct\n");fprintf(fp," Correct\n");
}else{
iWrong++;
printf(" Wrong\n");fprintf(fp," Wrong\n");
}
printf("--------------------------------------------------------------------\n");
fprintf(fp,"--------------------------------------------------------------------\n");
}else{
printf(" w2");fprintf(fp," w2");
if(h>10||h==10){
iCorrect++;
printf(" Correct\n");fprintf(fp," Correct\n");
}else{
iWrong++;
printf(" Wrong\n");
fprintf(fp," Wrong\n");
}
printf("--------------------------------------------------------------------\n");
fprintf(fp,"--------------------------------------------------------------------\n");
}
}//end for h=20
printf("Sample=20 Correct=%d Wrong=%d The empirical training error=%5.2f\n",iCorrect,iWrong,(double)iWrong/20);
printf("--------------------------------------------------------------------\n");
fprintf(fp,"Sample=20 Correct=%d Wrong=%d The empirical training error=%5.2f\n",iCorrect,iWrong,(double)iWrong/20);
fprintf(fp,"--------------------------------------------------------------------\n");
//Use the Bhattacharyya bound to bound the error
//P(e)<= ;k(1/2)=
setaplus[0][0]=(seta1[0][0]+seta2[0][0])/2;
setaplus[0][1]=(seta1[0][1]+seta2[0][1])/2;
setaplus[0][2]=(seta1[0][2]+seta2[0][2])/2;
setaplus[1][0]=(seta1[1][0]+seta2[1][0])/2;
setaplus[1][1]=(seta1[1][1]+seta2[1][1])/2;
setaplus[1][2]=(seta1[1][2]+seta2[1][2])/2;
setaplus[2][0]=(seta1[2][0]+seta2[2][0])/2;
setaplus[2][1]=(seta1[2][1]+seta2[2][1])/2;
setaplus[2][2]=(seta1[2][2]+seta2[2][2])/2;
setaplus_hl=hanglieshi(setaplus);
a11=setaplus[0][0];
a12=setaplus[0][1];
a13=setaplus[0][2];
a21=setaplus[1][0];
a22=setaplus[1][1];
a23=setaplus[1][2];
a31=setaplus[2][0];
a32=setaplus[2][1];
a33=setaplus[2][2];
b11=a22*a33-a23*a32; inver_setaplus[0][0]=b11;
b12=-(a21*a33-a23*a31); inver_setaplus[1][0]=b12;
b13=a21*a32-a22*a31; inver_setaplus[2][0]=b13;
b21=-(a12*a33-a13*a32); inver_setaplus[0][1]=b21;
b22=a11*a33-a13*a31; inver_setaplus[1][1]=b22;
b23=-(a11*a32-a12*a31); inver_setaplus[2][1]=b23;
b31=a12*a23-a13*a22; inver_setaplus[0][2]=b31;
b32=-(a11*a23-a13*a22); inver_setaplus[1][2]=b32;
b33=a11*a22-a12*a21; inver_setaplus[2][2]=b33;
JBln=log(setaplus_hl/(sqrt(seta_w1_hl*seta_w2_hl)))/2;
u2_u1[0]=u2[0]-u1[0];
u2_u1[1]=u2[1]-u1[1];
u2_u1[2]=u2[2]-u1[2];
for(i=0;i<3;i++)
{
dmatrixw1[i]=u2_u1[0]*inver_setaplus[0][i]+u2_u1[1]*inver_setaplus[1][i]+u2_u1[2]*inver_setaplus[2][i];
}
JBUT=(dmatrixw1[0]*u2_u1[0]+dmatrixw1[1]*u2_u1[1]+dmatrixw1[2]*u2_u1[2])/setaplus_hl;
JB=JBUT/8+JBln;
Pe=sqrt(0.25)*exp(-JB);
printf("Bhattachryya Bound P(error)<=%4.2lf k(1/2)=%6.2f\n",Pe,JB);
fprintf(fp,"Bhattachryya Bound P(error)<=%4.2lf k(1/2)=%6.2f\n",Pe,JB);
fclose(fp);
int quit1=0;
do{
printf("input 0 to quit program!:\n");
scanf("%d",&quit1);
}while(quit1!=0);
}//end main
long double hanglieshi(long double a[3][3])
{
int i,j,k=1,m,n,s,t;
long double sn,f=1,c,x;
for (i=0,j=0;i<3&&j<3;i++,j++)
{
if (a[i][j]==0)
{
for (m=i;a[m][j]==0;m++);
if (m==3)
{
sn=0; exit(0);
}
else
for (n=j;n<3;n++)
{
c=a[i][n];
a[i][n]=a[m][n];
a[m][n]=c;
}
k*=(-1);
}
for (s=3-1;s>i;s--)
{
x=a[s][j];
for (t=j;t<3;t++)
a[s][t]-=a[i][t]*(x/a[i][j]);
}
}
for (i=0;i<3;i++)
f*=a[i][i];
sn=k*f;
//printf("detA=%4.2lf\n",sn);
return(sn);
}
//MatrixOpp
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -