StepRegressions.cpp

//---------------------------------------------------------------------------

#include <vcl.h>

#pragma hdrstop

//---------------------------------------------------------------------------

USEFORM("StepRegressionsUnit1.cpp", FormMain);

//---------------------------------------------------------------------------

WINAPI WinMain(HINSTANCE, HINSTANCE, LPSTR, int)

{

try

{

Application->Initialize();

Application->CreateForm(__classid(TFormMain), &FormMain);

Application->Run();

}

catch (Exception &exception)

{

Application->ShowException(&exception);

}

catch (...)

{

try

{

throw Exception("");

}

catch (Exception &exception)

{

Application->ShowException(&exception);

}

}

return 0;

}

//---------------------------------------------------------------------------

StepRegressionsUnit1.h

//---------------------------------------------------------------------------

#ifndef StepRegressionsUnit1H

#define StepRegressionsUnit1H

//---------------------------------------------------------------------------

#include <Classes.hpp>

#include <Controls.hpp>

#include <StdCtrls.hpp>

#include <Forms.hpp>

#include <Menus.hpp>

#include <Dialogs.hpp>

#include <ExtCtrls.hpp>

#include <ComCtrls.hpp>

//---------------------------------------------------------------------------

class TFormMain : public TForm

{

__published: // IDE-managed Components

TMainMenu *MainMenu1;

TMenuItem *N1;

TMenuItem *N2;

TMenuItem *N3;

TMenuItem *N4;

TMenuItem *N5;

TMenuItem *N6;

TMenuItem *N7;

TMenuItem *N8;

TOpenDialog *OpenDialog1;

TSaveDialog *SaveDialog1;

TRadioGroup *RadioGroup1;

TRadioButton *RadioButton1;

TRadioButton *RadioButton2;

TRadioButton *RadioButton3;

TMemo *Memo1;

TStatusBar *StatusBar1;

TComboBox *ComboBoxFontSize;

TButton *ButtonMemoClear;

TGroupBox *GroupBox1;

TCheckBox *CheckBox1;

TCheckBox *CheckBoxMaxCorr;

TCheckBox *CheckBoxSvob;

TCheckBox *CheckBoxCorrel;

TGroupBox *GroupBox2;

TEdit *EditFInp;

TLabel *LabelFinp;

TLabel *LabelFout;

TEdit *EditFOut;

void __fastcall FormDestroy(TObject *Sender);

void __fastcall N1Click(TObject *Sender);

void __fastcall N3Click(TObject *Sender);

void __fastcall N4Click(TObject *Sender);

void __fastcall N2Click(TObject *Sender);

void __fastcall ComboBoxFontSizeChange(TObject *Sender);

void __fastcall ButtonMemoClearClick(TObject *Sender);

private:

AnsiString StrMess; // строка сообщения для вывода в Memo1

long m; // количество наблюдений

long n; // количество исследуемых регрессоров

double * y; // вектор отклика размера m

double ** XT; // транспонированная матрица X(m*n) наблюдений размера n*m

double ** XTX; // произведение матриц XT на X для оценки вектора параметров

bool * jInModel; // вектор флажков о включении в модель соответствующего регрессора

bool * jInTemporaryModel; // --//-- о включении в пробную модель --//--

long nCurrent; // текущее количество регрессоров в модели

long * J; // массив размера nCurrent номеров узлов, включенных в модель

double * f; // призведение матрицы XT на Y - есть правая часть СЛАУ для df-ия theta

double * theta;

double RSS;

double MinCould_Max_Fj;

double Fin;

double Fout; // вектор параметров размера n

long MaxCorrellation(double & Max_rXjY, long & IndexOfMax_rXjY);

int WhichAddRegressor(double & _Max_Fj, long & _IndexOfMax_Fj);

int WhichExpelRegressor(double & _Bad_Fj, long & _IndexOfBad_Fj);

int ThetaForInModel(); // вычисление theta текущей модели

int ThetaForInModelSlaeResolve(); // решение СЛАУ в определении theta

int OutPutResult(AnsiString Mes); // вывод уравнения регрессии

double RSSForInModel(); // остаточная сумма квадратов для текущей модели

int AddRegressor(long NumOfAddRegressor);

int ExpelRegressor(long NumOfExpelRegressor);

public: // User declarations

__fastcall TFormMain(TComponent* Owner);

};

//---------------------------------------------------------------------------

extern PACKAGE TFormMain *FormMain;

//---------------------------------------------------------------------------

#endif

//---------------------------------------------------------------------------

StepRegressionsUnit1.cpp

//---------------------------------------------------------------------------

#include <stdio.h>

#include <math.h>

#include <vcl.h>

//---------------------------------------------------------------------------

#pragma hdrstop

//---------------------------------------------------------------------------

#include "StepRegressionsUnit1.h"

//---------------------------------------------------------------------------

#pragma package(smart_init)

#pragma resource "*.dfm"

//---------------------------------------------------------------------------

TFormMain *FormMain;

//-Конструктор формы---------------------------------------------------------

__fastcall TFormMain::TFormMain(TComponent* Owner)

: TForm(Owner)

{

nCurrent=0;

Memo1->Font->Size=ComboBoxFontSize->Text.ToInt();

MinCould_Max_Fj=3.28;

Fin=EditFInp->Text.ToDouble();

Fout=EditFOut->Text.ToDouble();

}

//-Деструктор формы----------------------------------------------------------

void __fastcall TFormMain::FormDestroy(TObject *Sender)

{

for(long j=0; j<n; j++)

delete [] XT[j];

delete [] XT;

delete [] y;

delete [] theta;

delete [] jInModel;

delete [] J;

}

//-Загрузка входных данных---------------------------------------------------

void __fastcall TFormMain::N1Click(TObject *Sender)

{

do

if(!OpenDialog1->Execute()) return;

while (OpenDialog1->FileName=="");

FILE *fi=NULL;

AnsiString FileName;

FileName=OpenDialog1->FileName;

if((fi=fopen(FileName.c_str(),"r"))==NULL)

{ ShowMessage("Файл\n"+FileName+"\nне открыт.\nЗагрузка провалена!!!");

return;

}

fscanf(fi,"%d %d",&m,&n);

n+=1;

XT=new double*[n];

for(long i=0; i<n; i++)

XT[i]=new double[m];

y=new double[m];

for(long i=0; i<m; i++)

{ XT[0][i]=1.0;

for(long j=1; j<n; j++)

fscanf(fi,"%lf",&XT[j][i]);

fscanf(fi,"%lf",&y[i]);

}

fclose(fi);

ShowMessage("Загрузка из файла\n"+FileName+"\nзавершена.");

N2->Enabled=true;

return;

}

//-Сохранение результатов в файл---------------------------------------------

void __fastcall TFormMain::N3Click(TObject *Sender)

{

do

if(!SaveDialog1->Execute()) return;

while (SaveDialog1->FileName=="");

AnsiString FileName;

FileName=SaveDialog1->FileName;

Memo1->Lines->SaveToFile(FileName);

ShowMessage("Сохранение в файл\n"+FileName+"\nзавершено.");

return;

}

//-ВЫХОД---------------------------------------------------------------------

void __fastcall TFormMain::N4Click(TObject *Sender)

{ // здесь зделать проверку на наличие изменений и запрос на сохранение

exit(0);

}

//-Вычисления----------------------------------------------------------------

void __fastcall TFormMain::N2Click(TObject *Sender)

{

Fin=EditFInp->Text.ToDouble();

Fout=EditFOut->Text.ToDouble();

theta=new double[n];

J=new long[n]; // на всякий случай берем по полной а не nCurrent

for(long j=0; j<n; j++)

J[j]=-1;

jInModel=new bool[n]; // текущая модель

for(long j=0; j<n; j++)

jInModel[j]=false;

if(CheckBoxSvob->Checked) // вкл. СЧ изначально

{

AddRegressor(0);

Memo1->Lines->Add(

"Оценка параметров нулевого уравнения регрессии.");

ThetaForInModel();

if(OutPutResult("Нулевое")==0)

{

StrMess.sprintf("RSS = %G",RSSForInModel());

Memo1->Lines->Add(StrMess);

}

Memo1->Lines->Add("");

}

// добавление регрессора, максимально коррелированного с откликом

if(CheckBoxMaxCorr->Checked)

{

double Max_rXjY;

long IndexOfMax_rXjY;

MaxCorrellation(Max_rXjY,IndexOfMax_rXjY);

AddRegressor(IndexOfMax_rXjY); // добавление "первого" регрессора

StrMess.sprintf("В модель добавлен регрессор номер %d, максимально "

"коррелированный с откликом (rXjY=%G).",IndexOfMax_rXjY,Max_rXjY);

Memo1->Lines->Add(StrMess);

Memo1->Lines->Add(

"Оценка параметров первого уравнения регрессии.");

ThetaForInModel();

if(OutPutResult("Первое")==0)

{

StrMess.sprintf("RSS = %G",RSSForInModel());

Memo1->Lines->Add(StrMess);

}

Memo1->Lines->Add("");

}

bool FlagEnd1=false, FlagEnd2=false;

while(1)

{ // определение очередного регрессора на добавление

double Max_Fj=0.0;

long IndexOfMax_Fj=-1;

// если больше нет возможности добавить хотя бы один

// регрессор, то вычисления закончились.

if(WhichAddRegressor(Max_Fj,IndexOfMax_Fj)==-1)

FlagEnd1=true; // нет для добавления сейчас

else

{

AddRegressor(IndexOfMax_Fj);

StrMess.sprintf("В уравнение добавлен регрессор № %d.",

IndexOfMax_Fj);

Memo1->Lines->Add(StrMess);

Memo1->Lines->Add(

"Оценка параметров нового уравнения регрессии.");

ThetaForInModel();

if(OutPutResult("Новое")==0)

{ StrMess.sprintf("RSS = %G",RSSForInModel());

Memo1->Lines->Add(StrMess);

}

Memo1->Lines->Add("");

}

// определение регрессоров на исключение

while(1)

{ double Bad_Fj;

long IndexOfBad_Fj;

WhichExpelRegressor(Bad_Fj,IndexOfBad_Fj);

if(IndexOfBad_Fj==-1)

{ FlagEnd2=true;

break; // исключать нечего

}

// иначе исключить

ExpelRegressor(IndexOfBad_Fj);

StrMess.sprintf("Из уравнения исключен регрессор № %d.",

IndexOfBad_Fj);

Memo1->Lines->Add(StrMess);

Memo1->Lines->Add(

"Оценка параметров нового уравнения регрессии.");

ThetaForInModel();

if(OutPutResult("Новое")==0)

{ StrMess.sprintf("RSS = %G",RSSForInModel());

Memo1->Lines->Add(StrMess);

}

Memo1->Lines->Add("");

}

if(FlagEnd1 && FlagEnd2) break;

FlagEnd1=false;

FlagEnd2=false;

/* // если больше не будет возможности добавить хоть один

// регрессор, то выйти из цикла.

if(NotCouldAddRegressor()) break;*/

}

if(CheckBox1->Checked) // вкл все регрессоры

{

for(long i=0; i<n; i++)

AddRegressor(i);

Memo1->Lines->Add("Принудительно включены все регрессоры.");

Memo1->Lines->Add("");

}

ThetaForInModel(); // оценка параметров определенного регрессионного уравнения

if(OutPutResult("Полученное")==0)

{

RSS=RSSForInModel();

StrMess.sprintf("RSS = %G",RSS);

Memo1->Lines->Add(StrMess);

}

Memo1->Lines->Add("=====================================================================");

N3->Enabled=true;

}

//-Определение фактора, максимально коррелированного с откликом--------------

long TFormMain::MaxCorrellation(double & _Max_rXjY, long & _IndexOfMax_rXjY)

{

// определение средних значений:

// отклика

double MiddleY=0.0;

for(long i=0;i<m;i++)

MiddleY+=y[i];

MiddleY/=m;

// столбца Xj в цикле df-ия max корреляции:

double Max_rXjY=0.0;

double IndexOfMax_rXjY=-1;

if(CheckBoxCorrel->Checked)

Memo1->Lines->Add("Вектор корреляций с откликом:");

for(long j=1; j<n; j++)

{ // if(j!=0) // для своб.чл.

double MiddleXj=0;

for(long i=0;i<m;i++)

MiddleXj+=XT[j][i];

MiddleXj/=m;

double S1=0.0, S2=0.0, S3=0.0;

for(long i=0; i<m; i++)

{ S1+=(XT[j][i]-MiddleXj)*(y[i]-MiddleY);

S2+=pow(XT[j][i]-MiddleXj,2);

S3+=pow(y[i]-MiddleY,2);

}

AnsiString StrMess;

if(S1==0.0 || S2==0.0 || S3==0.0)

{ StrMess=StrMess.sprintf("Не посчитан %d", j);

Memo1->Lines->Add(StrMess);

continue;

}

double rXjY=S1/sqrt(S2*S3);

if(CheckBoxCorrel->Checked)

{

StrMess.cat_sprintf("j=%d\t%f\t", j,rXjY);

Memo1->Lines->Add(StrMess);

}

rXjY=(rXjY<0)?-rXjY:rXjY;

if(rXjY>Max_rXjY)

{ Max_rXjY=rXjY;

IndexOfMax_rXjY=j;

}

}

// Максимальный коэффициент корреляции Max_rXjY

// Номер регрессора с максимальным коэффициентом корреляции

// относительно отклика

_Max_rXjY=Max_rXjY;

_IndexOfMax_rXjY=IndexOfMax_rXjY;

return 0;

}

//---------------------------------------------------------------------------

int TFormMain::WhichAddRegressor(double & _Max_Fj, long & _IndexOfMax_Fj)

{

double Max_Fj=0.0;

long IndexOfMax_Fj=-1;

ThetaForInModel();

RSS=RSSForInModel();

double RSSj;

double Fj;

Memo1->Lines->Add("НА ДОБАВЛЕНИЕ");

// по тем, которых нет в модели

for(long j=0; j<n; j++)

if(!jInModel[j])

{

AddRegressor(j); // временно

ThetaForInModel();

RSSj=RSSForInModel();

Fj=(RSS-RSSj)/(RSSj)*(m-(nCurrent-1));

if(Fj>Max_Fj)

{ Max_Fj=Fj;

IndexOfMax_Fj=j;

}

StrMess.sprintf("RSS=%G\tRSSj=%G\tFj=%G\tj=%d"

,RSS,RSSj,Fj,j);

Memo1->Lines->Add(StrMess);

ExpelRegressor(j); // временно

ThetaForInModel();

}

if(Max_Fj<Fin)

{ _Max_Fj=0.0;

_IndexOfMax_Fj=-1;

Memo1->Lines->Add("добавления не будет...");

Memo1->Lines->Add("");

return -1;

}

_Max_Fj=Max_Fj;

_IndexOfMax_Fj=IndexOfMax_Fj;

Memo1->Lines->Add("будет добавление...");

Memo1->Lines->Add("");

return 0;

}

//---------------------------------------------------------------------------

int TFormMain::WhichExpelRegressor(double & _Bad_Fj, long & _IndexOfBad_Fj)

{

double Bad_Fj=0.0;

double IndexOfBad_Fj=-1;

ThetaForInModel();

RSS=RSSForInModel();

double RSSj;

double Fj;

Memo1->Lines->Add("НА ИСКЛЮЧЕНИЕ");

// по тем, которые в модели

for(long j=0; j<n; j++)

if(jInModel[j])

{

ExpelRegressor(j); // временно

ThetaForInModel();

RSSj=RSSForInModel();

Fj=(RSSj-RSS)/(RSSj)*(m-nCurrent);

if(Fj<Bad_Fj)

{ Bad_Fj=Fj;

IndexOfBad_Fj=j;

}

StrMess.sprintf("RSS=%G\tRSSj=%G\tFj=%G\tj=%d",RSS,RSSj,Fj,j);

Memo1->Lines->Add(StrMess);

AddRegressor(j); // временно

ThetaForInModel();

}

if(Bad_Fj>Fout)

{ _Bad_Fj=0.0;

_IndexOfBad_Fj=-1;

Memo1->Lines->Add("исключения не будет...");

Memo1->Lines->Add("");

return -1;

}

_Bad_Fj=Bad_Fj;

_IndexOfBad_Fj=IndexOfBad_Fj;

Memo1->Lines->Add("исключение будет...");

Memo1->Lines->Add("");

return 0;

}

//---------------------------------------------------------------------------

int TFormMain::ThetaForInModel()

{

// Выделение памяти для матрицы СЛАУ

XTX=new double*[nCurrent]; // столбец указателей

for(long i=0; i<nCurrent; i++) // создание строк

XTX[i]=new double[nCurrent];

// Получение XTX=XT*X (перепроверить!!!)

for(long i=0; i<nCurrent; i++)

for(long j=0; j<nCurrent; j++)

{

XTX[i][j]=0.0;

for(long k=0; k<m; k++)

XTX[i][j]+=XT[J[i]][k]*XT[J[j]][k];

}

// создание правой части СЛАУ как произведение XT*Y

f=new double[nCurrent];

for(long i=0; i<nCurrent; i++)

{ f[i]=0.0;

for(long j=0; j<m; j++)

f[i]+=XT[J[i]][j]*y[j];

}

ThetaForInModelSlaeResolve();

for(long j=0; j<n; j++)

theta[j]=0.0;

for(long j=0; j<nCurrent; j++)

theta[J[j]]=f[j];

for(long j=0; j<nCurrent; j++)

delete [] XTX[j];

delete [] XTX;

delete [] f;

return 0;

}

//---------------------------------------------------------------------------

int TFormMain::ThetaForInModelSlaeResolve()

{ // LDU

for(long np=0; np<nCurrent; np++)

{ // диагональный:

for(long k=0; k<np; k++)

XTX[np][np]-=XTX[np][k]*XTX[k][k]*XTX[k][np];

// нижняя часть

long j=np;

for(long i=np+1; i<nCurrent; i++)

{

for(long k=0; k<j; k++)

XTX[i][j]-=XTX[i][k]*XTX[k][k]*XTX[k][j];

XTX[i][j]/=XTX[j][j];

}

// верхняя часть

long i=np;

for(long j=np+1; j<nCurrent; j++)

{

for(long k=0; k<i; k++)

XTX[i][j]-=XTX[i][k]*XTX[k][k]*XTX[k][j];

XTX[i][j]/=XTX[i][i];

}

}

// прямой ход:

for(long i=0; i<nCurrent; i++)

for(long j=0; j<i; j++)

f[i]-=XTX[i][j]*f[j];

// диагональный ход:

for(long i=0; i<nCurrent; i++)

f[i]/=XTX[i][i];

// обратный ход:

for(long i=nCurrent-1; i>=0; i--)

for(long j=i+1; j<nCurrent; j++)

f[i]-=XTX[i][j]*f[j];

return 0;

}

//---------------------------------------------------------------------------

int TFormMain::OutPutResult(AnsiString Mes)

{

AnsiString Str0;

if(Mes=="")

Str0="Уравнение регрессии: ";

else Str0=Mes+" уравнение регрессии: ";

AnsiString Str1="y ";

if(nCurrent>1 || !jInModel[0])

Str1.cat_sprintf("(");

AnsiString Str2;

bool FlagFirstPrintArgument=true;

for(long j=1; true; j++)

{ if(jInModel[j])

{ if(!FlagFirstPrintArgument)

Str2+=StrMess.sprintf(", x%d",j);

else

{ Str2+=StrMess.sprintf("x%d",j);

FlagFirstPrintArgument=false;

}

}

if(j<n-1);

else break;

}

AnsiString Str3;

if(nCurrent>1 || !jInModel[0])

Str3.cat_sprintf(") = ");

else

Str3="= ";

AnsiString Str4;

bool FlagFirstPrintSlagaemoe=true;

if(jInModel[0])

{ Str4+=StrMess.sprintf("%G",theta[0]);

FlagFirstPrintSlagaemoe=false;

}

for(long j=1; true; j++)

{ if(jInModel[j])

{

if(!FlagFirstPrintSlagaemoe)

Str4+=StrMess.sprintf(" %+G*x%d",theta[j],j);

else

{ Str4+=StrMess.sprintf("%G*x%d",theta[j],j);

FlagFirstPrintSlagaemoe=false;

}

}

if(j<n-1);

else break;

}

if(FlagFirstPrintArgument && FlagFirstPrintSlagaemoe)

{ StrMess="Уравнение неопределено.";

Memo1->Lines->Add(StrMess);

return 1;

}

else

{

Memo1->Lines->Add(Str0);

StrMess=Str1+Str2+Str3+Str4;

Memo1->Lines->Add(StrMess);

return 0;

}

}

//---------------------------------------------------------------------------

void __fastcall TFormMain::ComboBoxFontSizeChange(TObject *Sender)

{

Memo1->Font->Size=ComboBoxFontSize->Text.ToInt();

Memo1->Refresh();

}

//---------------------------------------------------------------------------

double TFormMain::RSSForInModel()

{

double *raznost;

raznost=new double[m];

// умножение матрицы X на вектор параметров theta

double rss=0;

for(long i=0; i<m; i++)

{ raznost[i]=0;

for(long j=0; j<nCurrent; j++)

raznost[i]-=XT[J[j]][i]*theta[J[j]];

// вычитание полученного из отклика

raznost[i]+=y[i];

rss+=pow(raznost[i],2);

}

delete [] raznost;

return rss;

}

//---------------------------------------------------------------------------

int TFormMain::AddRegressor(long NumOfAddRegressor)

{

jInModel[NumOfAddRegressor]=true;

// Подсчет количества регрессоров

long nCurr=0;

for(long j=0; j<n; j++)

if(jInModel[j]) nCurr+=1;

// Перепись массива номеров включенных регрессоров

long j1=0;

for(long j=0; j<n; j++)

if(jInModel[j])

{ J[j1]=j;

j1+=1;

}

// заполнение оставшейся части массива J

for(;j1<n; j1++)

J[j1]=-1;

nCurrent=nCurr; // ~nCurrent+=1;

return 0;

}

//---------------------------------------------------------------------------

int TFormMain::ExpelRegressor(long NumOfExpelRegressor)

{

jInModel[NumOfExpelRegressor]=false;

// Подсчет количества регрессоров

long nCurr=0;

for(long j=0; j<n; j++)

if(jInModel[j]) nCurr+=1;

// Перепись массива номеров включенных регрессоров

long j1=0;

for(long j=0; j<n; j++)

if(jInModel[j])

{ J[j1]=j;

j1+=1;

}

// заполнение оставшейся части массива J

for(;j1<n; j1++)

J[j1]=-1;

nCurrent=nCurr; // ~nCurrent-=1;

return 0;

}

//---------------------------------------------------------------------------

void __fastcall TFormMain::ButtonMemoClearClick(TObject *Sender)

{

Memo1->Clear();

}

//---------------------------------------------------------------------------

Размещено на Allbest.ru