我按照教程,将C语言程序包装成了.mexw文件,但接下来不知道如何使用,在simulink中创建S-funtion组件之后,不会进一步设置。下面是C语言程序(多元线性回归,读取.csv文件中的表格,输入x,y。在所输入的最近区间求z=(x,y)的回归)求大佬教如何使用
#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:方程元数;
//Answer[count]:求解数组。返回0,求解成功。-1无解或无穷解;
int LinearEquations(double *data,int count,double *Answer)
{
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--)
d[m] += Answer[j] * dat[m][j];
Answer[m] = (dat[m][count]-d[m]) / dat[m][m];
}
FreeData(dat, d, count);
return 0;
}
//求多元 回归方程:Y=B0+B1X1+B2X2+......+BnXn
//data[rows*cols]二维数组:X1i,X2i......Xni,Yi(i=0 to rows-1)
//rows:数据行数;cols数据列表;Answer[cols]:返回回归系数数组(B0,B1......Bn)
//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)
}
n=LinearEquations(dat, cols, Answer); //计算方程式
//方差分析
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 += (*p * Answer[i]);
//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';
printf("回归方程式: Z= %.5lf", Answer[0]);
for(i=1; i<cols;i++)
printf("+%.5lf*%c",Answer[i], ch+i-1);
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 ++)
{
v= Answer[0];
for(j = 1; j < cols; j ++, p ++)
v += *p * Answer[j];
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;
//若符合矩阵格式,则进行矩阵的多元线性回归方程运算,求得Answer[0](常数),Answer[1](x的k值),Answer[2](y的k值);
double Answer[5],SquarePoor[4];
if(MultipleRegression((double*)data,4,3,Answer,SquarePoor)==0)
Display((double*)data, Answer, SquarePoor, 4, 3);
z=Answer[0]+x*Answer[1]+y*Answer[2]; //将x,y代入到求出的回归方程
printf("Z=%.5lf",z); //输出z的值
return 0;//结束
}