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
|
// Test GLM log-likelihood and score functions using numeric derivatives. The
// tests confirm that the analytic score function agrees with the numeric
// derivative of the log-likelihood function.
package glm
import (
"fmt"
"testing"
"github.com/kshedden/statmodel/statmodel"
"gonum.org/v1/gonum/diff/fd"
"gonum.org/v1/gonum/floats"
)
// A test problem
type difftestprob struct {
title string
family *Family
data statmodel.Dataset
xnames []string
weight bool
offset bool
params [][]float64
scale float64
l2wgt map[string]float64
}
var diffTests []difftestprob = []difftestprob{
{
title: "Gaussian 1",
family: NewFamily(GaussianFamily),
data: data1(),
xnames: []string{"x1", "x2"},
weight: false,
scale: 2,
params: [][]float64{{1, 0}, {0, 1}, {1, 1}, {-1, 1}},
},
{
title: "Gaussian 2",
family: NewFamily(GaussianFamily),
data: data1(),
xnames: []string{"x1", "x2"},
weight: true,
scale: 2,
params: [][]float64{{1, 0}, {0, 1}, {1, 1}, {-1, 1}},
},
{
title: "Poisson 1",
family: NewFamily(PoissonFamily),
data: data1(),
xnames: []string{"x1", "x2"},
weight: false,
scale: 1,
params: [][]float64{{1, 0}, {0, 1}, {1, 1}, {-1, 1}},
},
{
title: "Poisson 2",
family: NewFamily(PoissonFamily),
data: data1(),
xnames: []string{"x1", "x2"},
weight: true,
scale: 1,
params: [][]float64{{1, 0}, {0, 1}, {1, 1}, {-1, 1}},
},
{
title: "Binomial 1",
family: NewFamily(BinomialFamily),
data: data2(),
xnames: []string{"x1", "x2", "x3"},
weight: true,
params: [][]float64{{1, 0, 0}, {0, 1, 0}, {1, 1, 1}, {-1, 0, 1}},
scale: 1,
},
{
title: "Gamma 1",
family: NewFamily(GammaFamily),
data: data4(),
xnames: []string{"x1", "x2", "x3"},
weight: true,
params: [][]float64{{1, 0, 0}, {1, 1, 1}, {1, 0, -0.1}},
scale: 2,
},
{
title: "Inverse Gaussian 1",
family: NewFamily(InvGaussianFamily),
data: data4(),
xnames: []string{"x1", "x2", "x3"},
params: [][]float64{{1, 0, 0}, {1, 1, 1}, {1, 0, -0.1}},
weight: true,
scale: 0.5,
},
{
title: "Tweedie 1",
family: NewTweedieFamily(1.5, NewLink(LogLink)),
data: data1(),
xnames: []string{"x1", "x2"},
params: [][]float64{{1, 0}, {0, 1}, {1, 1}, {-1, 1}},
weight: false,
scale: 1.2,
},
}
func TestGrad(t *testing.T) {
for _, dt := range diffTests {
config := DefaultConfig()
if dt.weight {
config.WeightVar = "w"
}
glm, err := NewGLM(dt.data, "y", dt.xnames, config)
if err != nil {
panic(err)
}
p := len(dt.params[0])
ngrad := make([]float64, p)
score := make([]float64, p)
loglike := func(x []float64) float64 {
return glm.LogLike(&GLMParams{x, dt.scale}, true)
}
for _, params := range dt.params {
fd.Gradient(ngrad, loglike, params, nil)
glm.Score(&GLMParams{params, dt.scale}, score)
if !floats.EqualApprox(score, ngrad, 1e-5) {
fmt.Printf("%s\n", dt.title)
fmt.Printf("Numerical: %v\n", ngrad)
fmt.Printf("Analytical: %v\n", score)
t.Fail()
}
}
}
}
|
