1// Copyright ©2014 The Gonum 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 5package stat 6 7import ( 8 "fmt" 9 "math" 10 "reflect" 11 "strconv" 12 "testing" 13 14 "golang.org/x/exp/rand" 15 16 "gonum.org/v1/gonum/floats" 17 "gonum.org/v1/gonum/floats/scalar" 18) 19 20func ExampleCircularMean() { 21 x := []float64{0, 0.25 * math.Pi, 0.75 * math.Pi} 22 weights := []float64{1, 2, 2.5} 23 cmean := CircularMean(x, weights) 24 25 fmt.Printf("The circular mean is %.5f.\n", cmean) 26 // Output: 27 // The circular mean is 1.37037. 28} 29 30func TestCircularMean(t *testing.T) { 31 for i, test := range []struct { 32 x []float64 33 wts []float64 34 ans float64 35 }{ 36 // Values compared against scipy. 37 { 38 x: []float64{0, 2 * math.Pi}, 39 ans: 0, 40 }, 41 { 42 x: []float64{0, 0.5 * math.Pi}, 43 ans: 0.78539816339744, 44 }, 45 { 46 x: []float64{-1.5 * math.Pi, 0.5 * math.Pi, 2.5 * math.Pi}, 47 wts: []float64{1, 2, 3}, 48 ans: 0.5 * math.Pi, 49 }, 50 { 51 x: []float64{0, 0.5 * math.Pi}, 52 wts: []float64{1, 2}, 53 ans: 1.10714871779409, 54 }, 55 } { 56 c := CircularMean(test.x, test.wts) 57 if math.Abs(c-test.ans) > 1e-14 { 58 t.Errorf("Circular mean mismatch case %d: Expected %v, Found %v", i, test.ans, c) 59 } 60 } 61 if !panics(func() { CircularMean(make([]float64, 3), make([]float64, 2)) }) { 62 t.Errorf("CircularMean did not panic with x, wts length mismatch") 63 } 64} 65 66func ExampleCorrelation() { 67 x := []float64{8, -3, 7, 8, -4} 68 y := []float64{10, 5, 6, 3, -1} 69 w := []float64{2, 1.5, 3, 3, 2} 70 71 fmt.Println("Correlation computes the degree to which two datasets move together") 72 fmt.Println("about their mean. For example, x and y above move similarly.") 73 74 c := Correlation(x, y, w) 75 fmt.Printf("Correlation is %.5f\n", c) 76 77 // Output: 78 // Correlation computes the degree to which two datasets move together 79 // about their mean. For example, x and y above move similarly. 80 // Correlation is 0.59915 81} 82 83func TestCorrelation(t *testing.T) { 84 for i, test := range []struct { 85 x []float64 86 y []float64 87 w []float64 88 ans float64 89 }{ 90 { 91 x: []float64{8, -3, 7, 8, -4}, 92 y: []float64{8, -3, 7, 8, -4}, 93 w: nil, 94 ans: 1, 95 }, 96 { 97 x: []float64{8, -3, 7, 8, -4}, 98 y: []float64{8, -3, 7, 8, -4}, 99 w: []float64{1, 1, 1, 1, 1}, 100 ans: 1, 101 }, 102 { 103 x: []float64{8, -3, 7, 8, -4}, 104 y: []float64{8, -3, 7, 8, -4}, 105 w: []float64{1, 6, 7, 0.8, 2.1}, 106 ans: 1, 107 }, 108 { 109 x: []float64{8, -3, 7, 8, -4}, 110 y: []float64{10, 15, 4, 5, -1}, 111 w: nil, 112 ans: 0.0093334660769059, 113 }, 114 { 115 x: []float64{8, -3, 7, 8, -4}, 116 y: []float64{10, 15, 4, 5, -1}, 117 w: nil, 118 ans: 0.0093334660769059, 119 }, 120 { 121 x: []float64{8, -3, 7, 8, -4}, 122 y: []float64{10, 15, 4, 5, -1}, 123 w: []float64{1, 3, 1, 2, 2}, 124 ans: -0.13966633352689, 125 }, 126 } { 127 c := Correlation(test.x, test.y, test.w) 128 if math.Abs(test.ans-c) > 1e-14 { 129 t.Errorf("Correlation mismatch case %d. Expected %v, Found %v", i, test.ans, c) 130 } 131 } 132 if !panics(func() { Correlation(make([]float64, 2), make([]float64, 3), make([]float64, 3)) }) { 133 t.Errorf("Correlation did not panic with length mismatch") 134 } 135 if !panics(func() { Correlation(make([]float64, 2), make([]float64, 3), nil) }) { 136 t.Errorf("Correlation did not panic with length mismatch") 137 } 138 if !panics(func() { Correlation(make([]float64, 3), make([]float64, 3), make([]float64, 2)) }) { 139 t.Errorf("Correlation did not panic with weights length mismatch") 140 } 141} 142 143func ExampleKendall() { 144 x := []float64{8, -3, 7, 8, -4} 145 y := []float64{10, 5, 6, 3, -1} 146 w := []float64{2, 1.5, 3, 3, 2} 147 148 fmt.Println("Kendall correlation computes the number of ordered pairs") 149 fmt.Println("between two datasets.") 150 151 c := Kendall(x, y, w) 152 fmt.Printf("Kendall correlation is %.5f\n", c) 153 154 // Output: 155 // Kendall correlation computes the number of ordered pairs 156 // between two datasets. 157 // Kendall correlation is 0.25000 158} 159 160func TestKendall(t *testing.T) { 161 for i, test := range []struct { 162 x []float64 163 y []float64 164 weights []float64 165 ans float64 166 }{ 167 { 168 x: []float64{0, 1, 2, 3}, 169 y: []float64{0, 1, 2, 3}, 170 weights: nil, 171 ans: 1, 172 }, 173 { 174 x: []float64{0, 1}, 175 y: []float64{1, 0}, 176 weights: nil, 177 ans: -1, 178 }, 179 { 180 x: []float64{8, -3, 7, 8, -4}, 181 y: []float64{10, 15, 4, 5, -1}, 182 weights: nil, 183 ans: 0.2, 184 }, 185 { 186 x: []float64{8, -3, 7, 8, -4}, 187 y: []float64{10, 5, 6, 3, -1}, 188 weights: nil, 189 ans: 0.4, 190 }, 191 { 192 x: []float64{1, 2, 3, 4, 5}, 193 y: []float64{2, 3, 4, 5, 6}, 194 weights: []float64{1, 1, 1, 1, 1}, 195 ans: 1, 196 }, 197 { 198 x: []float64{1, 2, 3, 2, 1}, 199 y: []float64{2, 3, 2, 1, 0}, 200 weights: []float64{1, 1, 0, 0, 0}, 201 ans: 1, 202 }, 203 } { 204 c := Kendall(test.x, test.y, test.weights) 205 if math.Abs(test.ans-c) > 1e-14 { 206 t.Errorf("Correlation mismatch case %d. Expected %v, Found %v", i, test.ans, c) 207 } 208 } 209 if !panics(func() { Kendall(make([]float64, 2), make([]float64, 3), make([]float64, 3)) }) { 210 t.Errorf("Kendall did not panic with length mismatch") 211 } 212 if !panics(func() { Kendall(make([]float64, 2), make([]float64, 3), nil) }) { 213 t.Errorf("Kendall did not panic with length mismatch") 214 } 215 if !panics(func() { Kendall(make([]float64, 3), make([]float64, 3), make([]float64, 2)) }) { 216 t.Errorf("Kendall did not panic with weights length mismatch") 217 } 218} 219 220func ExampleCovariance() { 221 fmt.Println("Covariance computes the degree to which datasets move together") 222 fmt.Println("about their mean.") 223 x := []float64{8, -3, 7, 8, -4} 224 y := []float64{10, 2, 2, 4, 1} 225 cov := Covariance(x, y, nil) 226 fmt.Printf("Cov = %.4f\n", cov) 227 fmt.Println("If datasets move perfectly together, the variance equals the covariance") 228 y2 := []float64{12, 1, 11, 12, 0} 229 cov2 := Covariance(x, y2, nil) 230 varX := Variance(x, nil) 231 fmt.Printf("Cov2 is %.4f, VarX is %.4f", cov2, varX) 232 // Output: 233 // Covariance computes the degree to which datasets move together 234 // about their mean. 235 // Cov = 13.8000 236 // If datasets move perfectly together, the variance equals the covariance 237 // Cov2 is 37.7000, VarX is 37.7000 238} 239 240func TestCovariance(t *testing.T) { 241 for i, test := range []struct { 242 p []float64 243 q []float64 244 weights []float64 245 ans float64 246 }{ 247 { 248 p: []float64{0.75, 0.1, 0.05}, 249 q: []float64{0.5, 0.25, 0.25}, 250 ans: 0.05625, 251 }, 252 { 253 p: []float64{1, 2, 3}, 254 q: []float64{2, 4, 6}, 255 ans: 2, 256 }, 257 { 258 p: []float64{1, 2, 3}, 259 q: []float64{1, 4, 9}, 260 ans: 4, 261 }, 262 { 263 p: []float64{1, 2, 3}, 264 q: []float64{1, 4, 9}, 265 weights: []float64{1, 1.5, 1}, 266 ans: 3.2, 267 }, 268 { 269 p: []float64{1, 4, 9}, 270 q: []float64{1, 4, 9}, 271 weights: []float64{1, 1.5, 1}, 272 ans: 13.142857142857146, 273 }, 274 } { 275 c := Covariance(test.p, test.q, test.weights) 276 if math.Abs(c-test.ans) > 1e-14 { 277 t.Errorf("Covariance mismatch case %d: Expected %v, Found %v", i, test.ans, c) 278 } 279 } 280 281 // test the panic states 282 if !panics(func() { Covariance(make([]float64, 2), make([]float64, 3), nil) }) { 283 t.Errorf("Covariance did not panic with x, y length mismatch") 284 } 285 if !panics(func() { Covariance(make([]float64, 3), make([]float64, 3), make([]float64, 2)) }) { 286 t.Errorf("Covariance did not panic with x, weights length mismatch") 287 } 288 289} 290 291func TestCrossEntropy(t *testing.T) { 292 for i, test := range []struct { 293 p []float64 294 q []float64 295 ans float64 296 }{ 297 { 298 p: []float64{0.75, 0.1, 0.05}, 299 q: []float64{0.5, 0.25, 0.25}, 300 ans: 0.7278045395879426, 301 }, 302 { 303 p: []float64{0.75, 0.1, 0.05, 0, 0, 0}, 304 q: []float64{0.5, 0.25, 0.25, 0, 0, 0}, 305 ans: 0.7278045395879426, 306 }, 307 { 308 p: []float64{0.75, 0.1, 0.05, 0, 0, 0.1}, 309 q: []float64{0.5, 0.25, 0.25, 0, 0, 0}, 310 ans: math.Inf(1), 311 }, 312 { 313 p: nil, 314 q: nil, 315 ans: 0, 316 }, 317 } { 318 c := CrossEntropy(test.p, test.q) 319 if math.Abs(c-test.ans) > 1e-14 { 320 t.Errorf("Cross entropy mismatch case %d: Expected %v, Found %v", i, test.ans, c) 321 } 322 } 323 if !panics(func() { CrossEntropy(make([]float64, 3), make([]float64, 2)) }) { 324 t.Errorf("CrossEntropy did not panic with p, q length mismatch") 325 } 326} 327 328func ExampleEntropy() { 329 330 p := []float64{0.05, 0.1, 0.9, 0.05} 331 entP := Entropy(p) 332 333 q := []float64{0.2, 0.4, 0.25, 0.15} 334 entQ := Entropy(q) 335 336 r := []float64{0.2, 0, 0, 0.5, 0, 0.2, 0.1, 0, 0, 0} 337 entR := Entropy(r) 338 339 s := []float64{0, 0, 1, 0} 340 entS := Entropy(s) 341 342 fmt.Println("Entropy is a measure of the amount of uncertainty in a distribution") 343 fmt.Printf("The second bin of p is very likely to occur. It's entropy is %.4f\n", entP) 344 fmt.Printf("The distribution of q is more spread out. It's entropy is %.4f\n", entQ) 345 fmt.Println("Adding buckets with zero probability does not change the entropy.") 346 fmt.Printf("The entropy of r is: %.4f\n", entR) 347 fmt.Printf("A distribution with no uncertainty has entropy %.4f\n", entS) 348 // Output: 349 // Entropy is a measure of the amount of uncertainty in a distribution 350 // The second bin of p is very likely to occur. It's entropy is 0.6247 351 // The distribution of q is more spread out. It's entropy is 1.3195 352 // Adding buckets with zero probability does not change the entropy. 353 // The entropy of r is: 1.2206 354 // A distribution with no uncertainty has entropy 0.0000 355} 356 357func ExampleExKurtosis() { 358 fmt.Println(`Kurtosis is a measure of the 'peakedness' of a distribution, and the 359excess kurtosis is the kurtosis above or below that of the standard normal 360distribution`) 361 x := []float64{5, 4, -3, -2} 362 kurt := ExKurtosis(x, nil) 363 fmt.Printf("ExKurtosis = %.5f\n", kurt) 364 weights := []float64{1, 2, 3, 5} 365 wKurt := ExKurtosis(x, weights) 366 fmt.Printf("Weighted ExKurtosis is %.4f", wKurt) 367 // Output: 368 // Kurtosis is a measure of the 'peakedness' of a distribution, and the 369 // excess kurtosis is the kurtosis above or below that of the standard normal 370 // distribution 371 // ExKurtosis = -5.41200 372 // Weighted ExKurtosis is -0.6779 373} 374 375func TestExKurtosis(t *testing.T) { 376 // the example does a good job, this just has to cover the panic 377 if !panics(func() { ExKurtosis(make([]float64, 3), make([]float64, 2)) }) { 378 t.Errorf("ExKurtosis did not panic with x, weights length mismatch") 379 } 380} 381 382func ExampleGeometricMean() { 383 x := []float64{8, 2, 9, 15, 4} 384 weights := []float64{2, 2, 6, 7, 1} 385 mean := Mean(x, weights) 386 gmean := GeometricMean(x, weights) 387 388 logx := make([]float64, len(x)) 389 for i, v := range x { 390 logx[i] = math.Log(v) 391 } 392 expMeanLog := math.Exp(Mean(logx, weights)) 393 fmt.Printf("The arithmetic mean is %.4f, but the geometric mean is %.4f.\n", mean, gmean) 394 fmt.Printf("The exponential of the mean of the logs is %.4f\n", expMeanLog) 395 // Output: 396 // The arithmetic mean is 10.1667, but the geometric mean is 8.7637. 397 // The exponential of the mean of the logs is 8.7637 398} 399 400func TestGeometricMean(t *testing.T) { 401 for i, test := range []struct { 402 x []float64 403 wts []float64 404 ans float64 405 }{ 406 { 407 x: []float64{2, 8}, 408 ans: 4, 409 }, 410 { 411 x: []float64{3, 81}, 412 wts: []float64{2, 1}, 413 ans: 9, 414 }, 415 } { 416 c := GeometricMean(test.x, test.wts) 417 if math.Abs(c-test.ans) > 1e-14 { 418 t.Errorf("Geometric mean mismatch case %d: Expected %v, Found %v", i, test.ans, c) 419 } 420 } 421 if !panics(func() { GeometricMean(make([]float64, 3), make([]float64, 2)) }) { 422 t.Errorf("GeometricMean did not panic with x, wts length mismatch") 423 } 424} 425 426func ExampleHarmonicMean() { 427 x := []float64{8, 2, 9, 15, 4} 428 weights := []float64{2, 2, 6, 7, 1} 429 mean := Mean(x, weights) 430 hmean := HarmonicMean(x, weights) 431 432 fmt.Printf("The arithmetic mean is %.5f, but the harmonic mean is %.4f.\n", mean, hmean) 433 // Output: 434 // The arithmetic mean is 10.16667, but the harmonic mean is 6.8354. 435} 436 437func TestHarmonicMean(t *testing.T) { 438 for i, test := range []struct { 439 x []float64 440 wts []float64 441 ans float64 442 }{ 443 { 444 x: []float64{.5, .125}, 445 ans: .2, 446 }, 447 { 448 x: []float64{.5, .125}, 449 wts: []float64{2, 1}, 450 ans: .25, 451 }, 452 } { 453 c := HarmonicMean(test.x, test.wts) 454 if math.Abs(c-test.ans) > 1e-14 { 455 t.Errorf("Harmonic mean mismatch case %d: Expected %v, Found %v", i, test.ans, c) 456 } 457 } 458 if !panics(func() { HarmonicMean(make([]float64, 3), make([]float64, 2)) }) { 459 t.Errorf("HarmonicMean did not panic with x, wts length mismatch") 460 } 461} 462 463func TestHistogram(t *testing.T) { 464 for i, test := range []struct { 465 x []float64 466 weights []float64 467 dividers []float64 468 ans []float64 469 }{ 470 { 471 x: []float64{1, 3, 5, 6, 7, 8}, 472 dividers: []float64{0, 2, 4, 6, 7, 9}, 473 ans: []float64{1, 1, 1, 1, 2}, 474 }, 475 { 476 x: []float64{1, 3, 5, 6, 7, 8}, 477 dividers: []float64{1, 2, 4, 6, 7, 9}, 478 weights: []float64{1, 2, 1, 1, 1, 2}, 479 ans: []float64{1, 2, 1, 1, 3}, 480 }, 481 { 482 x: []float64{1, 8}, 483 dividers: []float64{0, 2, 4, 6, 7, 9}, 484 weights: []float64{1, 2}, 485 ans: []float64{1, 0, 0, 0, 2}, 486 }, 487 { 488 x: []float64{1, 8}, 489 dividers: []float64{0, 2, 4, 6, 7, 9}, 490 ans: []float64{1, 0, 0, 0, 1}, 491 }, 492 { 493 x: []float64{}, 494 dividers: []float64{1, 3}, 495 ans: []float64{0}, 496 }, 497 } { 498 hist := Histogram(nil, test.dividers, test.x, test.weights) 499 if !floats.Equal(hist, test.ans) { 500 t.Errorf("Hist mismatch case %d. Expected %v, Found %v", i, test.ans, hist) 501 } 502 // Test with non-zero values 503 Histogram(hist, test.dividers, test.x, test.weights) 504 if !floats.Equal(hist, test.ans) { 505 t.Errorf("Hist mismatch case %d. Expected %v, Found %v", i, test.ans, hist) 506 } 507 } 508 // panic cases 509 for _, test := range []struct { 510 name string 511 x []float64 512 weights []float64 513 dividers []float64 514 count []float64 515 }{ 516 { 517 name: "len(x) != len(weights)", 518 x: []float64{1, 3, 5, 6, 7, 8}, 519 weights: []float64{1, 1, 1, 1}, 520 }, 521 { 522 name: "len(count) != len(dividers) - 1", 523 x: []float64{1, 3, 5, 6, 7, 8}, 524 dividers: []float64{1, 4, 9}, 525 count: make([]float64, 6), 526 }, 527 { 528 name: "dividers not sorted", 529 x: []float64{1, 3, 5, 6, 7, 8}, 530 dividers: []float64{0, -1, 0}, 531 }, 532 { 533 name: "x not sorted", 534 x: []float64{1, 5, 2, 9, 7, 8}, 535 dividers: []float64{1, 4, 9}, 536 }, 537 { 538 name: "fewer than 2 dividers", 539 x: []float64{1, 2, 3}, 540 dividers: []float64{5}, 541 }, 542 { 543 name: "x too large", 544 x: []float64{1, 2, 3}, 545 dividers: []float64{1, 3}, 546 }, 547 { 548 name: "x too small", 549 x: []float64{1, 2, 3}, 550 dividers: []float64{2, 3}, 551 }, 552 } { 553 if !panics(func() { Histogram(test.count, test.dividers, test.x, test.weights) }) { 554 t.Errorf("Histogram did not panic when %s", test.name) 555 } 556 } 557} 558 559func ExampleHistogram() { 560 x := make([]float64, 101) 561 for i := range x { 562 x[i] = 1.1 * float64(i) // x data ranges from 0 to 110 563 } 564 dividers := []float64{0, 7, 20, 100, 1000} 565 fmt.Println(`Histogram counts the amount of data in the bins specified by 566the dividers. In this data set, there are 7 data points less than 7 (between dividers[0] 567and dividers[1]), 12 data points between 7 and 20 (dividers[1] and dividers[2]), 568and 0 data points above 1000. Since dividers has length 5, there will be 4 bins.`) 569 hist := Histogram(nil, dividers, x, nil) 570 fmt.Printf("Hist = %v\n", hist) 571 572 fmt.Println() 573 fmt.Println("For ease, the floats Span function can be used to set the dividers") 574 nBins := 10 575 dividers = make([]float64, nBins+1) 576 min := floats.Min(x) 577 max := floats.Max(x) 578 // Increase the maximum divider so that the maximum value of x is contained 579 // within the last bucket. 580 max++ 581 floats.Span(dividers, min, max) 582 // Span includes the min and the max. Trim the dividers to create 10 buckets 583 hist = Histogram(nil, dividers, x, nil) 584 fmt.Printf("Hist = %v\n", hist) 585 fmt.Println() 586 fmt.Println(`Histogram also works with weighted data, and allows reusing of 587the count field in order to avoid extra garbage`) 588 weights := make([]float64, len(x)) 589 for i := range weights { 590 weights[i] = float64(i + 1) 591 } 592 Histogram(hist, dividers, x, weights) 593 fmt.Printf("Weighted Hist = %v\n", hist) 594 595 // Output: 596 // Histogram counts the amount of data in the bins specified by 597 // the dividers. In this data set, there are 7 data points less than 7 (between dividers[0] 598 // and dividers[1]), 12 data points between 7 and 20 (dividers[1] and dividers[2]), 599 // and 0 data points above 1000. Since dividers has length 5, there will be 4 bins. 600 // Hist = [7 12 72 10] 601 // 602 // For ease, the floats Span function can be used to set the dividers 603 // Hist = [11 10 10 10 10 10 10 10 10 10] 604 // 605 // Histogram also works with weighted data, and allows reusing of 606 // the count field in order to avoid extra garbage 607 // Weighted Hist = [66 165 265 365 465 565 665 765 865 965] 608} 609 610func TestJensenShannon(t *testing.T) { 611 for i, test := range []struct { 612 p []float64 613 q []float64 614 }{ 615 { 616 p: []float64{0.5, 0.1, 0.3, 0.1}, 617 q: []float64{0.1, 0.4, 0.25, 0.25}, 618 }, 619 { 620 p: []float64{0.4, 0.6, 0.0}, 621 q: []float64{0.2, 0.2, 0.6}, 622 }, 623 { 624 p: []float64{0.1, 0.1, 0.0, 0.8}, 625 q: []float64{0.6, 0.3, 0.0, 0.1}, 626 }, 627 { 628 p: []float64{0.5, 0.1, 0.3, 0.1}, 629 q: []float64{0.5, 0, 0.25, 0.25}, 630 }, 631 { 632 p: []float64{0.5, 0.1, 0, 0.4}, 633 q: []float64{0.1, 0.4, 0.25, 0.25}, 634 }, 635 } { 636 637 m := make([]float64, len(test.p)) 638 p := test.p 639 q := test.q 640 floats.Add(m, p) 641 floats.Add(m, q) 642 floats.Scale(0.5, m) 643 644 js1 := 0.5*KullbackLeibler(p, m) + 0.5*KullbackLeibler(q, m) 645 js2 := JensenShannon(p, q) 646 647 if math.IsNaN(js2) { 648 t.Errorf("In case %v, JS distance is NaN", i) 649 } 650 651 if math.Abs(js1-js2) > 1e-14 { 652 t.Errorf("JS mismatch case %v. Expected %v, found %v.", i, js1, js2) 653 } 654 } 655 if !panics(func() { JensenShannon(make([]float64, 3), make([]float64, 2)) }) { 656 t.Errorf("JensenShannon did not panic with p, q length mismatch") 657 } 658} 659 660func TestKolmogorovSmirnov(t *testing.T) { 661 for i, test := range []struct { 662 x []float64 663 xWeights []float64 664 y []float64 665 yWeights []float64 666 dist float64 667 }{ 668 669 { 670 dist: 0, 671 }, 672 { 673 x: []float64{1}, 674 dist: 1, 675 }, 676 { 677 y: []float64{1}, 678 dist: 1, 679 }, 680 { 681 x: []float64{1}, 682 xWeights: []float64{8}, 683 dist: 1, 684 }, 685 { 686 y: []float64{1}, 687 yWeights: []float64{8}, 688 dist: 1, 689 }, 690 { 691 x: []float64{1}, 692 xWeights: []float64{8}, 693 y: []float64{1}, 694 yWeights: []float64{8}, 695 dist: 0, 696 }, 697 { 698 x: []float64{1, 1, 1}, 699 xWeights: []float64{2, 3, 7}, 700 y: []float64{1}, 701 yWeights: []float64{8}, 702 dist: 0, 703 }, 704 { 705 x: []float64{1, 1, 1, 1, 1}, 706 y: []float64{1, 1, 1}, 707 yWeights: []float64{2, 5, 2}, 708 dist: 0, 709 }, 710 711 { 712 x: []float64{1, 2, 3}, 713 y: []float64{1, 2, 3}, 714 dist: 0, 715 }, 716 { 717 x: []float64{1, 2, 3}, 718 y: []float64{1, 2, 3}, 719 yWeights: []float64{1, 1, 1}, 720 dist: 0, 721 }, 722 723 { 724 x: []float64{1, 2, 3}, 725 xWeights: []float64{1, 1, 1}, 726 y: []float64{1, 2, 3}, 727 yWeights: []float64{1, 1, 1}, 728 dist: 0, 729 }, 730 { 731 x: []float64{1, 2}, 732 xWeights: []float64{2, 5}, 733 y: []float64{1, 1, 2, 2, 2, 2, 2}, 734 dist: 0, 735 }, 736 { 737 x: []float64{1, 1, 2, 2, 2, 2, 2}, 738 y: []float64{1, 2}, 739 yWeights: []float64{2, 5}, 740 dist: 0, 741 }, 742 { 743 x: []float64{1, 1, 2, 2, 2}, 744 xWeights: []float64{0.5, 1.5, 1, 2, 2}, 745 y: []float64{1, 2}, 746 yWeights: []float64{2, 5}, 747 dist: 0, 748 }, 749 { 750 x: []float64{1, 2, 3, 4}, 751 y: []float64{5, 6}, 752 dist: 1, 753 }, 754 { 755 x: []float64{5, 6}, 756 y: []float64{1, 2, 3, 4}, 757 dist: 1, 758 }, 759 { 760 x: []float64{5, 6}, 761 xWeights: []float64{8, 7}, 762 y: []float64{1, 2, 3, 4}, 763 dist: 1, 764 }, 765 { 766 x: []float64{5, 6}, 767 xWeights: []float64{8, 7}, 768 y: []float64{1, 2, 3, 4}, 769 yWeights: []float64{9, 2, 1, 6}, 770 dist: 1, 771 }, 772 { 773 x: []float64{-4, 5, 6}, 774 xWeights: []float64{0, 8, 7}, 775 y: []float64{1, 2, 3, 4}, 776 yWeights: []float64{9, 2, 1, 6}, 777 dist: 1, 778 }, 779 { 780 x: []float64{-4, -2, -2, 5, 6}, 781 xWeights: []float64{0, 0, 0, 8, 7}, 782 y: []float64{1, 2, 3, 4}, 783 yWeights: []float64{9, 2, 1, 6}, 784 dist: 1, 785 }, 786 { 787 x: []float64{1, 2, 3}, 788 y: []float64{1, 1, 3}, 789 dist: 1.0 / 3.0, 790 }, 791 { 792 x: []float64{1, 2, 3}, 793 y: []float64{1, 3}, 794 yWeights: []float64{2, 1}, 795 dist: 1.0 / 3.0, 796 }, 797 { 798 x: []float64{1, 2, 3}, 799 xWeights: []float64{2, 2, 2}, 800 y: []float64{1, 3}, 801 yWeights: []float64{2, 1}, 802 dist: 1.0 / 3.0, 803 }, 804 { 805 x: []float64{2, 3, 4}, 806 y: []float64{1, 5}, 807 dist: 1.0 / 2.0, 808 }, 809 { 810 x: []float64{1, 2, math.NaN()}, 811 y: []float64{1, 1, 3}, 812 dist: math.NaN(), 813 }, 814 { 815 x: []float64{1, 2, 3}, 816 y: []float64{1, 1, math.NaN()}, 817 dist: math.NaN(), 818 }, 819 } { 820 dist := KolmogorovSmirnov(test.x, test.xWeights, test.y, test.yWeights) 821 if math.Abs(dist-test.dist) > 1e-14 && !(math.IsNaN(test.dist) && math.IsNaN(dist)) { 822 t.Errorf("Distance mismatch case %v: Expected: %v, Found: %v", i, test.dist, dist) 823 } 824 } 825 // panic cases 826 for _, test := range []struct { 827 name string 828 x []float64 829 xWeights []float64 830 y []float64 831 yWeights []float64 832 }{ 833 { 834 name: "len(x) != len(xWeights)", 835 x: []float64{1, 3, 5, 6, 7, 8}, 836 xWeights: []float64{1, 1, 1, 1}, 837 }, 838 { 839 name: "len(y) != len(yWeights)", 840 x: []float64{1, 3, 5, 6, 7, 8}, 841 y: []float64{1, 3, 5, 6, 7, 8}, 842 yWeights: []float64{1, 1, 1, 1}, 843 }, 844 { 845 name: "x not sorted", 846 x: []float64{10, 3, 5, 6, 7, 8}, 847 y: []float64{1, 3, 5, 6, 7, 8}, 848 }, 849 { 850 name: "y not sorted", 851 x: []float64{1, 3, 5, 6, 7, 8}, 852 y: []float64{10, 3, 5, 6, 7, 8}, 853 }, 854 } { 855 if !panics(func() { KolmogorovSmirnov(test.x, test.xWeights, test.y, test.yWeights) }) { 856 t.Errorf("KolmogorovSmirnov did not panic when %s", test.name) 857 } 858 } 859} 860 861func ExampleKullbackLeibler() { 862 863 p := []float64{0.05, 0.1, 0.9, 0.05} 864 q := []float64{0.2, 0.4, 0.25, 0.15} 865 s := []float64{0, 0, 1, 0} 866 867 klPQ := KullbackLeibler(p, q) 868 klPS := KullbackLeibler(p, s) 869 klPP := KullbackLeibler(p, p) 870 871 fmt.Println("Kullback-Leibler is one measure of the difference between two distributions") 872 fmt.Printf("The K-L distance between p and q is %.4f\n", klPQ) 873 fmt.Println("It is impossible for s and p to be the same distribution, because") 874 fmt.Println("the first bucket has zero probability in s and non-zero in p. Thus,") 875 fmt.Printf("the K-L distance between them is %.4f\n", klPS) 876 fmt.Printf("The K-L distance between identical distributions is %.4f\n", klPP) 877 878 // Kullback-Leibler is one measure of the difference between two distributions 879 // The K-L distance between p and q is 0.8900 880 // It is impossible for s and p to be the same distribution, because 881 // the first bucket has zero probability in s and non-zero in p. Thus, 882 // the K-L distance between them is +Inf 883 // The K-L distance between identical distributions is 0.0000 884} 885 886func TestKullbackLeibler(t *testing.T) { 887 if !panics(func() { KullbackLeibler(make([]float64, 3), make([]float64, 2)) }) { 888 t.Errorf("KullbackLeibler did not panic with p, q length mismatch") 889 } 890} 891 892var linearRegressionTests = []struct { 893 name string 894 895 x, y []float64 896 weights []float64 897 origin bool 898 899 alpha float64 900 beta float64 901 r float64 902 903 tol float64 904}{ 905 { 906 name: "faithful", 907 908 x: faithful.waiting, 909 y: faithful.eruptions, 910 911 // Values calculated by R using lm(eruptions ~ waiting, data=faithful). 912 alpha: -1.87402, 913 beta: 0.07563, 914 r: 0.8114608, 915 916 tol: 1e-5, 917 }, 918 { 919 name: "faithful through origin", 920 921 x: faithful.waiting, 922 y: faithful.eruptions, 923 origin: true, 924 925 // Values calculated by R using lm(eruptions ~ waiting - 1, data=faithful). 926 alpha: 0, 927 beta: 0.05013, 928 r: 0.9726036, 929 930 tol: 1e-5, 931 }, 932 { 933 name: "faithful explicit weights", 934 935 x: faithful.waiting, 936 y: faithful.eruptions, 937 weights: func() []float64 { 938 w := make([]float64, len(faithful.eruptions)) 939 for i := range w { 940 w[i] = 1 941 } 942 return w 943 }(), 944 945 // Values calculated by R using lm(eruptions ~ waiting, data=faithful). 946 alpha: -1.87402, 947 beta: 0.07563, 948 r: 0.8114608, 949 950 tol: 1e-5, 951 }, 952 { 953 name: "faithful non-uniform weights", 954 955 x: faithful.waiting, 956 y: faithful.eruptions, 957 weights: faithful.waiting, // Just an arbitrary set of non-uniform weights. 958 959 // Values calculated by R using lm(eruptions ~ waiting, data=faithful, weights=faithful$waiting). 960 alpha: -1.79268, 961 beta: 0.07452, 962 r: 0.7840372, 963 964 tol: 1e-5, 965 }, 966} 967 968func TestLinearRegression(t *testing.T) { 969 for _, test := range linearRegressionTests { 970 alpha, beta := LinearRegression(test.x, test.y, test.weights, test.origin) 971 var r float64 972 if test.origin { 973 r = RNoughtSquared(test.x, test.y, test.weights, beta) 974 } else { 975 r = RSquared(test.x, test.y, test.weights, alpha, beta) 976 ests := make([]float64, len(test.y)) 977 for i, x := range test.x { 978 ests[i] = alpha + beta*x 979 } 980 rvals := RSquaredFrom(ests, test.y, test.weights) 981 if r != rvals { 982 t.Errorf("%s: RSquared and RSquaredFrom mismatch: %v != %v", test.name, r, rvals) 983 } 984 } 985 if !scalar.EqualWithinAbsOrRel(alpha, test.alpha, test.tol, test.tol) { 986 t.Errorf("%s: unexpected alpha estimate: want:%v got:%v", test.name, test.alpha, alpha) 987 } 988 if !scalar.EqualWithinAbsOrRel(beta, test.beta, test.tol, test.tol) { 989 t.Errorf("%s: unexpected beta estimate: want:%v got:%v", test.name, test.beta, beta) 990 } 991 if !scalar.EqualWithinAbsOrRel(r, test.r, test.tol, test.tol) { 992 t.Errorf("%s: unexpected r estimate: want:%v got:%v", test.name, test.r, r) 993 } 994 } 995} 996 997func BenchmarkLinearRegression(b *testing.B) { 998 rnd := rand.New(rand.NewSource(1)) 999 slope, offset := 2.0, 3.0 1000 1001 maxn := 10000 1002 xs := make([]float64, maxn) 1003 ys := make([]float64, maxn) 1004 weights := make([]float64, maxn) 1005 for i := range xs { 1006 x := rnd.Float64() 1007 xs[i] = x 1008 ys[i] = slope*x + offset 1009 weights[i] = rnd.Float64() 1010 } 1011 1012 for _, n := range []int{10, 100, 1000, maxn} { 1013 for _, weighted := range []bool{true, false} { 1014 for _, origin := range []bool{true, false} { 1015 name := "n" + strconv.Itoa(n) 1016 if weighted { 1017 name += "wt" 1018 } else { 1019 name += "wf" 1020 } 1021 if origin { 1022 name += "ot" 1023 } else { 1024 name += "of" 1025 } 1026 b.Run(name, func(b *testing.B) { 1027 for i := 0; i < b.N; i++ { 1028 var ws []float64 1029 if weighted { 1030 ws = weights[:n] 1031 } 1032 LinearRegression(xs[:n], ys[:n], ws, origin) 1033 } 1034 }) 1035 } 1036 } 1037 } 1038} 1039 1040func TestChiSquare(t *testing.T) { 1041 for i, test := range []struct { 1042 p []float64 1043 q []float64 1044 res float64 1045 }{ 1046 { 1047 p: []float64{16, 18, 16, 14, 12, 12}, 1048 q: []float64{16, 16, 16, 16, 16, 8}, 1049 res: 3.5, 1050 }, 1051 { 1052 p: []float64{16, 18, 16, 14, 12, 12}, 1053 q: []float64{8, 20, 20, 16, 12, 12}, 1054 res: 9.25, 1055 }, 1056 { 1057 p: []float64{40, 60, 30, 45}, 1058 q: []float64{50, 50, 50, 50}, 1059 res: 12.5, 1060 }, 1061 { 1062 p: []float64{40, 60, 30, 45, 0, 0}, 1063 q: []float64{50, 50, 50, 50, 0, 0}, 1064 res: 12.5, 1065 }, 1066 } { 1067 resultpq := ChiSquare(test.p, test.q) 1068 1069 if math.Abs(resultpq-test.res) > 1e-10 { 1070 t.Errorf("ChiSquare distance mismatch in case %d. Expected %v, Found %v", i, test.res, resultpq) 1071 } 1072 } 1073 if !panics(func() { ChiSquare(make([]float64, 2), make([]float64, 3)) }) { 1074 t.Errorf("ChiSquare did not panic with length mismatch") 1075 } 1076} 1077 1078// panics returns true if the called function panics during evaluation. 1079func panics(fun func()) (b bool) { 1080 defer func() { 1081 err := recover() 1082 if err != nil { 1083 b = true 1084 } 1085 }() 1086 fun() 1087 return 1088} 1089 1090func TestBhattacharyya(t *testing.T) { 1091 for i, test := range []struct { 1092 p []float64 1093 q []float64 1094 res float64 1095 }{ 1096 { 1097 p: []float64{0.5, 0.1, 0.3, 0.1}, 1098 q: []float64{0.1, 0.4, 0.25, 0.25}, 1099 res: 0.15597338718671386, 1100 }, 1101 { 1102 p: []float64{0.4, 0.6, 0.0}, 1103 q: []float64{0.2, 0.2, 0.6}, 1104 res: 0.46322207765351153, 1105 }, 1106 { 1107 p: []float64{0.1, 0.1, 0.0, 0.8}, 1108 q: []float64{0.6, 0.3, 0.0, 0.1}, 1109 res: 0.3552520032137785, 1110 }, 1111 } { 1112 resultpq := Bhattacharyya(test.p, test.q) 1113 resultqp := Bhattacharyya(test.q, test.p) 1114 1115 if math.Abs(resultpq-test.res) > 1e-10 { 1116 t.Errorf("Bhattacharyya distance mismatch in case %d. Expected %v, Found %v", i, test.res, resultpq) 1117 } 1118 if math.Abs(resultpq-resultqp) > 1e-10 { 1119 t.Errorf("Bhattacharyya distance is assymmetric in case %d.", i) 1120 } 1121 } 1122 // Bhattacharyya should panic if the inputs have different length 1123 if !panics(func() { Bhattacharyya(make([]float64, 2), make([]float64, 3)) }) { 1124 t.Errorf("Bhattacharyya did not panic with length mismatch") 1125 } 1126} 1127 1128func TestHellinger(t *testing.T) { 1129 for i, test := range []struct { 1130 p []float64 1131 q []float64 1132 res float64 1133 }{ 1134 { 1135 p: []float64{0.5, 0.1, 0.3, 0.1}, 1136 q: []float64{0.1, 0.4, 0.25, 0.25}, 1137 res: 0.3800237367441919, 1138 }, 1139 { 1140 p: []float64{0.4, 0.6, 0.0}, 1141 q: []float64{0.2, 0.2, 0.6}, 1142 res: 0.6088900771170487, 1143 }, 1144 { 1145 p: []float64{0.1, 0.1, 0.0, 0.8}, 1146 q: []float64{0.6, 0.3, 0.0, 0.1}, 1147 res: 0.5468118803484205, 1148 }, 1149 } { 1150 resultpq := Hellinger(test.p, test.q) 1151 resultqp := Hellinger(test.q, test.p) 1152 1153 if math.Abs(resultpq-test.res) > 1e-10 { 1154 t.Errorf("Hellinger distance mismatch in case %d. Expected %v, Found %v", i, test.res, resultpq) 1155 } 1156 if math.Abs(resultpq-resultqp) > 1e-10 { 1157 t.Errorf("Hellinger distance is assymmetric in case %d.", i) 1158 } 1159 } 1160 if !panics(func() { Hellinger(make([]float64, 2), make([]float64, 3)) }) { 1161 t.Errorf("Hellinger did not panic with length mismatch") 1162 } 1163} 1164 1165func ExampleMean() { 1166 x := []float64{8.2, -6, 5, 7} 1167 mean := Mean(x, nil) 1168 fmt.Printf("The mean of the samples is %.4f\n", mean) 1169 w := []float64{2, 6, 3, 5} 1170 weightedMean := Mean(x, w) 1171 fmt.Printf("The weighted mean of the samples is %.4f\n", weightedMean) 1172 x2 := []float64{8.2, 8.2, -6, -6, -6, -6, -6, -6, 5, 5, 5, 7, 7, 7, 7, 7} 1173 mean2 := Mean(x2, nil) 1174 fmt.Printf("The mean of x2 is %.4f\n", mean2) 1175 fmt.Println("The weights act as if there were more samples of that number") 1176 // Output: 1177 // The mean of the samples is 3.5500 1178 // The weighted mean of the samples is 1.9000 1179 // The mean of x2 is 1.9000 1180 // The weights act as if there were more samples of that number 1181} 1182func TestMean(t *testing.T) { 1183 if !panics(func() { Mean(make([]float64, 3), make([]float64, 2)) }) { 1184 t.Errorf("Mean did not panic with x, weights length mismatch") 1185 } 1186} 1187 1188func TestMode(t *testing.T) { 1189 for i, test := range []struct { 1190 x []float64 1191 weights []float64 1192 ans float64 1193 count float64 1194 }{ 1195 {}, 1196 { 1197 x: []float64{1, 6, 1, 9, -2}, 1198 ans: 1, 1199 count: 2, 1200 }, 1201 { 1202 x: []float64{1, 6, 1, 9, -2}, 1203 weights: []float64{1, 7, 3, 5, 0}, 1204 ans: 6, 1205 count: 7, 1206 }, 1207 } { 1208 m, count := Mode(test.x, test.weights) 1209 if test.ans != m { 1210 t.Errorf("Mode mismatch case %d. Expected %v, found %v", i, test.ans, m) 1211 } 1212 if test.count != count { 1213 t.Errorf("Mode count mismatch case %d. Expected %v, found %v", i, test.count, count) 1214 } 1215 } 1216 if !panics(func() { Mode(make([]float64, 3), make([]float64, 2)) }) { 1217 t.Errorf("Mode did not panic with x, weights length mismatch") 1218 } 1219} 1220 1221func TestMixedMoment(t *testing.T) { 1222 for i, test := range []struct { 1223 x, y, weights []float64 1224 r, s float64 1225 ans float64 1226 }{ 1227 { 1228 x: []float64{10, 2, 1, 8, 5}, 1229 y: []float64{8, 15, 1, 6, 3}, 1230 r: 1, 1231 s: 1, 1232 ans: 0.48, 1233 }, 1234 { 1235 x: []float64{10, 2, 1, 8, 5}, 1236 y: []float64{8, 15, 1, 6, 3}, 1237 weights: []float64{1, 1, 1, 1, 1}, 1238 r: 1, 1239 s: 1, 1240 ans: 0.48, 1241 }, 1242 { 1243 x: []float64{10, 2, 1, 8, 5}, 1244 y: []float64{8, 15, 1, 6, 3}, 1245 weights: []float64{2, 3, 0.2, 8, 4}, 1246 r: 1, 1247 s: 1, 1248 ans: -4.786371011357490, 1249 }, 1250 { 1251 x: []float64{10, 2, 1, 8, 5}, 1252 y: []float64{8, 15, 1, 6, 3}, 1253 weights: []float64{2, 3, 0.2, 8, 4}, 1254 r: 2, 1255 s: 3, 1256 ans: 1.598600579313326e+03, 1257 }, 1258 } { 1259 m := BivariateMoment(test.r, test.s, test.x, test.y, test.weights) 1260 if math.Abs(test.ans-m) > 1e-14 { 1261 t.Errorf("Moment mismatch case %d. Expected %v, found %v", i, test.ans, m) 1262 } 1263 } 1264 if !panics(func() { BivariateMoment(1, 1, make([]float64, 3), make([]float64, 2), nil) }) { 1265 t.Errorf("Moment did not panic with x, y length mismatch") 1266 } 1267 if !panics(func() { BivariateMoment(1, 1, make([]float64, 2), make([]float64, 3), nil) }) { 1268 t.Errorf("Moment did not panic with x, y length mismatch") 1269 } 1270 if !panics(func() { BivariateMoment(1, 1, make([]float64, 2), make([]float64, 2), make([]float64, 3)) }) { 1271 t.Errorf("Moment did not panic with x, weights length mismatch") 1272 } 1273 if !panics(func() { BivariateMoment(1, 1, make([]float64, 2), make([]float64, 2), make([]float64, 1)) }) { 1274 t.Errorf("Moment did not panic with x, weights length mismatch") 1275 } 1276} 1277 1278func TestMoment(t *testing.T) { 1279 for i, test := range []struct { 1280 x []float64 1281 weights []float64 1282 moment float64 1283 ans float64 1284 }{ 1285 { 1286 x: []float64{6, 2, 4, 8, 10}, 1287 moment: 5, 1288 ans: 0, 1289 }, 1290 { 1291 x: []float64{6, 2, 4, 8, 10}, 1292 weights: []float64{1, 2, 2, 2, 1}, 1293 moment: 5, 1294 ans: 121.875, 1295 }, 1296 } { 1297 m := Moment(test.moment, test.x, test.weights) 1298 if math.Abs(test.ans-m) > 1e-14 { 1299 t.Errorf("Moment mismatch case %d. Expected %v, found %v", i, test.ans, m) 1300 } 1301 } 1302 if !panics(func() { Moment(1, make([]float64, 3), make([]float64, 2)) }) { 1303 t.Errorf("Moment did not panic with x, weights length mismatch") 1304 } 1305 if !panics(func() { Moment(1, make([]float64, 2), make([]float64, 3)) }) { 1306 t.Errorf("Moment did not panic with x, weights length mismatch") 1307 } 1308} 1309 1310func TestMomentAbout(t *testing.T) { 1311 for i, test := range []struct { 1312 x []float64 1313 weights []float64 1314 moment float64 1315 mean float64 1316 ans float64 1317 }{ 1318 { 1319 x: []float64{6, 2, 4, 8, 9}, 1320 mean: 3, 1321 moment: 5, 1322 ans: 2.2288e3, 1323 }, 1324 { 1325 x: []float64{6, 2, 4, 8, 9}, 1326 weights: []float64{1, 2, 2, 2, 1}, 1327 mean: 3, 1328 moment: 5, 1329 ans: 1.783625e3, 1330 }, 1331 } { 1332 m := MomentAbout(test.moment, test.x, test.mean, test.weights) 1333 if math.Abs(test.ans-m) > 1e-14 { 1334 t.Errorf("MomentAbout mismatch case %d. Expected %v, found %v", i, test.ans, m) 1335 } 1336 } 1337 if !panics(func() { MomentAbout(1, make([]float64, 3), 0, make([]float64, 2)) }) { 1338 t.Errorf("MomentAbout did not panic with x, weights length mismatch") 1339 } 1340} 1341 1342func TestCDF(t *testing.T) { 1343 cumulantKinds := []CumulantKind{Empirical} 1344 for i, test := range []struct { 1345 q []float64 1346 x []float64 1347 weights []float64 1348 ans [][]float64 1349 }{ 1350 {}, 1351 { 1352 q: []float64{0, 0.9, 1, 1.1, 2.9, 3, 3.1, 4.9, 5, 5.1}, 1353 x: []float64{1, 2, 3, 4, 5}, 1354 ans: [][]float64{{0, 0, 0.2, 0.2, 0.4, 0.6, 0.6, 0.8, 1, 1}}, 1355 }, 1356 { 1357 q: []float64{0, 0.9, 1, 1.1, 2.9, 3, 3.1, 4.9, 5, 5.1}, 1358 x: []float64{1, 2, 3, 4, 5}, 1359 weights: []float64{1, 1, 1, 1, 1}, 1360 ans: [][]float64{{0, 0, 0.2, 0.2, 0.4, 0.6, 0.6, 0.8, 1, 1}}, 1361 }, 1362 { 1363 q: []float64{0, 0.9, 1}, 1364 x: []float64{math.NaN()}, 1365 ans: [][]float64{{math.NaN(), math.NaN(), math.NaN()}}, 1366 }, 1367 } { 1368 copyX := make([]float64, len(test.x)) 1369 copy(copyX, test.x) 1370 var copyW []float64 1371 if test.weights != nil { 1372 copyW = make([]float64, len(test.weights)) 1373 copy(copyW, test.weights) 1374 } 1375 for j, q := range test.q { 1376 for k, kind := range cumulantKinds { 1377 v := CDF(q, kind, test.x, test.weights) 1378 if !floats.Equal(copyX, test.x) && !math.IsNaN(v) { 1379 t.Errorf("x changed for case %d kind %d percentile %v", i, k, q) 1380 } 1381 if !floats.Equal(copyW, test.weights) { 1382 t.Errorf("x changed for case %d kind %d percentile %v", i, k, q) 1383 } 1384 if v != test.ans[k][j] && !(math.IsNaN(v) && math.IsNaN(test.ans[k][j])) { 1385 t.Errorf("mismatch case %d kind %d percentile %v. Expected: %v, found: %v", i, k, q, test.ans[k][j], v) 1386 } 1387 } 1388 } 1389 } 1390 1391 // these test cases should all result in a panic 1392 for i, test := range []struct { 1393 name string 1394 q float64 1395 kind CumulantKind 1396 x []float64 1397 weights []float64 1398 }{ 1399 { 1400 name: "len(x) != len(weights)", 1401 q: 1.5, 1402 kind: Empirical, 1403 x: []float64{1, 2, 3, 4, 5}, 1404 weights: []float64{1, 2, 3}, 1405 }, 1406 { 1407 name: "unsorted x", 1408 q: 1.5, 1409 kind: Empirical, 1410 x: []float64{3, 2, 1}, 1411 }, 1412 { 1413 name: "unknown CumulantKind", 1414 q: 1.5, 1415 kind: CumulantKind(1000), // bogus 1416 x: []float64{1, 2, 3}, 1417 }, 1418 } { 1419 if !panics(func() { CDF(test.q, test.kind, test.x, test.weights) }) { 1420 t.Errorf("did not panic as expected with %s for case %d kind %d percentile %v x %v weights %v", test.name, i, test.kind, test.q, test.x, test.weights) 1421 } 1422 } 1423 1424} 1425 1426func TestQuantile(t *testing.T) { 1427 cumulantKinds := []CumulantKind{ 1428 Empirical, 1429 LinInterp, 1430 } 1431 for i, test := range []struct { 1432 p []float64 1433 x []float64 1434 w []float64 1435 ans [][]float64 1436 }{ 1437 { 1438 p: []float64{0, 0.05, 0.1, 0.15, 0.45, 0.5, 0.55, 0.85, 0.9, 0.95, 1}, 1439 x: []float64{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, 1440 w: nil, 1441 ans: [][]float64{ 1442 {1, 1, 1, 2, 5, 5, 6, 9, 9, 10, 10}, 1443 {1, 1, 1, 1.5, 4.5, 5, 5.5, 8.5, 9, 9.5, 10}, 1444 }, 1445 }, 1446 { 1447 p: []float64{0, 0.05, 0.1, 0.15, 0.45, 0.5, 0.55, 0.85, 0.9, 0.95, 1}, 1448 x: []float64{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, 1449 w: []float64{3, 3, 3, 3, 3, 3, 3, 3, 3, 3}, 1450 ans: [][]float64{ 1451 {1, 1, 1, 2, 5, 5, 6, 9, 9, 10, 10}, 1452 {1, 1, 1, 1.5, 4.5, 5, 5.5, 8.5, 9, 9.5, 10}, 1453 }, 1454 }, 1455 { 1456 p: []float64{0, 0.05, 0.1, 0.15, 0.45, 0.5, 0.55, 0.85, 0.9, 0.95, 1}, 1457 x: []float64{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, 1458 w: []float64{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0}, 1459 ans: [][]float64{ 1460 {1, 2, 3, 4, 7, 7, 8, 10, 10, 10, 10}, 1461 {1, 1.875, 2.833333333333333, 3.5625, 6.535714285714286, 6.928571428571429, 7.281250000000001, 9.175, 9.45, 9.725, 10}, 1462 }, 1463 }, 1464 { 1465 p: []float64{0.5}, 1466 x: []float64{1, 2, 3, 4, 5, 6, 7, 8, math.NaN(), 10}, 1467 ans: [][]float64{ 1468 {math.NaN()}, 1469 {math.NaN()}, 1470 }, 1471 }, 1472 } { 1473 copyX := make([]float64, len(test.x)) 1474 copy(copyX, test.x) 1475 var copyW []float64 1476 if test.w != nil { 1477 copyW = make([]float64, len(test.w)) 1478 copy(copyW, test.w) 1479 } 1480 for j, p := range test.p { 1481 for k, kind := range cumulantKinds { 1482 v := Quantile(p, kind, test.x, test.w) 1483 if !floats.Same(copyX, test.x) { 1484 t.Errorf("x changed for case %d kind %d percentile %v", i, k, p) 1485 } 1486 if !floats.Same(copyW, test.w) { 1487 t.Errorf("x changed for case %d kind %d percentile %v", i, k, p) 1488 } 1489 if v != test.ans[k][j] && !(math.IsNaN(v) && math.IsNaN(test.ans[k][j])) { 1490 t.Errorf("mismatch case %d kind %d percentile %v. Expected: %v, found: %v", i, k, p, test.ans[k][j], v) 1491 } 1492 } 1493 } 1494 } 1495} 1496 1497func TestQuantileInvalidInput(t *testing.T) { 1498 cumulantKinds := []CumulantKind{ 1499 Empirical, 1500 LinInterp, 1501 } 1502 for _, test := range []struct { 1503 name string 1504 p float64 1505 x []float64 1506 w []float64 1507 }{ 1508 { 1509 name: "p < 0", 1510 p: -1, 1511 }, 1512 { 1513 name: "p > 1", 1514 p: 2, 1515 }, 1516 { 1517 name: "p is NaN", 1518 p: math.NaN(), 1519 }, 1520 { 1521 name: "len(x) != len(weights)", 1522 p: .5, 1523 x: make([]float64, 4), 1524 w: make([]float64, 2), 1525 }, 1526 { 1527 name: "x not sorted", 1528 p: .5, 1529 x: []float64{3, 2, 1}, 1530 }, 1531 } { 1532 for _, kind := range cumulantKinds { 1533 if !panics(func() { Quantile(test.p, kind, test.x, test.w) }) { 1534 t.Errorf("Quantile did not panic when %s", test.name) 1535 } 1536 } 1537 } 1538} 1539 1540func TestQuantileInvalidCumulantKind(t *testing.T) { 1541 if !panics(func() { Quantile(0.5, CumulantKind(1000), []float64{1, 2, 3}, nil) }) { 1542 t.Errorf("Quantile did not panic when CumulantKind is unknown") 1543 } 1544} 1545 1546func ExampleStdDev() { 1547 x := []float64{8, 2, -9, 15, 4} 1548 stdev := StdDev(x, nil) 1549 fmt.Printf("The standard deviation of the samples is %.4f\n", stdev) 1550 1551 weights := []float64{2, 2, 6, 7, 1} 1552 weightedStdev := StdDev(x, weights) 1553 fmt.Printf("The weighted standard deviation of the samples is %.4f\n", weightedStdev) 1554 // Output: 1555 // The standard deviation of the samples is 8.8034 1556 // The weighted standard deviation of the samples is 10.5733 1557} 1558 1559func ExamplePopStdDev() { 1560 x := []float64{8, 2, -9, 15, 4} 1561 stdev := PopStdDev(x, nil) 1562 fmt.Printf("The standard deviation of the population is %.4f\n", stdev) 1563 1564 weights := []float64{2, 2, 6, 7, 1} 1565 weightedStdev := PopStdDev(x, weights) 1566 fmt.Printf("The weighted standard deviation of the population is %.4f\n", weightedStdev) 1567 // Output: 1568 // The standard deviation of the population is 7.8740 1569 // The weighted standard deviation of the population is 10.2754 1570} 1571 1572func ExampleStdErr() { 1573 x := []float64{8, 2, -9, 15, 4} 1574 weights := []float64{2, 2, 6, 7, 1} 1575 mean := Mean(x, weights) 1576 stdev := StdDev(x, weights) 1577 nSamples := floats.Sum(weights) 1578 stdErr := StdErr(stdev, nSamples) 1579 fmt.Printf("The standard deviation is %.4f and there are %g samples, so the mean\nis likely %.4f ± %.4f.", stdev, nSamples, mean, stdErr) 1580 // Output: 1581 // The standard deviation is 10.5733 and there are 18 samples, so the mean 1582 // is likely 4.1667 ± 2.4921. 1583} 1584 1585func TestSkew(t *testing.T) { 1586 for i, test := range []struct { 1587 x []float64 1588 weights []float64 1589 ans float64 1590 }{ 1591 { 1592 x: []float64{8, 3, 7, 8, 4}, 1593 weights: nil, 1594 ans: -0.581456499151665, 1595 }, 1596 { 1597 x: []float64{8, 3, 7, 8, 4}, 1598 weights: []float64{1, 1, 1, 1, 1}, 1599 ans: -0.581456499151665, 1600 }, 1601 { 1602 x: []float64{8, 3, 7, 8, 4}, 1603 weights: []float64{2, 1, 2, 1, 1}, 1604 ans: -1.12066646837198, 1605 }, 1606 } { 1607 skew := Skew(test.x, test.weights) 1608 if math.Abs(skew-test.ans) > 1e-14 { 1609 t.Errorf("Skew mismatch case %d. Expected %v, Found %v", i, test.ans, skew) 1610 } 1611 } 1612 if !panics(func() { Skew(make([]float64, 3), make([]float64, 2)) }) { 1613 t.Errorf("Skew did not panic with x, weights length mismatch") 1614 } 1615} 1616 1617func TestSortWeighted(t *testing.T) { 1618 for i, test := range []struct { 1619 x []float64 1620 w []float64 1621 ansx []float64 1622 answ []float64 1623 }{ 1624 { 1625 x: []float64{8, 3, 7, 8, 4}, 1626 ansx: []float64{3, 4, 7, 8, 8}, 1627 }, 1628 { 1629 x: []float64{8, 3, 7, 8, 4}, 1630 w: []float64{.5, 1, 1, .5, 1}, 1631 ansx: []float64{3, 4, 7, 8, 8}, 1632 answ: []float64{1, 1, 1, .5, .5}, 1633 }, 1634 } { 1635 SortWeighted(test.x, test.w) 1636 if !floats.Same(test.x, test.ansx) { 1637 t.Errorf("SortWeighted mismatch case %d. Expected x %v, Found x %v", i, test.ansx, test.x) 1638 } 1639 if !(test.w == nil) && !floats.Same(test.w, test.answ) { 1640 t.Errorf("SortWeighted mismatch case %d. Expected w %v, Found w %v", i, test.answ, test.w) 1641 } 1642 } 1643 if !panics(func() { SortWeighted(make([]float64, 3), make([]float64, 2)) }) { 1644 t.Errorf("SortWeighted did not panic with x, weights length mismatch") 1645 } 1646} 1647 1648func TestSortWeightedLabeled(t *testing.T) { 1649 for i, test := range []struct { 1650 x []float64 1651 l []bool 1652 w []float64 1653 ansx []float64 1654 ansl []bool 1655 answ []float64 1656 }{ 1657 { 1658 x: []float64{8, 3, 7, 8, 4}, 1659 ansx: []float64{3, 4, 7, 8, 8}, 1660 }, 1661 { 1662 x: []float64{8, 3, 7, 8, 4}, 1663 w: []float64{.5, 1, 1, .5, 1}, 1664 ansx: []float64{3, 4, 7, 8, 8}, 1665 answ: []float64{1, 1, 1, .5, .5}, 1666 }, 1667 { 1668 x: []float64{8, 3, 7, 8, 4}, 1669 l: []bool{false, false, true, false, true}, 1670 ansx: []float64{3, 4, 7, 8, 8}, 1671 ansl: []bool{false, true, true, false, false}, 1672 }, 1673 { 1674 x: []float64{8, 3, 7, 8, 4}, 1675 l: []bool{false, false, true, false, true}, 1676 w: []float64{.5, 1, 1, .5, 1}, 1677 ansx: []float64{3, 4, 7, 8, 8}, 1678 ansl: []bool{false, true, true, false, false}, 1679 answ: []float64{1, 1, 1, .5, .5}, 1680 }, 1681 } { 1682 SortWeightedLabeled(test.x, test.l, test.w) 1683 if !floats.Same(test.x, test.ansx) { 1684 t.Errorf("SortWeightedLabelled mismatch case %d. Expected x %v, Found x %v", i, test.ansx, test.x) 1685 } 1686 if (test.l != nil) && !reflect.DeepEqual(test.l, test.ansl) { 1687 t.Errorf("SortWeightedLabelled mismatch case %d. Expected l %v, Found l %v", i, test.ansl, test.l) 1688 } 1689 if (test.w != nil) && !floats.Same(test.w, test.answ) { 1690 t.Errorf("SortWeightedLabelled mismatch case %d. Expected w %v, Found w %v", i, test.answ, test.w) 1691 } 1692 } 1693 if !panics(func() { SortWeightedLabeled(make([]float64, 3), make([]bool, 2), make([]float64, 3)) }) { 1694 t.Errorf("SortWeighted did not panic with x, labels length mismatch") 1695 } 1696 if !panics(func() { SortWeightedLabeled(make([]float64, 3), make([]bool, 2), nil) }) { 1697 t.Errorf("SortWeighted did not panic with x, labels length mismatch") 1698 } 1699 if !panics(func() { SortWeightedLabeled(make([]float64, 3), make([]bool, 3), make([]float64, 2)) }) { 1700 t.Errorf("SortWeighted did not panic with x, weights length mismatch") 1701 } 1702 if !panics(func() { SortWeightedLabeled(make([]float64, 3), nil, make([]float64, 2)) }) { 1703 t.Errorf("SortWeighted did not panic with x, weights length mismatch") 1704 } 1705} 1706 1707func TestVariance(t *testing.T) { 1708 for i, test := range []struct { 1709 x []float64 1710 weights []float64 1711 ans float64 1712 }{ 1713 { 1714 x: []float64{8, -3, 7, 8, -4}, 1715 weights: nil, 1716 ans: 37.7, 1717 }, 1718 { 1719 x: []float64{8, -3, 7, 8, -4}, 1720 weights: []float64{1, 1, 1, 1, 1}, 1721 ans: 37.7, 1722 }, 1723 { 1724 x: []float64{8, 3, 7, 8, 4}, 1725 weights: []float64{2, 1, 2, 1, 1}, 1726 ans: 4.2857142857142865, 1727 }, 1728 { 1729 x: []float64{1, 4, 9}, 1730 weights: []float64{1, 1.5, 1}, 1731 ans: 13.142857142857146, 1732 }, 1733 { 1734 x: []float64{1, 2, 3}, 1735 weights: []float64{1, 1.5, 1}, 1736 ans: .8, 1737 }, 1738 } { 1739 variance := Variance(test.x, test.weights) 1740 if math.Abs(variance-test.ans) > 1e-14 { 1741 t.Errorf("Variance mismatch case %d. Expected %v, Found %v", i, test.ans, variance) 1742 } 1743 } 1744 if !panics(func() { Variance(make([]float64, 3), make([]float64, 2)) }) { 1745 t.Errorf("Variance did not panic with x, weights length mismatch") 1746 } 1747} 1748 1749func TestPopVariance(t *testing.T) { 1750 for i, test := range []struct { 1751 x []float64 1752 weights []float64 1753 ans float64 1754 }{ 1755 { 1756 x: []float64{8, -3, 7, 8, -4}, 1757 weights: nil, 1758 ans: 30.16, 1759 }, 1760 { 1761 x: []float64{8, -3, 7, 8, -4}, 1762 weights: []float64{1, 1, 1, 1, 1}, 1763 ans: 30.16, 1764 }, 1765 { 1766 x: []float64{8, 3, 7, 8, 4}, 1767 weights: []float64{2, 1, 2, 1, 1}, 1768 ans: 3.6734693877551026, 1769 }, 1770 { 1771 x: []float64{1, 4, 9}, 1772 weights: []float64{1, 1.5, 1}, 1773 ans: 9.387755102040817, 1774 }, 1775 { 1776 x: []float64{1, 2, 3}, 1777 weights: []float64{1, 1.5, 1}, 1778 ans: 0.5714285714285714, 1779 }, 1780 { 1781 x: []float64{2}, 1782 weights: nil, 1783 ans: 0, 1784 }, 1785 { 1786 x: []float64{2}, 1787 weights: []float64{2}, 1788 ans: 0, 1789 }, 1790 } { 1791 variance := PopVariance(test.x, test.weights) 1792 if math.Abs(variance-test.ans) > 1e-14 { 1793 t.Errorf("PopVariance mismatch case %d. Expected %v, Found %v", i, test.ans, variance) 1794 } 1795 } 1796 if !panics(func() { PopVariance(make([]float64, 3), make([]float64, 2)) }) { 1797 t.Errorf("PopVariance did not panic with x, weights length mismatch") 1798 } 1799} 1800 1801func ExampleVariance() { 1802 x := []float64{8, 2, -9, 15, 4} 1803 variance := Variance(x, nil) 1804 fmt.Printf("The variance of the samples is %.4f\n", variance) 1805 1806 weights := []float64{2, 2, 6, 7, 1} 1807 weightedVariance := Variance(x, weights) 1808 fmt.Printf("The weighted variance of the samples is %.4f\n", weightedVariance) 1809 // Output: 1810 // The variance of the samples is 77.5000 1811 // The weighted variance of the samples is 111.7941 1812} 1813 1814func ExamplePopVariance() { 1815 x := []float64{8, 2, -9, 15, 4} 1816 variance := PopVariance(x, nil) 1817 fmt.Printf("The biased variance of the samples is %.4f\n", variance) 1818 1819 weights := []float64{2, 2, 6, 7, 1} 1820 weightedVariance := PopVariance(x, weights) 1821 fmt.Printf("The weighted biased variance of the samples is %.4f\n", weightedVariance) 1822 // Output: 1823 // The biased variance of the samples is 62.0000 1824 // The weighted biased variance of the samples is 105.5833 1825} 1826 1827func TestStdScore(t *testing.T) { 1828 for i, test := range []struct { 1829 x float64 1830 u float64 1831 s float64 1832 z float64 1833 }{ 1834 { 1835 x: 4, 1836 u: -6, 1837 s: 5, 1838 z: 2, 1839 }, 1840 { 1841 x: 1, 1842 u: 0, 1843 s: 1, 1844 z: 1, 1845 }, 1846 } { 1847 z := StdScore(test.x, test.u, test.s) 1848 if math.Abs(z-test.z) > 1e-14 { 1849 t.Errorf("StdScore mismatch case %d. Expected %v, Found %v", i, test.z, z) 1850 } 1851 } 1852} 1853