/*************************************************** * Linear Least Squares Demonstration Program * * ------------------------------------------------ * * Reference: BASIC Scientific Subroutines, Vol. II * * By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].* * * * C++ version by J-P Moreau, Paris * * with possibility of a graph. * * (to be used with lstsqr1.mak and graph2d.cpp) * * (www.jpmoreau.fr) * * ------------------------------------------------ * * This program calculates a linear least squares * * fit to a given data set. * * * * INSTRUCTIONS * * ------------ * * * * The number of data coordinates provided must be * * greater than two. Otherwise, a divide by zero * * error may result. * * * * SAMPLE DATA: * * * * Number of data points : 5 * * * * X[0]=1 Y[0]=1.15 * * X[1]=2 Y[1]=2.22 * * X[2]=3 Y[2]=2.94 * * X[3]=4 Y[3]=3.85 * * X[4]=5 Y[4]=5.10 * * * * Fitted equation is: * * Y = 0.193000 + 0.953000 X * * * * Standard deviation of fit: 0.154262 * * * * (A graph is displayed). * ***************************************************/ #include #include #include #include HDC hdc; RECT rect; HPEN hpen, hpenOld; //Functions used here of Graph2D.cpp void Conversion(double,double,int *,int *); void MinMax(int,float *,double *,double *); void Initwindow(int,double,double,double,double); void MoveXY(double,double); void LineXY(double,double); void Legendes (int,char *,char *,char *,int,int,char *,int,int,char *); //"home made" graphic commands for hdc environment used by above functions void DrawPixel(int ix,int iy) { //sorry, no other available command found Rectangle(hdc,rect.left+ix,rect.top+iy, rect.left+ix+2,rect.top+iy+1); } void Swap(int *i1,int *i2) { int it; it=*i1; *i1=*i2; *i2=it; } void DrawLine(int ix1,int iy1,int ix2,int iy2) { int i,il,ix,iy; double dx,dy; if (ix2 2). The returned * * parameters are: a,b, coefficients of equation * * Y = a + b X, and d, standard deviation of fit. * **************************************************/ void Least_Square() { double a1,a2,b0,b1,d1; int m; a1 = 0; a2 = 0; b0 = 0; b1 = 0; for (m = 0; m < n; m++) { a1 = a1 + X[m]; a2 = a2 + X[m] * X[m]; b0 = b0 + Y[m]; b1 = b1 + Y[m] * X[m]; } a1 = a1 / n; a2 = a2 / n; b0 = b0 / n; b1 = b1 / n; d = a1 * a1 - a2; a = a1 * b1 - a2 * b0; a = a / d; b = a1 * b0 - b1; b = b / d; // Evaluation of standard deviation d (unbiased estimate) d = 0; for (m = 0; m < n; m++) { d1 = Y[m] - a - b * X[m]; d += d1 * d1; } d = sqrt(d / (n - 2)); } // Draw a cross at physical point (x,y) void CrossXY(double x,double y) { int jx,jy; Conversion(x,y,&jx,&jy); DrawLine(jx-5,jy,jx+5,jy); DrawLine(jx,jy-5,jx,jy+5); } void Draw_graph() { double xmn,xmx,ymn,ymx; extern MaxX, MaxY; char buf[40],titre[40]; int m; //client drawing zone in pixels MaxX=(rect.right-rect.left); MaxY=(rect.bottom-rect.top); //seek minimum, maximum of table Y MinMax(n,X,&xmn,&xmx); MinMax(n,Y,&ymn,&ymx); Initwindow(10,xmn,xmx,ymn,ymx); //Crosses represent input data points for (m = 0; m0) sprintf(titre," Y = %lf + %lf X ",a,b); else sprintf(titre," Y = %lf %lf X ",a,b); //Display graph caption, names of axes and //standard deviation, d sprintf(buf," Standard deviation = %lf",d); Legendes(10,titre,"X","Y",50,50,buf,0,0,""); } void Lstsqr1() { Data(); Least_Square(); Draw_graph(); } //Handle window Paint and Destroy messages long FAR PASCAL /*_export*/ WndProc(HWND hwnd, UINT message, UINT wParam, LONG lParam) { PAINTSTRUCT ps; switch (message) { case WM_PAINT: hdc = BeginPaint(hwnd, &ps); GetClientRect(hwnd, &rect); hpen = CreatePen(PS_SOLID, 1, RGB(0, 0, 255)); hpenOld = SelectObject(hdc,hpen); Lstsqr1(); SelectObject(hdc,hpenOld); DeleteObject(hpen); EndPaint(hwnd, &ps); return 0; case WM_DESTROY: PostQuitMessage(0); return 0; } return DefWindowProc(hwnd, message, wParam, lParam); } // main program int PASCAL WinMain(HANDLE hInstance, HANDLE hPrevInstance, LPSTR lpszCmdParam, int nCmdShow) { static char szAppName[] = "Lstsqr1"; HWND hwnd; MSG msg; WNDCLASS wndclass; if (!hPrevInstance) { wndclass.style = CS_HREDRAW | CS_VREDRAW; wndclass.lpfnWndProc = WndProc; wndclass.cbClsExtra = 0; wndclass.cbWndExtra = 0; wndclass.hInstance = hInstance; wndclass.hIcon = LoadIcon(NULL, IDI_APPLICATION); wndclass.hCursor = LoadCursor(NULL, IDC_ARROW); wndclass.hbrBackground = GetStockObject(WHITE_BRUSH); wndclass.lpszMenuName = NULL; wndclass.lpszClassName = szAppName; RegisterClass(&wndclass); } hwnd = CreateWindow(szAppName, // window class name " LINEAR LEAST SQUARES", // window caption WS_OVERLAPPEDWINDOW, // window style CW_USEDEFAULT, // initial x position CW_USEDEFAULT, // initial y position CW_USEDEFAULT, // initial x size CW_USEDEFAULT, // initial y size NULL, // parent window handle NULL, // window menu handle hInstance, // program instance handle NULL); // creation parameters ShowWindow(hwnd, nCmdShow); UpdateWindow(hwnd); while (GetMessage(&msg, NULL, 0, 0)) { TranslateMessage(&msg); DispatchMessage(&msg); } return msg.wParam; } //end of file lstsqr1.cpp