Documentation ¶
Overview ¶
Package fit provides functions to fit data.
Index ¶
Examples ¶
Constants ¶
This section is empty.
Variables ¶
This section is empty.
Functions ¶
func Curve1D ¶
Curve1D returns the result of a non-linear least squares to fit a function f to the underlying data with method m.
Example (Exponential) ¶
package main import ( "bufio" "image/color" "log" "math" "os" "strconv" "strings" "go-hep.org/x/hep/fit" "go-hep.org/x/hep/hplot" "gonum.org/v1/gonum/floats" "gonum.org/v1/gonum/optimize" "gonum.org/v1/plot/plotter" "gonum.org/v1/plot/vg" ) func main() { const ( a = 0.3 b = 0.1 ndf = 2.0 ) xdata, ydata, err := readXY("testdata/exp-data.txt") if err != nil { log.Fatal(err) } exp := func(x, a, b float64) float64 { return math.Exp(a*x + b) } res, err := fit.Curve1D( fit.Func1D{ F: func(x float64, ps []float64) float64 { return exp(x, ps[0], ps[1]) }, X: xdata, Y: ydata, N: 2, }, nil, &optimize.NelderMead{}, ) if err != nil { log.Fatal(err) } if err := res.Status.Err(); err != nil { log.Fatal(err) } if got, want := res.X, []float64{a, b}; !floats.EqualApprox(got, want, 0.1) { log.Fatalf("got= %v\nwant=%v\n", got, want) } { p := hplot.New() p.X.Label.Text = "exp(a*x+b)" p.Y.Label.Text = "y-data" p.Y.Min = 0 p.Y.Max = 5 p.X.Min = 0 p.X.Max = 5 s := hplot.NewS2D(hplot.ZipXY(xdata, ydata)) s.Color = color.RGBA{0, 0, 255, 255} p.Add(s) f := plotter.NewFunction(func(x float64) float64 { return exp(x, res.X[0], res.X[1]) }) f.Color = color.RGBA{255, 0, 0, 255} f.Samples = 1000 p.Add(f) p.Add(plotter.NewGrid()) err := p.Save(20*vg.Centimeter, -1, "testdata/exp-plot.png") if err != nil { log.Fatal(err) } } } func readXY(fname string) (xs, ys []float64, err error) { f, err := os.Open(fname) if err != nil { return xs, ys, err } defer f.Close() scan := bufio.NewScanner(f) for scan.Scan() { line := scan.Text() toks := strings.Split(line, " ") x, err := strconv.ParseFloat(toks[0], 64) if err != nil { return xs, ys, err } xs = append(xs, x) y, err := strconv.ParseFloat(toks[1], 64) if err != nil { return xs, ys, err } ys = append(ys, y) } return }
Output:
Example (Gaussian) ¶
package main import ( "bufio" "image/color" "log" "math" "os" "strconv" "strings" "go-hep.org/x/hep/fit" "go-hep.org/x/hep/hplot" "gonum.org/v1/gonum/floats" "gonum.org/v1/gonum/optimize" "gonum.org/v1/plot/plotter" "gonum.org/v1/plot/vg" ) func main() { var ( cst = 3.0 mean = 30.0 sigma = 20.0 want = []float64{cst, mean, sigma} ) xdata, ydata, err := readXY("testdata/gauss-data.txt") if err != nil { log.Fatal(err) } gauss := func(x, cst, mu, sigma float64) float64 { v := (x - mu) return cst * math.Exp(-v*v/sigma) } res, err := fit.Curve1D( fit.Func1D{ F: func(x float64, ps []float64) float64 { return gauss(x, ps[0], ps[1], ps[2]) }, X: xdata, Y: ydata, Ps: []float64{10, 10, 10}, }, nil, &optimize.NelderMead{}, ) if err != nil { log.Fatal(err) } if err := res.Status.Err(); err != nil { log.Fatal(err) } if got := res.X; !floats.EqualApprox(got, want, 1e-3) { log.Fatalf("got= %v\nwant=%v\n", got, want) } { p := hplot.New() p.X.Label.Text = "Gauss" p.Y.Label.Text = "y-data" s := hplot.NewS2D(hplot.ZipXY(xdata, ydata)) s.Color = color.RGBA{0, 0, 255, 255} p.Add(s) f := plotter.NewFunction(func(x float64) float64 { return gauss(x, res.X[0], res.X[1], res.X[2]) }) f.Color = color.RGBA{255, 0, 0, 255} f.Samples = 1000 p.Add(f) p.Add(plotter.NewGrid()) err := p.Save(20*vg.Centimeter, -1, "testdata/gauss-plot.png") if err != nil { log.Fatal(err) } } } func readXY(fname string) (xs, ys []float64, err error) { f, err := os.Open(fname) if err != nil { return xs, ys, err } defer f.Close() scan := bufio.NewScanner(f) for scan.Scan() { line := scan.Text() toks := strings.Split(line, " ") x, err := strconv.ParseFloat(toks[0], 64) if err != nil { return xs, ys, err } xs = append(xs, x) y, err := strconv.ParseFloat(toks[1], 64) if err != nil { return xs, ys, err } ys = append(ys, y) } return }
Output:
Example (Poly) ¶
package main import ( "image/color" "log" "math/rand" "go-hep.org/x/hep/fit" "go-hep.org/x/hep/hplot" "gonum.org/v1/gonum/floats" "gonum.org/v1/gonum/optimize" "gonum.org/v1/plot/plotter" "gonum.org/v1/plot/vg" ) func main() { var ( a = 1.0 b = 2.0 ps = []float64{a, b} want = []float64{1.38592513, 1.98485122} // from scipy.curve_fit ) poly := func(x float64, ps []float64) float64 { return ps[0] + ps[1]*x*x } xdata, ydata := genXY(100, poly, ps, -10, 10) res, err := fit.Curve1D( fit.Func1D{ F: poly, X: xdata, Y: ydata, Ps: []float64{1, 1}, }, nil, &optimize.NelderMead{}, ) if err != nil { log.Fatal(err) } if err := res.Status.Err(); err != nil { log.Fatal(err) } if got := res.X; !floats.EqualApprox(got, want, 1e-6) { log.Fatalf("got= %v\nwant=%v\n", got, want) } { p := hplot.New() p.X.Label.Text = "f(x) = a + b*x*x" p.Y.Label.Text = "y-data" p.X.Min = -10 p.X.Max = +10 p.Y.Min = 0 p.Y.Max = 220 s := hplot.NewS2D(hplot.ZipXY(xdata, ydata)) s.Color = color.RGBA{0, 0, 255, 255} p.Add(s) f := plotter.NewFunction(func(x float64) float64 { return poly(x, res.X) }) f.Color = color.RGBA{255, 0, 0, 255} f.Samples = 1000 p.Add(f) p.Add(plotter.NewGrid()) err := p.Save(20*vg.Centimeter, -1, "testdata/poly-plot.png") if err != nil { log.Fatal(err) } } } func genXY(n int, f func(x float64, ps []float64) float64, ps []float64, xmin, xmax float64) ([]float64, []float64) { xdata := make([]float64, n) ydata := make([]float64, n) rnd := rand.New(rand.NewSource(1234)) xstep := (xmax - xmin) / float64(n) p := make([]float64, len(ps)) for i := 0; i < n; i++ { x := xmin + xstep*float64(i) for j := range p { v := rnd.NormFloat64() p[j] = ps[j] + v*0.2 } xdata[i] = x ydata[i] = f(x, p) } return xdata, ydata }
Output:
Example (Powerlaw) ¶
package main import ( "bufio" "image/color" "log" "math" "os" "strconv" "strings" "go-hep.org/x/hep/fit" "go-hep.org/x/hep/hbook" "go-hep.org/x/hep/hplot" "gonum.org/v1/gonum/floats" "gonum.org/v1/gonum/optimize" "gonum.org/v1/plot/plotter" "gonum.org/v1/plot/vg" ) func main() { var ( amp = 11.021171432949746 index = -2.027389113217428 want = []float64{amp, index} ) xdata, ydata, yerrs, err := readXYerr("testdata/powerlaw-data.txt") if err != nil { log.Fatal(err) } plaw := func(x, amp, index float64) float64 { return amp * math.Pow(x, index) } res, err := fit.Curve1D( fit.Func1D{ F: func(x float64, ps []float64) float64 { return plaw(x, ps[0], ps[1]) }, X: xdata, Y: ydata, Err: yerrs, Ps: []float64{1, 1}, }, nil, &optimize.NelderMead{}, ) if err != nil { log.Fatal(err) } if err := res.Status.Err(); err != nil { log.Fatal(err) } if got := res.X; !floats.EqualApprox(got, want, 1e-3) { log.Fatalf("got= %v\nwant=%v\n", got, want) } { p := hplot.New() p.X.Label.Text = "f(x) = a * x^b" p.Y.Label.Text = "y-data" p.X.Min = 0 p.X.Max = 10 p.Y.Min = 0 p.Y.Max = 10 pts := make([]hbook.Point2D, len(xdata)) for i := range pts { pts[i].X = xdata[i] pts[i].Y = ydata[i] pts[i].ErrY.Min = 0.5 * yerrs[i] pts[i].ErrY.Max = 0.5 * yerrs[i] } s := hplot.NewS2D(hbook.NewS2D(pts...), hplot.WithYErrBars) s.Color = color.RGBA{0, 0, 255, 255} p.Add(s) f := plotter.NewFunction(func(x float64) float64 { return plaw(x, res.X[0], res.X[1]) }) f.Color = color.RGBA{255, 0, 0, 255} f.Samples = 1000 p.Add(f) p.Add(plotter.NewGrid()) err := p.Save(20*vg.Centimeter, -1, "testdata/powerlaw-plot.png") if err != nil { log.Fatal(err) } } } func readXYerr(fname string) (xs, ys, yerrs []float64, err error) { f, err := os.Open(fname) if err != nil { return xs, ys, yerrs, err } defer f.Close() scan := bufio.NewScanner(f) for scan.Scan() { line := scan.Text() toks := strings.Split(line, " ") x, err := strconv.ParseFloat(toks[0], 64) if err != nil { return xs, ys, yerrs, err } xs = append(xs, x) y, err := strconv.ParseFloat(toks[1], 64) if err != nil { return xs, ys, yerrs, err } ys = append(ys, y) yerr, err := strconv.ParseFloat(toks[2], 64) if err != nil { return xs, ys, yerrs, err } yerrs = append(yerrs, yerr) } return }
Output:
func H1D ¶
func H1D(h *hbook.H1D, f Func1D, settings *optimize.Settings, m optimize.Method) (*optimize.Result, error)
H1D returns the fit of histogram h with function f and optimization method m.
Only bins with at least an entry are considered for the fit. In case settings is nil, the optimize.DefaultSettingsLocal is used. In case m is nil, the same default optimization method than for Curve1D is used.
Example (Gaussian) ¶
package main import ( "image/color" "log" "math" "go-hep.org/x/hep/fit" "go-hep.org/x/hep/hbook" "go-hep.org/x/hep/hplot" "golang.org/x/exp/rand" "gonum.org/v1/gonum/floats" "gonum.org/v1/gonum/optimize" "gonum.org/v1/gonum/stat/distuv" "gonum.org/v1/plot/plotter" "gonum.org/v1/plot/vg" ) func main() { var ( mean = 2.0 sigma = 4.0 // Values from gonum/optimize: want = []float64{447.0483517081991, 2.02127773281178, 3.9965893891862687} // Values from ROOT: // want = []float64{4.53720e+02, 1.93218e+00, 3.93188e+00} ) const npoints = 10000 // Create a normal distribution. dist := distuv.Normal{ Mu: mean, Sigma: sigma, Src: rand.New(rand.NewSource(0)), } // Draw some random values from the standard // normal distribution. hist := hbook.NewH1D(100, -20, +25) for i := 0; i < npoints; i++ { v := dist.Rand() hist.Fill(v, 1) } gauss := func(x, cst, mu, sigma float64) float64 { v := (x - mu) / sigma return cst * math.Exp(-0.5*v*v) } res, err := fit.H1D( hist, fit.Func1D{ F: func(x float64, ps []float64) float64 { return gauss(x, ps[0], ps[1], ps[2]) }, N: len(want), }, nil, &optimize.NelderMead{}, ) if err != nil { log.Fatal(err) } if err := res.Status.Err(); err != nil { log.Fatal(err) } if got := res.X; !floats.EqualApprox(got, want, 1e-3) { log.Fatalf("got= %v\nwant=%v\n", got, want) } { p := hplot.New() p.X.Label.Text = "f(x) = cst * exp(-0.5 * ((x-mu)/sigma)^2)" p.Y.Label.Text = "y-data" p.Y.Min = 0 h := hplot.NewH1D(hist) h.Color = color.RGBA{0, 0, 255, 255} p.Add(h) f := plotter.NewFunction(func(x float64) float64 { return gauss(x, res.X[0], res.X[1], res.X[2]) }) f.Color = color.RGBA{255, 0, 0, 255} f.Samples = 1000 p.Add(f) p.Add(plotter.NewGrid()) err := p.Save(20*vg.Centimeter, -1, "testdata/h1d-gauss-plot.png") if err != nil { log.Fatal(err) } } }
Output:
Types ¶
type Func1D ¶
type Func1D struct { // F is the function to minimize. // ps is the slice of parameters to optimize during the fit. F func(x float64, ps []float64) float64 // N is the number of parameters to optimize during the fit. // If N is 0, Ps must not be nil. N int // Ps is the initial values for the parameters. // If Ps is nil, the set of initial parameters values is a slice of // length N filled with zeros. Ps []float64 X []float64 Y []float64 Err []float64 // contains filtered or unexported fields }
Func1D describes a 1D function to fit some data.
Click to show internal directories.
Click to hide internal directories.