# Tech Tips

1. Uncategorized
2. 49 view

# Least Square Method

20140706 zuqqhi2-lsm-v1 from Hidetomo Suzuki

# Program

```#!/bin/bash

for num in `seq 2 31`
do
go run lsm-polynomial.go \$num
done
```
```package main

import (
"image/color"
"math/rand"
"math"
"strconv"
"os"
"fmt"
)

// Make random value with normal distribution by Box-Muller Transform
func normalRand(mu, sigma float64) float64 {
z := math.Sqrt(-2.0 * math.Log(rand.Float64())) * math.Sin(2.0 * math.Pi * rand.Float64())
return sigma*z + mu
}

// Make sequence of numbers with common difference
func linspace(start, end float64, n int, x plotter.XYs) {
for i := 0; i &amp;amp;lt; n; i++ {
t := float64(i) / float64(n-1)
x[i].X = (1.0 - t) * start + t * end
}
}

// Elementary matrix with Gauss-Jordan Method
func gaussElimination(mat [][]float64, numParams int) []float64 {
// Pivoting
for x := 0; x &amp;amp;lt; numParams; x++ {
pivot := x
maxVal := 0.0
for y := x; y &amp;amp;lt; numParams; y++ {
if math.Abs(mat[y][x]) &amp;amp;gt; maxVal {
maxVal = math.Abs(mat[y][x])
pivot = y
}
}

// switch column when pivot != x
if pivot != x {
for k := 0; k &amp;amp;lt; numParams+1; k++ {
tmp           := mat[x][k]
mat[x][k]     =  mat[pivot][k]
mat[pivot][k] =  tmp
}
}
}

// Make upper triangular matrix with basic deformation of matrix
for k := 0; k &amp;amp;lt; numParams; k++ {
// Save diagonal element
diagElem := mat[k][k]
// Diagonal element must be 1.0
mat[k][k] = 1.0
// Basic deformation
for x := k+1; x &amp;amp;lt; numParams+1; x++ { mat[k][x] /= diagElem  }
for y := k+1; y &amp;amp;lt; numParams; y++ {
tmp := mat[y][k]
for x := k+1; x &amp;amp;lt; numParams+1; x++ {
mat[y][x] -= tmp*mat[k][x]
}
mat[y][k] = 0.0
}
}

// Solve
result := make([]float64, numParams, numParams)
for y := numParams-1; y &amp;amp;gt;= 0; y-- {
result[y] = mat[y][numParams]
for x := numParams-1; x &amp;amp;gt; y; x-- {
result[y] -= mat[y][x] * result[x]
}
}

return result
}

func main() {
//===================================================
// Check Command Line Arguments

numParams, err := strconv.Atoi(os.Args[1])
if err != nil {
panic(err)
}

//===================================================
// Make observed points

rand.Seed(int64(0))

// Prepare X axis of observed points
n := 50

// make observed points
pix := make([]float64, n)
for i := 0; i &amp;amp;lt; n; i++ {
}
for i := 0; i &amp;amp;lt; n; i++ {
answer[i].Y = math.Sin(pix[i]) / pix[i] + 0.1 * answer[i].X + normalRand(1.0, 0.05)
}

//====================================================
// LSM

// Make a matrix to solve
mat := make([][]float64, numParams, numParams)
fmt.Printf("[Matrix for Gausss Elimination]\n")
for y := 0; y &amp;amp;lt; numParams; y++ {
mat[y] = make([]float64, numParams+1, numParams+1)
for x := 0; x &amp;amp;lt; numParams+1; x++ {
mat[y][x] = 0.0
if x &amp;amp;lt; numParams {
for k := 0; k &amp;amp;lt; n; k++ {
}
} else {
for k := 0; k &amp;amp;lt; n; k++ {
}
}
fmt.Printf("%f\t", mat[y][x])
}
fmt.Printf("\n")
}
fmt.Printf("\n")

// Solve y=phi*theta
params := gaussElimination(mat, numParams)

// Output
fmt.Printf("[Optimized Params]\n")
for i := 0; i &amp;amp;lt; numParams; i++ {
fmt.Printf("param %d : %f\n", i, params[i])
}

// Calculation Error
sumError := 0.0
for i := 0; i &amp;amp;lt; n; i++ {
y := 0.0
for k := 0; k &amp;amp;lt; numParams; k++ {
}

sumError += math.Pow(answer[i].Y - y, 2.0)
}

sumError = 0.0
for i := 0; i &amp;amp;lt; n; i++ {
// New Data
x := rand.Float64() * 6.0 - 3.0
ansY := math.Sin(math.Pi * x) / (math.Pi * x) + 0.1 * x + normalRand(1.0, 0.05)

// Estimation
y := 0.0
for k := 0; k &amp;amp;lt; numParams; k++ {
y += params[k]*math.Pow(x, float64(k))
}

// Error
sumError += math.Pow(ansY - y, 2.0)
}

//====================================================
// Graph Setting

// Result graph
N := 1000
result := make(plotter.XYs, N)
linspace(-3, 3, N, result)
for i := 0; i &amp;amp;lt; N; i++ {
result[i].Y = 0.0
for k := 0; k &amp;amp;lt; numParams; k++ {
result[i].Y += params[k]*math.Pow(result[i].X, float64(k))
}
}

// Create a new plot, set its title and axis labels
p, err := plot.New()
if err != nil {
panic(err)
}
p.Title.Text = "Least Square Method"
p.X.Label.Text = "X"
p.Y.Label.Text = "Y"

// Make a scatter plotter and set its style
// Make a line plotter with points and set its style.
if err != nil {
panic(err)
}
lpLine.Color = color.RGBA{G: 255, A: 255}
lpPoints.Shape = plot.PyramidGlyph{}
lpPoints.Color = color.RGBA{R: 255, A: 255}

// Result line
lpResultLine, _, err := plotter.NewLinePoints(result)
if err != nil {
panic(err)
}
lpResultLine.Color = color.RGBA{B: 255, A: 255}

// Save the plot to a PNG file.
filename := "lsm-polynomial" + strconv.Itoa(numParams) + ".png"
if err := p.Save(4, 4, filename); err != nil {
panic(err)
}
}
```

# Result

In this time, looks n=10(go run lsm-polynomial.go 11) is best because generalization error is minimum in this simulation.