Least Squares-Based Linear Regression

"Linear equation that best fits the data"

    We use opencv library to plot the results of the algorithm. This shows how to implement these codes in this page.

#include <iostream> #include "opencv4/opencv2/opencv.hpp" #include <math.h> using namespace std; using namespace cv;
import cv2 from random import randint import numpy as np from math import atan

    The variables below are used in the algorithm, as well as a canvas image for plotting the results. The "canvas.png" can be downloaded here.

const int nSample = 100; const int maxSize = 1000; int particles_x[nSample]; int particles_y[nSample]; Mat canvas = imread("canvas.png"); resize(canvas, canvas, Size(maxSize, maxSize)); float sumX = 0.0; float sumY = 0.0; float sumXY = 0.0; float sumX2 = 0.0;
nSample = 100 maxSize = 1000 particles_x = np.array([]) particles_y = np.array([]) canvas = cv2.imread("canvas.png") canvas = cv2.resize(canvas, (maxSize, maxSize)) sumX = 0.0 sumY = 0.0 sumXY = 0.0 sumX2 = 0.0

    Then, we need to generate random samples or particles as the input for the algorithm. This can be changed to a specific problem that we want to solve.

srand((unsigned int)time(NULL)); for (int i=0; i<nSample; i++) { particles_x[i] = round((float)rand()/RAND_MAX*maxSize); particles_y[i] = 300 + round((float)rand()/RAND_MAX*400); circle(canvas, Point(particles_x[i], particles_y[i]), 10, Scalar(128,128,128), -1); }
for i in range(nSample): rnd_x = randint(0, maxSize) rnd_y = randint(300, 700) particles_x = np.append(particles_x, rnd_x) particles_y = np.append(particles_y, rnd_y) canvas = cv2.circle(canvas, (int(particles_x[i]), int(particles_y[i])), 10, (128,128,128), -1)

    In the main algorithm, we use following equations. Previously, we need to get the sum value (Σ) of x, yxy, and x^2 of all the samples, as the general rule of Least Squares method.

            y = mx + b

            m(n.Σxy - Σx.Σy) / (n.Σx^2 - (Σx)^2)

            b = (Σy - m.Σx) / n

            θ = atan (m)

    Where, m is a slope or gradient of the line, b is y-intercept or y coordinate when x = 0, n is the number of samples, θ is the line angle, x and y are the values that we want to play with. 

for (int i=0; i<nSample; i++) { sumX += particles_x[i]; sumY += particles_y[i]; sumXY += particles_x[i] * particles_y[i]; sumX2 += particles_x[i] * particles_x[i]; } float m = ((nSample * sumXY) - (sumX * sumY)) / ((nSample * sumX2) - (sumX * sumX)); float b = (sumY - (m * sumX)) / nSample; float theta = -atan(m) * 180 / 3.142857; int new_x = 0, new_y = 0; for (int i=0; i<nSample; i++) { if (fabs(theta) <= 45.0) { //horizontal line, find y new_x = particles_x[i]; new_y = (m * new_x) + b; } else { //vertical line, find x new_y = particles_y[i]; new_x = (new_y - b) / m; } line(canvas, Point(new_x, new_y), Point(particles_x[i], particles_y[i]), Scalar(0,0,0), 1); }
for i in range(nSample): sumX += particles_x[i] sumY += particles_y[i] sumXY += particles_x[i] * particles_y[i] sumX2 += particles_x[i] * particles_x[i] m = ((nSample * sumXY) - (sumX * sumY)) / ((nSample * sumX2) - (sumX * sumX)) b = (sumY - (m * sumX)) / nSample theta = -atan(m) * 180 / 3.142857 new_x = 0; new_y = 0 for i in range(nSample): if abs(theta) <= 45.0: #horizontal line, find y new_x = particles_x[i] new_y = (m * new_x) + b else: #vertical line, find x new_y = particles_y[i] new_x = (new_y - b) / m canvas = cv2.line(canvas, (int(new_x), int(new_y)), (int(particles_x[i]), int(particles_y[i])), (0,0,0), 1) #canvas = cv2.circle(canvas, (int(new_x), int(new_y)), 10, (255,0,0), -1)

    The last, plot the final result. It is similar to RANSAC in terms of results. However, each of them has its own pros and cons depending on the use case. Therefore, figuring out the most suitable for our implementation is important.

int min_x = 0; int max_x = maxSize; int min_y = round((m * min_x) + b); int max_y = round((m * max_x) + b); line(canvas, Point(min_x, min_y), Point(max_x, max_y), Scalar(0,0,255), 5); string theta_str = to_string((int)round(theta)); putText(canvas, "Theta: "+to_string((int)round(theta))+"deg. m: "+to_string(m)+" b: "+to_string(b), Point(10,maxSize-10), FONT_HERSHEY_COMPLEX, 1, Scalar(0,0,255), 2); while (1) { imshow("Least Squares-Based Linear Regression @The JPID Coder", canvas); waitKey(1); }
min_x = 0 max_x = maxSize min_y = (m * min_x) + b max_y = (m * max_x) + b m_str = str(round(m * 100.0) / 100.0) b_str = str(round(b * 100.0) / 100.0) canvas = cv2.line(canvas, (int(min_x), int(min_y)), (int(max_x), int(max_y)), (0,0,255), 5) canvas = cv2.putText(canvas, "Theta: "+str(round(theta))+"deg. m: "+m_str+" b: "+b_str, (10,maxSize-10), cv2.FONT_HERSHEY_SIMPLEX, 1, (0,0,255), 2) while True: cv2.imshow("Least Squares-Based Linear Regression @The JPID Coder", canvas) cv2.waitKey(1)

    Please look at for the implementation of this algorithm.

References
[1] Wikipedia, Linear Least Squares.


Home     Algorithms     Applications     Live Simulations     YouTube     About us

Comments

© 2022 The JPID Coder