# fit

package
v0.34.1 Latest Latest

Go to latest
Published: Oct 12, 2023 License: BSD-3-Clause

### fit

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

#### H1D

##### Fit a gaussian

``````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),
},
)
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}

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

err := p.Save(20*vg.Centimeter, -1, "testdata/h1d-gauss-plot.png")
if err != nil {
log.Fatal(err)
}
}
}
``````

#### Curve1D

##### Fit a gaussian

``````func ExampleCurve1D_gaussian() {
var (
cst   = 3.0
mean  = 30.0
sigma = 20.0
want  = []float64{cst, mean, sigma}
)

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},
},
)
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}

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

err := p.Save(20*vg.Centimeter, -1, "testdata/gauss-plot.png")
if err != nil {
log.Fatal(err)
}
}
}
``````
##### Fit a powerlaw (with Y-errors)

``````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},
},
)
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}

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

err := p.Save(20*vg.Centimeter, -1, "testdata/powerlaw-plot.png")
if err != nil {
log.Fatal(err)
}
}
}
``````
##### Fit an exponential

``````func ExampleCurve1D_exponential() {
const (
a   = 0.3
b   = 0.1
ndf = 2.0
)

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,
},
)
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}

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

err := p.Save(20*vg.Centimeter, -1, "testdata/exp-plot.png")
if err != nil {
log.Fatal(err)
}
}
}
``````
##### Fit a polynomial

``````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},
},
)
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}

f := plotter.NewFunction(func(x float64) float64 {
return poly(x, res.X)
})
f.Color = color.RGBA{255, 0, 0, 255}
f.Samples = 1000

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

``````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,
},
)
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}

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

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.

### 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
)

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,
},
)
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}

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

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}
)

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},
},
)
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}

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

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},
},
)
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}

f := plotter.NewFunction(func(x float64) float64 {
return poly(x, res.X)
})
f.Color = color.RGBA{255, 0, 0, 255}
f.Samples = 1000

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},
},
)
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}

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

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,
},
)
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}

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

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),
},
)
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}

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

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.