/
fit.go
169 lines (136 loc) · 3.9 KB
/
fit.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
// Copyright ©2017 The go-hep Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
// Package fit provides functions to fit data.
package fit // import "go-hep.org/x/hep/fit"
import (
"gonum.org/v1/gonum/diff/fd"
"gonum.org/v1/gonum/mat"
)
//go:generate go get github.com/campoy/embedmd
//go:generate embedmd -w README.md
// Func1D describes a 1D function to fit some data.
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
sig2 []float64 // inverse of squares of measurement errors along Y.
fct func(ps []float64) float64 // cost function (objective function)
grad func(grad, ps []float64)
hess func(hess *mat.SymDense, x []float64)
}
func (f *Func1D) init() {
f.sig2 = make([]float64, len(f.Y))
switch {
default:
for i := range f.Y {
f.sig2[i] = 1
}
case f.Err != nil:
for i, v := range f.Err {
f.sig2[i] = 1 / (v * v)
}
}
if f.Ps == nil {
f.Ps = make([]float64, f.N)
}
if len(f.Ps) == 0 {
panic("fit: invalid number of initial parameters")
}
if len(f.X) != len(f.Y) {
panic("fit: mismatch length")
}
if len(f.sig2) != len(f.Y) {
panic("fit: mismatch length")
}
f.fct = func(ps []float64) float64 {
var chi2 float64
for i := range f.X {
res := f.F(f.X[i], ps) - f.Y[i]
chi2 += res * res * f.sig2[i]
}
return 0.5 * chi2
}
f.grad = func(grad, ps []float64) {
fd.Gradient(grad, f.fct, ps, nil)
}
f.hess = func(hess *mat.SymDense, x []float64) {
fd.Hessian(hess, f.fct, x, nil)
}
}
// FuncND describes a multivariate function F(x0, x1... xn; p0, p1... pn)
// for which the parameters ps can be found with a fit.
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
sig2 []float64 // inverse of squares of measurement errors along Y.
fct func(ps []float64) float64 // cost function (objective function)
grad func(grad, ps []float64)
hess func(hess *mat.SymDense, x []float64) // hessian matrix
}
func (f *FuncND) init() {
f.sig2 = make([]float64, len(f.Y))
switch {
default:
for i := range f.Y {
f.sig2[i] = 1
}
case f.Err != nil:
for i, v := range f.Err {
f.sig2[i] = 1 / (v * v)
}
}
if f.Ps == nil {
f.Ps = make([]float64, f.N)
}
if len(f.Ps) == 0 {
panic("fit: invalid number of initial parameters")
}
if len(f.X) != len(f.Y) {
panic("fit: mismatch length")
}
if len(f.sig2) != len(f.Y) {
panic("fit: mismatch length")
}
f.fct = func(ps []float64) float64 {
var chi2 float64
for i := range f.X {
res := f.F(f.X[i], ps) - f.Y[i]
chi2 += res * res * f.sig2[i]
}
return 0.5 * chi2
}
f.grad = func(grad []float64, ps []float64) {
fd.Gradient(grad, f.fct, ps, nil)
}
f.hess = func(hess *mat.SymDense, x []float64) {
fd.Hessian(hess, f.fct, x, nil)
}
}