机器学习-学习笔记(三)第三章 习题3.3

用C语言实现了一下3.3中的对率回归问题,C语言求逆矩阵是个很麻烦的事情,这里参考了一些网上的博客,偷了懒。

#include<stdio.h>
#include<string.h>
#include<math.h> 
typedef struct parameter{
	double w[10];
	int w_lenght;
	double b;
}parameter;
typedef struct input{
	double x[10];
	int x_lenght;
	double const_num;
}input;
double mix[3][3];
void matrix_a(double temp[3][3])
{
	int k, k2;
	for (k = 0; k < 3; k++){
		for(k2 = 0; k2 < 3; k2++){
			mix[k][k2]+=temp[k][k2]; 
		}
	}
}
int  matrix_inv(int ndimen)
{
	double tmp, tmp2, b_tmp[20], c_tmp[20];
	int k, k1, k2, k3, j, i, j2, i2, kme[20], kmf[20];
	i2 = j2 = 0;
	for (k = 0; k < ndimen; k++){
		tmp2 = 0.0;
		for (i = k; i < ndimen; i++){
			for (j = k; j < ndimen; j++){
				if (fabs(mix[i][j] ) <= fabs(tmp2)) 
					continue;
				tmp2 = mix[i][j];
				i2 = i;
				j2 = j;
			}  
		}
		if (i2 != k){
			for (j = 0; j < ndimen; j++){
				tmp = mix[i2][j];
				mix[i2][j] = mix[k][j];
				mix[k][j] = tmp;
			}
		}
		if (j2 != k){
			for (i = 0; i < ndimen; i++){
				tmp = mix[i][j2];
				mix[i][j2] = mix[i][k];
				mix[i][k] = tmp;
			}    
		}
		kme[k] = i2;
		kmf[k] = j2;
		for (j = 0; j < ndimen; j++){
			if (j == k){
				b_tmp[j] = 1.0 / tmp2;
				c_tmp[j] = 1.0;
			}
			else{
				b_tmp[j] = -mix[k][j] / tmp2;
				c_tmp[j] = mix[j][k];
			}
			mix[k][j] = 0.0;
			mix[j][k] = 0.0;
		}
		for (i = 0; i < ndimen; i++){
			for (j = 0; j < ndimen; j++){
				mix[i][j] = mix[i][j] + c_tmp[i] * b_tmp[j];
			}  
		}
	}
	for (k3 = 0; k3 < ndimen;  k3++){
		k  = ndimen - k3 - 1;
		k1 = kme[k];
		k2 = kmf[k];
		if (k1 != k){
			for (i = 0; i < ndimen; i++){
				tmp = mix[i][k1];
				mix[i][k1] = mix[i][k];
				mix[i][k] = tmp;
			}  
		}
		if (k2 != k){
			for(j = 0; j < ndimen; j++)  
			{
				tmp = mix[k2][j];
				mix[k2][j] = mix[k][j];
				mix[k][j] = tmp;
			}
		}
	}
	return (0);
}

double p1(input x,parameter p){
	double answer;
	answer=exp(x.x[0]*p.w[0]+x.x[1]*p.w[1]+p.b)/(exp(x.x[0]*p.w[0]+x.x[1]*p.w[1]+p.b)+1);
	//printf("p1==%lf\n",answer);
	return answer;
}

parameter OneDerivative(input inputs[17],parameter p){
	parameter answer;
	answer.b=0;
	memset(answer.w,0,sizeof(answer.w));
	
	for(int i=0;i<17;i++){
		double temp=inputs[i].const_num-p1(inputs[i],p);
		answer.w[0]+=(-inputs[i].x[0]*temp);
		answer.w[1]+=(-inputs[i].x[1]*temp);
		answer.b+=(-temp);
	}
	return answer;
}

double TwoDerivative(input inputs[17],parameter p){
	double answer=0; 
	double temp[3][3],nm[3],x;
	for(int i=1;i<17;i++){
		nm[0]=inputs[i].x[0];
		nm[1]=inputs[i].x[1];
		nm[2]=1;
		x=p1(inputs[i],p);
		for(int j=0;j<3;j++){
			for(int z=0;z<3;z++){
				temp[j][z]=nm[j]*nm[z]*x*(1-x);
			}
		}
		matrix_a(temp);
	}
	return answer;
}

void prin(parameter p){
	printf("w1=%lf w2=%lf b=%lf\n",p.w[0],p.w[1],p.b);
}
double lps(input inputs[17],parameter p){
	double lp=0;
	for(int i=0;i<17;i++){
		lp+=(log(1+exp(inputs[i].x[0]*p.w[0]+inputs[i].x[1]*p.w[1]+p.b))-(inputs[i].const_num*(inputs[i].x[0]*p.w[0]+inputs[i].x[1]*p.w[1]+p.b)));
	}
	return lp;
}
int main(){
	double old_lp=1,new_lp=0;
	parameter par;
	input inputs[17],tests[4];
	memset(inputs,0,sizeof(inputs));
	memset(tests,0,sizeof(tests));
	int N,flag=0;
	par.b=1;
	memset(par.w,0,sizeof(par.w));
	par.w[0]=0;
	par.w[1]=0;
	
	freopen("input.txt","r",stdin);
	scanf("%d",&N);
	for(int i=0;i<N;i++){
		scanf("%lf %lf %lf",&inputs[i].x[0],&inputs[i].x[1],&inputs[i].const_num);
	}
	while(fabs(old_lp-new_lp)>0.001){
		if(flag==1){
			old_lp=new_lp;
		}	
		parameter temp=OneDerivative(inputs,par);
		memset(mix,0,sizeof(mix));
		TwoDerivative(inputs,par);
		
		matrix_inv(3);
		
		double temps[3];
		for(int i=0;i<3;i++){
			temps[i]=mix[i][0]*temp.w[0]+mix[i][1]*temp.w[1]+mix[i][2]*temp.b;
		}
		par.w[0]=par.w[0]-temps[0];
		par.w[1]=par.w[1]-temps[1];
		par.b=par.b-temps[2];
		new_lp=lps(inputs,par);
		
		if(flag==0){
			flag=1;
		}
		
	}
	prin(par);
	//freopen("test.txt","r",stdin);
	scanf("%d",&N);
	for(int i=0;i<N;i++){
		scanf("%lf %lf",&tests[i].x[0],&tests[i].x[1]);
	}
	for(int i=0;i<N;i++){
		double z=(tests[i].x[0]*par.w[0]+tests[i].x[1]*par.w[1]+par.b);
		z=0-z;
		printf("y==%lf\n",1/(1+exp(z)));
	}
	return 0;
}

输出回归后的参数和对样例的测试结果。

点赞