Skip to content

Instantly share code, notes, and snippets.

@MadFaill
Forked from sukamenev/least_squares.go
Created April 21, 2024 06:54
Show Gist options
  • Save MadFaill/a24983638c79d9b6e319460fab401620 to your computer and use it in GitHub Desktop.
Save MadFaill/a24983638c79d9b6e319460fab401620 to your computer and use it in GitHub Desktop.
Golang least squares method: traditional and minimizing sum of squares of perpendicular lengths approach.
/*
Код взят отсюда и доработан
https://gist.github.com/ota42y/db4ff0298d9c945cd261
Код применим только для y = ax + b
Для вертикальных прямых неприменим.
Least Squares Fitting: Perpendicular Offsets
Мы используем продвинутую версию метода наименьших квадратов, где меряется расстояние перпендикуляра до прямой
https://mathworld.wolfram.com/LeastSquaresFittingPerpendicularOffsets.html
http://web.archive.org/web/20160322074618/http://cgm.computergraphics.ru/content/view/41
Наивная реализация метода наименьших квадратов
http://www.cleverstudents.ru/articles/mnk.html
С другой стороны, если нам нужно делать предсказания, то метод наименьших квадратов вроде как подходит,
так как он минимизирует разницу по оси Y.
https://habr.com/ru/post/514818/
*/
package main
import (
"fmt"
"math"
)
type Point struct {
x float64
y float64
}
// Traditional least squares method
// y = a*x + b
// Это японская версия
// Speed 48.7M/sec
func LeastSquaresMethod(points *[]Point) (a, b float64) {
n := float64(len(*points))
sumX := 0.0
sumY := 0.0
sumXY := 0.0
sumXX := 0.0
for _, p := range *points{
sumX += p.x
sumY += p.y
sumXY += p.x * p.y
sumXX += p.x * p.x
}
/*
По методу Крамера находим корни:
base := (n * sumXX - sumX * sumX)
a = (n * sumXY - sumX * sumY) / base
b = (sumXX * sumY - sumXY * sumX) / base
*/
// Метод подстановки. Более быстрый. На одну операцию умножения меньше.
a = (n * sumXY - sumY * sumX) / (n * sumXX - sumX * sumX)
b = (sumY - a * sumX) / n
return a, b
}
// Least squares - minimum of squares distances of squares of lens perpendiculars to line.
// Least Squares Fitting: Perpendicular Offsets
// Технику взял отсюда
// https://mathworld.wolfram.com/LeastSquaresFittingPerpendicularOffsets.html
// Скорость 15.8М в секунду
// Скорость 13.7М в секунду c вызовом SumPerpDistance(points, a1, b1)
func LeastSquaresDist(points *[]Point) (float64, float64) {
n := float64(len(*points))
sumX := 0.0
sumY := 0.0
sumXY := 0.0
sumXX := 0.0
sumYY := 0.0
// Для проверки на фиксированность Y
countYConst := 0
firstY := (*points)[0].y
for _, p := range *points{
if p.y == firstY {
countYConst++
}
sumX += p.x
sumY += p.y
sumXY += p.x * p.y
sumXX += p.x * p.x
sumYY += p.y * p.y
}
// Если все Y одинаковы, то мы знаем ответ
if countYConst == len(*points) {
return 0.0, firstY
}
sumXn := sumX / n
sumYn := sumY / n
numerator := (sumYY - sumY * sumYn) - (sumXX - sumX * sumXn)
denominator := sumX * sumYn - sumXY
// Для тех случаев, когда уравнение не решается.
// Если все X равны константе, то метод неприменим.
if denominator == 0.0 {
return LeastSquaresMethod(points)
}
B := 0.5 * numerator / denominator
D := math.Sqrt(B*B + 1)
a1 := -B + D
a2 := -B - D
// a находим используя подстановку
b1 := sumYn - a1 * sumXn
b2 := sumYn - a2 * sumXn
/*
if SumPerpDistance(points, a1, b1) < SumPerpDistance(points, a2, b2) {
return a1, b1
}*/
// Определяем абсолютный минимум из двух локальных
sum1 := 0.0
sum2 := 0.0
for _, p := range *points{
// Квадрат нельзя использовать. в 7 раз медленнее. В проде, наверное, лучше обойтись просто длиной...
sum1 += math.Pow(p.y - b1 - a1 * p.x, 2)
sum2 += math.Pow(p.y - b2 - a2 * p.x, 2)
}
sum1 = sum1 / (a1*a1 + 1)
sum2 = sum2 / (a2*a2 + 1)
if sum1 < sum2 {
return a1, b1
}
return a2, b2
}
func SumPerpDistance(points *[]Point, a, b float64) float64 {
sum := 0.0
// more fast
// Enough to choose the best extremum
/*
for _, p := range *points{
sum += math.Abs(p.y - b - a * p.x)
}
*/
// more correctness
// lower speed
for _, p := range *points{
sum += math.Pow(p.y - b - a * p.x, 2)
}
sum = sum / (a*a + 1)
return sum
}
func main() {
points := make([]Point, 0)
// Example for testing equality traditional and perpendicular offset approach
points = append(points, Point{x:0.0, y:1.0,})
points = append(points, Point{x:0.1, y:1.5,})
points = append(points, Point{x:0.2, y:2.0,})
points = append(points, Point{x:0.3, y:2.5,})
points = append(points, Point{x:0.4, y:3.0,})
points = append(points, Point{x:0.5, y:3.5,})
points = append(points, Point{x:0.6, y:4.0,})
/* Example with more accuracy in Least Squares Perpendicular Offset approach
// See plot on https://www.desmos.com/calculator/gmlhcrwg05?lang=ru
// Blue (LST), Red (LSPO)
points = append(points, Point{x:0.0, y:1.0,})
points = append(points, Point{x:1.0, y:0.0,})
points = append(points, Point{x:2.0, y:2.50,})
points = append(points, Point{x:3.0, y:5.0,})
points = append(points, Point{x:4.0, y:4.0,})
points = append(points, Point{x:5.0, y:6.0,})
points = append(points, Point{x:6.0, y:9.0,})
*/
a, b := LeastSquaresMethod(&points)
fmt.Println("Least Squares traditional")
fmt.Println(a)
fmt.Println(b)
fmt.Println("Distance (sum of squares perp) ", SumPerpDistance(&points, a, b))
fmt.Println("\nLeast Squares Fitting: Perpendicular Offsets")
a1, b1 := LeastSquaresDist(&points)
fmt.Println("\nBest root")
fmt.Println(a1)
fmt.Println(b1)
fmt.Println("Distance (sum of squares perp) ", SumPerpDistance(&points, a1, b1))
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment