Васильев N1

//============================================================================ // Name  : MNK.cpp // Author  : // Version  : // Copyright  : Your copyright notice // Description : Hello World in C++, Ansi-style //============================================================================

1. include <iostream>
2. include <fstream>
3. include <cstdlib>
4. include <cmath>
5. include <cstring>

using namespace std;

void PrintArr(double *x, int N) { for (int i = 0; i < N; i++) { cout << x[i] << " "; } cout << endl; }

void AllocationMatrix(double **x, int size_1, int size_2) { x = new double*[size_1]; for (int i = 0; i < size_1; i++) { x[i] = new double[size_2]; } }

void TransposeMatrix (double **x, double **xT, int size_1, int size_2) { for (int i = 0; i < size_1; i++) { for (int j = 0; j < size_2; j++) { xT[j][i] = x[i][j]; } } }

void MultMatrix(double **x, double **y, double **c, int size_11, int size_22, int size_12) { for (int i = 0; i < size_11; i++) { for (int j = 0; j < size_22; j++) { for (int k = 0; k < size_12; k++) { c[i][j] += x[i][k] * y[k][j]; } } } }

void reverseMatrix (double **x, int size) { double **e; AllocationMatrix(e, size, size); double a = 0;

for (int i = 0; i < size; i++) for (int j = 0; j < size; j++) i == j ? e[i][j] = 1 : e[i][j] = 0;

for(int i = 0; i < size; i++){

a = x[i][i]; for(int j = 0 ; j < size; j++){ x[i][j] /= a; e[i][j] /= a; }

for (int k = i + 1; k < size; k++ ){

a = x[k][i]; for(int j = 0; j < size; j++){ x[k][j] -= x[i][j] * a; e[k][j] -= e[i][j] * a; } } }

for(int i = size - 1; i >0; i--){ for(int j = i - 1; j >= 0; j--){

a = x[j][i]; for(int k = 0; k < size; k++){ x[j][k] -= x[i][k] * a; e[j][k] -= e[i][k] * a; } } }

for (int i = 0; i < size; i++) for (int j = 0; j < size; j++) x[i][j] = e[i][j]; }

void ReadFile(char *filename, double *x1, double *x2, double *y, int N) { char buff[10]; ifstream fin(filename); for (int i = 0; i < N; i++) { for (int j = 0; j < 4; j++) { fin >> buff; if (j == 1) { x1[i] = atof(buff); } else if (j == 2) { x2[i] = atof(buff); } else if (j == 3) { y[i] = atof(buff); } } } }

double* MNK_first_method(double *x1_matrix, double *x2_matrix, double *y_matrix, int N) { double x1 = 0, x2 = 0, y = 0, yx1 = 0, x12 = 0, x11 = 0, yx2 = 0, x22 = 0; double *B = new double[3];

for (int i = 0; i < N; i++) { x1 += x1_matrix[i]; x2 += x2_matrix[i]; y += y_matrix[i]; yx1 += x1_matrix[i]*y_matrix[i]; x12 += x1_matrix[i]*x2_matrix[i]; x11 += x1_matrix[i]*x1_matrix[i]; yx2 += x2_matrix[i]*y_matrix[i]; x22 += x2_matrix[i]*x2_matrix[i]; } x1 /= N, x2 /= N, y /= N, yx1 /= N, x12 /= N, x11 /= N, yx2 /= N, x22 /= N;

B[2] = ((yx2 - y*x2)*(x1*x1 - x11) + (y*x1 - yx1)*(x1*x2 - x12))/(-2*x2*x1*x12+x12*x12 + x2*x2*x11 + x22*x1*x1 - x22*x11); B[1] = (y*x1 - yx1 - B[2]*(x2*x1 - x12))/(x1*x1 - x11); B[0] = y - B[1]*x1 - B[2]*x2;

return B;

}

double* MNK_second_method(double *x1, double *x2, double *y, int size) {

double **x; AllocationMatrix(x, size, 3); double **yMatrix; AllocationMatrix(yMatrix, size, 1); double **xT; AllocationMatrix(xT, 3, size); double **x_xT; AllocationMatrix(x_xT, 3, 3); double **x_xT_xT; AllocationMatrix(x_xT_xT, 3, size); double **B; AllocationMatrix(B, 3, 1);

for (int i = 0; i < size; i++) { x[i][0] = 1; x[i][1] = x1[i]; x[i][2] = x2[i]; yMatrix[i][0] = y[i]; }

// (xT*x)^-1 * xT * y = B TransposeMatrix(x, xT, size, 3); MultMatrix(xT, x, x_xT, 3, 3, size); reverseMatrix(x_xT, 3); MultMatrix(x_xT, xT, x_xT_xT, 3, size, 3); MultMatrix(x_xT_xT, yMatrix, B, 3, 1, size);

double *B_final = new double[3]; for (int i = 0; i < 3; i++) { B_final[i] = B[i][0]; }

return B_final; }

double* Regres(double *x1, double *x2, double *y, int N) { double m1x1 = 0, m2x1 = 0, sigmaX1 = 0; double m1x2 = 0, m2x2 = 0, sigmaX2 = 0;

double m1y = 0, m2y = 0, sigmaY = 0; double rx1y = 0, rx2y = 0;

double *b = new double[3];

for (int i = 0; i < N; i++) { m1x1 += x1[i]; m1x2 += x2[i]; m1y += y[i]; }

m1x1 /= N, m1x2 /= N, m1y /= N;

for (int i = 0; i < N; i++) { m2x1 += (x1[i] - m1x1)*(x1[i] - m1x1); m2x2 += (x2[i] - m1x2)*(x2[i] - m1x2); m2y += (y[i] - m1y)*(y[i] - m1y); }

m2x1 /= (N - 1), m2x2 /= (N - 1), m2y /= (N - 1);

sigmaX1 = pow(m2x1, 0.5), sigmaX2 = pow(m2x2, 0.5), sigmaY = pow(m2y, 0.5);

for (int i = 0; i < N; i++) { rx1y += (x1[i] - m1x1)*(y[i] - m1y); rx2y += (x2[i] - m1x2)*(y[i] - m1y); }

rx1y /= ((N-1)*sigmaX1*sigmaY); rx2y /= ((N-1)*sigmaX2*sigmaY);

b[1] = rx1y*sigmaY/sigmaX1; b[2] = rx2y*sigmaY/sigmaX2;

double b1x1 = 0, b2x2 = 0; for (int i = 0; i < N; i++) { b1x1 += b[1]*x1[i]; b2x2 += b[2]*x2[i]; }

b[0] = m1y - (b[1]*m1x1 + b[2]*m1x2);

cout << "Print: " << endl; cout << "sigmaX1" << sigmaX1 << " sigmaX2" << sigmaX2 << " sigmaY" << sigmaY << endl; cout << "rx1y=" << rx1y << " rx2y=" << rx2y << endl;

return b; }

void Fisher(double *B1, double *B2, double *B3, double *x1, double *x2, double *y, double *y_parall, int N, int N_parall) {

double My1 = 0, My2 = 0, My3 = 0, Myp = 0;

for (int i = 0; i < N; i++) { My1 += (B1[0] + B1[1]*x1[i] + B1[2]*x2[i] - y[i]); My2 += (B2[0] + B2[1]*x1[i] + B2[2]*x2[i] - y[i]); My3 += (B3[0] + B3[1]*x1[i] + B3[2]*x2[i] - y[i]);

if (i < N_parall) Myp += y_parall[i]; }

My1 /= N, My2 /= N, My3 /= N, Myp /= N_parall;

cout << My1 << " " << My2 << endl;

double Dy1 = 0, Dy2 = 0, Dy3 = 0, Dyp = 0;

for (int i = 0; i < N; i++) { Dy1 += pow((B1[0] + B1[1]*x1[i] + B1[2]*x2[i] - y[i] - My1), 2); Dy2 += pow((B2[0] + B2[1]*x1[i] + B2[2]*x2[i] - y[i] - My2), 2); Dy3 += pow((B3[0] + B3[1]*x1[i] + B3[2]*x2[i] - y[i] - My3), 2);

if (i < N_parall) Dyp += (y_parall[i] - Myp)*(y_parall[i] - Myp); }

Dy1 /= (N - 3), Dy2 /= (N - 3), Dy3 /= (N - 3), Dyp /= (N_parall - 1);

cout << Dy1/Dyp << " " << Dy2/Dyp << " " << Dy3/Dyp << endl; }

int main() {

int N = 49; int N_parall = 6;

double *x1 = new double[N]; double *x2 = new double[N]; double *y = new double[N];

double *x1_parall = new double[N_parall]; double *x2_parall = new double[N_parall]; double *y_parall = new double[N_parall];

ReadFile("./data.mnk", x1, x2, y, N);

ReadFile("./parall.mnk", x1_parall, x2_parall, y_parall, N_parall);

double * B1 = new double[3]; double * B2 = new double[3]; double * B3 = new double[3];

B1 = MNK_first_method(x1, x2, y, N); B2 = MNK_first_method(x1, x2, y, N); B3 = Regres(x1, x2, y, N);

cout << "MNK_first_method: " << B1[0] << " " << B1[1] << " " << B1[2] << endl; cout << "MNK_second_method: " << B2[0] << " " << B2[1] << " " << B2[2] << endl; cout << "Regres: " << B3[0] << " " << B3[1] << " " << B3[2] << endl;

Fisher(B1, B2, B3, x1, x2, y, y_parall, N, N_parall);

return 0; }