go-hep.org/x/hep@v0.38.1/fit/hist_example_test.go (about)

     1  // Copyright ©2017 The go-hep Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package fit_test
     6  
     7  import (
     8  	"image/color"
     9  	"log"
    10  	"math"
    11  	"math/rand/v2"
    12  
    13  	"go-hep.org/x/hep/fit"
    14  	"go-hep.org/x/hep/hbook"
    15  	"go-hep.org/x/hep/hplot"
    16  	"gonum.org/v1/gonum/floats"
    17  	"gonum.org/v1/gonum/optimize"
    18  	"gonum.org/v1/gonum/stat/distuv"
    19  	"gonum.org/v1/plot/plotter"
    20  	"gonum.org/v1/plot/vg"
    21  )
    22  
    23  func ExampleH1D_gaussian() {
    24  	var (
    25  		mean  = 2.0
    26  		sigma = 4.0
    27  		// Values from gonum/optimize:
    28  		want = []float64{450.56454241860934, 2.0146898541006277, 3.9613086671267466}
    29  		// Values from ROOT:
    30  		// want  = []float64{4.53720e+02, 1.93218e+00, 3.93188e+00}
    31  	)
    32  
    33  	const npoints = 10000
    34  
    35  	// Create a normal distribution.
    36  	dist := distuv.Normal{
    37  		Mu:    mean,
    38  		Sigma: sigma,
    39  		Src:   rand.New(rand.NewPCG(0, 0)),
    40  	}
    41  
    42  	// Draw some random values from the standard
    43  	// normal distribution.
    44  	hist := hbook.NewH1D(100, -20, +25)
    45  	for range npoints {
    46  		v := dist.Rand()
    47  		hist.Fill(v, 1)
    48  	}
    49  
    50  	gauss := func(x, cst, mu, sigma float64) float64 {
    51  		v := (x - mu) / sigma
    52  		return cst * math.Exp(-0.5*v*v)
    53  	}
    54  
    55  	res, err := fit.H1D(
    56  		hist,
    57  		fit.Func1D{
    58  			F: func(x float64, ps []float64) float64 {
    59  				return gauss(x, ps[0], ps[1], ps[2])
    60  			},
    61  			N: len(want),
    62  		},
    63  		nil, &optimize.NelderMead{},
    64  	)
    65  	if err != nil {
    66  		log.Fatal(err)
    67  	}
    68  
    69  	if err := res.Status.Err(); err != nil {
    70  		log.Fatal(err)
    71  	}
    72  	if got := res.X; !floats.EqualApprox(got, want, 1e-3) {
    73  		log.Fatalf("got= %v\nwant=%v\n", got, want)
    74  	}
    75  
    76  	{
    77  		p := hplot.New()
    78  		p.X.Label.Text = "f(x) = cst * exp(-0.5 * ((x-mu)/sigma)^2)"
    79  		p.Y.Label.Text = "y-data"
    80  		p.Y.Min = 0
    81  
    82  		h := hplot.NewH1D(hist)
    83  		h.Color = color.RGBA{0, 0, 255, 255}
    84  		p.Add(h)
    85  
    86  		f := plotter.NewFunction(func(x float64) float64 {
    87  			return gauss(x, res.X[0], res.X[1], res.X[2])
    88  		})
    89  		f.Color = color.RGBA{255, 0, 0, 255}
    90  		f.Samples = 1000
    91  		p.Add(f)
    92  
    93  		p.Add(plotter.NewGrid())
    94  
    95  		err := p.Save(20*vg.Centimeter, -1, "testdata/h1d-gauss-plot.png")
    96  		if err != nil {
    97  			log.Fatal(err)
    98  		}
    99  	}
   100  }