Golang least squares method: traditional and minimizing sum of squares of perpendicular lengths approach.

least_squares.go

/*
 Код взят отсюда и доработан
 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))
}