Васильев N1

Материал из Wiki
(Различия между версиями)
Перейти к: навигация, поиск
(Новая страница: « hjvbhbhbh»)
 
 
(не показаны 6 промежуточных версий 2 участников)
Строка 1: Строка 1:
 +
Print:
  
−
[[hjvbhbhbh]]
+
sigmaX127.263 sigmaX226.5292 sigmaY205.134
 +
 
 +
rx1y=0.831126 rx2y=0.566734
 +
 
 +
MNK_first_method: 15.0154 6.19976 4.30021
 +
 
 +
MNK_second_method: 15.0154 6.19976 4.30021
 +
 
 +
Regres: 7.3335 6.25361 4.3822
 +
 
 +
5.80035e-16 5.80035e-16
 +
 
 +
1.93219 1.93219 1300.63
 +
 
 +
//============================================================================
 +
// Name        : MNK.cpp
 +
// Author      :
 +
// Version    :
 +
// Copyright  : Your copyright notice
 +
// Description : Hello World in C++, Ansi-style
 +
//============================================================================
 +
 
 +
#include <iostream>
 +
#include <fstream>
 +
#include <cstdlib>
 +
#include <cmath>
 +
#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;
 +
}

Текущая версия на 12:20, 24 мая 2019

Print:

sigmaX127.263 sigmaX226.5292 sigmaY205.134

rx1y=0.831126 rx2y=0.566734

MNK_first_method: 15.0154 6.19976 4.30021

MNK_second_method: 15.0154 6.19976 4.30021

Regres: 7.3335 6.25361 4.3822

5.80035e-16 5.80035e-16

1.93219 1.93219 1300.63

//============================================================================ // 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; }

Источник — «http://v.michm.ru/index.php?title=%D0%92%D0%B0%D1%81%D0%B8%D0%BB%D1%8C%D0%B5%D0%B2_N1&oldid=20144»
Персональные инструменты
Пространства имён

Варианты
Действия
Навигация
Инструменты