fit

package
Version: v0.32.1 Latest Latest
Warning

This package is not in the latest version of its module.

 Go to latest Published: Sep 16, 2022 License: BSD-3-Clause Imports: 4 Imported by: 2

README

fit

GoDoc

fit is a WIP package to provide easy fitting models and curve fitting functions.

H1D

Fit a gaussian

h1d-gaussian-example

func ExampleH1D_gaussian() {
	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)
		}
	}
}

Curve1D

Fit a gaussian

func1d-gaussian-example

func ExampleCurve1D_gaussian() {
	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)
		}
	}
}
Fit a powerlaw (with Y-errors)

func1d-powerlaw-example

func ExampleCurve1D_powerlaw() {
	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(true))
		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)
		}
	}
}
Fit an exponential

func1d-exp-example

func ExampleCurve1D_exponential() {
	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)
		}
	}
}
Fit a polynomial

func1d-poly-example

func ExampleCurve1D_poly() {
	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)
		}
	}
}

Fitting with more than one independent variable (x has more than one dimension)

Fit a flat plane

2d-example

func ExampleCurveND_plane() {
	var (
		m1    = 0.3
		m2    = 0.1
		c     = 0.2
		ps    = []float64{m1, m2, c}
		n0    = 10
		n1    = 10
		x0min = -1.
		x0max = 1.
		x1min = -1.
		x1max = 1.
	)

	plane := func(x, ps []float64) float64 {
		return ps[0]*x[0] + ps[1]*x[1] + ps[2]
	}

	xData, yData := genData2D(n0, n1, plane, ps, x0min, x0max, x1min, x1max)

	res, err := fit.CurveND(
		fit.FuncND{
			F: func(x []float64, ps []float64) float64 {
				return plane(x, ps)
			},
			X: xData,
			Y: yData,
			N: 3,
		},
		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{m1, m2, c}; !floats.EqualApprox(got, want, 0.1) {
		log.Fatalf("got= %v\nwant=%v\n", got, want)
	}

	{
		// slicing for a particular x0 value to plot y as a function of x1,
		// to visualise how well the fit is working for a given x0.
		x0Selection := 8
		if 0 > x0Selection || x0Selection > n0 {
			log.Fatalf("x0 slice, %d, is not in valid range [0 - %d]", x0Selection, n0)
		}
		x0SlicePos := x0min + ((x0max-x0min)/float64(n0))*float64(x0Selection)

		var x1Slice []float64
		var ySlice []float64

		for i := range xData {
			if xData[i][0] == x0SlicePos {
				x1Slice = append(x1Slice, xData[i][1])
				ySlice = append(ySlice, yData[i])
			}
		}

		p := hplot.New()
		p.Title.Text = fmt.Sprintf("Slice of plane at x0 = %.2f", x0SlicePos)
		p.X.Label.Text = "x1"
		p.Y.Label.Text = "y"
		p.Y.Min = x1min
		p.Y.Max = x1max
		p.X.Min = x0min
		p.X.Max = x0max

		s := hplot.NewS2D(hplot.ZipXY(x1Slice, ySlice))
		s.Color = color.RGBA{B: 255, A: 255}
		p.Add(s)

		shiftLine := func(x, m, c, mxOtherAxis float64) float64 {
			return m*x + c + mxOtherAxis
		}

		f := plotter.NewFunction(func(x float64) float64 {
			return shiftLine(x, res.X[1], res.X[2], res.X[0]*x0SlicePos)
		})
		f.Color = color.RGBA{R: 255, A: 255}
		f.Samples = 1000
		p.Add(f)

		p.Add(plotter.NewGrid())
		err := p.Save(20*vg.Centimeter, -1, "testdata/2d-plane-plot.png")
		if err != nil {
			log.Fatal(err)
		}
	}
}

Documentation

Overview

Package fit provides functions to fit data.

Index

Examples

Constants

This section is empty.

Variables

This section is empty.

Functions

func Curve1D 

func Curve1D(f Func1D, settings *optimize.Settings, m optimize.Method) (*optimize.Result, error)

Curve1D returns the result of a non-linear least squares to fit a function f to the underlying data with method m.

Example (Exponential) 
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)
	}
}

Output:
Example (Gaussian) 
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)
	}
}

Output:
Example (Poly) 
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)
	}
}

Output:
Example (Powerlaw) 
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(true))
	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)
	}
}

Output:

func CurveND added in v0.27.0 

func CurveND(f FuncND, settings *optimize.Settings, m optimize.Method) (*optimize.Result, error)

CurveND returns the result of a non-linear least squares to fit a function f to the underlying data with method m, where there is more than one independent variable.

Example (Plane) 
var (
	m1    = 0.3
	m2    = 0.1
	c     = 0.2
	ps    = []float64{m1, m2, c}
	n0    = 10
	n1    = 10
	x0min = -1.
	x0max = 1.
	x1min = -1.
	x1max = 1.
)

plane := func(x, ps []float64) float64 {
	return ps[0]*x[0] + ps[1]*x[1] + ps[2]
}

xData, yData := genData2D(n0, n1, plane, ps, x0min, x0max, x1min, x1max)

res, err := fit.CurveND(
	fit.FuncND{
		F: func(x []float64, ps []float64) float64 {
			return plane(x, ps)
		},
		X: xData,
		Y: yData,
		N: 3,
	},
	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{m1, m2, c}; !floats.EqualApprox(got, want, 0.1) {
	log.Fatalf("got= %v\nwant=%v\n", got, want)
}

{
	// slicing for a particular x0 value to plot y as a function of x1,
	// to visualise how well the fit is working for a given x0.
	x0Selection := 8
	if 0 > x0Selection || x0Selection > n0 {
		log.Fatalf("x0 slice, %d, is not in valid range [0 - %d]", x0Selection, n0)
	}
	x0SlicePos := x0min + ((x0max-x0min)/float64(n0))*float64(x0Selection)

	var x1Slice []float64
	var ySlice []float64

	for i := range xData {
		if xData[i][0] == x0SlicePos {
			x1Slice = append(x1Slice, xData[i][1])
			ySlice = append(ySlice, yData[i])
		}
	}

	p := hplot.New()
	p.Title.Text = fmt.Sprintf("Slice of plane at x0 = %.2f", x0SlicePos)
	p.X.Label.Text = "x1"
	p.Y.Label.Text = "y"
	p.Y.Min = x1min
	p.Y.Max = x1max
	p.X.Min = x0min
	p.X.Max = x0max

	s := hplot.NewS2D(hplot.ZipXY(x1Slice, ySlice))
	s.Color = color.RGBA{B: 255, A: 255}
	p.Add(s)

	shiftLine := func(x, m, c, mxOtherAxis float64) float64 {
		return m*x + c + mxOtherAxis
	}

	f := plotter.NewFunction(func(x float64) float64 {
		return shiftLine(x, res.X[1], res.X[2], res.X[0]*x0SlicePos)
	})
	f.Color = color.RGBA{R: 255, A: 255}
	f.Samples = 1000
	p.Add(f)

	p.Add(plotter.NewGrid())
	err := p.Save(20*vg.Centimeter, -1, "testdata/2d-plane-plot.png")
	if err != nil {
		log.Fatal(err)
	}
}

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.

type FuncND added in v0.27.0 

type FuncND struct {
	// F is the function to minimize.
	// ps is the slice of parameters to optimize during the fit.
	// x is the slice of independent variables.
	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 is the multidimensional slice of the independent variables,
	// it must be structured so that the X[i] is a list of values for the
	// independent variables that corresponds to a single Y value.
	// In other words, the sequence of rows must correspond to the sequence
	// of independent variable values.
	X   [][]float64
	Y   []float64
	Err []float64
	// contains filtered or unexported fields
}

FuncND describes a multivariate function F(x0, x1... xn; p0, p1... pn) for which the parameters ps can be found with a fit.

Source Files

View all

Jump to

Keyboard shortcuts

? : This menu
/ : Search site
f or F : Jump to
y or Y : Canonical URL