莓红的十字架

2019-11-16 12:37

``````#include<stdio.h>
#include"math.h"
#include <stdlib.h>
#include<string.h>
void FreeData(double **dat, double *d, int count)
{
int i, j;
free(d);
for(i = 0;i < count;  i ++)
free(dat[i]);
free(dat);
}
//解线性方程。data[count*[count+1])矩阵数组；count：方程元数；
{
int j ,m ,n;
double tmp, **dat, *d=data;
dat = (double**)malloc(count * sizeof(double*));
for (m=0;m<count;m++,d+=(count +1))
{
dat[m] = (double*)malloc((count+1) * sizeof(double));
memcpy(dat[m],d, (count+1) * sizeof(double));
}
d = (double*)malloc((count + 1) * sizeof(double));
for(m = 0; m < count - 1;m ++)
{
//如果主对角线元素为0，行交换；
for(n = m + 1;n < count && dat[m][m] == 0.0;n ++)
{
if( dat[n][m] != 0.0)
{
memcpy(d, dat[m], (count + 1)*sizeof(double));
memcpy(dat[m], dat[n], (count + 1) * sizeof(double));
memcpy(dat[n], d, (count + 1) * sizeof(double));
}
}
//行交换后，主对角线元素仍然为0，无解，返回-1；
if ( dat[m][m] == 0.0)
{
FreeData(dat, d, count);
return -1;
}
//消元
for(n = m + 1; n < count; n++)
{
tmp= dat[n][m] / dat[m][m];
for(j=m; j <= count; j++)
dat[n][j] -= tmp * dat[m][j];
}
}
for(j=0; j<count; j++)
d[j] = 0.0;
//求得count - 1 的元
Answer[count - 1]= dat[count - 1][count] / dat[count - 1][count - 1];
//逐行代入求各元
for (m = count - 2;m >= 0; m --)
{
for(j=count-1; j>m ; j--)
}
FreeData(dat, d, count);
return 0;

}
//求多元 回归方程：Y=B0+B1X1+B2X2+......+BnXn
//data[rows*cols]二维数组：X1i,X2i......Xni,Yi(i=0 to rows-1)
//SquarePoor[4]:返回方差分析指标：回归平方和，剩余平方和，回归平方差，剩余平方差
//返回值：0求解成功，-1错误；
int MultipleRegression(double *data, int rows, int cols, double *Answer, double *SquarePoor)
{
int m, n, i, count = cols - 1;
double *dat, *p, a, b;
if(data == 0 || Answer == 0 || rows<2 || cols<2)
return -1;
dat = (double*)malloc(cols * (cols + 1) * sizeof(double));
dat[0] = (double)rows;
for(n=0;n<count;n++)                                                           //n=0 to cols-2
{
a = b = 0.0;
for(p = data + n, m = 0; m < rows; m ++, p += cols)
{
a += *p;
b += (*p * *p);
}
dat[n + 1] = a;                                                //dat[0,n+1]=Sum(Xn)
dat[(n + 1) * (cols + 1)] = a;                            //dat[n+1,0]=Sum(Xn)
dat[(n + 1) * (cols + 1) + n + 1] = b;                         //dat[n+1,n+1]=Sum(Xn*Xn)
for(i = n + 1; i < count; i++)                                             //i=n+1 to cols-2
{
for(a = 0.0, p = data, m = 0; m < rows; m ++, p += cols)
a += (p[n] * p[i]);
dat[(n+1) * (cols + 1) + i + 1] = a;            //dat[n+1,i+1]=Sum(Xn*Xi)
dat[(i+1) * (cols + 1) + n + 1] = a;           //dat[i+1,n+1]=Sum(Xn*Xi)
}
}
for(b = 0.0, m = 0, p = data + n; m < rows; m++, p += cols)
b += *p;
dat[cols]= b;                                                                //dat[0,cols]=Sum(Y)
for(n = 0;n < count; n++)
{
for(a = 0.0,p = data, m = 0; m < rows; m ++,p += cols)
a += (p[n] * p[count]);
dat[(n+1) * (cols + 1) + cols] = a;               //dat[n+1,cols]=Sum(Xn*Y)
}
//方差分析
if(n == 0 && SquarePoor)
{
b = b / rows;                                               //b=Y的平均值
SquarePoor[0] = SquarePoor[1] = 0.0;
p = data;
for(m = 0; m < rows; m ++, p ++)
{
for( i=1, a = Answer[0]; i < cols;i ++,p ++)

//a=Ym的估计值
SquarePoor[0] += ((a - b) * (a - b));
//U(回归平方和）
SquarePoor[1] += ((*p - a)*(*p - a));
//Q(剩余平方和）（*p=Ym）
}
SquarePoor[2] = SquarePoor[0] / count;
//回归方差
if(rows - cols > 0.0)
SquarePoor[3] = SquarePoor[1] / (rows - cols);//剩余方差
else
SquarePoor[3] = 0.0;
}
free(dat);
return n;
}
//输出回归方程，并输出误差估计
void Display(double *dat, double *Answer, double *SquarePoor, int rows, int cols)
{
double v, *p;
int i, j;
char ch='X';
for(i=1; i<cols;i++)
printf(" \n");
printf("回归显著性检验：");
printf("回归平方和:   %12.4lf  \n  回归方差：%12.4lf\n", SquarePoor[0], SquarePoor[2]);
printf("剩余平方和：%12.4lf  \n  剩余方差：%12.4lf\n", SquarePoor[1], SquarePoor[3]);
printf("离差平方和：%12.4lf  \n  标准误差：%12.4lf\n", SquarePoor[1], SquarePoor[3]);
printf("离差平方和：%12.4lf  \n  标准误差：%12.4lf\n", SquarePoor[0] + SquarePoor[1], sqrt(SquarePoor[3]));
printf("F      检  验 :  %12.4lf  \n  相关系数:  %12.4lf\n" ,SquarePoor[2] / SquarePoor[3], sqrt(SquarePoor[0] / (SquarePoor[0] + SquarePoor[1])));
printf("剩余分析: \n");
printf("     观察值      估计值       剩余值     剩余平方   \n");
for(i = 0, p = dat; i < rows; i ++, p ++)
{
for(j = 1; j < cols; j ++, p ++)
printf("%12.2lf%12.2lf%12.2lf%12.2lf\n", *p, v, *p - v, (*p - v) * (*p - v));
}
system("pause");
}
//主程序
int main()
{
double data[4][3];//定义矩阵，4列3行，4列为临近的四个点，3行为X，Y，Z；

FILE *fp = fopen("C://BK.csv", "r");//打开文件（对应的文件名和路径）
if (fp == NULL)  //如果文件打开失败则结束
{
printf("file open error\n");
return -1;
}
//定义Y的数组，Y[0]为Y的值，Y[1-1000]为该Y对应的Z值
double A[1000];
double B[1000];
double C[1000];
double D[1000];
double E[1000];
double F[1000];
double G[1000];
double H[1000];
double I[1000];
double J[1000];
double K[1000];
//运用循环语句，将文件中的数字矩阵存入到数组
for (int i = 0;i<255; i++)
{
fscanf(fp, "%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf", &A[i], &B[i], &C[i], &D[i], &E[i], &F[i], &G[i], &H[i], &I[i], &J[i], &K[i]);
}
double x,y,z;
scanf("%lf%lf",&x,&y);                      //输入已知的x，y
double X1,X2,Y1,Y2,Z1,Z2,Z3,Z4;
int j;
for(j=0;j<255;j++)
{
if(A[j]<x && x<A[j+1])                //判断x处于表格的哪个X值区间
{
//判断y处于表格的哪个Y值区间，并将锁定位置最近的四个数据记为（X1,Y1,Z1)(X2,Y1,Z2)(X1,Y2,Z3)(X2,Y2,Z4)
if(B[0]<y && y<C[0])
{
X1=A[j];
X2=A[j+1];
Y1=B[0];
Y2=C[0];
Z1=B[j];
Z2=B[j+1];
Z3=C[j];
Z4=C[j+1];
}
else if(C[0]<y && y<D[0])
{
X1=A[j];
X2=A[j+1];
Y1=C[0];
Y2=D[0];
Z1=C[j];
Z2=C[j+1];
Z3=D[j];
Z4=D[j+1];
}
else if(D[0]<y && y<E[0])
{
X1=A[j];
X2=A[j+1];
Y1=D[0];
Y2=E[0];
Z1=D[j];
Z2=D[j+1];
Z3=E[j];
Z4=E[j+1];
}
else if(E[0]<y && y<F[0])
{
X1=A[j];
X2=A[j+1];
Y1=E[0];
Y2=F[0];
Z1=E[j];
Z2=E[j+1];
Z3=F[j];
Z4=F[j+1];
}
else if(F[0]<y && y<G[0])
{
X1=A[j];
X2=A[j+1];
Y1=F[0];
Y2=G[0];
Z1=F[j];
Z2=F[j+1];
Z3=G[j];
Z4=G[j+1];
}
else if(G[0]<y && y<H[0])
{
X1=A[j];
X2=A[j+1];
Y1=G[0];
Y2=H[0];
Z1=G[j];
Z2=G[j+1];
Z3=H[j];
Z4=H[j+1];
}
else if(H[0]<y && y<I[0])
{
X1=A[j];
X2=A[j+1];
Y1=H[0];
Y2=I[0];
Z1=H[j];
Z2=H[j+1];
Z3=I[j];
Z4=I[j+1];
}
else if(I[0]<y && y<J[0])
{
X1=A[j];
X2=A[j+1];
Y1=I[0];
Y2=J[0];
Z1=I[j];
Z2=I[j+1];
Z3=J[j];
Z4=J[j+1];
}
else if(J[0]<y && y<K[0])
{
X1=A[j];
X2=A[j+1];
Y1=J[0];
Y2=K[0];
Z1=J[j];
Z2=J[j+1];
Z3=K[j];
Z4=K[j+1];
}
}

}
fclose(fp);             //结束文件读取
system("pause");           //关闭文件
//将（X1,Y1,Z1)(X2,Y1,Z2)(X1,Y2,Z3)(X2,Y2,Z4)，输入到data矩阵
data[0][0]=X1;
data[0][1]=Y1;
data[0][2]=Z1;
data[1][0]=X2;
data[1][1]=Y1;
data[1][2]=Z2;
data[2][0]=X1;
data[2][1]=Y2;
data[2][2]=Z3;
data[3][0]=X2;
data[3][1]=Y2;
data[3][2]=Z4;
printf("Z=%.5lf",z);    //输出z的值
return 0;//结束

}
``````
• 点赞
• 写回答
• 关注问题
• 收藏
• 复制链接分享
• 邀请回答