用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;
}
输出回归后的参数和对样例的测试结果。