-
-
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.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* | |
Код взят отсюда и доработан | |
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