diff --git a/analysis/app/app.go b/analysis/app/app.go
index 5ab9459..c4d3029 100644
--- a/analysis/app/app.go
+++ b/analysis/app/app.go
@@ -23,6 +23,7 @@
 	mux.HandleFunc("/", a.index)
 	mux.HandleFunc("/search", a.search)
 	mux.HandleFunc("/compare", a.compare)
+	mux.HandleFunc("/trend", a.trend)
 }
 
 // search handles /search.
diff --git a/analysis/app/kza.go b/analysis/app/kza.go
new file mode 100644
index 0000000..9ad6fa6
--- /dev/null
+++ b/analysis/app/kza.go
@@ -0,0 +1,125 @@
+// Copyright 2015 The Go 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 app
+
+import "math"
+
+// TODO: This all assumes that data is sampled at a regular interval
+// and there are no missing values. It could be generalized to accept
+// missing values (perhaps represented by NaN), or generalized much
+// further by accepting (t, x) pairs and a vector of times at which to
+// evaluate the filter (and an arbitrary window size). I would have to
+// figure out how that affects the difference array in KZA.
+
+// TODO: These can generate a lot of garbage. Perhaps the caller
+// should pass in the target slice? Or these should just overwrite the
+// input array and leave it to the caller to copy if necessary?
+
+// MovingAverage performs a moving average (MA) filter of xs with
+// window size m. m must be a positive odd integer.
+//
+// Note that this is filter is often described in terms of the half
+// length of the window (m-1)/2.
+func MovingAverage(xs []float64, m int) []float64 {
+	if m <= 0 || m%2 != 1 {
+		panic("m must be a positive, odd integer")
+	}
+	ys := make([]float64, len(xs))
+	sum, n := 0.0, 0
+	for l, i, r := -m, -(m-1)/2, 0; i < len(ys); l, i, r = l+1, i+1, r+1 {
+		if l >= 0 {
+			sum -= xs[l]
+			n--
+		}
+		if r < len(xs) {
+			sum += xs[r]
+			n++
+		}
+		if i >= 0 {
+			ys[i] = sum / float64(n)
+		}
+	}
+	return ys
+}
+
+// KolmogorovZurbenko performs a Kolmogorov-Zurbenko (KZ) filter of xs
+// with window size m and k iterations. m must be a positive odd
+// integer. k must be positive.
+func KolmogorovZurbenko(xs []float64, m, k int) []float64 {
+	// k is typically small, and MA is quite efficient, so just do
+	// the iterated moving average rather than bothering to
+	// compute the binomial coefficient kernel.
+	for i := 0; i < k; i++ {
+		// TODO: Generate less garbage.
+		xs = MovingAverage(xs, m)
+	}
+	return xs
+}
+
+// AdaptiveKolmogorovZurbenko performs an adaptive Kolmogorov-Zurbenko
+// (KZA) filter of xs using an initial window size m and k iterations.
+// m must be a positive odd integer. k must be positive.
+//
+// See Zurbenko, et al. 1996: Detecting discontinuities in time series
+// of upper air data: Demonstration of an adaptive filter technique.
+// Journal of Climate, 9, 3548–3560.
+func AdaptiveKolmogorovZurbenko(xs []float64, m, k int) []float64 {
+	// Perform initial KZ filter.
+	z := KolmogorovZurbenko(xs, m, k)
+
+	// Compute differenced values.
+	q := (m - 1) / 2
+	d := make([]float64, len(z)+1)
+	maxD := 0.0
+	for i := q; i < len(z)-q; i++ {
+		d[i] = math.Abs(z[i+q] - z[i-q])
+		if d[i] > maxD {
+			maxD = d[i]
+		}
+	}
+
+	if maxD == 0 {
+		// xs is constant, so no amount of filtering will do
+		// anything. Avoid dividing 0/0 below.
+		return xs
+	}
+
+	// Compute adaptive filter.
+	ys := make([]float64, len(xs))
+	for t := range ys {
+		dPrime := d[t+1] - d[t]
+		f := 1 - d[t]/maxD
+
+		qt := q
+		if dPrime <= 0 {
+			// Zurbenko doesn't specify what to do with
+			// the fractional part of qt and qh, so we
+			// interpret this as summing all points of xs
+			// between qt and qh.
+			qt = int(math.Ceil(float64(q) * f))
+		}
+		if t-qt < 0 {
+			qt = t
+		}
+
+		qh := q
+		if dPrime >= 0 {
+			qh = int(math.Floor(float64(q) * f))
+		}
+		if t+qh >= len(xs) {
+			qh = len(xs) - t - 1
+		}
+
+		sum := 0.0
+		for i := t - qt; i <= t+qh; i++ {
+			sum += xs[i]
+		}
+		// Zurbenko divides by qh+qt, but this undercounts the
+		// number of terms in the sum by 1.
+		ys[t] = sum / float64(qh+qt+1)
+	}
+
+	return ys
+}
diff --git a/analysis/app/kza_test.go b/analysis/app/kza_test.go
new file mode 100644
index 0000000..df65223
--- /dev/null
+++ b/analysis/app/kza_test.go
@@ -0,0 +1,54 @@
+// Copyright 2015 The Go 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 app
+
+import (
+	"math/rand"
+	"testing"
+)
+
+// Aeq returns true if expect and got are equal to 8 significant
+// figures (1 part in 100 million).
+func Aeq(expect, got float64) bool {
+	if expect < 0 && got < 0 {
+		expect, got = -expect, -got
+	}
+	return expect*0.99999999 <= got && got*0.99999999 <= expect
+}
+
+func TestMovingAverage(t *testing.T) {
+	// Test MovingAverage against the obvious (but slow)
+	// implementation.
+	xs := make([]float64, 100)
+	for iter := 0; iter < 10; iter++ {
+		for i := range xs {
+			xs[i] = rand.Float64()
+		}
+		m := 1 + 2*rand.Intn(100)
+		ys1, ys2 := MovingAverage(xs, m), slowMovingAverage(xs, m)
+
+		// TODO: Use stuff from mathtest.
+		for i, y1 := range ys1 {
+			if !Aeq(y1, ys2[i]) {
+				t.Fatalf("want %v, got %v", ys2, ys1)
+			}
+		}
+	}
+}
+
+func slowMovingAverage(xs []float64, m int) []float64 {
+	ys := make([]float64, len(xs))
+	for i := range ys {
+		psum, n := 0.0, 0
+		for j := i - (m-1)/2; j <= i+(m-1)/2; j++ {
+			if 0 <= j && j < len(xs) {
+				psum += xs[j]
+				n++
+			}
+		}
+		ys[i] = psum / float64(n)
+	}
+	return ys
+}
diff --git a/analysis/app/trend.go b/analysis/app/trend.go
new file mode 100644
index 0000000..cf4837f
--- /dev/null
+++ b/analysis/app/trend.go
@@ -0,0 +1,410 @@
+// Copyright 2017 The Go Authors.  All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Loosely based on github.com/aclements/go-misc/benchplot
+
+package app
+
+import (
+	"bytes"
+	"encoding/json"
+	"fmt"
+	"html/template"
+	"io/ioutil"
+	"math"
+	"net/http"
+	"sort"
+	"strconv"
+	"strings"
+
+	"github.com/aclements/go-gg/generic/slice"
+	"github.com/aclements/go-gg/ggstat"
+	"github.com/aclements/go-gg/table"
+	"golang.org/x/net/context"
+	"golang.org/x/perf/storage"
+)
+
+// trend handles /trend.
+// With no query, it prints the list of recent uploads containg a "trend" key.
+// With a query, it shows a graph of the matching benchmark results.
+func (a *App) trend(w http.ResponseWriter, r *http.Request) {
+	ctx := requestContext(r)
+
+	if err := r.ParseForm(); err != nil {
+		http.Error(w, err.Error(), 500)
+		return
+	}
+
+	q := r.Form.Get("q")
+
+	tmpl, err := ioutil.ReadFile("template/trend.html")
+	if err != nil {
+		http.Error(w, err.Error(), 500)
+		return
+	}
+
+	t, err := template.New("main").Parse(string(tmpl))
+	if err != nil {
+		http.Error(w, err.Error(), 500)
+		return
+	}
+
+	data := a.trendQuery(ctx, q)
+
+	w.Header().Set("Content-Type", "text/html; charset=utf-8")
+	if err := t.Execute(w, data); err != nil {
+		http.Error(w, err.Error(), 500)
+		return
+	}
+}
+
+// trendData is the struct passed to the trend.html template.
+type trendData struct {
+	Q            string
+	Error        string
+	TrendUploads []storage.UploadInfo
+	PlotData     template.JS
+}
+
+// trendData computes the values for the template and returns a trendData for display.
+func (a *App) trendQuery(ctx context.Context, q string) *trendData {
+	d := &trendData{Q: q}
+	if q == "" {
+		ul := a.StorageClient.ListUploads(`trend>`, []string{"by", "upload-time", "trend"}, 16)
+		defer ul.Close()
+		for ul.Next() {
+			d.TrendUploads = append(d.TrendUploads, ul.Info())
+		}
+		if err := ul.Err(); err != nil {
+			errorf(ctx, "failed to fetch recent trend uploads: %v", err)
+		}
+		return d
+	}
+
+	// TODO(quentin): Chunk query based on matching upload IDs.
+	res := a.StorageClient.Query(q)
+	defer res.Close()
+	t, resultCols := queryToTable(res)
+	for _, col := range []string{"commit", "commit-time", "branch"} {
+		if t.Column(col) == nil {
+			d.Error = fmt.Sprintf("results missing %q label", col)
+			return d
+		}
+	}
+	data := plot(t, resultCols)
+
+	// TODO(quentin): Give the user control over across vs. plotting in separate graphs, instead of only showing one graph with ns/op for each benchmark.
+
+	// Pivot all of the benchmarks into columns of a single table.
+	ar := &aggResults{
+		Across: "name",
+		Values: []string{"filtered normalized mean result", "normalized mean result", "normalized min result", "normalized max result"},
+	}
+	data = ggstat.Agg("commit", "branch", "commit-index")(ar.agg).F(data)
+
+	tables := data.Tables()
+	infof(ctx, "tables: %v", tables)
+	columns := []column{{Name: "commit-index"}}
+	for _, prefix := range ar.Prefixes {
+		if len(ar.Prefixes) == 1 {
+			columns = append(columns,
+				column{Name: prefix + "/normalized mean result"},
+				column{Name: prefix + "/normalized min result", Role: "interval"},
+				column{Name: prefix + "/normalized max result", Role: "interval"},
+			)
+		}
+		columns = append(columns,
+			column{Name: prefix + "/filtered normalized mean result"},
+		)
+	}
+	d.PlotData = tableToJS(data.Table(tables[0]), columns)
+	return d
+}
+
+// queryToTable converts the result of a Query into a Table for later processing.
+// Each label is placed in a column named after the key.
+// Each metric is placed in a separate result column named after the unit.
+func queryToTable(q *storage.Query) (t *table.Table, resultCols []string) {
+	var names []string
+	labels := make(map[string][]string)
+	results := make(map[string][]float64)
+	i := 0
+	for q.Next() {
+		res := q.Result()
+		// TODO(quentin): Handle multiple results with the same name but different NameLabels.
+		names = append(names, res.NameLabels["name"])
+		for k := range res.Labels {
+			if labels[k] == nil {
+				labels[k] = make([]string, i)
+			}
+		}
+		for k := range labels {
+			labels[k] = append(labels[k], res.Labels[k])
+		}
+		f := strings.Fields(res.Content)
+		metrics := make(map[string]float64)
+		for j := 2; j+2 <= len(f); j += 2 {
+			val, err := strconv.ParseFloat(f[j], 64)
+			if err != nil {
+				continue
+			}
+			unit := f[j+1]
+			if results[unit] == nil {
+				results[unit] = make([]float64, i)
+			}
+			metrics[unit] = val
+		}
+		for k := range results {
+			results[k] = append(results[k], metrics[k])
+		}
+		i++
+	}
+
+	tab := new(table.Builder).Add("name", names)
+
+	for k, v := range labels {
+		tab.Add(k, v)
+	}
+	for k, v := range results {
+		tab.Add(k, v)
+		resultCols = append(resultCols, k)
+	}
+
+	sort.Strings(resultCols)
+
+	return tab.Done(), resultCols
+}
+
+// plot takes raw benchmark data in t and produces a Grouping object containing filtered, normalized metric results for a graph.
+// t must contain the string columns "commit", "commit-time", "branch". resultCols specifies the names of float64 columns containing metric results.
+// The returned grouping has columns "commit", "commit-time", "commit-index", "branch", "metric", "normalized min result", "normalized max result", "normalized mean result", "filtered normalized mean result".
+// This is roughly the algorithm from github.com/aclements/go-misc/benchplot
+func plot(t table.Grouping, resultCols []string) table.Grouping {
+	nrows := len(table.GroupBy(t, "name").Tables())
+
+	// Turn ordered commit-time into a "commit-index" column.
+	t = table.SortBy(t, "commit-time", "commit")
+	t = commitIndex{}.F(t)
+
+	// Unpivot all of the metrics into one column.
+	t = table.Unpivot(t, "metric", "result", resultCols...)
+
+	// TODO(quentin): Let user choose which metric(s) to keep.
+	t = table.FilterEq(t, "metric", "ns/op")
+
+	// Average each result at each commit (but keep columns names
+	// the same to keep things easier to read).
+	t = ggstat.Agg("commit", "name", "metric", "branch", "commit-index")(ggstat.AggMean("result"), ggstat.AggMin("result"), ggstat.AggMax("result")).F(t)
+	y := "mean result"
+
+	// Normalize to earliest commit on master. It's important to
+	// do this before the geomean if there are commits missing.
+	// Unfortunately, that also means we have to *temporarily*
+	// group by name and metric, since the geomean needs to be
+	// done on a different grouping.
+	t = table.GroupBy(t, "name", "metric")
+	t = ggstat.Normalize{X: "branch", By: firstMasterIndex, Cols: []string{"mean result", "max result", "min result"}, DenomCols: []string{"mean result", "mean result", "mean result"}}.F(t)
+	y = "normalized " + y
+	for _, col := range []string{"mean result", "max result", "min result"} {
+		t = table.Remove(t, col)
+	}
+	t = table.Ungroup(table.Ungroup(t))
+
+	// Compute geomean for each metric at each commit if there's
+	// more than one benchmark.
+	if len(table.GroupBy(t, "name").Tables()) > 1 {
+		gt := removeNaNs(t, y)
+		gt = ggstat.Agg("commit", "metric", "branch", "commit-index")(ggstat.AggGeoMean(y), ggstat.AggMin("normalized min result"), ggstat.AggMax("normalized max result")).F(gt)
+		gt = table.MapTables(gt, func(_ table.GroupID, t *table.Table) *table.Table {
+			return table.NewBuilder(t).AddConst("name", " geomean").Done()
+		})
+		gt = table.Rename(gt, "geomean "+y, y)
+		gt = table.Rename(gt, "min normalized min result", "normalized min result")
+		gt = table.Rename(gt, "max normalized max result", "normalized max result")
+		t = table.Concat(t, gt)
+		nrows++
+	}
+
+	// Filter the data to reduce noise.
+	t = table.GroupBy(t, "name", "metric")
+	t = kza{y, 15, 3}.F(t)
+	y = "filtered " + y
+	t = table.Ungroup(table.Ungroup(t))
+
+	return t
+}
+
+// aggResults pivots the table, taking the columns in Values and making a new column for each distinct value in Across.
+// aggResults("in", []string{"value1", "value2"} will reshape a table like
+//   in		value1	value2
+//   one	1	2
+//   two	3	4
+// and will turn in into a table like
+//   one/value1	one/value2	two/value1	two/value2
+//   1		2		3		4
+// across columns must be []string, and value columns must be []float64.
+type aggResults struct {
+	// Across is the name of the column whose values are the column prefix.
+	Across string
+	// Values is the name of the columns to split.
+	Values []string
+	// Prefixes is filled in after calling agg with the name of each prefix that was found.
+	Prefixes []string
+}
+
+// agg implements ggstat.Aggregator and allows using a with ggstat.Agg.
+func (a *aggResults) agg(input table.Grouping, output *table.Builder) {
+	var prefixes []string
+	rows := len(input.Tables())
+	columns := make(map[string][]float64)
+	for i, gid := range input.Tables() {
+		var vs [][]float64
+		for _, col := range a.Values {
+			vs = append(vs, input.Table(gid).MustColumn(col).([]float64))
+		}
+		as := input.Table(gid).MustColumn(a.Across).([]string)
+		for j, prefix := range as {
+			for k, col := range a.Values {
+				key := prefix + "/" + col
+				if columns[key] == nil {
+					if k == 0 {
+						// First time we've seen this prefix, track it.
+						prefixes = append(prefixes, prefix)
+					}
+					columns[key] = make([]float64, rows)
+					for i := range columns[key] {
+						columns[key][i] = math.NaN()
+					}
+				}
+				columns[key][i] = vs[k][j]
+			}
+		}
+	}
+	sort.Strings(prefixes)
+	a.Prefixes = prefixes
+	for _, prefix := range prefixes {
+		for _, col := range a.Values {
+			key := prefix + "/" + col
+			output.Add(key, columns[key])
+		}
+	}
+}
+
+// firstMasterIndex returns the index of the first commit on master.
+// This is used to find the value to normalize against.
+func firstMasterIndex(bs []string) int {
+	return slice.Index(bs, "master")
+}
+
+// commitIndex is a gg.Stat that adds a column called "commit-index" sequentially counting unique values of the column "commit".
+type commitIndex struct{}
+
+func (commitIndex) F(g table.Grouping) table.Grouping {
+	return table.MapTables(g, func(_ table.GroupID, t *table.Table) *table.Table {
+		idxs := make([]int, t.Len())
+		last, idx := "", -1
+		for i, hash := range t.MustColumn("commit").([]string) {
+			if hash != last {
+				idx++
+				last = hash
+			}
+			idxs[i] = idx
+		}
+		t = table.NewBuilder(t).Add("commit-index", idxs).Done()
+
+		return t
+	})
+}
+
+// removeNaNs returns a new Grouping with rows containg NaN in col removed.
+func removeNaNs(g table.Grouping, col string) table.Grouping {
+	return table.Filter(g, func(result float64) bool {
+		return !math.IsNaN(result)
+	}, col)
+}
+
+// kza implements adaptive Kolmogorov-Zurbenko filtering on the data in X.
+type kza struct {
+	X    string
+	M, K int
+}
+
+func (k kza) F(g table.Grouping) table.Grouping {
+	return table.MapTables(g, func(_ table.GroupID, t *table.Table) *table.Table {
+		var xs []float64
+		slice.Convert(&xs, t.MustColumn(k.X))
+		nxs := AdaptiveKolmogorovZurbenko(xs, k.M, k.K)
+		return table.NewBuilder(t).Add("filtered "+k.X, nxs).Done()
+	})
+}
+
+// column represents a column in a google.visualization.DataTable
+type column struct {
+	Name string `json:"id"`
+	Role string `json:"role,omitempty"`
+	// These fields are filled in by tableToJS if unspecified.
+	Type  string `json:"type"`
+	Label string `json:"label"`
+}
+
+// tableToJS converts a Table to a javascript literal which can be passed to "new google.visualization.DataTable".
+func tableToJS(t *table.Table, columns []column) template.JS {
+	var out bytes.Buffer
+	fmt.Fprint(&out, "{cols: [")
+	var slices []table.Slice
+	for i, c := range columns {
+		if i > 0 {
+			fmt.Fprint(&out, ",\n")
+		}
+		col := t.Column(c.Name)
+		slices = append(slices, col)
+		if c.Type == "" {
+			switch col.(type) {
+			case []string:
+				c.Type = "string"
+			case []int, []float64:
+				c.Type = "number"
+			}
+		}
+		if c.Label == "" {
+			c.Label = c.Name
+		}
+		data, err := json.Marshal(c)
+		if err != nil {
+			panic(err)
+		}
+		out.Write(data)
+	}
+	fmt.Fprint(&out, "],\nrows: [")
+	for i := 0; i < t.Len(); i++ {
+		if i > 0 {
+			fmt.Fprint(&out, ",\n")
+		}
+		fmt.Fprint(&out, "{c:[")
+		for j := range columns {
+			if j > 0 {
+				fmt.Fprint(&out, ", ")
+			}
+			fmt.Fprint(&out, "{v: ")
+			var value []byte
+			var err error
+			switch column := slices[j].(type) {
+			case []string:
+				value, err = json.Marshal(column[i])
+			case []int:
+				value, err = json.Marshal(column[i])
+			case []float64:
+				value, err = json.Marshal(column[i])
+			}
+			if err != nil {
+				panic(err)
+			}
+			out.Write(value)
+			fmt.Fprint(&out, "}")
+		}
+		fmt.Fprint(&out, "]}")
+	}
+	fmt.Fprint(&out, "]}")
+	return template.JS(out.String())
+}
diff --git a/analysis/app/trend_test.go b/analysis/app/trend_test.go
new file mode 100644
index 0000000..8bf2157
--- /dev/null
+++ b/analysis/app/trend_test.go
@@ -0,0 +1,29 @@
+// Copyright 2017 The Go 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 app
+
+import (
+	"testing"
+
+	"github.com/aclements/go-gg/table"
+	"golang.org/x/perf/internal/diff"
+)
+
+func TestTableToJS(t *testing.T) {
+	in := table.TableFromStrings(
+		[]string{"text", "num"},
+		[][]string{
+			{"hello", "15.1"},
+			{"world", "20"},
+		}, true)
+	have := tableToJS(in, []column{{Name: "text"}, {Name: "num"}})
+	want := `{cols: [{"id":"text","type":"string","label":"text"},
+{"id":"num","type":"number","label":"num"}],
+rows: [{c:[{v: "hello"}, {v: 15.1}]},
+{c:[{v: "world"}, {v: 20}]}]}`
+	if d := diff.Diff(string(have), want); d != "" {
+		t.Errorf("tableToJS returned wrong JS (- have, + want):\n%s", d)
+	}
+}
diff --git a/analysis/appengine/template/trend.html b/analysis/appengine/template/trend.html
new file mode 100644
index 0000000..cd1d8eb
--- /dev/null
+++ b/analysis/appengine/template/trend.html
@@ -0,0 +1,73 @@
+<!DOCTYPE html>
+<html>
+  <head>
+    <title>Performance Result Comparison</title>
+    <script type="text/javascript" src="https://www.gstatic.com/charts/loader.js"></script>
+    <style type="text/css">
+#header h1 {
+  display: inline;
+}
+#search {
+  padding: 1em .5em;
+  width: 100%;
+}
+input[type="text"] {
+  font-size: 100%;
+}
+#results {
+  border-top: 1px solid black;
+}
+    </style>
+  </head>
+  <body>
+    <div id="header">
+      <h1>Go Performance Dashboard</h1>
+      <a href="/">about</a>
+    </div>
+    <div id="search">
+      <form action="/trend">
+        <input type="text" name="q" value="{{.Q}}" size="120">
+        <input type="submit" value="Search">
+      </form>
+    </div>
+    <div id="results">
+      {{if not .Q}}
+        <h2>Recent Uploads</h2>
+        <table>
+          <tr><th>Upload ID</th><th>trend</th></tr>
+          {{range .TrendUploads}}
+          <tr><td><a href="/trend?q=upload:{{.UploadID}}">{{.UploadID}}</a></td><td>{{.LabelValues.trend}}</td></tr>
+          {{end}}
+        </table>
+      {{else}}
+        {{with .Error}}
+        <p>{{.}}</p>
+        {{else}}
+          <div id="chart" style="height: 600px"></div>
+          <script type="text/javascript">
+            google.charts.load('current', {'packages':['corechart']});
+            google.charts.setOnLoadCallback(draw);
+            function draw() {
+              var dt = new google.visualization.DataTable({{.PlotData}});
+              var options = {
+                title: 'Benchmark Trend',
+                hAxis: {
+                  title: 'commit index',
+                },
+                vAxis: {
+                  title: 'normalized ns/op',
+                },
+                explorer: {
+                  actions: ['dragToZoom', 'rightClickToReset'],
+                  maxZoomIn: 0.05,
+                },
+              };
+              var chart = new google.visualization.LineChart(document.getElementById('chart'));
+              chart.draw(dt, options);
+            }
+          </script>
+        {{end}}
+      {{end}}
+    </div>
+  </body>
+</html>
diff --git a/vendor/github.com/aclements/go-gg/LICENSE b/vendor/github.com/aclements/go-gg/LICENSE
new file mode 100644
index 0000000..a5389da
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/LICENSE
@@ -0,0 +1,28 @@
+Copyright (c) 2016 The Go Authors. All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+   * Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+   * Redistributions in binary form must reproduce the above
+copyright notice, this list of conditions and the following disclaimer
+in the documentation and/or other materials provided with the
+distribution.
+   * Neither the name of Google Inc. nor the names of its
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
diff --git a/vendor/github.com/aclements/go-gg/generic/doc.go b/vendor/github.com/aclements/go-gg/generic/doc.go
new file mode 100644
index 0000000..0df3e8a
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/doc.go
@@ -0,0 +1,6 @@
+// Copyright 2016 The Go 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 generic provides type-generic functions.
+package generic
diff --git a/vendor/github.com/aclements/go-gg/generic/error.go b/vendor/github.com/aclements/go-gg/generic/error.go
new file mode 100644
index 0000000..a04e5e5
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/error.go
@@ -0,0 +1,21 @@
+// Copyright 2016 The Go 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 generic
+
+import "reflect"
+
+type TypeError struct {
+	Type1, Type2 reflect.Type
+	Extra        string
+}
+
+func (e TypeError) Error() string {
+	msg := e.Type1.String()
+	if e.Type2 != nil {
+		msg += " and " + e.Type2.String()
+	}
+	msg += " " + e.Extra
+	return msg
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/order.go b/vendor/github.com/aclements/go-gg/generic/order.go
new file mode 100644
index 0000000..551373d
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/order.go
@@ -0,0 +1,96 @@
+// Copyright 2016 The Go 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 generic
+
+import "reflect"
+
+// CanOrder returns whether the values a and b are orderable according
+// to the Go language specification.
+func CanOrder(a, b interface{}) bool {
+	ak, bk := reflect.ValueOf(a).Kind(), reflect.ValueOf(b).Kind()
+	if ak != bk {
+		return false
+	}
+	return CanOrderR(ak)
+}
+
+var orderable = map[reflect.Kind]bool{
+	reflect.Int:     true,
+	reflect.Int8:    true,
+	reflect.Int16:   true,
+	reflect.Int32:   true,
+	reflect.Int64:   true,
+	reflect.Uint:    true,
+	reflect.Uintptr: true,
+	reflect.Uint8:   true,
+	reflect.Uint16:  true,
+	reflect.Uint32:  true,
+	reflect.Uint64:  true,
+	reflect.Float32: true,
+	reflect.Float64: true,
+	reflect.String:  true,
+}
+
+// CanOrderR returns whether two values of kind k are orderable
+// according to the Go language specification.
+func CanOrderR(k reflect.Kind) bool {
+	return orderable[k]
+}
+
+// Order returns the order of values a and b: -1 if a < b, 0 if a ==
+// b, 1 if a > b. The results are undefined if either a or b is NaN.
+//
+// Order panics if a and b are not orderable according to the Go
+// language specification.
+func Order(a, b interface{}) int {
+	return OrderR(reflect.ValueOf(a), reflect.ValueOf(b))
+}
+
+// OrderR is equivalent to Order, but operates on reflect.Values.
+func OrderR(a, b reflect.Value) int {
+	if a.Kind() != b.Kind() {
+		panic(&TypeError{a.Type(), b.Type(), "are not orderable because they are different kinds"})
+	}
+
+	switch a.Kind() {
+	case reflect.Float32, reflect.Float64:
+		a, b := a.Float(), b.Float()
+		if a < b {
+			return -1
+		} else if a > b {
+			return 1
+		}
+		return 0
+
+	case reflect.Int, reflect.Int8, reflect.Int16, reflect.Int32, reflect.Int64:
+		a, b := a.Int(), b.Int()
+		if a < b {
+			return -1
+		} else if a > b {
+			return 1
+		}
+		return 0
+
+	case reflect.Uint, reflect.Uintptr, reflect.Uint8, reflect.Uint16, reflect.Uint32, reflect.Uint64:
+		a, b := a.Uint(), b.Uint()
+		if a < b {
+			return -1
+		} else if a > b {
+			return 1
+		}
+		return 0
+
+	case reflect.String:
+		a, b := a.String(), b.String()
+		if a < b {
+			return -1
+		} else if a > b {
+			return 1
+		}
+		return 0
+	}
+
+	panic(&TypeError{a.Type(), nil, "is not orderable"})
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/concat.go b/vendor/github.com/aclements/go-gg/generic/slice/concat.go
new file mode 100644
index 0000000..e538790
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/concat.go
@@ -0,0 +1,41 @@
+// Copyright 2016 The Go 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 slice
+
+import (
+	"reflect"
+
+	"github.com/aclements/go-gg/generic"
+)
+
+// Concat returns the concatenation of all of ss. The types of all of
+// the arguments must be identical or Concat will panic with a
+// *generic.TypeError. The returned slice will have the same type as the
+// inputs. If there are 0 arguments, Concat returns nil. Concat does
+// not modify any of the input slices.
+func Concat(ss ...T) T {
+	if len(ss) == 0 {
+		return nil
+	}
+
+	rvs := make([]reflect.Value, len(ss))
+	total := 0
+	var typ reflect.Type
+	for i, s := range ss {
+		rvs[i] = reflectSlice(s)
+		total += rvs[i].Len()
+		if i == 0 {
+			typ = rvs[i].Type()
+		} else if rvs[i].Type() != typ {
+			panic(&generic.TypeError{typ, rvs[i].Type(), "have different types"})
+		}
+	}
+
+	out := reflect.MakeSlice(typ, 0, total)
+	for _, rv := range rvs {
+		out = reflect.AppendSlice(out, rv)
+	}
+	return out.Interface()
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/convert.go b/vendor/github.com/aclements/go-gg/generic/slice/convert.go
new file mode 100644
index 0000000..c3d46a3
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/convert.go
@@ -0,0 +1,55 @@
+// Copyright 2016 The Go 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 slice
+
+import (
+	"reflect"
+
+	"github.com/aclements/go-gg/generic"
+)
+
+// Convert converts each element in from and assigns it to *to. to
+// must be a pointer to a slice. Convert slices or extends *to to
+// len(from) and then assigns to[i] = T(from[i]) where T is the type
+// of *to's elements. If from and *to have the same element type, it
+// simply assigns *to = from.
+func Convert(to interface{}, from T) {
+	fv := reflectSlice(from)
+	tv := reflect.ValueOf(to)
+	if tv.Kind() != reflect.Ptr {
+		panic(&generic.TypeError{tv.Type(), nil, "is not a *[]T"})
+	}
+	tst := tv.Type().Elem()
+	if tst.Kind() != reflect.Slice {
+		panic(&generic.TypeError{tv.Type(), nil, "is not a *[]T"})
+	}
+
+	if fv.Type().AssignableTo(tst) {
+		tv.Elem().Set(fv)
+		return
+	}
+
+	eltt := tst.Elem()
+	if !fv.Type().Elem().ConvertibleTo(eltt) {
+		panic(&generic.TypeError{fv.Type(), tst, "cannot be converted"})
+	}
+
+	switch to := to.(type) {
+	case *[]float64:
+		// This is extremely common.
+		*to = (*to)[:0]
+		for i, len := 0, fv.Len(); i < len; i++ {
+			*to = append(*to, fv.Index(i).Convert(eltt).Float())
+		}
+
+	default:
+		tsv := tv.Elem()
+		tsv.SetLen(0)
+		for i, len := 0, fv.Len(); i < len; i++ {
+			tsv = reflect.Append(tsv, fv.Index(i).Convert(eltt))
+		}
+		tv.Elem().Set(tsv)
+	}
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/cycle.go b/vendor/github.com/aclements/go-gg/generic/slice/cycle.go
new file mode 100644
index 0000000..a94e5e2
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/cycle.go
@@ -0,0 +1,45 @@
+// Copyright 2016 The Go 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 slice
+
+import "reflect"
+
+// Cycle constructs a slice of length length by repeatedly
+// concatenating s to itself. If len(s) >= length, it returns
+// s[:length]. Otherwise, it allocates a new slice. If len(s) == 0 and
+// length != 0, Cycle panics.
+func Cycle(s T, length int) T {
+	rv := reflectSlice(s)
+	if rv.Len() >= length {
+		return rv.Slice(0, length).Interface()
+	}
+
+	if rv.Len() == 0 {
+		panic("empty slice")
+	}
+
+	// Allocate a new slice of the appropriate length.
+	out := reflect.MakeSlice(rv.Type(), length, length)
+
+	// Copy elements to out.
+	for pos := 0; pos < length; {
+		pos += reflect.Copy(out.Slice(pos, length), rv)
+	}
+
+	return out.Interface()
+}
+
+// Repeat returns a slice consisting of length copies of v.
+func Repeat(v interface{}, length int) T {
+	if length < 0 {
+		length = 0
+	}
+	rv := reflect.ValueOf(v)
+	out := reflect.MakeSlice(reflect.SliceOf(rv.Type()), length, length)
+	for i := 0; i < length; i++ {
+		out.Index(i).Set(rv)
+	}
+	return out.Interface()
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/doc.go b/vendor/github.com/aclements/go-gg/generic/slice/doc.go
new file mode 100644
index 0000000..2ebb862
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/doc.go
@@ -0,0 +1,6 @@
+// Copyright 2016 The Go 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 slice provides generic slice functions.
+package slice
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/find.go b/vendor/github.com/aclements/go-gg/generic/slice/find.go
new file mode 100644
index 0000000..09c7d53
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/find.go
@@ -0,0 +1,51 @@
+// Copyright 2016 The Go 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 slice
+
+import (
+	"reflect"
+
+	"github.com/aclements/go-gg/generic"
+)
+
+// Index returns the index of the first instance of val in s, or -1 if
+// val is not present in s. val's type must be s's element type.
+func Index(s T, val interface{}) int {
+	rs := reflectSlice(s)
+	if vt := reflect.TypeOf(val); rs.Type().Elem() != vt {
+		// TODO: Better "<seq> is not a sequence of <val>".
+		panic(&generic.TypeError{rs.Type(), vt, "cannot find"})
+	}
+
+	for i, l := 0, rs.Len(); i < l; i++ {
+		if rs.Index(i).Interface() == val {
+			return i
+		}
+	}
+	return -1
+}
+
+// LastIndex returns the index of the last instance of val in s, or -1
+// if val is not present in s. val's type must be s's element type.
+func LastIndex(s T, val interface{}) int {
+	rs := reflectSlice(s)
+	if vt := reflect.TypeOf(val); rs.Type().Elem() != vt {
+		// TODO: Better "<seq> is not a sequence of <val>".
+		panic(&generic.TypeError{rs.Type(), vt, "cannot find"})
+	}
+
+	for i := rs.Len() - 1; i >= 0; i-- {
+		if rs.Index(i).Interface() == val {
+			return i
+		}
+	}
+	return -1
+}
+
+// Contains reports whether val is within s. val's type must be s's
+// element type.
+func Contains(s T, val interface{}) bool {
+	return Index(s, val) >= 0
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/index.go b/vendor/github.com/aclements/go-gg/generic/slice/index.go
new file mode 100644
index 0000000..780b4ac
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/index.go
@@ -0,0 +1,82 @@
+// Copyright 2016 The Go 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 slice
+
+import (
+	"reflect"
+
+	"github.com/aclements/go-gg/generic"
+)
+
+// Select returns a slice w such that w[i] = v[indexes[i]].
+func Select(v T, indexes []int) T {
+	switch v := v.(type) {
+	case []int:
+		res := make([]int, len(indexes))
+		for i, x := range indexes {
+			res[i] = v[x]
+		}
+		return res
+
+	case []float64:
+		res := make([]float64, len(indexes))
+		for i, x := range indexes {
+			res[i] = v[x]
+		}
+		return res
+
+	case []string:
+		res := make([]string, len(indexes))
+		for i, x := range indexes {
+			res[i] = v[x]
+		}
+		return res
+	}
+
+	rv := reflectSlice(v)
+	res := reflect.MakeSlice(rv.Type(), len(indexes), len(indexes))
+	for i, x := range indexes {
+		res.Index(i).Set(rv.Index(x))
+	}
+	return res.Interface()
+}
+
+// SelectInto assigns out[i] = in[indexes[i]]. in and out must have
+// the same types and len(out) must be >= len(indexes). If in and out
+// overlap, the results are undefined.
+func SelectInto(out, in T, indexes []int) {
+	// TODO: Maybe they should only have to be assignable?
+	if it, ot := reflect.TypeOf(in), reflect.TypeOf(out); it != ot {
+		panic(&generic.TypeError{it, ot, "must be the same type"})
+	}
+
+	switch in := in.(type) {
+	case []int:
+		out := out.([]int)
+		for i, x := range indexes {
+			out[i] = in[x]
+		}
+		return
+
+	case []float64:
+		out := out.([]float64)
+		for i, x := range indexes {
+			out[i] = in[x]
+		}
+		return
+
+	case []string:
+		out := out.([]string)
+		for i, x := range indexes {
+			out[i] = in[x]
+		}
+		return
+	}
+
+	inv, outv := reflectSlice(in), reflectSlice(out)
+	for i, x := range indexes {
+		outv.Index(i).Set(inv.Index(x))
+	}
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/min.go b/vendor/github.com/aclements/go-gg/generic/slice/min.go
new file mode 100644
index 0000000..7cda3fc
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/min.go
@@ -0,0 +1,99 @@
+// Copyright 2016 The Go 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 slice
+
+import (
+	"reflect"
+	"sort"
+
+	"github.com/aclements/go-gg/generic"
+)
+
+// Min returns the minimum value in v. v must either implement
+// sort.Interface or its elements must be orderable. Min panics if v
+// is empty.
+func Min(v T) interface{} {
+	x, _ := minmax(v, -1, true)
+	return x.Interface()
+}
+
+// ArgMin returns the index of the minimum value in v. If there are
+// multiple indexes equal to the minimum value, ArgMin returns the
+// lowest of them. v must be a slice whose elements are orderable, or
+// must implement sort.Interface. ArgMin panics if v is empty.
+func ArgMin(v interface{}) int {
+	_, i := minmax(v, -1, false)
+	return i
+}
+
+// Max returns the maximum value in v. v must either implement
+// sort.Interface or its elements must be orderable. Max panics if v
+// is empty.
+func Max(v T) interface{} {
+	x, _ := minmax(v, 1, true)
+	return x.Interface()
+}
+
+// ArgMax returns the index of the maximum value in v. If there are
+// multiple indexes equal to the maximum value, ArgMax returns the
+// lowest of them. v must be a slice whose elements are orderable, or
+// must implement sort.Interface. ArgMax panics if v is empty.
+func ArgMax(v interface{}) int {
+	_, i := minmax(v, 1, false)
+	return i
+}
+
+func minmax(v interface{}, keep int, val bool) (reflect.Value, int) {
+	switch v := v.(type) {
+	case sort.Interface:
+		if v.Len() == 0 {
+			if keep < 0 {
+				panic("zero-length sequence has no minimum")
+			} else {
+				panic("zero-length sequence has no maximum")
+			}
+		}
+		maxi := 0
+		if keep < 0 {
+			for i, len := 0, v.Len(); i < len; i++ {
+				if v.Less(i, maxi) {
+					maxi = i
+				}
+			}
+		} else {
+			for i, len := 0, v.Len(); i < len; i++ {
+				if v.Less(maxi, i) {
+					maxi = i
+				}
+			}
+		}
+
+		if !val {
+			return reflect.Value{}, maxi
+		}
+
+		rv := reflectSlice(v)
+		return rv.Index(maxi), maxi
+	}
+
+	rv := reflectSlice(v)
+	if !generic.CanOrderR(rv.Type().Elem().Kind()) {
+		panic(&generic.TypeError{rv.Type().Elem(), nil, "is not orderable"})
+	}
+	if rv.Len() == 0 {
+		if keep < 0 {
+			panic("zero-length slice has no minimum")
+		} else {
+			panic("zero-length slice has no maximum")
+		}
+	}
+	max, maxi := rv.Index(0), 0
+	for i, len := 1, rv.Len(); i < len; i++ {
+		if elt := rv.Index(i); generic.OrderR(elt, max) == keep {
+			max, maxi = elt, i
+		}
+	}
+	return max, maxi
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/nub.go b/vendor/github.com/aclements/go-gg/generic/slice/nub.go
new file mode 100644
index 0000000..0c5f69c
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/nub.go
@@ -0,0 +1,52 @@
+// Copyright 2016 The Go 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 slice
+
+import "reflect"
+
+var trueVal = reflect.ValueOf(true)
+
+// Nub returns v with duplicates removed. It keeps the first instance
+// of each distinct value and preserves their order.
+func Nub(v T) T {
+	rv := reflectSlice(v)
+	indexes := make([]int, 0)
+	set := reflect.MakeMap(reflect.MapOf(rv.Type().Elem(), trueVal.Type()))
+	for i, l := 0, rv.Len(); i < l; i++ {
+		x := rv.Index(i)
+		if set.MapIndex(x).IsValid() {
+			continue
+		}
+		set.SetMapIndex(x, trueVal)
+		indexes = append(indexes, i)
+	}
+	return Select(v, indexes)
+}
+
+// NubAppend is equivalent to appending all of the slices in vs and
+// then calling Nub on the result, but more efficient.
+func NubAppend(vs ...T) T {
+	if len(vs) == 0 {
+		return nil
+	}
+
+	rv := reflectSlice(vs[0])
+	set := reflect.MakeMap(reflect.MapOf(rv.Type().Elem(), trueVal.Type()))
+	out := reflect.MakeSlice(rv.Type(), 0, 0)
+
+	for _, v := range vs {
+		rv := reflectSlice(v)
+		for i, l := 0, rv.Len(); i < l; i++ {
+			x := rv.Index(i)
+			if set.MapIndex(x).IsValid() {
+				continue
+			}
+			set.SetMapIndex(x, trueVal)
+			out = reflect.Append(out, x)
+		}
+	}
+
+	return out.Interface()
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/seq.go b/vendor/github.com/aclements/go-gg/generic/slice/seq.go
new file mode 100644
index 0000000..47812e4
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/seq.go
@@ -0,0 +1,28 @@
+// Copyright 2016 The Go 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 slice
+
+import (
+	"reflect"
+
+	"github.com/aclements/go-gg/generic"
+)
+
+// T is a Go slice value of type []U.
+//
+// This is primarily for documentation. There is no way to statically
+// enforce this in Go; however, functions that expect a slice will
+// panic with a *generic.TypeError if passed a non-slice value.
+type T interface{}
+
+// reflectSlice checks that s is a slice and returns its
+// reflect.Value. It panics with a *generic.TypeError if s is not a slice.
+func reflectSlice(s T) reflect.Value {
+	rv := reflect.ValueOf(s)
+	if rv.Kind() != reflect.Slice {
+		panic(&generic.TypeError{rv.Type(), nil, "is not a slice"})
+	}
+	return rv
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/sort.go b/vendor/github.com/aclements/go-gg/generic/slice/sort.go
new file mode 100644
index 0000000..e5ef8b6
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/sort.go
@@ -0,0 +1,137 @@
+// Copyright 2016 The Go 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 slice
+
+import (
+	"reflect"
+	"sort"
+	"time"
+
+	"github.com/aclements/go-gg/generic"
+)
+
+// CanSort returns whether the value v can be sorted.
+func CanSort(v interface{}) bool {
+	switch v.(type) {
+	case sort.Interface, []time.Time:
+		return true
+	}
+	return generic.CanOrderR(reflect.TypeOf(v).Elem().Kind())
+}
+
+// Sort sorts v in increasing order. v must implement sort.Interface,
+// be a slice whose elements are orderable, or be a []time.Time.
+func Sort(v interface{}) {
+	sort.Sort(Sorter(v))
+}
+
+// Sorter returns a sort.Interface for sorting v. v must implement
+// sort.Interface, be a slice whose elements are orderable, or be a
+// []time.Time.
+func Sorter(v interface{}) sort.Interface {
+	switch v := v.(type) {
+	case []int:
+		return sort.IntSlice(v)
+	case []float64:
+		return sort.Float64Slice(v)
+	case []string:
+		return sort.StringSlice(v)
+	case []time.Time:
+		return sortTimeSlice(v)
+	case sort.Interface:
+		return v
+	}
+
+	rv := reflectSlice(v)
+	switch rv.Type().Elem().Kind() {
+	case reflect.Int, reflect.Int8, reflect.Int16, reflect.Int32, reflect.Int64:
+		return sortIntSlice{rv}
+	case reflect.Uint, reflect.Uint8, reflect.Uint16, reflect.Uint32, reflect.Uint64, reflect.Uintptr:
+		return sortUintSlice{rv}
+	case reflect.Float32, reflect.Float64:
+		return sortFloatSlice{rv}
+	case reflect.String:
+		return sortStringSlice{rv}
+	}
+	panic(&generic.TypeError{rv.Type().Elem(), nil, "is not orderable"})
+}
+
+type sortIntSlice struct {
+	reflect.Value
+}
+
+func (s sortIntSlice) Len() int {
+	return s.Value.Len()
+}
+
+func (s sortIntSlice) Less(i, j int) bool {
+	return s.Index(i).Int() < s.Index(j).Int()
+}
+
+func (s sortIntSlice) Swap(i, j int) {
+	a, b := s.Index(i).Int(), s.Index(j).Int()
+	s.Index(i).SetInt(b)
+	s.Index(j).SetInt(a)
+}
+
+type sortUintSlice struct {
+	reflect.Value
+}
+
+func (s sortUintSlice) Len() int {
+	return s.Value.Len()
+}
+
+func (s sortUintSlice) Less(i, j int) bool {
+	return s.Index(i).Uint() < s.Index(j).Uint()
+}
+
+func (s sortUintSlice) Swap(i, j int) {
+	a, b := s.Index(i).Uint(), s.Index(j).Uint()
+	s.Index(i).SetUint(b)
+	s.Index(j).SetUint(a)
+}
+
+type sortFloatSlice struct {
+	reflect.Value
+}
+
+func (s sortFloatSlice) Len() int {
+	return s.Value.Len()
+}
+
+func (s sortFloatSlice) Less(i, j int) bool {
+	return s.Index(i).Float() < s.Index(j).Float()
+}
+
+func (s sortFloatSlice) Swap(i, j int) {
+	a, b := s.Index(i).Float(), s.Index(j).Float()
+	s.Index(i).SetFloat(b)
+	s.Index(j).SetFloat(a)
+}
+
+type sortStringSlice struct {
+	reflect.Value
+}
+
+func (s sortStringSlice) Len() int {
+	return s.Value.Len()
+}
+
+func (s sortStringSlice) Less(i, j int) bool {
+	return s.Index(i).String() < s.Index(j).String()
+}
+
+func (s sortStringSlice) Swap(i, j int) {
+	a, b := s.Index(i).String(), s.Index(j).String()
+	s.Index(i).SetString(b)
+	s.Index(j).SetString(a)
+}
+
+type sortTimeSlice []time.Time
+
+func (s sortTimeSlice) Len() int           { return len(s) }
+func (s sortTimeSlice) Less(i, j int) bool { return s[i].Before(s[j]) }
+func (s sortTimeSlice) Swap(i, j int)      { s[i], s[j] = s[j], s[i] }
diff --git a/vendor/github.com/aclements/go-gg/ggstat/agg.go b/vendor/github.com/aclements/go-gg/ggstat/agg.go
new file mode 100644
index 0000000..072d621
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/agg.go
@@ -0,0 +1,297 @@
+// Copyright 2016 The Go 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 ggstat
+
+import (
+	"fmt"
+	"reflect"
+
+	"github.com/aclements/go-gg/generic/slice"
+	"github.com/aclements/go-gg/table"
+	"github.com/aclements/go-moremath/stats"
+	"github.com/aclements/go-moremath/vec"
+)
+
+// TODO: AggFirst, AggTukey. StdDev?
+
+// Agg constructs an Aggregate transform from a grouping column and a
+// set of Aggregators.
+//
+// TODO: Does this belong in ggstat? The specific aggregator functions
+// probably do, but the concept could go in package table.
+func Agg(xs ...string) func(aggs ...Aggregator) Aggregate {
+	return func(aggs ...Aggregator) Aggregate {
+		return Aggregate{xs, aggs}
+	}
+}
+
+// Aggregate computes aggregate functions of a table grouped by
+// distinct values of a column or set of columns.
+//
+// Aggregate first groups the table by the Xs columns. Each of these
+// groups produces a single row in the output table, where the unique
+// value of each of the Xs columns appears in the output row, along
+// with constant columns from the input, as well as any columns that
+// have a unique value within every group (they're "effectively"
+// constant). Additional columns in the output row are produced by
+// applying the Aggregator functions to the group.
+type Aggregate struct {
+	// Xs is the list column names to group values by before
+	// computing aggregate functions.
+	Xs []string
+
+	// Aggregators is the set of Aggregator functions to apply to
+	// each group of values.
+	Aggregators []Aggregator
+}
+
+// An Aggregator is a function that aggregates each group of input
+// into one row and adds it to output. It may be based on multiple
+// columns from input and may add multiple columns to output.
+type Aggregator func(input table.Grouping, output *table.Builder)
+
+func (s Aggregate) F(g table.Grouping) table.Grouping {
+	isConst := make([]bool, len(g.Columns()))
+	for i := range isConst {
+		isConst[i] = true
+	}
+
+	subgroups := map[table.GroupID]table.Grouping{}
+	for _, gid := range g.Tables() {
+		g := table.GroupBy(g.Table(gid), s.Xs...)
+		subgroups[gid] = g
+
+		for i, col := range g.Columns() {
+			if !isConst[i] {
+				continue
+			}
+			// Can this column be promoted to constant?
+			for _, gid2 := range g.Tables() {
+				t := g.Table(gid2)
+				isConst[i] = isConst[i] && checkConst(t, col)
+			}
+		}
+	}
+
+	return table.MapTables(g, func(_ table.GroupID, t *table.Table) *table.Table {
+		g := table.GroupBy(t, s.Xs...)
+		var nt table.Builder
+
+		// Construct X columns.
+		rows := len(g.Tables())
+		for colidx, xcol := range s.Xs {
+			xs := reflect.MakeSlice(table.ColType(t, xcol), rows, rows)
+			for i, gid := range g.Tables() {
+				for j := 0; j < len(s.Xs)-colidx-1; j++ {
+					gid = gid.Parent()
+				}
+				xs.Index(i).Set(reflect.ValueOf(gid.Label()))
+			}
+
+			nt.Add(xcol, xs.Interface())
+		}
+
+		// Apply Aggregators.
+		for _, agg := range s.Aggregators {
+			agg(g, &nt)
+		}
+
+		// Keep constant and effectively constant columns.
+		for i := range isConst {
+			col := t.Columns()[i]
+			if !isConst[i] || nt.Has(col) {
+				continue
+			}
+			if cv, ok := t.Const(col); ok {
+				nt.AddConst(col, cv)
+				continue
+			}
+
+			ncol := reflect.MakeSlice(table.ColType(t, col), len(g.Tables()), len(g.Tables()))
+			for i, gid := range g.Tables() {
+				v := reflect.ValueOf(g.Table(gid).Column(col))
+				ncol.Index(i).Set(v.Index(0))
+			}
+			nt.Add(col, ncol.Interface())
+		}
+		return nt.Done()
+	})
+}
+
+func checkConst(t *table.Table, col string) bool {
+	if _, ok := t.Const(col); ok {
+		return true
+	}
+	v := reflect.ValueOf(t.Column(col))
+	if v.Len() <= 1 {
+		return true
+	}
+	if !v.Type().Elem().Comparable() {
+		return false
+	}
+	elem := v.Index(0).Interface()
+	for i, l := 1, v.Len(); i < l; i++ {
+		if elem != v.Index(i).Interface() {
+			return false
+		}
+	}
+	return true
+}
+
+// AggCount returns an aggregate function that computes the number of
+// rows in each group. The resulting column will be named label, or
+// "count" if label is "".
+func AggCount(label string) Aggregator {
+	if label == "" {
+		label = "count"
+	}
+
+	return func(input table.Grouping, b *table.Builder) {
+		counts := make([]int, 0, len(input.Tables()))
+		for _, gid := range input.Tables() {
+			counts = append(counts, input.Table(gid).Len())
+		}
+		b.Add(label, counts)
+	}
+}
+
+// AggMean returns an aggregate function that computes the mean of
+// each of cols. The resulting columns will be named "mean <col>" and
+// will have the same type as <col>.
+func AggMean(cols ...string) Aggregator {
+	return aggFn(stats.Mean, "mean ", cols...)
+}
+
+// AggGeoMean returns an aggregate function that computes the
+// geometric mean of each of cols. The resulting columns will be named
+// "geomean <col>" and will have the same type as <col>.
+func AggGeoMean(cols ...string) Aggregator {
+	return aggFn(stats.GeoMean, "geomean ", cols...)
+}
+
+// AggMin returns an aggregate function that computes the minimum of
+// each of cols. The resulting columns will be named "min <col>" and
+// will have the same type as <col>.
+func AggMin(cols ...string) Aggregator {
+	min := func(xs []float64) float64 {
+		x, _ := stats.Bounds(xs)
+		return x
+	}
+	return aggFn(min, "min ", cols...)
+}
+
+// AggMax returns an aggregate function that computes the maximum of
+// each of cols. The resulting columns will be named "max <col>" and
+// will have the same type as <col>.
+func AggMax(cols ...string) Aggregator {
+	max := func(xs []float64) float64 {
+		_, x := stats.Bounds(xs)
+		return x
+	}
+	return aggFn(max, "max ", cols...)
+}
+
+// AggSum returns an aggregate function that computes the sum of each
+// of cols. The resulting columns will be named "sum <col>" and will
+// have the same type as <col>.
+func AggSum(cols ...string) Aggregator {
+	return aggFn(vec.Sum, "sum ", cols...)
+}
+
+// AggQuantile returns an aggregate function that computes a quantile
+// of each of cols. quantile has a range of [0,1]. The resulting
+// columns will be named "<prefix> <col>" and will have the same type
+// as <col>.
+func AggQuantile(prefix string, quantile float64, cols ...string) Aggregator {
+	// "prefix" could be autogenerated (e.g. fmt.Sprintf("p%g ",
+	// quantile * 100)), but then the caller would need to do the
+	// same fmt.Sprintf to compute the column name they had just
+	// created. Perhaps Aggregator should provide a way to find
+	// the generated column names.
+	return aggFn(func(data []float64) float64 {
+		return stats.Sample{Xs: data}.Quantile(quantile)
+	}, prefix+" ", cols...)
+}
+
+func aggFn(f func([]float64) float64, prefix string, cols ...string) Aggregator {
+	ocols := make([]string, len(cols))
+	for i, col := range cols {
+		ocols[i] = prefix + col
+	}
+
+	return func(input table.Grouping, b *table.Builder) {
+		for coli, col := range cols {
+			means := make([]float64, 0, len(input.Tables()))
+
+			var xs []float64
+			var ct reflect.Type
+			for i, gid := range input.Tables() {
+				v := input.Table(gid).MustColumn(col)
+				if i == 0 {
+					ct = reflect.TypeOf(v)
+				}
+				slice.Convert(&xs, v)
+				means = append(means, f(xs))
+			}
+
+			if ct == float64SliceType {
+				b.Add(ocols[coli], means)
+			} else {
+				// Convert means back to the type of col.
+				outptr := reflect.New(ct)
+				slice.Convert(outptr.Interface(), means)
+				b.Add(ocols[coli], outptr.Elem().Interface())
+			}
+		}
+	}
+}
+
+// AggUnique returns an aggregate function retains the unique value of
+// each of cols within each aggregate group, or panics if some group
+// contains more than one value for one of these columns.
+//
+// Note that Aggregate will automatically retain columns that happen
+// to be unique. AggUnique can be used to enforce at aggregation time
+// that certain columns *must* be unique (and get a nice error if they
+// are not).
+func AggUnique(cols ...string) Aggregator {
+	return func(input table.Grouping, b *table.Builder) {
+		if len(cols) == 0 {
+			return
+		}
+		if len(input.Tables()) == 0 {
+			panic(fmt.Sprintf("unknown column: %q", cols[0]))
+		}
+
+		for _, col := range cols {
+			ctype := table.ColType(input, col)
+			rows := len(input.Tables())
+			vs := reflect.MakeSlice(ctype, rows, rows)
+			for i, gid := range input.Tables() {
+				// Get values in this column.
+				xs := reflect.ValueOf(input.Table(gid).MustColumn(col))
+
+				// Check for uniqueness.
+				if xs.Len() == 0 {
+					panic(fmt.Sprintf("cannot AggUnique empty column %q", col))
+				}
+				uniquev := xs.Index(0)
+				unique := uniquev.Interface()
+				for i, len := 1, xs.Len(); i < len; i++ {
+					other := xs.Index(i).Interface()
+					if unique != other {
+						panic(fmt.Sprintf("column %q is not unique; contains at least %v and %v", col, unique, other))
+					}
+				}
+
+				// Store unique value.
+				vs.Index(i).Set(uniquev)
+			}
+
+			// Add unique values slice to output table.
+			b.Add(col, vs.Interface())
+		}
+	}
+}
diff --git a/vendor/github.com/aclements/go-gg/ggstat/bin.go b/vendor/github.com/aclements/go-gg/ggstat/bin.go
new file mode 100644
index 0000000..2602004
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/bin.go
@@ -0,0 +1,196 @@
+// Copyright 2016 The Go 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 ggstat
+
+import (
+	"math"
+	"reflect"
+	"sort"
+
+	"github.com/aclements/go-gg/generic"
+	"github.com/aclements/go-gg/generic/slice"
+	"github.com/aclements/go-gg/table"
+	"github.com/aclements/go-moremath/vec"
+)
+
+// XXX If this is just based on the number of bins, it can come up
+// with really ugly boundary numbers. If the bin width is specified,
+// then you could also specify the left edge and bins will be placed
+// at [align+width*N, align+width*(N+1)]. ggplot2 also lets you
+// specify the center alignment.
+//
+// XXX In Matlab and NumPy, bins are open on the right *except* for
+// the last bin, which is closed on both.
+//
+// XXX Number of bins/bin width/specify boundaries, same bins across
+// all groups/separate for each group/based on shared scales (don't
+// have that information here), relative or absolute histogram (Matlab
+// has lots more).
+//
+// XXX Scale transform.
+//
+// The result of Bin has two columns in addition to constant columns from the input:
+//
+// - Column X is the left edge of the bin.
+//
+// - Column W is the sum of the rows' weights, or column "count" is
+//   the number of rows in the bin.
+type Bin struct {
+	// X is the name of the column to use for samples.
+	X string
+
+	// W is the optional name of the column to use for sample
+	// weights. It may be "" to weight each sample as 1.
+	W string
+
+	// Width controls how wide each bin should be. If not provided
+	// or 0, a width will be chosen to produce 30 bins. If X is an
+	// integer column, this width will be treated as an integer as
+	// well.
+	Width float64
+
+	// Center controls the center point of each bin. To center on
+	// integers, for example, you could use {Width: 1, Center:
+	// 0}.
+	// XXX What does center mean for integers? Should an unspecified center yield an autochosen one, or 0?
+	//Center float64
+
+	// Breaks is the set of break points to use as boundaries
+	// between bins. The interval of each bin is [Breaks[i],
+	// Breaks[i+1]). Data points before the first break are
+	// dropped. If provided, Width and Center are ignored.
+	Breaks table.Slice
+
+	// SplitGroups indicates that each group in the table should
+	// have separate bounds based on the data in that group alone.
+	// The default, false, indicates that the binning function
+	// should use the bounds of all of the data combined. This
+	// makes it easier to compare bins across groups.
+	SplitGroups bool
+}
+
+func (b Bin) F(g table.Grouping) table.Grouping {
+	breaks := reflect.ValueOf(b.Breaks)
+	agg := AggCount("count")
+	if b.W != "" {
+		agg = aggFn(vec.Sum, "", b.W)
+	}
+	if !breaks.IsValid() && !b.SplitGroups {
+		breaks = b.computeBreaks(g)
+	}
+	// Change b.X to the start of the bin.
+	g = table.MapTables(g, func(_ table.GroupID, t *table.Table) *table.Table {
+		breaks := breaks
+		if !breaks.IsValid() {
+			breaks = b.computeBreaks(t)
+		}
+		nbreaks := breaks.Len()
+
+		in := reflect.ValueOf(t.MustColumn(b.X))
+		nin := in.Len()
+
+		out := reflect.MakeSlice(breaks.Type(), nin, nin)
+		var found []int
+		for i := 0; i < nin; i++ {
+			elt := in.Index(i)
+			bin := sort.Search(nbreaks, func(j int) bool {
+				return generic.OrderR(elt, breaks.Index(j)) < 0
+			})
+			// 0 means the row doesn't fit on the front
+			// XXX Allow configuring the first and last bin as infinite or not.
+			bin = bin - 1
+			if bin >= 0 {
+				found = append(found, i)
+				out.Index(i).Set(breaks.Index(bin))
+			}
+		}
+		var nt table.Builder
+		for _, col := range t.Columns() {
+			if col == b.X {
+				nt.Add(col, slice.Select(out.Interface(), found))
+			} else if c, ok := t.Const(col); ok {
+				nt.AddConst(col, c)
+			} else {
+				nt.Add(col, slice.Select(t.Column(col), found))
+			}
+		}
+		return nt.Done()
+	})
+	// Group by the found bin
+	return Agg(b.X)(agg).F(g)
+}
+
+func (b Bin) computeBreaks(g table.Grouping) reflect.Value {
+	var cols []slice.T
+	for _, gid := range g.Tables() {
+		cols = append(cols, g.Table(gid).MustColumn(b.X))
+	}
+	data := slice.Concat(cols...)
+
+	min := slice.Min(data)
+	max := slice.Max(data)
+
+	rv := reflect.ValueOf(min)
+	switch rv.Type().Kind() {
+	case reflect.Int, reflect.Int8, reflect.Int16, reflect.Int32, reflect.Int64:
+		min, max := rv.Int(), reflect.ValueOf(max).Int()
+		width := int64(b.Width)
+		if width == 0 {
+			width = (max - min) / 30
+			if width < 1 {
+				width = 1
+			}
+		}
+		// XXX: This assumes boundaries should be aligned with
+		// 0. We should support explicit Center or Boundary
+		// requests.
+		min -= (min % width)
+		var breaks []int64
+		for i := min; i < max; i += width {
+			breaks = append(breaks, i)
+		}
+		outs := reflect.New(reflect.ValueOf(cols[0]).Type())
+		slice.Convert(outs.Interface(), breaks)
+		return outs.Elem()
+	case reflect.Uint, reflect.Uint8, reflect.Uint16, reflect.Uint32, reflect.Uint64, reflect.Uintptr:
+		min, max := rv.Uint(), reflect.ValueOf(max).Uint()
+		width := uint64(b.Width)
+		if width == 0 {
+			width = (max - min) / 30
+			if width < 1 {
+				width = 1
+			}
+		}
+		min -= (min % width)
+		var breaks []uint64
+		for i := min; i < max; i += width {
+			breaks = append(breaks, i)
+		}
+		outs := reflect.New(reflect.ValueOf(cols[0]).Type())
+		slice.Convert(outs.Interface(), breaks)
+		return outs.Elem()
+	case reflect.Float32, reflect.Float64:
+		min, max := rv.Float(), reflect.ValueOf(max).Float()
+		width := b.Width
+		if width == 0 {
+			width = (max - min) / 30
+			if width == 0 {
+				width = 1
+			}
+		}
+		min -= math.Mod(min, width)
+		var breaks []float64
+		for i := min; i < max; i += width {
+			breaks = append(breaks, i)
+		}
+		outs := reflect.New(reflect.ValueOf(cols[0]).Type())
+		slice.Convert(outs.Interface(), breaks)
+		return outs.Elem()
+	default:
+		panic("can't compute breaks for unknown type")
+	}
+}
+
+// TODO: Count for categorical data.
diff --git a/vendor/github.com/aclements/go-gg/ggstat/common.go b/vendor/github.com/aclements/go-gg/ggstat/common.go
new file mode 100644
index 0000000..db0d7a5
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/common.go
@@ -0,0 +1,10 @@
+// Copyright 2016 The Go 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 ggstat
+
+import "reflect"
+
+var float64Type = reflect.TypeOf(float64(0))
+var float64SliceType = reflect.TypeOf([]float64(nil))
diff --git a/vendor/github.com/aclements/go-gg/ggstat/density.go b/vendor/github.com/aclements/go-gg/ggstat/density.go
new file mode 100644
index 0000000..73b8a60
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/density.go
@@ -0,0 +1,123 @@
+// Copyright 2016 The Go 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 ggstat
+
+import (
+	"github.com/aclements/go-gg/generic/slice"
+	"github.com/aclements/go-gg/table"
+	"github.com/aclements/go-moremath/stats"
+	"github.com/aclements/go-moremath/vec"
+)
+
+// TODO: Default to first (and second) column for X (and Y)?
+
+// Density constructs a probability density estimate from a set of
+// samples using kernel density estimation.
+//
+// X is the only required field. All other fields have reasonable
+// default zero values.
+//
+// The result of Density has three columns in addition to constant
+// columns from the input:
+//
+// - Column X is the points at which the density estimate is sampled.
+//
+// - Column "probability density" is the density estimate.
+//
+// - Column "cumulative density" is the cumulative density estimate.
+type Density struct {
+	// X is the name of the column to use for samples.
+	X string
+
+	// W is the optional name of the column to use for sample
+	// weights. It may be "" to uniformly weight samples.
+	W string
+
+	// N is the number of points to sample the KDE at. If N is 0,
+	// a reasonable default is used.
+	//
+	// TODO: This is particularly sensitive to the scale
+	// transform.
+	//
+	// TODO: Base the default on the bandwidth. If the bandwidth
+	// is really narrow, we may need a lot of samples to exceed
+	// the Nyquist rate.
+	N int
+
+	// Domain specifies the domain at which to sample this function.
+	// If Domain is nil, it defaults to DomainData{}.
+	Domain FunctionDomainer
+
+	// Kernel is the kernel to use for the KDE.
+	Kernel stats.KDEKernel
+
+	// Bandwidth is the bandwidth to use for the KDE.
+	//
+	// If this is zero, the bandwidth is computed from the data
+	// using a default bandwidth estimator (currently
+	// stats.BandwidthScott).
+	Bandwidth float64
+
+	// BoundaryMethod is the boundary correction method to use for
+	// the KDE. The default value is BoundaryReflect; however, the
+	// default bounds are effectively +/-inf, which is equivalent
+	// to performing no boundary correction.
+	BoundaryMethod stats.KDEBoundaryMethod
+
+	// [BoundaryMin, BoundaryMax) specify a bounded support for
+	// the KDE. If both are 0 (their default values), they are
+	// treated as +/-inf.
+	//
+	// To specify a half-bounded support, set Min to math.Inf(-1)
+	// or Max to math.Inf(1).
+	BoundaryMin float64
+	BoundaryMax float64
+}
+
+func (d Density) F(g table.Grouping) table.Grouping {
+	kde := stats.KDE{
+		Kernel:         d.Kernel,
+		Bandwidth:      d.Bandwidth,
+		BoundaryMethod: d.BoundaryMethod,
+		BoundaryMin:    d.BoundaryMin,
+		BoundaryMax:    d.BoundaryMax,
+	}
+	dname, cname := "probability density", "cumulative density"
+
+	addEmpty := func(out *table.Builder) {
+		out.Add(dname, []float64{})
+		out.Add(cname, []float64{})
+	}
+
+	return Function{
+		X: d.X, N: d.N, Domain: d.Domain,
+		Fn: func(gid table.GroupID, in *table.Table, sampleAt []float64, out *table.Builder) {
+			if len(sampleAt) == 0 {
+				addEmpty(out)
+				return
+			}
+
+			// Get input sample.
+			var sample stats.Sample
+			slice.Convert(&sample.Xs, in.MustColumn(d.X))
+			if d.W != "" {
+				slice.Convert(&sample.Weights, in.MustColumn(d.W))
+				if sample.Weight() == 0 {
+					addEmpty(out)
+					return
+				}
+			}
+
+			// Compute KDE.
+			kde.Sample = sample
+			if d.Bandwidth == 0 {
+				kde.Bandwidth = stats.BandwidthScott(sample)
+			}
+
+			out.Add(dname, vec.Map(kde.PDF, sampleAt))
+			out.Add(cname, vec.Map(kde.CDF, sampleAt))
+		},
+	}.F(g)
+}
diff --git a/vendor/github.com/aclements/go-gg/ggstat/domain.go b/vendor/github.com/aclements/go-gg/ggstat/domain.go
new file mode 100644
index 0000000..77bb78e
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/domain.go
@@ -0,0 +1,109 @@
+// Copyright 2016 The Go 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 ggstat
+
+import (
+	"math"
+
+	"github.com/aclements/go-gg/generic/slice"
+	"github.com/aclements/go-gg/table"
+	"github.com/aclements/go-moremath/stats"
+)
+
+// A FunctionDomainer computes the domain over which to evaluate a
+// statistical function.
+type FunctionDomainer interface {
+	// FunctionDomain computes the domain of a particular column
+	// within a table. It takes a Grouping and a column in that
+	// Grouping to compute the domain of and returns a function
+	// that computes the domain for a specific group in the
+	// Grouping. This makes it possible for FunctionDomain to
+	// easily compute either Grouping-wide domains, or per-Table
+	// domains.
+	//
+	// The returned domain may be (NaN, NaN) to indicate that
+	// there is no data and the domain is vacuous.
+	FunctionDomain(g table.Grouping, col string) func(gid table.GroupID) (min, max float64)
+}
+
+// DomainFixed is a FunctionDomainer that returns a fixed domain.
+type DomainFixed struct {
+	Min, Max float64
+}
+
+var _ FunctionDomainer = DomainFixed{}
+
+func (r DomainFixed) FunctionDomain(g table.Grouping, col string) func(gid table.GroupID) (min, max float64) {
+	return func(table.GroupID) (min, max float64) {
+		return r.Min, r.Max
+	}
+}
+
+// DomainData is a FunctionDomainer that computes domains based on the
+// bounds of the data.
+type DomainData struct {
+	// Widen expands the domain by Widen times the span of the
+	// data.
+	//
+	// A value of 1.0 means to use exactly the bounds of the data.
+	// If Widen is 0, it is treated as 1.1 (that is, widen the
+	// domain by 10%, or 5% on the left and 5% on the right).
+	Widen float64
+
+	// SplitGroups indicates that each group in the table should
+	// have a separate domain based on the data in that group
+	// alone. The default, false, indicates that the domain should
+	// be based on all of the data in the table combined. This
+	// makes it possible to stack functions and easier to compare
+	// them across groups.
+	SplitGroups bool
+}
+
+var _ FunctionDomainer = DomainData{}
+
+const defaultWiden = 1.1
+
+func (r DomainData) FunctionDomain(g table.Grouping, col string) func(gid table.GroupID) (min, max float64) {
+	widen := r.Widen
+	if widen <= 0 {
+		widen = defaultWiden
+	}
+
+	var xs []float64
+	if !r.SplitGroups {
+		// Compute combined bounds.
+		gmin, gmax := math.NaN(), math.NaN()
+		for _, gid := range g.Tables() {
+			t := g.Table(gid)
+			slice.Convert(&xs, t.MustColumn(col))
+			xmin, xmax := stats.Bounds(xs)
+			if xmin < gmin || math.IsNaN(gmin) {
+				gmin = xmin
+			}
+			if xmax > gmax || math.IsNaN(gmax) {
+				gmax = xmax
+			}
+		}
+
+		// Widen bounds.
+		span := gmax - gmin
+		gmin, gmax = gmin-span*(widen-1)/2, gmax+span*(widen-1)/2
+
+		return func(table.GroupID) (min, max float64) {
+			return gmin, gmax
+		}
+	}
+
+	return func(gid table.GroupID) (min, max float64) {
+		// Compute bounds.
+		slice.Convert(&xs, g.Table(gid).MustColumn(col))
+		min, max = stats.Bounds(xs)
+
+		// Widen bounds.
+		span := max - min
+		min, max = min-span*(widen-1)/2, max+span*(widen-1)/2
+		return
+	}
+}
diff --git a/vendor/github.com/aclements/go-gg/ggstat/ecdf.go b/vendor/github.com/aclements/go-gg/ggstat/ecdf.go
new file mode 100644
index 0000000..46364a8
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/ecdf.go
@@ -0,0 +1,138 @@
+// Copyright 2016 The Go 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 ggstat
+
+import (
+	"github.com/aclements/go-gg/generic/slice"
+	"github.com/aclements/go-gg/table"
+	"github.com/aclements/go-moremath/vec"
+)
+
+// ECDF constructs an empirical CDF from a set of samples.
+//
+// X is the only required field. All other fields have reasonable
+// default zero values.
+//
+// The result of ECDF has three columns in addition to constant
+// columns from the input. The names of the columns depend on whether
+// Label is "".
+//
+// - Column X is the points at which the CDF changes (a subset of the
+// samples).
+//
+// - Column "cumulative density" or "cumulative density of <label>" is
+// the cumulative density estimate.
+//
+// - Column "cumulative count" (if W and Label are ""), "cumulative
+// weight" (if W is not "", but Label is "") or "cumulative <label>"
+// (if Label is not "") is the cumulative count or weight of samples.
+// That is, cumulative density times the total weight of the samples.
+type ECDF struct {
+	// X is the name of the column to use for samples.
+	X string
+
+	// W is the optional name of the column to use for sample
+	// weights. It may be "" to uniformly weight samples.
+	W string
+
+	// Label, if not "", gives a label for the samples. It is used
+	// to construct more specific names for the output columns. It
+	// should be a plural noun.
+	Label string
+
+	// Domain specifies the domain of the returned ECDF. If the
+	// domain is wider than the bounds of the data in a group,
+	// ECDF will add a point below the smallest sample and above
+	// the largest sample to make the 0 and 1 levels clear. If
+	// Domain is nil, it defaults to DomainData{}.
+	Domain FunctionDomainer
+}
+
+func (s ECDF) F(g table.Grouping) table.Grouping {
+	// Set defaults.
+	if s.Domain == nil {
+		s.Domain = DomainData{}
+	}
+
+	// Construct output column names.
+	dname, cname := "cumulative density", "cumulative count"
+	if s.Label != "" {
+		dname += " of " + s.Label
+		cname = "cumulative " + s.Label
+	} else if s.W != "" {
+		cname = "cumulative weight"
+	}
+
+	g = table.SortBy(g, s.X)
+	domain := s.Domain.FunctionDomain(g, s.X)
+
+	return table.MapTables(g, func(gid table.GroupID, t *table.Table) *table.Table {
+		// Get input columns.
+		var xs, ws []float64
+		slice.Convert(&xs, t.MustColumn(s.X))
+		if s.W != "" {
+			slice.Convert(&ws, t.MustColumn(s.W))
+		}
+
+		// Ignore empty tables.
+		if len(xs) == 0 {
+			nt := new(table.Builder).Add(s.X, []float64{}).Add(cname, []float64{}).Add(dname, []float64{})
+			preserveConsts(nt, t)
+			return nt.Done()
+		}
+
+		// Get domain.
+		min, max := domain(gid)
+
+		// Create output columns.
+		xo, do, co := make([]float64, 0), make([]float64, 0), make([]float64, 0)
+		if min < xs[0] {
+			// Extend to the left.
+			xo = append(xo, min)
+			do = append(do, 0)
+			co = append(co, 0)
+		}
+
+		// Compute total weight.
+		var total float64
+		if ws == nil {
+			total = float64(t.Len())
+		} else {
+			total = vec.Sum(ws)
+		}
+
+		// Create ECDF.
+		cum := 0.0
+		for i := 0; i < len(xs); {
+			j := i
+			for j < len(xs) && xs[i] == xs[j] {
+				if ws == nil {
+					cum += 1
+				} else {
+					cum += ws[j]
+				}
+				j++
+			}
+
+			xo = append(xo, xs[i])
+			do = append(do, cum/total)
+			co = append(co, cum)
+
+			i = j
+		}
+
+		if xs[len(xs)-1] < max {
+			// Extend to the right.
+			xo = append(xo, max)
+			do = append(do, 1)
+			co = append(co, cum)
+		}
+
+		// Construct results table.
+		nt := new(table.Builder).Add(s.X, xo).Add(dname, do).Add(cname, co)
+		preserveConsts(nt, t)
+		return nt.Done()
+	})
+}
diff --git a/vendor/github.com/aclements/go-gg/ggstat/fn.go b/vendor/github.com/aclements/go-gg/ggstat/fn.go
new file mode 100644
index 0000000..0a4de33
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/fn.go
@@ -0,0 +1,123 @@
+// Copyright 2016 The Go 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 ggstat
+
+import (
+	"math"
+	"reflect"
+
+	"github.com/aclements/go-gg/generic/slice"
+	"github.com/aclements/go-gg/table"
+	"github.com/aclements/go-moremath/vec"
+)
+
+// Function samples a continuous univariate function at N points in
+// the domain computed by Domain.
+//
+// The result of Function binds column X to the X values at which the
+// function is sampled and retains constant columns from the input.
+// The computed function can add arbitrary columns for its output.
+type Function struct {
+	// X is the name of the column to use for input domain of this
+	// function.
+	X string
+
+	// N is the number of points to sample the function at. If N
+	// is 0, a reasonable default is used.
+	N int
+
+	// Domain specifies the domain of which to sample this function.
+	// If Domain is nil, it defaults to DomainData{}.
+	Domain FunctionDomainer
+
+	// Fn is the continuous univariate function to sample. Fn will
+	// be called with each table in the grouping and the X values
+	// at which it should be sampled. Fn must add its output
+	// columns to out. The output table will already contain the
+	// sample points bound to the X column.
+	Fn func(gid table.GroupID, in *table.Table, sampleAt []float64, out *table.Builder)
+}
+
+const defaultFunctionSamples = 200
+
+func (f Function) F(g table.Grouping) table.Grouping {
+	// Set defaults.
+	if f.N <= 0 {
+		f.N = defaultFunctionSamples
+	}
+	if f.Domain == nil {
+		f.Domain = DomainData{}
+	}
+
+	domain := f.Domain.FunctionDomain(g, f.X)
+	return table.MapTables(g, func(gid table.GroupID, t *table.Table) *table.Table {
+		min, max := domain(gid)
+
+		// Compute sample points. If there's no data, there
+		// are no sample points, but we still have to run the
+		// function to get the right output columns.
+		var ss []float64
+		if math.IsNaN(min) {
+			ss = []float64{}
+		} else {
+			ss = vec.Linspace(min, max, f.N)
+		}
+
+		var nt table.Builder
+		ctype := table.ColType(t, f.X)
+		if ctype == float64Type {
+			// Bind output X column.
+			nt.Add(f.X, ss)
+		} else {
+			// Convert to the column type.
+			vsp := reflect.New(ctype)
+			slice.Convert(vsp.Interface(), ss)
+			vs := vsp.Elem()
+			// This may have produced duplicate values.
+			// Eliminate those.
+			if vs.Len() > 0 {
+				prev, i := vs.Index(0).Interface(), 1
+				for j := 1; j < vs.Len(); j++ {
+					next := vs.Index(j).Interface()
+					if prev == next {
+						// Skip duplicate.
+						continue
+					}
+
+					if i != j {
+						vs.Index(i).Set(vs.Index(j))
+					}
+					i++
+					prev = next
+				}
+				vs.SetLen(i)
+			}
+			// Bind column-typed values to output X.
+			nt.Add(f.X, vs.Interface())
+			// And convert back to []float64 so we can
+			// apply the function.
+			slice.Convert(&ss, vs.Interface())
+		}
+
+		// Apply the function to the sample points.
+		f.Fn(gid, t, ss, &nt)
+
+		preserveConsts(&nt, t)
+		return nt.Done()
+	})
+}
+
+// preserveConsts copies the constant columns from t into nt.
+func preserveConsts(nt *table.Builder, t *table.Table) {
+	for _, col := range t.Columns() {
+		if nt.Has(col) {
+			// Don't overwrite existing columns in nt.
+			continue
+		}
+		if cv, ok := t.Const(col); ok {
+			nt.AddConst(col, cv)
+		}
+	}
+}
diff --git a/vendor/github.com/aclements/go-gg/ggstat/normalize.go b/vendor/github.com/aclements/go-gg/ggstat/normalize.go
new file mode 100644
index 0000000..4644246
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/normalize.go
@@ -0,0 +1,175 @@
+// Copyright 2016 The Go 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 ggstat
+
+import (
+	"reflect"
+
+	"github.com/aclements/go-gg/generic/slice"
+	"github.com/aclements/go-gg/table"
+)
+
+// Normalize normalizes each group such that some data point is 1.
+//
+// Either X or Index is required (though 0 is a reasonable value of
+// Index).
+//
+// The result of Normalize is the same as the input table, plus
+// additional columns for each normalized column. These columns will
+// be named "normalized <col>" where <col> is the name of the original
+// column and will have type []float64.
+type Normalize struct {
+	// X is the name of the column to use to find the denominator
+	// row. If X is "", Index is used instead.
+	X string
+
+	// Index is the row index of the denominator row if X is ""
+	// (otherwise it is ignored). Index may be negative, in which
+	// case it is added to the number of rows (e.g., -1 is the
+	// last row).
+	Index int
+
+	// By is a function func([]T) int that returns the index of
+	// the denominator row given column X. By may be nil, in which
+	// case it defaults to generic.ArgMin.
+	By interface{}
+
+	// Cols is a slice of the names of columns to normalize
+	// relative to the corresponding DenomCols value in the
+	// denominator row. Cols may be nil, in which case it defaults
+	// to all integral and floating point columns.
+	Cols []string
+
+	// DenomCols is a slice of the names of columns used as the
+	// demoninator. DenomCols may be nil, in which case it
+	// defaults to Cols (i.e. each column will be normalized to
+	// the value from that column in the denominator row.)
+	// Otherwise, DenomCols must be the same length as Cols.
+	DenomCols []string
+}
+
+func (s Normalize) F(g table.Grouping) table.Grouping {
+	// Find the columns to normalize.
+	if s.Cols == nil {
+		cols := []string{}
+		for i, ct := range colTypes(g) {
+			if canNormalize(ct.Elem().Kind()) {
+				cols = append(cols, g.Columns()[i])
+			}
+		}
+		s.Cols = cols
+	}
+	if len(s.Cols) == 0 {
+		return g
+	}
+
+	// Construct new column names.
+	newcols := make([]string, len(s.Cols))
+	for i, col := range s.Cols {
+		newcols[i] = "normalized " + col
+	}
+
+	// Get "by" function.
+	var byv reflect.Value
+	byargs := make([]reflect.Value, 1)
+	if s.By != nil {
+		byv = reflect.ValueOf(s.By)
+		// TODO: Type check byv better.
+	}
+
+	return table.MapTables(g, func(_ table.GroupID, t *table.Table) *table.Table {
+		if t.Len() == 0 {
+			return t
+		}
+
+		// Find the denominator row.
+		var drow int
+		if s.X == "" {
+			drow = s.Index
+			if drow < 0 {
+				drow += t.Len()
+			}
+		} else {
+			xs := t.MustColumn(s.X)
+			if s.By == nil {
+				drow = slice.ArgMin(xs)
+			} else {
+				byargs[0] = reflect.ValueOf(xs)
+				byout := byv.Call(byargs)
+				drow = int(byout[0].Int())
+			}
+		}
+
+		// Normalize columns.
+		newt := table.NewBuilder(t)
+		denomCols := s.DenomCols
+		if denomCols == nil {
+			denomCols = s.Cols
+		}
+		for coli, col := range s.Cols {
+			denom := denomValue(t.MustColumn(denomCols[coli]), drow)
+			out := normalizeTo(t.MustColumn(col), denom)
+			newt.Add(newcols[coli], out)
+		}
+
+		return newt.Done()
+	})
+}
+
+func colTypes(g table.Grouping) []reflect.Type {
+	cts := make([]reflect.Type, len(g.Columns()))
+	for i, col := range g.Columns() {
+		cts[i] = table.ColType(g, col)
+	}
+	return cts
+}
+
+var canNormalizeKinds = map[reflect.Kind]bool{
+	reflect.Float32: true,
+	reflect.Float64: true,
+	reflect.Int:     true,
+	reflect.Int8:    true,
+	reflect.Int16:   true,
+	reflect.Int32:   true,
+	reflect.Int64:   true,
+	reflect.Uint:    true,
+	reflect.Uintptr: true,
+	reflect.Uint8:   true,
+	reflect.Uint16:  true,
+	reflect.Uint32:  true,
+	reflect.Uint64:  true,
+}
+
+func canNormalize(k reflect.Kind) bool {
+	return canNormalizeKinds[k]
+}
+
+func denomValue(s interface{}, index int) float64 {
+	switch s := s.(type) {
+	case []float64:
+		return s[index]
+	}
+	return reflect.ValueOf(s).Index(index).Convert(float64Type).Float()
+}
+
+func normalizeTo(s interface{}, denom float64) interface{} {
+	switch s := s.(type) {
+	case []float64:
+		out := make([]float64, len(s))
+		for i, numer := range s {
+			out[i] = numer / denom
+		}
+		return out
+	}
+
+	sv := reflect.ValueOf(s)
+
+	out := reflect.MakeSlice(float64SliceType, sv.Len(), sv.Len())
+	for i, len := 0, sv.Len(); i < len; i++ {
+		numer := sv.Index(i).Convert(float64Type).Float()
+		out.Index(i).SetFloat(numer / denom)
+	}
+	return out.Interface()
+}
diff --git a/vendor/github.com/aclements/go-gg/table/concat.go b/vendor/github.com/aclements/go-gg/table/concat.go
new file mode 100644
index 0000000..14a8f4f
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/concat.go
@@ -0,0 +1,104 @@
+// Copyright 2016 The Go 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 table
+
+import (
+	"fmt"
+
+	"github.com/aclements/go-gg/generic/slice"
+)
+
+// Concat returns the concatenation of the rows in each matching group
+// across gs. All Groupings in gs must have the same set of columns
+// (though they need not be in the same order; the column order from
+// gs[0] will be used). The GroupIDs in the returned Grouping will be
+// the union of the GroupIDs in gs.
+func Concat(gs ...Grouping) Grouping {
+	if len(gs) == 0 {
+		return new(Table)
+	}
+
+	// Check that all Groupings have the same set of columns. They
+	// can be in different orders.
+	colSet := map[string]bool{}
+	for _, col := range gs[0].Columns() {
+		colSet[col] = true
+	}
+	for i, g2 := range gs[1:] {
+		diff := len(g2.Columns()) != len(colSet)
+		if !diff {
+			for _, col := range g2.Columns() {
+				if !colSet[col] {
+					diff = true
+					break
+				}
+			}
+		}
+		if diff {
+			panic(fmt.Sprintf("columns in Groupings 0 and %d differ: %q vs %q", i+1, gs[0].Columns(), g2.Columns()))
+		}
+	}
+
+	// Collect group IDs.
+	haveGID := map[GroupID]bool{}
+	gids := []GroupID{}
+	for _, g := range gs {
+		for _, gid := range g.Tables() {
+			if haveGID[gid] {
+				continue
+			}
+			haveGID[gid] = true
+			gids = append(gids, gid)
+		}
+	}
+
+	// Build output groups.
+	var ng GroupingBuilder
+	for _, gid := range gids {
+		// Build output table.
+		var nt Builder
+		var cols []slice.T
+		for _, col := range gs[0].Columns() {
+			// Is it constant?
+			isConst := false
+			var cv interface{}
+			for _, g := range gs {
+				t := g.Table(gid)
+				if t == nil {
+					continue
+				}
+				if cv1, ok := t.Const(col); ok {
+					if !isConst {
+						isConst = true
+						cv = cv1
+					} else if cv != cv1 {
+						isConst = false
+						break
+					}
+				} else {
+					isConst = false
+					break
+				}
+			}
+			if isConst {
+				nt.AddConst(col, cv)
+				continue
+			}
+
+			// Not a constant. Collect slices.
+			for _, g := range gs {
+				t := g.Table(gid)
+				if t == nil {
+					continue
+				}
+				cols = append(cols, t.Column(col))
+			}
+			nt.Add(col, slice.Concat(cols...))
+			cols = cols[:0]
+		}
+		ng.Add(gid, nt.Done())
+	}
+	return ng.Done()
+}
diff --git a/vendor/github.com/aclements/go-gg/table/filter.go b/vendor/github.com/aclements/go-gg/table/filter.go
new file mode 100644
index 0000000..5008138
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/filter.go
@@ -0,0 +1,94 @@
+// Copyright 2016 The Go 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 table
+
+import (
+	"fmt"
+	"reflect"
+
+	"github.com/aclements/go-gg/generic/slice"
+)
+
+var boolType = reflect.TypeOf(false)
+
+// Filter filters g to only rows where pred returns true. pred must be
+// a function that returns bool and takes len(cols) arguments where
+// the type of col[i] is assignable to argument i.
+//
+// TODO: Create a faster batch variant where pred takes slices.
+func Filter(g Grouping, pred interface{}, cols ...string) Grouping {
+	// TODO: Use generic.TypeError.
+	predv := reflect.ValueOf(pred)
+	predt := predv.Type()
+	if predt.Kind() != reflect.Func || predt.NumIn() != len(cols) || predt.NumOut() != 1 || predt.Out(0) != boolType {
+		panic("predicate function must be func(col[0], col[1], ...) bool")
+	}
+	if len(cols) == 0 {
+		return g
+	}
+	if len(g.Tables()) == 0 {
+		panic(fmt.Sprintf("unknown column %q", cols[0]))
+	}
+	// Type check arguments.
+	for i, col := range cols {
+		colt := ColType(g, col)
+		if !colt.Elem().AssignableTo(predt.In(i)) {
+			panic(fmt.Sprintf("column %d (type %s) is not assignable to predicate argument %d (type %s)", i, colt.Elem(), i, predt.In(i)))
+		}
+	}
+
+	args := make([]reflect.Value, len(cols))
+	colvs := make([]reflect.Value, len(cols))
+	match := make([]int, 0)
+	return MapTables(g, func(_ GroupID, t *Table) *Table {
+		// Get columns.
+		for i, col := range cols {
+			colvs[i] = reflect.ValueOf(t.MustColumn(col))
+		}
+
+		// Find the set of row indexes that satisfy pred.
+		match = match[:0]
+		for r, len := 0, t.Len(); r < len; r++ {
+			for c, colv := range colvs {
+				args[c] = colv.Index(r)
+			}
+			if predv.Call(args)[0].Bool() {
+				match = append(match, r)
+			}
+		}
+
+		// Create the new table.
+		if len(match) == t.Len() {
+			return t
+		}
+		var nt Builder
+		for _, col := range t.Columns() {
+			nt.Add(col, slice.Select(t.Column(col), match))
+		}
+		return nt.Done()
+	})
+}
+
+// FilterEq filters g to only rows where the value in col equals val.
+func FilterEq(g Grouping, col string, val interface{}) Grouping {
+	match := make([]int, 0)
+	return MapTables(g, func(_ GroupID, t *Table) *Table {
+		// Find the set of row indexes that match val.
+		seq := t.MustColumn(col)
+		match = match[:0]
+		rv := reflect.ValueOf(seq)
+		for i, len := 0, rv.Len(); i < len; i++ {
+			if rv.Index(i).Interface() == val {
+				match = append(match, i)
+			}
+		}
+
+		var nt Builder
+		for _, col := range t.Columns() {
+			nt.Add(col, slice.Select(t.Column(col), match))
+		}
+		return nt.Done()
+	})
+}
diff --git a/vendor/github.com/aclements/go-gg/table/group.go b/vendor/github.com/aclements/go-gg/table/group.go
new file mode 100644
index 0000000..afabba2
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/group.go
@@ -0,0 +1,243 @@
+// Copyright 2016 The Go 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 table
+
+import (
+	"fmt"
+	"reflect"
+	"strings"
+
+	"github.com/aclements/go-gg/generic/slice"
+)
+
+// GroupID identifies a group. GroupIDs form a tree, rooted at
+// RootGroupID (which is also the zero GroupID).
+type GroupID struct {
+	*groupNode
+}
+
+// RootGroupID is the root of the GroupID tree.
+var RootGroupID = GroupID{}
+
+type groupNode struct {
+	parent GroupID
+	label  interface{}
+}
+
+// String returns the path to GroupID g in the form "/l1/l2/l3". If g
+// is RootGroupID, it returns "/". Each level in the group is formed
+// by formatting the label using fmt's "%v" verb. Note that this is
+// purely diagnostic; this string may not uniquely identify g.
+func (g GroupID) String() string {
+	if g == RootGroupID {
+		return "/"
+	}
+	parts := []string{}
+	for p := g; p != RootGroupID; p = p.parent {
+		part := fmt.Sprintf("/%v", p.label)
+		parts = append(parts, part)
+	}
+	for i, j := 0, len(parts)-1; i < j; i, j = i+1, j-1 {
+		parts[i], parts[j] = parts[j], parts[i]
+	}
+	return strings.Join(parts, "")
+}
+
+// Extend returns a new GroupID that is a child of GroupID g. The
+// returned GroupID will not be equal to any existing GroupID (even if
+// label is not unique among g's children). The label is primarily
+// diagnostic; the table package uses it only when printing tables,
+// but callers may store semantic information in group labels.
+func (g GroupID) Extend(label interface{}) GroupID {
+	return GroupID{&groupNode{g, label}}
+}
+
+// Parent returns the parent of g. The parent of RootGroupID is
+// RootGroupID.
+func (g GroupID) Parent() GroupID {
+	if g == RootGroupID {
+		return RootGroupID
+	}
+	return g.parent
+}
+
+// Label returns the label of g.
+func (g GroupID) Label() interface{} {
+	return g.label
+}
+
+// GroupBy sub-divides all groups such that all of the rows in each
+// group have equal values for all of the named columns. The relative
+// order of rows with equal values for the named columns is
+// maintained. Grouped-by columns become constant columns within each
+// group.
+func GroupBy(g Grouping, cols ...string) Grouping {
+	// TODO: This would generate much less garbage if we grouped
+	// all of cols in one pass.
+	//
+	// TODO: This constructs one slice per column per input group,
+	// but it would be even better if it constructed just one
+	// slice per column.
+
+	if len(cols) == 0 {
+		return g
+	}
+
+	var out GroupingBuilder
+	for _, gid := range g.Tables() {
+		t := g.Table(gid)
+
+		if cv, ok := t.Const(cols[0]); ok {
+			// Grouping by a constant is trivial.
+			subgid := gid.Extend(cv)
+			out.Add(subgid, t)
+			continue
+		}
+
+		c := t.MustColumn(cols[0])
+
+		// Create an index on c.
+		type subgroupInfo struct {
+			key  interface{}
+			rows []int
+		}
+		subgroups := []subgroupInfo{}
+		keys := make(map[interface{}]int)
+		seq := reflect.ValueOf(c)
+		for i := 0; i < seq.Len(); i++ {
+			x := seq.Index(i).Interface()
+			sg, ok := keys[x]
+			if !ok {
+				sg = len(subgroups)
+				subgroups = append(subgroups, subgroupInfo{x, []int{}})
+				keys[x] = sg
+			}
+			subgroup := &subgroups[sg]
+			subgroup.rows = append(subgroup.rows, i)
+		}
+
+		// Count rows in each subgroup.
+		offsets := make([]int, 1+len(subgroups))
+		for i := range subgroups {
+			offsets[i+1] = offsets[i] + len(subgroups[i].rows)
+		}
+
+		// Split each column.
+		builders := make([]Builder, len(subgroups))
+		for _, name := range t.Columns() {
+			if name == cols[0] {
+				// Promote the group-by column to a
+				// constant.
+				for i := range subgroups {
+					builders[i].AddConst(name, subgroups[i].key)
+				}
+				continue
+			}
+
+			if cv, ok := t.Const(name); ok {
+				// Keep constants constant.
+				for i := range builders {
+					builders[i].AddConst(name, cv)
+				}
+				continue
+			}
+
+			// Create a slice for all of the values.
+			col := t.Column(name)
+			ncol := reflect.MakeSlice(reflect.TypeOf(col), t.Len(), t.Len())
+
+			// Shuffle each subgroup into ncol.
+			for i := range subgroups {
+				subcol := ncol.Slice(offsets[i], offsets[i+1]).Interface()
+				slice.SelectInto(subcol, col, subgroups[i].rows)
+				builders[i].Add(name, subcol)
+			}
+		}
+
+		// Add tables to output Grouping.
+		for i := range builders {
+			subgid := gid.Extend(subgroups[i].key)
+			out.Add(subgid, builders[i].Done())
+		}
+	}
+
+	return GroupBy(out.Done(), cols[1:]...)
+}
+
+// Ungroup concatenates adjacent Tables in g that share a group parent
+// into a Table identified by the parent, undoing the effects of the
+// most recent GroupBy operation.
+func Ungroup(g Grouping) Grouping {
+	groups := g.Tables()
+	if len(groups) == 0 || len(groups) == 1 && groups[0] == RootGroupID {
+		return g
+	}
+
+	var out GroupingBuilder
+	runGid := groups[0].Parent()
+	runTabs := []*Table{}
+	for _, gid := range groups {
+		if gid.Parent() != runGid {
+			// Flush the run.
+			out.Add(runGid, concatRows(runTabs...))
+
+			runGid = gid.Parent()
+			runTabs = runTabs[:0]
+		}
+		runTabs = append(runTabs, g.Table(gid))
+	}
+	// Flush the last run.
+	out.Add(runGid, concatRows(runTabs...))
+
+	return out.Done()
+}
+
+// Flatten concatenates all of the groups in g into a single Table.
+// This is equivalent to repeatedly Ungrouping g.
+func Flatten(g Grouping) *Table {
+	groups := g.Tables()
+	switch len(groups) {
+	case 0:
+		return new(Table)
+
+	case 1:
+		return g.Table(groups[0])
+	}
+
+	tabs := make([]*Table, len(groups))
+	for i, gid := range groups {
+		tabs[i] = g.Table(gid)
+	}
+
+	return concatRows(tabs...)
+}
+
+// concatRows concatenates the rows of tabs into a single Table. All
+// Tables in tabs must all have the same column set.
+func concatRows(tabs ...*Table) *Table {
+	// TODO: Consider making this public. It would have to check
+	// the columns, and we would probably also want a concatCols.
+
+	switch len(tabs) {
+	case 0:
+		return new(Table)
+
+	case 1:
+		return tabs[0]
+	}
+
+	// Construct each column.
+	var out Builder
+	seqs := make([]slice.T, len(tabs))
+	for _, col := range tabs[0].Columns() {
+		seqs = seqs[:0]
+		for _, tab := range tabs {
+			seqs = append(seqs, tab.Column(col))
+		}
+		out.Add(col, slice.Concat(seqs...))
+	}
+
+	return out.Done()
+}
diff --git a/vendor/github.com/aclements/go-gg/table/head.go b/vendor/github.com/aclements/go-gg/table/head.go
new file mode 100644
index 0000000..3c1b2ec
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/head.go
@@ -0,0 +1,69 @@
+// Copyright 2016 The Go 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 table
+
+import "reflect"
+
+// Head returns the first n rows in each Table of g.
+func Head(g Grouping, n int) Grouping {
+	return headTail(g, n, false)
+}
+
+// Tail returns the last n rows in each Table of g.
+func Tail(g Grouping, n int) Grouping {
+	return headTail(g, n, true)
+}
+
+func headTail(g Grouping, n int, tail bool) Grouping {
+	return MapTables(g, func(_ GroupID, t *Table) *Table {
+		if t.Len() <= n {
+			return t
+		}
+
+		var nt Builder
+		for _, col := range t.Columns() {
+			if cv, ok := t.Const(col); ok {
+				nt.AddConst(col, cv)
+				continue
+			}
+
+			cv := reflect.ValueOf(t.Column(col))
+			if tail {
+				cv = cv.Slice(t.Len()-n, t.Len())
+			} else {
+				cv = cv.Slice(0, n)
+			}
+			nt.Add(col, cv.Interface())
+		}
+		return nt.Done()
+	})
+}
+
+// HeadTables returns the first n tables in g.
+func HeadTables(g Grouping, n int) Grouping {
+	return headTailTables(g, n, false)
+}
+
+// TailTables returns the first n tables in g.
+func TailTables(g Grouping, n int) Grouping {
+	return headTailTables(g, n, true)
+}
+
+func headTailTables(g Grouping, n int, tail bool) Grouping {
+	tables := g.Tables()
+	if len(tables) <= n {
+		return g
+	} else if tail {
+		tables = tables[len(tables)-n:]
+	} else {
+		tables = tables[:n]
+	}
+
+	var ng GroupingBuilder
+	for _, gid := range tables {
+		ng.Add(gid, g.Table(gid))
+	}
+	return ng.Done()
+}
diff --git a/vendor/github.com/aclements/go-gg/table/join.go b/vendor/github.com/aclements/go-gg/table/join.go
new file mode 100644
index 0000000..03ce51b
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/join.go
@@ -0,0 +1,77 @@
+// Copyright 2016 The Go 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 table
+
+import (
+	"reflect"
+
+	"github.com/aclements/go-gg/generic/slice"
+)
+
+// Join joins g1 and g2 on tables with identical group IDs where col1
+// in g1 equals col2 in g2. It maintains the group order of g1, except
+// that groups that aren't in g2 are removed, and maintains the row
+// order of g1, followed by the row order of g2.
+//
+// TODO: Support join on more than one column.
+func Join(g1 Grouping, col1 string, g2 Grouping, col2 string) Grouping {
+	var ng GroupingBuilder
+	for _, gid := range g1.Tables() {
+		t1, t2 := g1.Table(gid), g2.Table(gid)
+		if t2 == nil {
+			continue
+		}
+
+		// TODO: Optimize for cases where col1 and/or col2 are
+		// constant.
+
+		// Index col2 in t2.
+		ridx := make(map[interface{}][]int)
+		rv := reflect.ValueOf(t2.MustColumn(col2))
+		for i, l := 0, rv.Len(); i < l; i++ {
+			v := rv.Index(i).Interface()
+			ridx[v] = append(ridx[v], i)
+		}
+
+		// For each row in t1, find the matching rows in col2
+		// and build up the row indexes for t1 and t2.
+		idx1, idx2 := []int{}, []int{}
+		lv := reflect.ValueOf(t1.MustColumn(col1))
+		for i, l := 0, lv.Len(); i < l; i++ {
+			r := ridx[lv.Index(i).Interface()]
+			for range r {
+				idx1 = append(idx1, i)
+			}
+			idx2 = append(idx2, r...)
+		}
+
+		// Build the joined table.
+		var nt Builder
+		for _, col := range t1.Columns() {
+			if cv, ok := t1.Const(col); ok {
+				nt.Add(col, cv)
+				continue
+			}
+			nt.Add(col, slice.Select(t1.Column(col), idx1))
+		}
+		for _, col := range t2.Columns() {
+			// Often the join column is the same in both
+			// and we can skip it because we added it from
+			// the first table.
+			if col == col1 && col == col2 {
+				continue
+			}
+
+			if cv, ok := t2.Const(col); ok {
+				nt.Add(col, cv)
+				continue
+			}
+			nt.Add(col, slice.Select(t2.Column(col), idx2))
+		}
+
+		ng.Add(gid, nt.Done())
+	}
+	return ng.Done()
+}
diff --git a/vendor/github.com/aclements/go-gg/table/map.go b/vendor/github.com/aclements/go-gg/table/map.go
new file mode 100644
index 0000000..87641a5
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/map.go
@@ -0,0 +1,170 @@
+// Copyright 2016 The Go 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 table
+
+import (
+	"fmt"
+	"reflect"
+
+	"github.com/aclements/go-gg/generic"
+)
+
+// MapTables applies f to each Table in g and returns a new Grouping
+// with the same group structure as g, but with the Tables returned by
+// f.
+func MapTables(g Grouping, f func(gid GroupID, table *Table) *Table) Grouping {
+	var out GroupingBuilder
+	for _, gid := range g.Tables() {
+		out.Add(gid, f(gid, g.Table(gid)))
+	}
+	return out.Done()
+}
+
+// MapCols applies f to a set of input columns to construct a set of
+// new output columns.
+//
+// For each Table in g, MapCols calls f(in[0], in[1], ..., out[0],
+// out[1], ...) where in[i] is column incols[i]. f should process the
+// values in the input column slices and fill output columns slices
+// out[j] accordingly. MapCols returns a new Grouping that adds each
+// outcols[j] bound to out[j].
+//
+// If all of the input columns are constant for a given table, MapCols
+// will call f with all slices of length 1. The input column slices
+// will contain the constant column values and MapCols will bind each
+// output column value out[i][0] as a constant.
+func MapCols(g Grouping, f interface{}, incols ...string) func(outcols ...string) Grouping {
+	return func(outcols ...string) Grouping {
+		fv := reflect.ValueOf(f)
+		if fv.Kind() != reflect.Func {
+			panic(&generic.TypeError{fv.Type(), nil, "must be a function"})
+		}
+		ft := fv.Type()
+		if ft.NumIn() != len(incols)+len(outcols) {
+			panic(&generic.TypeError{ft, nil, fmt.Sprintf("has the wrong number of arguments; expected %d", len(incols)+len(outcols))})
+		}
+		if ft.NumOut() != 0 {
+			panic(&generic.TypeError{ft, nil, "has the wrong number of results; expected 0"})
+		}
+
+		// Create output column slices.
+		totalRows := 0
+		for _, gid := range g.Tables() {
+			t := g.Table(gid)
+		colloop:
+			for _, incol := range incols {
+				if _, ok := t.Const(incol); !ok {
+					totalRows += g.Table(gid).Len()
+					break colloop
+				}
+			}
+		}
+		ocols := make([]reflect.Value, len(outcols))
+		for i := range ocols {
+			ocols[i] = reflect.MakeSlice(ft.In(i+len(incols)), totalRows, totalRows)
+		}
+
+		// Apply f to each group.
+		var out GroupingBuilder
+		args := make([]reflect.Value, len(incols)+len(outcols))
+		opos := 0
+		for _, gid := range g.Tables() {
+			t := g.Table(gid)
+
+			// Are all inputs are constants?
+			allConst := true
+			for _, incol := range incols {
+				if _, ok := t.Const(incol); !ok {
+					allConst = false
+					break
+				}
+			}
+			if allConst {
+				for i, incol := range incols {
+					cv, _ := t.Const(incol)
+					args[i] = reflect.MakeSlice(ColType(t, incol), 1, 1)
+					args[i].Index(0).Set(reflect.ValueOf(cv))
+				}
+				for i, ocol := range ocols {
+					args[i+len(incols)] = reflect.MakeSlice(ocol.Type(), 1, 1)
+				}
+
+				fv.Call(args)
+
+				tb := NewBuilder(t)
+				for i, outcol := range outcols {
+					tb.AddConst(outcol, args[i+len(incols)].Index(0).Interface())
+				}
+				out.Add(gid, tb.Done())
+				continue
+			}
+
+			// Prepare arguments.
+			for i, incol := range incols {
+				args[i] = reflect.ValueOf(t.MustColumn(incol))
+			}
+			for i, ocol := range ocols {
+				args[i+len(incols)] = ocol.Slice(opos, opos+t.Len())
+			}
+			opos += t.Len()
+
+			// Call f.
+			fv.Call(args)
+
+			// Add output columns.
+			tb := NewBuilder(t)
+			for i, outcol := range outcols {
+				tb.Add(outcol, args[i+len(incols)].Interface())
+			}
+			out.Add(gid, tb.Done())
+		}
+		return out.Done()
+	}
+}
+
+// Rename returns g with column 'from' renamed to 'to'. The column
+// retains its position.
+func Rename(g Grouping, from, to string) Grouping {
+	return MapTables(g, func(_ GroupID, t *Table) *Table {
+		t.MustColumn(from)
+		var nt Builder
+		for _, col := range t.Columns() {
+			if col == to {
+				continue
+			}
+
+			ncol := col
+			if col == from {
+				ncol = to
+			}
+
+			if cv, ok := t.Const(col); ok {
+				nt.AddConst(ncol, cv)
+			} else {
+				nt.Add(ncol, t.Column(col))
+			}
+		}
+		return nt.Done()
+	})
+}
+
+// Remove returns g with column 'col' removed.
+func Remove(g Grouping, col string) Grouping {
+	return MapTables(g, func(_ GroupID, t *Table) *Table {
+		t.MustColumn(col)
+		var nt Builder
+		for _, col2 := range t.Columns() {
+			if col == col2 {
+				continue
+			}
+			if cv, ok := t.Const(col2); ok {
+				nt.AddConst(col2, cv)
+			} else {
+				nt.Add(col2, t.Column(col2))
+			}
+		}
+		return nt.Done()
+	})
+}
diff --git a/vendor/github.com/aclements/go-gg/table/new.go b/vendor/github.com/aclements/go-gg/table/new.go
new file mode 100644
index 0000000..8094593
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/new.go
@@ -0,0 +1,105 @@
+// Copyright 2016 The Go 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 table
+
+import (
+	"reflect"
+	"strconv"
+
+	"github.com/aclements/go-gg/generic"
+)
+
+// TableFromStructs converts a []T where T is a struct to a Table
+// where the columns of the table correspond to T's exported fields.
+func TableFromStructs(structs Slice) *Table {
+	s := reflectSlice(structs)
+	st := s.Type()
+	if st.Elem().Kind() != reflect.Struct {
+		panic(&generic.TypeError{st, nil, "is not a slice of struct"})
+	}
+
+	var t Builder
+	rows := s.Len()
+	var rec func(reflect.Type, []int)
+	rec = func(typ reflect.Type, index []int) {
+		for fn := 0; fn < typ.NumField(); fn++ {
+			field := typ.Field(fn)
+			if field.PkgPath != "" {
+				continue
+			}
+			oldIndexLen := len(index)
+			index = append(index, field.Index...)
+			if field.Anonymous {
+				rec(field.Type, index)
+			} else {
+				col := reflect.MakeSlice(reflect.SliceOf(field.Type), rows, rows)
+				for i := 0; i < rows; i++ {
+					col.Index(i).Set(s.Index(i).FieldByIndex(index))
+				}
+				t.Add(field.Name, col.Interface())
+			}
+			index = index[:oldIndexLen]
+		}
+	}
+	rec(st.Elem(), []int{})
+	return t.Done()
+}
+
+// TableFromStrings converts a [][]string to a Table. This is intended
+// for processing external data, such as from CSV files. If coerce is
+// true, TableFromStrings will convert columns to []int or []float
+// when every string in that column is accepted by strconv.ParseInt or
+// strconv.ParseFloat, respectively.
+func TableFromStrings(cols []string, rows [][]string, coerce bool) *Table {
+	var t Builder
+	for i, col := range cols {
+		slice := make([]string, len(rows))
+		for j, row := range rows {
+			slice[j] = row[i]
+		}
+
+		var colData interface{} = slice
+		switch {
+		case coerce && len(slice) > 0:
+			// Try []int.
+			var err error
+			for _, str := range slice {
+				_, err = strconv.ParseInt(str, 10, 0)
+				if err != nil {
+					break
+				}
+			}
+			if err == nil {
+				nslice := make([]int, len(rows))
+				for i, str := range slice {
+					v, _ := strconv.ParseInt(str, 10, 0)
+					nslice[i] = int(v)
+				}
+				colData = nslice
+				break
+			}
+
+			// Try []float64. This must be done after
+			// []int. It's also more expensive.
+			for _, str := range slice {
+				_, err = strconv.ParseFloat(str, 64)
+				if err != nil {
+					break
+				}
+			}
+			if err == nil {
+				nslice := make([]float64, len(rows))
+				for i, str := range slice {
+					nslice[i], _ = strconv.ParseFloat(str, 64)
+				}
+				colData = nslice
+				break
+			}
+		}
+
+		t.Add(col, colData)
+	}
+	return t.Done()
+}
diff --git a/vendor/github.com/aclements/go-gg/table/pivot.go b/vendor/github.com/aclements/go-gg/table/pivot.go
new file mode 100644
index 0000000..6f488f9
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/pivot.go
@@ -0,0 +1,153 @@
+// Copyright 2016 The Go 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 table
+
+import (
+	"reflect"
+
+	"github.com/aclements/go-gg/generic"
+)
+
+// Pivot converts rows of g into columns. label and value must name
+// columns in g, and the label column must have type []string. Pivot
+// returns a Grouping with a new column named after each distinct
+// value in the label column, where the values in that column
+// correspond to the values from the value column. All other columns
+// (besides label and value) are copied to the output. If, for a given
+// column in an output row, no input row has that column in the label
+// column, the output cell will have the zero value for its type.
+func Pivot(g Grouping, label, value string) Grouping {
+	// Find all unique values of label. These are the new columns.
+	labels := []string{}
+	lset := map[string]int{}
+	for _, gid := range g.Tables() {
+		for _, l := range g.Table(gid).MustColumn(label).([]string) {
+			if _, ok := lset[l]; !ok {
+				lset[l] = len(lset)
+				labels = append(labels, l)
+			}
+		}
+	}
+
+	// Get all columns that are not label or value.
+	groupCols := []string{}
+	for _, col := range g.Columns() {
+		if col != label && col != value {
+			groupCols = append(groupCols, col)
+		}
+	}
+
+	return MapTables(g, func(_ GroupID, t *Table) *Table {
+		var nt Builder
+
+		// Group by all other columns. Each group in gg
+		// becomes an output row.
+		gg := GroupBy(t, groupCols...)
+
+		// Copy grouped-by values.
+		for _, groupCol := range groupCols {
+			cv := reflect.MakeSlice(reflect.TypeOf(t.Column(groupCol)), len(gg.Tables()), len(gg.Tables()))
+			for i, gid := range gg.Tables() {
+				sub := gg.Table(gid)
+				cv.Index(i).Set(reflect.ValueOf(sub.Column(groupCol)).Index(0))
+			}
+			nt.Add(groupCol, cv.Interface())
+		}
+
+		// Initialize new columns.
+		newCols := make([]reflect.Value, len(lset))
+		vt := reflect.TypeOf(t.MustColumn(value))
+		for i := range newCols {
+			newCols[i] = reflect.MakeSlice(vt, len(gg.Tables()), len(gg.Tables()))
+		}
+
+		// Fill in new columns.
+		for i, gid := range gg.Tables() {
+			sub := gg.Table(gid)
+
+			vcol := reflect.ValueOf(sub.MustColumn(value))
+			for j, l := range sub.MustColumn(label).([]string) {
+				val := vcol.Index(j)
+				newCols[lset[l]].Index(i).Set(val)
+			}
+		}
+
+		// Add new columns to output table.
+		for i, newCol := range newCols {
+			nt.Add(labels[i], newCol.Interface())
+		}
+
+		return nt.Done()
+	})
+}
+
+// Unpivot converts columns of g into rows. The returned Grouping
+// consists of the columns of g *not* listed in cols, plus two columns
+// named by the label and value arguments. For each input row in g,
+// the returned Grouping will have len(cols) output rows. The i'th
+// such output row corresponds to column cols[i] in the input row. The
+// label column will contain the name of the unpivoted column,
+// cols[i], and the value column will contain that column's value from
+// the input row. The values of all other columns in the input row
+// will be repeated across the output rows. All columns in cols must
+// have the same type.
+func Unpivot(g Grouping, label, value string, cols ...string) Grouping {
+	if len(cols) == 0 {
+		panic("Unpivot requires at least 1 column")
+	}
+
+	colSet := map[string]bool{}
+	for _, col := range cols {
+		colSet[col] = true
+	}
+
+	return MapTables(g, func(_ GroupID, t *Table) *Table {
+		var nt Builder
+
+		// Repeat all other columns len(cols) times.
+		ntlen := t.Len() * len(cols)
+		for _, name := range t.Columns() {
+			if colSet[name] || name == label || name == value {
+				continue
+			}
+
+			col := reflect.ValueOf(t.Column(name))
+			ncol := reflect.MakeSlice(col.Type(), ntlen, ntlen)
+			for i, l := 0, col.Len(); i < l; i++ {
+				v := col.Index(i)
+				for j := range cols {
+					ncol.Index(i*len(cols) + j).Set(v)
+				}
+			}
+
+			nt.Add(name, ncol.Interface())
+		}
+
+		// Get input columns.
+		var vt reflect.Type
+		colvs := make([]reflect.Value, len(cols))
+		for i, col := range cols {
+			colvs[i] = reflect.ValueOf(t.MustColumn(col))
+			if i == 0 {
+				vt = colvs[i].Type()
+			} else if vt != colvs[i].Type() {
+				panic(&generic.TypeError{vt, colvs[i].Type(), "; cannot Unpivot columns with different types"})
+			}
+		}
+
+		// Create label and value columns.
+		lcol := make([]string, 0, ntlen)
+		vcol := reflect.MakeSlice(vt, ntlen, ntlen)
+		for i := 0; i < t.Len(); i++ {
+			lcol = append(lcol, cols...)
+			for j, colv := range colvs {
+				vcol.Index(i*len(cols) + j).Set(colv.Index(i))
+			}
+		}
+		nt.Add(label, lcol).Add(value, vcol.Interface())
+
+		return nt.Done()
+	})
+}
diff --git a/vendor/github.com/aclements/go-gg/table/print.go b/vendor/github.com/aclements/go-gg/table/print.go
new file mode 100644
index 0000000..1795865
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/print.go
@@ -0,0 +1,119 @@
+// Copyright 2016 The Go 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 table
+
+import (
+	"fmt"
+	"io"
+	"os"
+	"reflect"
+	"strings"
+)
+
+// TODO: Have a format struct with options for things like column
+// separator, and header separator. Provide some defaults ones for,
+// e.g., Markdown, CSV, TSV, and such. Make top-level Print and Fprint
+// call methods in some default format.
+
+// Print(...) is shorthand for Fprint(os.Stderr, ...).
+func Print(g Grouping, formats ...string) error {
+	return Fprint(os.Stdout, g, formats...)
+}
+
+// Fprint prints Grouping g to w. formats[i] specifies a fmt-style
+// format string for column i. If there are more columns than formats,
+// remaining columns are formatted with %v (in particular, formats may
+// be omitted entirely to use %v for all columns). Numeric columns are
+// right aligned; all other column types are left aligned.
+func Fprint(w io.Writer, g Grouping, formats ...string) error {
+	if g.Columns() == nil {
+		return nil
+	}
+
+	// Convert each column to strings.
+	ss := make([][]string, len(g.Columns()))
+	rowFmts := make([]string, len(g.Columns()))
+	for i, col := range g.Columns() {
+		format := "%v"
+		if i < len(formats) {
+			format = formats[i]
+		}
+
+		// Format column.
+		var valKind reflect.Kind
+		ss[i] = []string{col}
+		for _, gid := range g.Tables() {
+			seq := reflect.ValueOf(g.Table(gid).Column(col))
+			for row := 0; row < seq.Len(); row++ {
+				str := fmt.Sprintf(format, seq.Index(row).Interface())
+				ss[i] = append(ss[i], str)
+			}
+
+			if valKind == reflect.Invalid {
+				valKind = seq.Type().Elem().Kind()
+			}
+		}
+
+		// Find column width.
+		width := 0
+		for _, s := range ss[i] {
+			if len(s) > width {
+				width = len(s)
+			}
+		}
+
+		// If it's a numeric column, right align.
+		//
+		// TODO: Even better would be to decimal align, though
+		// that may require some understanding of the format;
+		// or we could only do it for the default format.
+		switch valKind {
+		case reflect.Int, reflect.Int8, reflect.Int16, reflect.Int32, reflect.Int64, reflect.Uint, reflect.Uint8, reflect.Uint16, reflect.Uint32, reflect.Uint64, reflect.Uintptr, reflect.Float32, reflect.Float64:
+			width = -width
+		}
+
+		if i == len(g.Columns())-1 && width > 0 {
+			// Don't pad the last column.
+			rowFmts[i] = "%s"
+		} else {
+			rowFmts[i] = fmt.Sprintf("%%%ds", -width)
+		}
+	}
+
+	// Compute group headers.
+	groups := []GroupID{}
+	groupPos := []int{}
+	lastPos := 1
+	for _, gid := range g.Tables() {
+		groups = append(groups, gid)
+		groupPos = append(groupPos, lastPos)
+		lastPos += g.Table(gid).Len()
+	}
+	if len(groups) == 1 && groups[0] == RootGroupID {
+		groups, groupPos = nil, nil
+	}
+
+	// Print rows.
+	rowFmt := strings.Join(rowFmts, "  ") + "\n"
+	rowBuf := make([]interface{}, len(rowFmts))
+	for row := 0; row < len(ss[0]); row++ {
+		if len(groupPos) > 0 && row == groupPos[0] {
+			_, err := fmt.Fprintf(w, "-- %s\n", groups[0])
+			if err != nil {
+				return err
+			}
+			groups, groupPos = groups[1:], groupPos[1:]
+		}
+
+		for col := range rowBuf {
+			rowBuf[col] = ss[col][row]
+		}
+		_, err := fmt.Fprintf(w, rowFmt, rowBuf...)
+		if err != nil {
+			return err
+		}
+	}
+	return nil
+}
diff --git a/vendor/github.com/aclements/go-gg/table/sort.go b/vendor/github.com/aclements/go-gg/table/sort.go
new file mode 100644
index 0000000..ec64e00
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/sort.go
@@ -0,0 +1,93 @@
+// Copyright 2016 The Go 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 table
+
+import (
+	"sort"
+
+	"github.com/aclements/go-gg/generic/slice"
+)
+
+// SortBy sorts each group of g by the named columns. If a column's
+// type implements sort.Interface, rows will be sorted according to
+// that order. Otherwise, the values in the column must be naturally
+// ordered (their types must be orderable by the Go specification). If
+// neither is true, SortBy panics with a *generic.TypeError. If more
+// than one column is given, SortBy sorts by the tuple of the columns;
+// that is, if two values in the first column are equal, they are
+// sorted by the second column, and so on.
+func SortBy(g Grouping, cols ...string) Grouping {
+	// Sort each group.
+	sorters := make([]sort.Interface, len(cols))
+	return MapTables(g, func(_ GroupID, t *Table) *Table {
+		// Create sorters for each column.
+		sorters = sorters[:0]
+		for _, col := range cols {
+			if _, ok := t.Const(col); ok {
+				continue
+			}
+			seq := t.MustColumn(col)
+			sorter := slice.Sorter(seq)
+			if sort.IsSorted(sorter) {
+				continue
+			}
+			sorters = append(sorters, sorter)
+		}
+
+		if len(sorters) == 0 {
+			// Avoid shuffling everything by the identity
+			// permutation.
+			return t
+		}
+
+		// Generate an initial permutation sequence.
+		perm := make([]int, t.Len())
+		for i := range perm {
+			perm[i] = i
+		}
+
+		// Sort the permutation sequence.
+		sort.Stable(&permSort{perm, sorters})
+
+		// Permute all columns.
+		var nt Builder
+		for _, name := range t.Columns() {
+			if cv, ok := t.Const(name); ok {
+				nt.AddConst(name, cv)
+				continue
+			}
+			seq := t.Column(name)
+			seq = slice.Select(seq, perm)
+			nt.Add(name, seq)
+		}
+		return nt.Done()
+	})
+}
+
+type permSort struct {
+	perm []int
+	keys []sort.Interface
+}
+
+func (s *permSort) Len() int {
+	return len(s.perm)
+}
+
+func (s *permSort) Less(i, j int) bool {
+	// Since there's no way to ask about equality, we have to do
+	// extra work for all of the keys except the last.
+	for _, key := range s.keys[:len(s.keys)-1] {
+		if key.Less(s.perm[i], s.perm[j]) {
+			return true
+		} else if key.Less(s.perm[j], s.perm[i]) {
+			return false
+		}
+	}
+	return s.keys[len(s.keys)-1].Less(s.perm[i], s.perm[j])
+}
+
+func (s *permSort) Swap(i, j int) {
+	s.perm[i], s.perm[j] = s.perm[j], s.perm[i]
+}
diff --git a/vendor/github.com/aclements/go-gg/table/table.go b/vendor/github.com/aclements/go-gg/table/table.go
new file mode 100644
index 0000000..5821752
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/table.go
@@ -0,0 +1,489 @@
+// Copyright 2016 The Go 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 table implements ordered, grouped two dimensional relations.
+//
+// There are two related abstractions: Table and Grouping.
+//
+// A Table is an ordered relation of rows and columns. Each column is
+// a Go slice and hence must be homogeneously typed, but different
+// columns may have different types. All columns in a Table have the
+// same number of rows.
+//
+// A Grouping generalizes a Table by grouping the Table's rows into
+// zero or more groups. A Table is itself a Grouping with zero or one
+// groups. Most operations take a Grouping and operate on each group
+// independently, though some operations sub-divide or combine groups.
+//
+// The structures of both Tables and Groupings are immutable. They are
+// constructed using a Builder or a GroupingBuilder, respectively, and
+// then "frozen" into their respective immutable data structures.
+package table
+
+import (
+	"fmt"
+	"reflect"
+
+	"github.com/aclements/go-gg/generic"
+	"github.com/aclements/go-gg/generic/slice"
+)
+
+// TODO
+//
+// Rename Table to T?
+//
+// Make Table an interface? Then columns could be constructed lazily.
+//
+// Do all transformation functions as func(g Grouping) Grouping? That
+// could be a "Transform" type that has easy methods for chaining. In
+// a lot of cases, transformation functions could just return the
+// Transform returned by another function (like MapTables).
+//
+// Make an error type for "unknown column".
+
+// A Table is an immutable, ordered two dimensional relation. It
+// consists of a set of named columns. Each column is a sequence of
+// values of a consistent type or a constant value. All (non-constant)
+// columns have the same length.
+//
+// The zero value of Table is the "empty table": it has no rows and no
+// columns. Note that a Table may have one or more columns, but no
+// rows; such a Table is *not* considered empty.
+//
+// A Table is also a trivial Grouping. If a Table is empty, it has no
+// groups and hence the zero value of Table is also the "empty group".
+// Otherwise, it consists only of the root group, RootGroupID.
+type Table struct {
+	cols     map[string]Slice
+	consts   map[string]interface{}
+	colNames []string
+	len      int
+}
+
+// A Builder constructs a Table one column at a time.
+//
+// The zero value of a Builder represents an empty Table.
+type Builder struct {
+	t Table
+}
+
+// A Grouping is an immutable set of tables with identical sets of
+// columns, each identified by a distinct GroupID.
+//
+// Visually, a Grouping can be thought of as follows:
+//
+//	   Col A  Col B  Col C
+//	------ group /a ------
+//	0   5.4    "x"     90
+//	1   -.2    "y"     30
+//	------ group /b ------
+//	0   9.3    "a"     10
+//
+// Like a Table, a Grouping's structure is immutable. To construct a
+// Grouping, use a GroupingBuilder.
+//
+// Despite the fact that GroupIDs form a hierarchy, a Grouping ignores
+// this hierarchy and simply operates on a flat map of distinct
+// GroupIDs to Tables.
+type Grouping interface {
+	// Columns returns the names of the columns in this Grouping,
+	// or nil if there are no Tables or the group consists solely
+	// of empty Tables. All Tables in this Grouping have the same
+	// set of columns.
+	Columns() []string
+
+	// Tables returns the group IDs of the tables in this
+	// Grouping.
+	Tables() []GroupID
+
+	// Table returns the Table in group gid, or nil if there is no
+	// such Table.
+	Table(gid GroupID) *Table
+}
+
+// A GroupingBuilder constructs a Grouping one table a time.
+//
+// The zero value of a GroupingBuilder represents an empty Grouping
+// with no tables and no columns.
+type GroupingBuilder struct {
+	g        groupedTable
+	colTypes []reflect.Type
+}
+
+type groupedTable struct {
+	tables   map[GroupID]*Table
+	groups   []GroupID
+	colNames []string
+}
+
+// A Slice is a Go slice value.
+//
+// This is primarily for documentation. There is no way to statically
+// enforce this in Go; however, functions that expect a Slice will
+// panic with a *generic.TypeError if passed a non-slice value.
+type Slice interface{}
+
+func reflectSlice(s Slice) reflect.Value {
+	rv := reflect.ValueOf(s)
+	if rv.Kind() != reflect.Slice {
+		panic(&generic.TypeError{rv.Type(), nil, "is not a slice"})
+	}
+	return rv
+}
+
+// NewBuilder returns a new Builder. If t is non-nil, it populates the
+// new Builder with the columns from t.
+func NewBuilder(t *Table) *Builder {
+	if t == nil {
+		return new(Builder)
+	}
+	b := Builder{Table{
+		cols:     make(map[string]Slice),
+		consts:   make(map[string]interface{}),
+		colNames: append([]string(nil), t.Columns()...),
+		len:      t.len,
+	}}
+	for k, v := range t.cols {
+		b.t.cols[k] = v
+	}
+	for k, v := range t.consts {
+		b.t.consts[k] = v
+	}
+	return &b
+}
+
+// Add adds a column to b, or removes the named column if data is nil.
+// If b already has a column with the given name, Add replaces it. If
+// data is non-nil, it must have the same length as any existing
+// columns or Add will panic.
+func (b *Builder) Add(name string, data Slice) *Builder {
+	if data == nil {
+		// Remove the column.
+		if _, ok := b.t.cols[name]; !ok {
+			if _, ok := b.t.consts[name]; !ok {
+				// Nothing to remove.
+				return b
+			}
+		}
+		delete(b.t.cols, name)
+		delete(b.t.consts, name)
+		for i, n := range b.t.colNames {
+			if n == name {
+				copy(b.t.colNames[i:], b.t.colNames[i+1:])
+				b.t.colNames = b.t.colNames[:len(b.t.colNames)-1]
+				break
+			}
+		}
+		return b
+	}
+
+	// Are we replacing an existing column?
+	_, replace := b.t.cols[name]
+	if !replace {
+		_, replace = b.t.consts[name]
+	}
+
+	// Check the column and add it.
+	rv := reflectSlice(data)
+	dataLen := rv.Len()
+	if len(b.t.cols) == 0 || (replace && len(b.t.cols) == 1) {
+		if b.t.cols == nil {
+			b.t.cols = make(map[string]Slice)
+		}
+		// First non-constant column (possibly replacing the
+		// only non-constant column).
+		b.t.cols[name] = data
+		b.t.len = dataLen
+	} else if b.t.len != dataLen {
+		panic(fmt.Sprintf("cannot add column %q with %d elements to table with %d rows", name, dataLen, b.t.len))
+	} else {
+		b.t.cols[name] = data
+	}
+
+	if replace {
+		// Make sure it's not in constants.
+		delete(b.t.consts, name)
+	} else {
+		b.t.colNames = append(b.t.colNames, name)
+	}
+
+	return b
+}
+
+// AddConst adds a constant column to b whose value is val. If b
+// already has a column with this name, AddConst replaces it.
+//
+// A constant column has the same value in every row of the Table. It
+// does not itself have an inherent length.
+func (b *Builder) AddConst(name string, val interface{}) *Builder {
+	// Are we replacing an existing column?
+	_, replace := b.t.cols[name]
+	if !replace {
+		_, replace = b.t.consts[name]
+	}
+
+	if b.t.consts == nil {
+		b.t.consts = make(map[string]interface{})
+	}
+	b.t.consts[name] = val
+
+	if replace {
+		// Make sure it's not in cols.
+		delete(b.t.cols, name)
+	} else {
+		b.t.colNames = append(b.t.colNames, name)
+	}
+
+	return b
+}
+
+// Has returns true if b has a column named "name".
+func (b *Builder) Has(name string) bool {
+	if _, ok := b.t.cols[name]; ok {
+		return true
+	}
+	if _, ok := b.t.consts[name]; ok {
+		return true
+	}
+	return false
+}
+
+// Done returns the constructed Table and resets b.
+func (b *Builder) Done() *Table {
+	if len(b.t.colNames) == 0 {
+		return new(Table)
+	}
+	t := b.t
+	b.t = Table{}
+	return &t
+}
+
+// Len returns the number of rows in Table t.
+func (t *Table) Len() int {
+	return t.len
+}
+
+// Columns returns the names of the columns in Table t, or nil if this
+// Table is empty.
+func (t *Table) Columns() []string {
+	return t.colNames
+}
+
+// Column returns the slice of data in column name of Table t, or nil
+// if there is no such column. If name is a constant column, this
+// returns a slice with the constant value repeated to the length of
+// the Table.
+func (t *Table) Column(name string) Slice {
+	if c, ok := t.cols[name]; ok {
+		// It's a regular column or a constant column with a
+		// cached expansion.
+		return c
+	}
+
+	if cv, ok := t.consts[name]; ok {
+		// Expand the constant column and cache the result.
+		expanded := slice.Repeat(cv, t.len)
+		t.cols[name] = expanded
+		return expanded
+	}
+
+	return nil
+}
+
+// MustColumn is like Column, but panics if there is no such column.
+func (t *Table) MustColumn(name string) Slice {
+	if c := t.Column(name); c != nil {
+		return c
+	}
+	panic(fmt.Sprintf("unknown column %q", name))
+}
+
+// Const returns the value of constant column name. If this column
+// does not exist or is not a constant column, Const returns nil,
+// false.
+func (t *Table) Const(name string) (val interface{}, ok bool) {
+	cv, ok := t.consts[name]
+	return cv, ok
+}
+
+// isEmpty returns true if t is an empty Table, meaning it has no rows
+// or columns.
+func (t *Table) isEmpty() bool {
+	return t.colNames == nil
+}
+
+// Tables returns the groups IDs in this Table. If t is empty, there
+// are no group IDs. Otherwise, there is only RootGroupID.
+func (t *Table) Tables() []GroupID {
+	if t.isEmpty() {
+		return []GroupID{}
+	}
+	return []GroupID{RootGroupID}
+}
+
+// Table returns t if gid is RootGroupID and t is not empty; otherwise
+// it returns nil.
+func (t *Table) Table(gid GroupID) *Table {
+	if gid == RootGroupID && !t.isEmpty() {
+		return t
+	}
+	return nil
+}
+
+// NewGroupingBuilder returns a new GroupingBuilder. If g is non-nil,
+// it populates the new GroupingBuilder with the tables from g.
+func NewGroupingBuilder(g Grouping) *GroupingBuilder {
+	if g == nil {
+		return new(GroupingBuilder)
+	}
+	b := GroupingBuilder{groupedTable{
+		tables:   make(map[GroupID]*Table),
+		groups:   append([]GroupID(nil), g.Tables()...),
+		colNames: append([]string(nil), g.Columns()...),
+	}, nil}
+	for _, gid := range g.Tables() {
+		t := g.Table(gid)
+		b.g.tables[gid] = t
+		if b.colTypes == nil {
+			b.colTypes = colTypes(t)
+		}
+	}
+	return &b
+}
+
+func colTypes(t *Table) []reflect.Type {
+	colTypes := make([]reflect.Type, len(t.colNames))
+	for i, col := range t.colNames {
+		if c, ok := t.cols[col]; ok {
+			colTypes[i] = reflect.TypeOf(c).Elem()
+		} else {
+			colTypes[i] = reflect.TypeOf(t.consts[col])
+		}
+	}
+	return colTypes
+}
+
+// Add adds a Table to b, or removes a table if t is nil. If t is the
+// empty Table, this is a no-op because the empty Table contains no
+// groups. If gid already exists, Add replaces it. Table t must have
+// the same columns as any existing Tables in this Grouping and they
+// must have identical types; otherwise, Add will panic.
+//
+// TODO This doesn't make it easy to combine two Groupings. It could
+// instead take a Grouping and reparent it.
+func (b *GroupingBuilder) Add(gid GroupID, t *Table) *GroupingBuilder {
+	if t == nil {
+		if _, ok := b.g.tables[gid]; !ok {
+			// Nothing to remove.
+			return b
+		}
+		delete(b.g.tables, gid)
+		for i, g2 := range b.g.groups {
+			if g2 == gid {
+				copy(b.g.groups[i:], b.g.groups[i+1:])
+				b.g.groups = b.g.groups[:len(b.g.groups)-1]
+				break
+			}
+		}
+		return b
+	}
+
+	if t != nil && t.isEmpty() {
+		// Adding an empty table has no effect.
+		return b
+	}
+
+	if len(b.g.groups) == 1 && b.g.groups[0] == gid {
+		// We're replacing the only group. This is allowed to
+		// change the shape of the Grouping.
+		b.g.tables[gid] = t
+		b.g.colNames = t.Columns()
+		b.colTypes = colTypes(t)
+		return b
+	} else if len(b.g.groups) == 0 {
+		b.g.tables = map[GroupID]*Table{gid: t}
+		b.g.groups = []GroupID{gid}
+		b.g.colNames = t.Columns()
+		b.colTypes = colTypes(t)
+		return b
+	}
+
+	// Check that t's column names match.
+	matches := true
+	if len(t.colNames) != len(b.g.colNames) {
+		matches = false
+	} else {
+		for i, n := range t.colNames {
+			if b.g.colNames[i] != n {
+				matches = false
+				break
+			}
+		}
+	}
+	if !matches {
+		panic(fmt.Sprintf("table columns %q do not match group columns %q", t.colNames, b.g.colNames))
+	}
+
+	// Check that t's column types match.
+	for i, col := range b.g.colNames {
+		t0 := b.colTypes[i]
+		var t1 reflect.Type
+		if c, ok := t.cols[col]; ok {
+			t1 = reflect.TypeOf(c).Elem()
+		} else if cv, ok := t.consts[col]; ok {
+			t1 = reflect.TypeOf(cv)
+		}
+		if t0 != t1 {
+			panic(&generic.TypeError{t0, t1, fmt.Sprintf("for column %q are not the same", col)})
+		}
+	}
+
+	// Add t.
+	if _, ok := b.g.tables[gid]; !ok {
+		b.g.groups = append(b.g.groups, gid)
+	}
+	b.g.tables[gid] = t
+
+	return b
+}
+
+// Done returns the constructed Grouping and resets b.
+func (b *GroupingBuilder) Done() Grouping {
+	if len(b.g.groups) == 0 {
+		return new(groupedTable)
+	}
+	g := b.g
+	b.g = groupedTable{}
+	return &g
+}
+
+func (g *groupedTable) Columns() []string {
+	return g.colNames
+}
+
+func (g *groupedTable) Tables() []GroupID {
+	return g.groups
+}
+
+func (g *groupedTable) Table(gid GroupID) *Table {
+	return g.tables[gid]
+}
+
+// ColType returns the type of column col in g. This will always be a
+// slice type, even if col is a constant column. ColType panics if col
+// is unknown.
+//
+// TODO: If I introduce a first-class representation for a grouped
+// column, this should probably be in that.
+func ColType(g Grouping, col string) reflect.Type {
+	tables := g.Tables()
+	if len(tables) == 0 {
+		panic(fmt.Sprintf("unknown column %q", col))
+	}
+	t0 := g.Table(tables[0])
+	if cv, ok := t0.Const(col); ok {
+		return reflect.SliceOf(reflect.TypeOf(cv))
+	}
+	return reflect.TypeOf(t0.MustColumn(col))
+}
diff --git a/vendor/github.com/aclements/go-moremath/LICENSE b/vendor/github.com/aclements/go-moremath/LICENSE
new file mode 100644
index 0000000..d29b372
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/LICENSE
@@ -0,0 +1,27 @@
+Copyright (c) 2015 The Go Authors. All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+   * Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+   * Redistributions in binary form must reproduce the above
+copyright notice, this list of conditions and the following disclaimer
+in the documentation and/or other materials provided with the
+distribution.
+   * Neither the name of Google Inc. nor the names of its
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/vendor/github.com/aclements/go-moremath/mathx/beta.go b/vendor/github.com/aclements/go-moremath/mathx/beta.go
new file mode 100644
index 0000000..49f8722
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/mathx/beta.go
@@ -0,0 +1,93 @@
+// Copyright 2015 The Go 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 mathx
+
+import "math"
+
+func lgamma(x float64) float64 {
+	y, _ := math.Lgamma(x)
+	return y
+}
+
+// Beta returns the value of the complete beta function B(a, b).
+func Beta(a, b float64) float64 {
+	// B(x,y) = Γ(x)Γ(y) / Γ(x+y)
+	return math.Exp(lgamma(a) + lgamma(b) - lgamma(a+b))
+}
+
+// BetaInc returns the value of the regularized incomplete beta
+// function Iₓ(a, b).
+//
+// This is not to be confused with the "incomplete beta function",
+// which can be computed as BetaInc(x, a, b)*Beta(a, b).
+//
+// If x < 0 or x > 1, returns NaN.
+func BetaInc(x, a, b float64) float64 {
+	// Based on Numerical Recipes in C, section 6.4. This uses the
+	// continued fraction definition of I:
+	//
+	//  (xᵃ*(1-x)ᵇ)/(a*B(a,b)) * (1/(1+(d₁/(1+(d₂/(1+...))))))
+	//
+	// where B(a,b) is the beta function and
+	//
+	//  d_{2m+1} = -(a+m)(a+b+m)x/((a+2m)(a+2m+1))
+	//  d_{2m}   = m(b-m)x/((a+2m-1)(a+2m))
+	if x < 0 || x > 1 {
+		return math.NaN()
+	}
+	bt := 0.0
+	if 0 < x && x < 1 {
+		// Compute the coefficient before the continued
+		// fraction.
+		bt = math.Exp(lgamma(a+b) - lgamma(a) - lgamma(b) +
+			a*math.Log(x) + b*math.Log(1-x))
+	}
+	if x < (a+1)/(a+b+2) {
+		// Compute continued fraction directly.
+		return bt * betacf(x, a, b) / a
+	} else {
+		// Compute continued fraction after symmetry transform.
+		return 1 - bt*betacf(1-x, b, a)/b
+	}
+}
+
+// betacf is the continued fraction component of the regularized
+// incomplete beta function Iₓ(a, b).
+func betacf(x, a, b float64) float64 {
+	const maxIterations = 200
+	const epsilon = 3e-14
+
+	raiseZero := func(z float64) float64 {
+		if math.Abs(z) < math.SmallestNonzeroFloat64 {
+			return math.SmallestNonzeroFloat64
+		}
+		return z
+	}
+
+	c := 1.0
+	d := 1 / raiseZero(1-(a+b)*x/(a+1))
+	h := d
+	for m := 1; m <= maxIterations; m++ {
+		mf := float64(m)
+
+		// Even step of the recurrence.
+		numer := mf * (b - mf) * x / ((a + 2*mf - 1) * (a + 2*mf))
+		d = 1 / raiseZero(1+numer*d)
+		c = raiseZero(1 + numer/c)
+		h *= d * c
+
+		// Odd step of the recurrence.
+		numer = -(a + mf) * (a + b + mf) * x / ((a + 2*mf) * (a + 2*mf + 1))
+		d = 1 / raiseZero(1+numer*d)
+		c = raiseZero(1 + numer/c)
+		hfac := d * c
+		h *= hfac
+
+		if math.Abs(hfac-1) < epsilon {
+			return h
+		}
+	}
+	panic("betainc: a or b too big; failed to converge")
+}
diff --git a/vendor/github.com/aclements/go-moremath/mathx/choose.go b/vendor/github.com/aclements/go-moremath/mathx/choose.go
new file mode 100644
index 0000000..54dc27c
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/mathx/choose.go
@@ -0,0 +1,62 @@
+// Copyright 2015 The Go 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 mathx
+
+import "math"
+
+const smallFactLimit = 20 // 20! => 62 bits
+var smallFact [smallFactLimit + 1]int64
+
+func init() {
+	smallFact[0] = 1
+	fact := int64(1)
+	for n := int64(1); n <= smallFactLimit; n++ {
+		fact *= n
+		smallFact[n] = fact
+	}
+}
+
+// Choose returns the binomial coefficient of n and k.
+func Choose(n, k int) float64 {
+	if k == 0 || k == n {
+		return 1
+	}
+	if k < 0 || n < k {
+		return 0
+	}
+	if n <= smallFactLimit { // Implies k <= smallFactLimit
+		// It's faster to do several integer multiplications
+		// than it is to do an extra integer division.
+		// Remarkably, this is also faster than pre-computing
+		// Pascal's triangle (presumably because this is very
+		// cache efficient).
+		numer := int64(1)
+		for n1 := int64(n - (k - 1)); n1 <= int64(n); n1++ {
+			numer *= n1
+		}
+		denom := smallFact[k]
+		return float64(numer / denom)
+	}
+
+	return math.Exp(lchoose(n, k))
+}
+
+// Lchoose returns math.Log(Choose(n, k)).
+func Lchoose(n, k int) float64 {
+	if k == 0 || k == n {
+		return 0
+	}
+	if k < 0 || n < k {
+		return math.NaN()
+	}
+	return lchoose(n, k)
+}
+
+func lchoose(n, k int) float64 {
+	a, _ := math.Lgamma(float64(n + 1))
+	b, _ := math.Lgamma(float64(k + 1))
+	c, _ := math.Lgamma(float64(n - k + 1))
+	return a - b - c
+}
diff --git a/vendor/github.com/aclements/go-moremath/mathx/gamma.go b/vendor/github.com/aclements/go-moremath/mathx/gamma.go
new file mode 100644
index 0000000..d11096e
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/mathx/gamma.go
@@ -0,0 +1,96 @@
+// Copyright 2015 The Go 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 mathx
+
+import "math"
+
+// GammaInc returns the value of the incomplete gamma function (also
+// known as the regularized gamma function):
+//
+//   P(a, x) = 1 / Γ(a) * ∫₀ˣ exp(-t) t**(a-1) dt
+func GammaInc(a, x float64) float64 {
+	// Based on Numerical Recipes in C, section 6.2.
+
+	if a <= 0 || x < 0 || math.IsNaN(a) || math.IsNaN(x) {
+		return math.NaN()
+	}
+
+	if x < a+1 {
+		// Use the series representation, which converges more
+		// rapidly in this range.
+		return gammaIncSeries(a, x)
+	} else {
+		// Use the continued fraction representation.
+		return 1 - gammaIncCF(a, x)
+	}
+}
+
+// GammaIncComp returns the complement of the incomplete gamma
+// function 1 - GammaInc(a, x). This is more numerically stable for
+// values near 0.
+func GammaIncComp(a, x float64) float64 {
+	if a <= 0 || x < 0 || math.IsNaN(a) || math.IsNaN(x) {
+		return math.NaN()
+	}
+
+	if x < a+1 {
+		return 1 - gammaIncSeries(a, x)
+	} else {
+		return gammaIncCF(a, x)
+	}
+}
+
+func gammaIncSeries(a, x float64) float64 {
+	const maxIterations = 200
+	const epsilon = 3e-14
+
+	if x == 0 {
+		return 0
+	}
+
+	ap := a
+	del := 1 / a
+	sum := del
+	for n := 0; n < maxIterations; n++ {
+		ap++
+		del *= x / ap
+		sum += del
+		if math.Abs(del) < math.Abs(sum)*epsilon {
+			return sum * math.Exp(-x+a*math.Log(x)-lgamma(a))
+		}
+	}
+	panic("a too large; failed to converge")
+}
+
+func gammaIncCF(a, x float64) float64 {
+	const maxIterations = 200
+	const epsilon = 3e-14
+
+	raiseZero := func(z float64) float64 {
+		if math.Abs(z) < math.SmallestNonzeroFloat64 {
+			return math.SmallestNonzeroFloat64
+		}
+		return z
+	}
+
+	b := x + 1 - a
+	c := math.MaxFloat64
+	d := 1 / b
+	h := d
+
+	for i := 1; i <= maxIterations; i++ {
+		an := -float64(i) * (float64(i) - a)
+		b += 2
+		d = raiseZero(an*d + b)
+		c = raiseZero(b + an/c)
+		d = 1 / d
+		del := d * c
+		h *= del
+		if math.Abs(del-1) < epsilon {
+			return math.Exp(-x+a*math.Log(x)-lgamma(a)) * h
+		}
+	}
+	panic("a too large; failed to converge")
+}
diff --git a/vendor/github.com/aclements/go-moremath/mathx/package.go b/vendor/github.com/aclements/go-moremath/mathx/package.go
new file mode 100644
index 0000000..9d5de0d
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/mathx/package.go
@@ -0,0 +1,11 @@
+// Copyright 2015 The Go 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 mathx implements special functions not provided by the
+// standard math package.
+package mathx // import "github.com/aclements/go-moremath/mathx"
+
+import "math"
+
+var nan = math.NaN()
diff --git a/vendor/github.com/aclements/go-moremath/mathx/sign.go b/vendor/github.com/aclements/go-moremath/mathx/sign.go
new file mode 100644
index 0000000..372e92f
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/mathx/sign.go
@@ -0,0 +1,18 @@
+// Copyright 2015 The Go 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 mathx
+
+// Sign returns the sign of x: -1 if x < 0, 0 if x == 0, 1 if x > 0.
+// If x is NaN, it returns NaN.
+func Sign(x float64) float64 {
+	if x == 0 {
+		return 0
+	} else if x < 0 {
+		return -1
+	} else if x > 0 {
+		return 1
+	}
+	return nan
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/alg.go b/vendor/github.com/aclements/go-moremath/stats/alg.go
new file mode 100644
index 0000000..f704cef
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/alg.go
@@ -0,0 +1,114 @@
+// Copyright 2015 The Go 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 stats
+
+// Miscellaneous helper algorithms
+
+import (
+	"fmt"
+
+	"github.com/aclements/go-moremath/mathx"
+)
+
+func maxint(a, b int) int {
+	if a > b {
+		return a
+	}
+	return b
+}
+
+func minint(a, b int) int {
+	if a < b {
+		return a
+	}
+	return b
+}
+
+func sumint(xs []int) int {
+	sum := 0
+	for _, x := range xs {
+		sum += x
+	}
+	return sum
+}
+
+// bisect returns an x in [low, high] such that |f(x)| <= tolerance
+// using the bisection method.
+//
+// f(low) and f(high) must have opposite signs.
+//
+// If f does not have a root in this interval (e.g., it is
+// discontiguous), this returns the X of the apparent discontinuity
+// and false.
+func bisect(f func(float64) float64, low, high, tolerance float64) (float64, bool) {
+	flow, fhigh := f(low), f(high)
+	if -tolerance <= flow && flow <= tolerance {
+		return low, true
+	}
+	if -tolerance <= fhigh && fhigh <= tolerance {
+		return high, true
+	}
+	if mathx.Sign(flow) == mathx.Sign(fhigh) {
+		panic(fmt.Sprintf("root of f is not bracketed by [low, high]; f(%g)=%g f(%g)=%g", low, flow, high, fhigh))
+	}
+	for {
+		mid := (high + low) / 2
+		fmid := f(mid)
+		if -tolerance <= fmid && fmid <= tolerance {
+			return mid, true
+		}
+		if mid == high || mid == low {
+			return mid, false
+		}
+		if mathx.Sign(fmid) == mathx.Sign(flow) {
+			low = mid
+			flow = fmid
+		} else {
+			high = mid
+			fhigh = fmid
+		}
+	}
+}
+
+// bisectBool implements the bisection method on a boolean function.
+// It returns x1, x2 ∈ [low, high], x1 < x2 such that f(x1) != f(x2)
+// and x2 - x1 <= xtol.
+//
+// If f(low) == f(high), it panics.
+func bisectBool(f func(float64) bool, low, high, xtol float64) (x1, x2 float64) {
+	flow, fhigh := f(low), f(high)
+	if flow == fhigh {
+		panic(fmt.Sprintf("root of f is not bracketed by [low, high]; f(%g)=%v f(%g)=%v", low, flow, high, fhigh))
+	}
+	for {
+		if high-low <= xtol {
+			return low, high
+		}
+		mid := (high + low) / 2
+		if mid == high || mid == low {
+			return low, high
+		}
+		fmid := f(mid)
+		if fmid == flow {
+			low = mid
+			flow = fmid
+		} else {
+			high = mid
+			fhigh = fmid
+		}
+	}
+}
+
+// series returns the sum of the series f(0), f(1), ...
+//
+// This implementation is fast, but subject to round-off error.
+func series(f func(float64) float64) float64 {
+	y, yp := 0.0, 1.0
+	for n := 0.0; y != yp; n++ {
+		yp = y
+		y += f(n)
+	}
+	return y
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/deltadist.go b/vendor/github.com/aclements/go-moremath/stats/deltadist.go
new file mode 100644
index 0000000..bb3ba3f
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/deltadist.go
@@ -0,0 +1,57 @@
+// Copyright 2015 The Go 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 stats
+
+// DeltaDist is the Dirac delta function, centered at T, with total
+// area 1.
+//
+// The CDF of the Dirac delta function is the Heaviside step function,
+// centered at T. Specifically, f(T) == 1.
+type DeltaDist struct {
+	T float64
+}
+
+func (d DeltaDist) PDF(x float64) float64 {
+	if x == d.T {
+		return inf
+	}
+	return 0
+}
+
+func (d DeltaDist) pdfEach(xs []float64) []float64 {
+	res := make([]float64, len(xs))
+	for i, x := range xs {
+		if x == d.T {
+			res[i] = inf
+		}
+	}
+	return res
+}
+
+func (d DeltaDist) CDF(x float64) float64 {
+	if x >= d.T {
+		return 1
+	}
+	return 0
+}
+
+func (d DeltaDist) cdfEach(xs []float64) []float64 {
+	res := make([]float64, len(xs))
+	for i, x := range xs {
+		res[i] = d.CDF(x)
+	}
+	return res
+}
+
+func (d DeltaDist) InvCDF(y float64) float64 {
+	if y < 0 || y > 1 {
+		return nan
+	}
+	return d.T
+}
+
+func (d DeltaDist) Bounds() (float64, float64) {
+	return d.T - 1, d.T + 1
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/dist.go b/vendor/github.com/aclements/go-moremath/stats/dist.go
new file mode 100644
index 0000000..048477d
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/dist.go
@@ -0,0 +1,210 @@
+// Copyright 2015 The Go 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 stats
+
+import "math/rand"
+
+// A DistCommon is a statistical distribution. DistCommon is a base
+// interface provided by both continuous and discrete distributions.
+type DistCommon interface {
+	// CDF returns the cumulative probability Pr[X <= x].
+	//
+	// For continuous distributions, the CDF is the integral of
+	// the PDF from -inf to x.
+	//
+	// For discrete distributions, the CDF is the sum of the PMF
+	// at all defined points from -inf to x, inclusive. Note that
+	// the CDF of a discrete distribution is defined for the whole
+	// real line (unlike the PMF) but has discontinuities where
+	// the PMF is non-zero.
+	//
+	// The CDF is a monotonically increasing function and has a
+	// domain of all real numbers. If the distribution has bounded
+	// support, it has a range of [0, 1]; otherwise it has a range
+	// of (0, 1). Finally, CDF(-inf)==0 and CDF(inf)==1.
+	CDF(x float64) float64
+
+	// Bounds returns reasonable bounds for this distribution's
+	// PDF/PMF and CDF. The total weight outside of these bounds
+	// should be approximately 0.
+	//
+	// For a discrete distribution, both bounds are integer
+	// multiples of Step().
+	//
+	// If this distribution has finite support, it returns exact
+	// bounds l, h such that CDF(l')=0 for all l' < l and
+	// CDF(h')=1 for all h' >= h.
+	Bounds() (float64, float64)
+}
+
+// A Dist is a continuous statistical distribution.
+type Dist interface {
+	DistCommon
+
+	// PDF returns the value of the probability density function
+	// of this distribution at x.
+	PDF(x float64) float64
+}
+
+// A DiscreteDist is a discrete statistical distribution.
+//
+// Most discrete distributions are defined only at integral values of
+// the random variable. However, some are defined at other intervals,
+// so this interface takes a float64 value for the random variable.
+// The probability mass function rounds down to the nearest defined
+// point. Note that float64 values can exactly represent integer
+// values between ±2**53, so this generally shouldn't be an issue for
+// integer-valued distributions (likewise, for half-integer-valued
+// distributions, float64 can exactly represent all values between
+// ±2**52).
+type DiscreteDist interface {
+	DistCommon
+
+	// PMF returns the value of the probability mass function
+	// Pr[X = x'], where x' is x rounded down to the nearest
+	// defined point on the distribution.
+	//
+	// Note for implementers: for integer-valued distributions,
+	// round x using int(math.Floor(x)). Do not use int(x), since
+	// that truncates toward zero (unless all x <= 0 are handled
+	// the same).
+	PMF(x float64) float64
+
+	// Step returns s, where the distribution is defined for sℕ.
+	Step() float64
+}
+
+// TODO: Add a Support method for finite support distributions? Or
+// maybe just another return value from Bounds indicating that the
+// bounds are exact?
+
+// TODO: Plot method to return a pre-configured Plot object with
+// reasonable bounds and an integral function? Have to distinguish
+// PDF/CDF/InvCDF. Three methods? Argument?
+//
+// Doesn't have to be a method of Dist. Could be just a function that
+// takes a Dist and uses Bounds.
+
+// InvCDF returns the inverse CDF function of the given distribution
+// (also known as the quantile function or the percent point
+// function). This is a function f such that f(dist.CDF(x)) == x. If
+// dist.CDF is only weakly monotonic (that it, there are intervals
+// over which it is constant) and y > 0, f returns the smallest x that
+// satisfies this condition. In general, the inverse CDF is not
+// well-defined for y==0, but for convenience if y==0, f returns the
+// largest x that satisfies this condition. For distributions with
+// infinite support both the largest and smallest x are -Inf; however,
+// for distributions with finite support, this is the lower bound of
+// the support.
+//
+// If y < 0 or y > 1, f returns NaN.
+//
+// If dist implements InvCDF(float64) float64, this returns that
+// method. Otherwise, it returns a function that uses a generic
+// numerical method to construct the inverse CDF at y by finding x
+// such that dist.CDF(x) == y. This may have poor precision around
+// points of discontinuity, including f(0) and f(1).
+func InvCDF(dist DistCommon) func(y float64) (x float64) {
+	type invCDF interface {
+		InvCDF(float64) float64
+	}
+	if dist, ok := dist.(invCDF); ok {
+		return dist.InvCDF
+	}
+
+	// Otherwise, use a numerical algorithm.
+	//
+	// TODO: For discrete distributions, use the step size to
+	// inform this computation.
+	return func(y float64) (x float64) {
+		const almostInf = 1e100
+		const xtol = 1e-16
+
+		if y < 0 || y > 1 {
+			return nan
+		} else if y == 0 {
+			l, _ := dist.Bounds()
+			if dist.CDF(l) == 0 {
+				// Finite support
+				return l
+			} else {
+				// Infinite support
+				return -inf
+			}
+		} else if y == 1 {
+			_, h := dist.Bounds()
+			if dist.CDF(h) == 1 {
+				// Finite support
+				return h
+			} else {
+				// Infinite support
+				return inf
+			}
+		}
+
+		// Find loX, hiX for which cdf(loX) < y <= cdf(hiX).
+		var loX, loY, hiX, hiY float64
+		x1, y1 := 0.0, dist.CDF(0)
+		xdelta := 1.0
+		if y1 < y {
+			hiX, hiY = x1, y1
+			for hiY < y && hiX != inf {
+				loX, loY, hiX = hiX, hiY, hiX+xdelta
+				hiY = dist.CDF(hiX)
+				xdelta *= 2
+			}
+		} else {
+			loX, loY = x1, y1
+			for y <= loY && loX != -inf {
+				hiX, hiY, loX = loX, loY, loX-xdelta
+				loY = dist.CDF(loX)
+				xdelta *= 2
+			}
+		}
+		if loX == -inf {
+			return loX
+		} else if hiX == inf {
+			return hiX
+		}
+
+		// Use bisection on the interval to find the smallest
+		// x at which cdf(x) <= y.
+		_, x = bisectBool(func(x float64) bool {
+			return dist.CDF(x) < y
+		}, loX, hiX, xtol)
+		return
+	}
+}
+
+// Rand returns a random number generator that draws from the given
+// distribution. The returned generator takes an optional source of
+// randomness; if this is nil, it uses the default global source.
+//
+// If dist implements Rand(*rand.Rand) float64, Rand returns that
+// method. Otherwise, it returns a generic generator based on dist's
+// inverse CDF (which may in turn use an efficient implementation or a
+// generic numerical implementation; see InvCDF).
+func Rand(dist DistCommon) func(*rand.Rand) float64 {
+	type distRand interface {
+		Rand(*rand.Rand) float64
+	}
+	if dist, ok := dist.(distRand); ok {
+		return dist.Rand
+	}
+
+	// Otherwise, use a generic algorithm.
+	inv := InvCDF(dist)
+	return func(r *rand.Rand) float64 {
+		var y float64
+		for y == 0 {
+			if r == nil {
+				y = rand.Float64()
+			} else {
+				y = r.Float64()
+			}
+		}
+		return inv(y)
+	}
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/hist.go b/vendor/github.com/aclements/go-moremath/stats/hist.go
new file mode 100644
index 0000000..8578c07
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/hist.go
@@ -0,0 +1,66 @@
+// Copyright 2015 The Go 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 stats
+
+import "math"
+
+// TODO: Implement histograms on top of scales.
+
+type Histogram interface {
+	// Add adds a sample with value x to histogram h.
+	Add(x float64)
+
+	// Counts returns the number of samples less than the lowest
+	// bin, a slice of the number of samples in each bin,
+	// and the number of samples greater than the highest bin.
+	Counts() (under uint, counts []uint, over uint)
+
+	// BinToValue returns the value that would appear at the given
+	// bin index.
+	//
+	// For integral values of bin, BinToValue returns the lower
+	// bound of bin.  That is, a sample value x will be in bin if
+	// bin is integral and
+	//
+	//    BinToValue(bin) <= x < BinToValue(bin + 1)
+	//
+	// For non-integral values of bin, BinToValue interpolates
+	// between the lower and upper bounds of math.Floor(bin).
+	//
+	// BinToValue is undefined if bin > 1 + the number of bins.
+	BinToValue(bin float64) float64
+}
+
+// HistogramQuantile returns the x such that n*q samples in hist are
+// <= x, assuming values are distibuted within each bin according to
+// hist's distribution.
+//
+// If the q'th sample falls below the lowest bin or above the highest
+// bin, returns NaN.
+func HistogramQuantile(hist Histogram, q float64) float64 {
+	under, counts, over := hist.Counts()
+	total := under + over
+	for _, count := range counts {
+		total += count
+	}
+
+	goal := uint(float64(total) * q)
+	if goal <= under || goal > total-over {
+		return math.NaN()
+	}
+	for bin, count := range counts {
+		if count > goal {
+			return hist.BinToValue(float64(bin) + float64(goal)/float64(count))
+		}
+		goal -= count
+	}
+	panic("goal count not reached")
+}
+
+// HistogramIQR returns the interquartile range of the samples in
+// hist.
+func HistogramIQR(hist Histogram) float64 {
+	return HistogramQuantile(hist, 0.75) - HistogramQuantile(hist, 0.25)
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/hypergdist.go b/vendor/github.com/aclements/go-moremath/stats/hypergdist.go
new file mode 100644
index 0000000..1ea05e2
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/hypergdist.go
@@ -0,0 +1,101 @@
+// Copyright 2015 The Go 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 stats
+
+import (
+	"math"
+
+	"github.com/aclements/go-moremath/mathx"
+)
+
+// HypergeometicDist is a hypergeometric distribution.
+type HypergeometicDist struct {
+	// N is the size of the population. N >= 0.
+	N int
+
+	// K is the number of successes in the population. 0 <= K <= N.
+	K int
+
+	// Draws is the number of draws from the population. This is
+	// usually written "n", but is called Draws here because of
+	// limitations on Go identifier naming. 0 <= Draws <= N.
+	Draws int
+}
+
+// PMF is the probability of getting exactly int(k) successes in
+// d.Draws draws with replacement from a population of size d.N that
+// contains exactly d.K successes.
+func (d HypergeometicDist) PMF(k float64) float64 {
+	ki := int(math.Floor(k))
+	l, h := d.bounds()
+	if ki < l || ki > h {
+		return 0
+	}
+	return d.pmf(ki)
+}
+
+func (d HypergeometicDist) pmf(k int) float64 {
+	return math.Exp(mathx.Lchoose(d.K, k) + mathx.Lchoose(d.N-d.K, d.Draws-k) - mathx.Lchoose(d.N, d.Draws))
+}
+
+// CDF is the probability of getting int(k) or fewer successes in
+// d.Draws draws with replacement from a population of size d.N that
+// contains exactly d.K successes.
+func (d HypergeometicDist) CDF(k float64) float64 {
+	// Based on Klotz, A Computational Approach to Statistics.
+	ki := int(math.Floor(k))
+	l, h := d.bounds()
+	if ki < l {
+		return 0
+	} else if ki >= h {
+		return 1
+	}
+	// Use symmetry to compute the smaller sum.
+	flip := false
+	if ki > (d.Draws+1)/(d.N+1)*(d.K+1) {
+		flip = true
+		ki = d.K - ki - 1
+		d.Draws = d.N - d.Draws
+	}
+	p := d.pmf(ki) * d.sum(ki)
+	if flip {
+		p = 1 - p
+	}
+	return p
+}
+
+func (d HypergeometicDist) sum(k int) float64 {
+	const epsilon = 1e-14
+	sum, ak := 1.0, 1.0
+	L := maxint(0, d.Draws+d.K-d.N)
+	for dk := 1; dk <= k-L && ak/sum > epsilon; dk++ {
+		ak *= float64(1+k-dk) / float64(d.Draws-k+dk)
+		ak *= float64(d.N-d.K-d.Draws+k+1-dk) / float64(d.K-k+dk)
+		sum += ak
+	}
+	return sum
+}
+
+func (d HypergeometicDist) bounds() (int, int) {
+	return maxint(0, d.Draws+d.K-d.N), minint(d.Draws, d.K)
+}
+
+func (d HypergeometicDist) Bounds() (float64, float64) {
+	l, h := d.bounds()
+	return float64(l), float64(h)
+}
+
+func (d HypergeometicDist) Step() float64 {
+	return 1
+}
+
+func (d HypergeometicDist) Mean() float64 {
+	return float64(d.Draws*d.K) / float64(d.N)
+}
+
+func (d HypergeometicDist) Variance() float64 {
+	return float64(d.Draws*d.K*(d.N-d.K)*(d.N-d.Draws)) /
+		float64(d.N*d.N*(d.N-1))
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/kde.go b/vendor/github.com/aclements/go-moremath/stats/kde.go
new file mode 100644
index 0000000..6bcbd49
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/kde.go
@@ -0,0 +1,350 @@
+// Copyright 2015 The Go 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 stats
+
+import (
+	"fmt"
+	"math"
+)
+
+// A KDE is a distribution that estimates the underlying distribution
+// of a Sample using kernel density estimation.
+//
+// Kernel density estimation is a method for constructing an estimate
+// ƒ̂(x) of a unknown distribution ƒ(x) given a sample from that
+// distribution. Unlike many techniques, kernel density estimation is
+// non-parametric: in general, it doesn't assume any particular true
+// distribution (note, however, that the resulting distribution
+// depends deeply on the selected bandwidth, and many bandwidth
+// estimation techniques assume normal reference rules).
+//
+// A kernel density estimate is similar to a histogram, except that it
+// is a smooth probability estimate and does not require choosing a
+// bin size and discretizing the data.
+//
+// Sample is the only required field. All others have reasonable
+// defaults.
+type KDE struct {
+	// Sample is the data sample underlying this KDE.
+	Sample Sample
+
+	// Kernel is the kernel to use for the KDE.
+	Kernel KDEKernel
+
+	// Bandwidth is the bandwidth to use for the KDE.
+	//
+	// If this is zero, the bandwidth is computed from the
+	// provided data using a default bandwidth estimator
+	// (currently BandwidthScott).
+	Bandwidth float64
+
+	// BoundaryMethod is the boundary correction method to use for
+	// the KDE. The default value is BoundaryReflect; however, the
+	// default bounds are effectively +/-inf, which is equivalent
+	// to performing no boundary correction.
+	BoundaryMethod KDEBoundaryMethod
+
+	// [BoundaryMin, BoundaryMax) specify a bounded support for
+	// the KDE. If both are 0 (their default values), they are
+	// treated as +/-inf.
+	//
+	// To specify a half-bounded support, set Min to math.Inf(-1)
+	// or Max to math.Inf(1).
+	BoundaryMin float64
+	BoundaryMax float64
+}
+
+// BandwidthSilverman is a bandwidth estimator implementing
+// Silverman's Rule of Thumb. It's fast, but not very robust to
+// outliers as it assumes data is approximately normal.
+//
+// Silverman, B. W. (1986) Density Estimation.
+func BandwidthSilverman(data interface {
+	StdDev() float64
+	Weight() float64
+}) float64 {
+	return 1.06 * data.StdDev() * math.Pow(data.Weight(), -1.0/5)
+}
+
+// BandwidthScott is a bandwidth estimator implementing Scott's Rule.
+// This is generally robust to outliers: it chooses the minimum
+// between the sample's standard deviation and an robust estimator of
+// a Gaussian distribution's standard deviation.
+//
+// Scott, D. W. (1992) Multivariate Density Estimation: Theory,
+// Practice, and Visualization.
+func BandwidthScott(data interface {
+	StdDev() float64
+	Weight() float64
+	Quantile(float64) float64
+}) float64 {
+	iqr := data.Quantile(0.75) - data.Quantile(0.25)
+	hScale := 1.06 * math.Pow(data.Weight(), -1.0/5)
+	stdDev := data.StdDev()
+	if stdDev < iqr/1.349 {
+		// Use Silverman's Rule of Thumb
+		return hScale * stdDev
+	} else {
+		// Use IQR/1.349 as a robust estimator of the standard
+		// deviation of a Gaussian distribution.
+		return hScale * (iqr / 1.349)
+	}
+}
+
+// TODO(austin) Implement bandwidth estimator from Botev, Grotowski,
+// Kroese. (2010) Kernel Density Estimation via Diffusion.
+
+// KDEKernel represents a kernel to use for a KDE.
+type KDEKernel int
+
+//go:generate stringer -type=KDEKernel
+
+const (
+	// An EpanechnikovKernel is a smooth kernel with bounded
+	// support. As a result, the KDE will also have bounded
+	// support. It is "optimal" in the sense that it minimizes the
+	// asymptotic mean integrated squared error (AMISE).
+	EpanechnikovKernel KDEKernel = iota
+
+	// A GaussianKernel is a Gaussian (normal) kernel.
+	GaussianKernel
+
+	// A DeltaKernel is a Dirac delta function. The PDF of such a
+	// KDE is not well-defined, but the CDF will represent each
+	// sample as an instantaneous increase. This kernel ignores
+	// bandwidth and never requires boundary correction.
+	DeltaKernel
+)
+
+// KDEBoundaryMethod represents a boundary correction method for
+// constructing a KDE with bounded support.
+type KDEBoundaryMethod int
+
+//go:generate stringer -type=KDEBoundaryMethod
+
+const (
+	// BoundaryReflect reflects the density estimate at the
+	// boundaries.  For example, for a KDE with support [0, inf),
+	// this is equivalent to ƒ̂ᵣ(x)=ƒ̂(x)+ƒ̂(-x) for x>=0.  This is a
+	// simple and fast technique, but enforces that ƒ̂ᵣ'(0)=0, so
+	// it may not be applicable to all distributions.
+	BoundaryReflect KDEBoundaryMethod = iota
+)
+
+type kdeKernel interface {
+	pdfEach(xs []float64) []float64
+	cdfEach(xs []float64) []float64
+}
+
+func (k *KDE) prepare() (kdeKernel, bool) {
+	// Compute bandwidth.
+	if k.Bandwidth == 0 {
+		k.Bandwidth = BandwidthScott(k.Sample)
+	}
+
+	// Construct kernel.
+	kernel := kdeKernel(nil)
+	switch k.Kernel {
+	default:
+		panic(fmt.Sprint("unknown kernel", k))
+	case EpanechnikovKernel:
+		kernel = epanechnikovKernel{k.Bandwidth}
+	case GaussianKernel:
+		kernel = NormalDist{0, k.Bandwidth}
+	case DeltaKernel:
+		kernel = DeltaDist{0}
+	}
+
+	// Use boundary correction?
+	bc := k.BoundaryMin != 0 || k.BoundaryMax != 0
+
+	return kernel, bc
+}
+
+// TODO: For KDEs of histograms, make histograms able to create a
+// weighted Sample and simply require the caller to provide a
+// good bandwidth from a StreamStats.
+
+// normalizedXs returns x - kde.Sample.Xs. Evaluating kernels shifted
+// by kde.Sample.Xs all at x is equivalent to evaluating one unshifted
+// kernel at x - kde.Sample.Xs.
+func (kde *KDE) normalizedXs(x float64) []float64 {
+	txs := make([]float64, len(kde.Sample.Xs))
+	for i, xi := range kde.Sample.Xs {
+		txs[i] = x - xi
+	}
+	return txs
+}
+
+func (kde *KDE) PDF(x float64) float64 {
+	kernel, bc := kde.prepare()
+
+	// Apply boundary
+	if bc && (x < kde.BoundaryMin || x >= kde.BoundaryMax) {
+		return 0
+	}
+
+	y := func(x float64) float64 {
+		// Shift kernel to each of kde.xs and evaluate at x
+		ys := kernel.pdfEach(kde.normalizedXs(x))
+
+		// Kernel samples are weighted according to the weights of xs
+		wys := Sample{Xs: ys, Weights: kde.Sample.Weights}
+
+		return wys.Sum() / wys.Weight()
+	}
+	if !bc {
+		return y(x)
+	}
+	switch kde.BoundaryMethod {
+	default:
+		panic("unknown boundary correction method")
+	case BoundaryReflect:
+		if math.IsInf(kde.BoundaryMax, 1) {
+			return y(x) + y(2*kde.BoundaryMin-x)
+		} else if math.IsInf(kde.BoundaryMin, -1) {
+			return y(x) + y(2*kde.BoundaryMax-x)
+		} else {
+			d := 2 * (kde.BoundaryMax - kde.BoundaryMin)
+			w := 2 * (x - kde.BoundaryMin)
+			return series(func(n float64) float64 {
+				// Points >= x
+				return y(x+n*d) + y(x+n*d-w)
+			}) + series(func(n float64) float64 {
+				// Points < x
+				return y(x-(n+1)*d+w) + y(x-(n+1)*d)
+			})
+		}
+	}
+}
+
+func (kde *KDE) CDF(x float64) float64 {
+	kernel, bc := kde.prepare()
+
+	// Apply boundary
+	if bc {
+		if x < kde.BoundaryMin {
+			return 0
+		} else if x >= kde.BoundaryMax {
+			return 1
+		}
+	}
+
+	y := func(x float64) float64 {
+		// Shift kernel integral to each of cdf.xs and evaluate at x
+		ys := kernel.cdfEach(kde.normalizedXs(x))
+
+		// Kernel samples are weighted according to the weights of xs
+		wys := Sample{Xs: ys, Weights: kde.Sample.Weights}
+
+		return wys.Sum() / wys.Weight()
+	}
+	if !bc {
+		return y(x)
+	}
+	switch kde.BoundaryMethod {
+	default:
+		panic("unknown boundary correction method")
+	case BoundaryReflect:
+		if math.IsInf(kde.BoundaryMax, 1) {
+			return y(x) - y(2*kde.BoundaryMin-x)
+		} else if math.IsInf(kde.BoundaryMin, -1) {
+			return y(x) + (1 - y(2*kde.BoundaryMax-x))
+		} else {
+			d := 2 * (kde.BoundaryMax - kde.BoundaryMin)
+			w := 2 * (x - kde.BoundaryMin)
+			return series(func(n float64) float64 {
+				// Windows >= x-w
+				return y(x+n*d) - y(x+n*d-w)
+			}) + series(func(n float64) float64 {
+				// Windows < x-w
+				return y(x-(n+1)*d) - y(x-(n+1)*d-w)
+			})
+		}
+	}
+}
+
+func (kde *KDE) Bounds() (low float64, high float64) {
+	_, bc := kde.prepare()
+
+	// TODO(austin) If this KDE came from a histogram, we'd better
+	// not sample at a significantly higher rate than the
+	// histogram.  Maybe we want to just return the bounds of the
+	// histogram?
+
+	// TODO(austin) It would be nice if this could be instructed
+	// to include all original data points, even if they are in
+	// the tail.  Probably that should just be up to the caller to
+	// pass an axis derived from the bounds of the original data.
+
+	// Use the lowest and highest samples as starting points
+	lowX, highX := kde.Sample.Bounds()
+	if lowX == highX {
+		lowX -= 1
+		highX += 1
+	}
+
+	// Find the end points that contain 99% of the CDF's weight.
+	// Since bisect requires that the root be bracketed, start by
+	// expanding our range if necessary.  TODO(austin) This can
+	// definitely be done faster.
+	const (
+		lowY      = 0.005
+		highY     = 0.995
+		tolerance = 0.001
+	)
+	for kde.CDF(lowX) > lowY {
+		lowX -= highX - lowX
+	}
+	for kde.CDF(highX) < highY {
+		highX += highX - lowX
+	}
+	// Explicitly accept discontinuities, since we may be using a
+	// discontiguous kernel.
+	low, _ = bisect(func(x float64) float64 { return kde.CDF(x) - lowY }, lowX, highX, tolerance)
+	high, _ = bisect(func(x float64) float64 { return kde.CDF(x) - highY }, lowX, highX, tolerance)
+
+	// Expand width by 20% to give some margins
+	width := high - low
+	low, high = low-0.1*width, high+0.1*width
+
+	// Limit to bounds
+	if bc {
+		low = math.Max(low, kde.BoundaryMin)
+		high = math.Min(high, kde.BoundaryMax)
+	}
+
+	return
+}
+
+type epanechnikovKernel struct {
+	h float64
+}
+
+func (d epanechnikovKernel) pdfEach(xs []float64) []float64 {
+	ys := make([]float64, len(xs))
+	a := 0.75 / d.h
+	invhh := 1 / (d.h * d.h)
+	for i, x := range xs {
+		if -d.h < x && x < d.h {
+			ys[i] = a * (1 - x*x*invhh)
+		}
+	}
+	return ys
+}
+
+func (d epanechnikovKernel) cdfEach(xs []float64) []float64 {
+	ys := make([]float64, len(xs))
+	invh := 1 / d.h
+	for i, x := range xs {
+		if x > d.h {
+			ys[i] = 1
+		} else if x > -d.h {
+			u := x * invh
+			ys[i] = 0.25 * (2 + 3*u - u*u*u)
+		}
+	}
+	return ys
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/kdeboundarymethod_string.go b/vendor/github.com/aclements/go-moremath/stats/kdeboundarymethod_string.go
new file mode 100644
index 0000000..e01d0e7
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/kdeboundarymethod_string.go
@@ -0,0 +1,16 @@
+// generated by stringer -type=KDEBoundaryMethod; DO NOT EDIT
+
+package stats
+
+import "fmt"
+
+const _KDEBoundaryMethod_name = "BoundaryReflect"
+
+var _KDEBoundaryMethod_index = [...]uint8{0, 15}
+
+func (i KDEBoundaryMethod) String() string {
+	if i < 0 || i+1 >= KDEBoundaryMethod(len(_KDEBoundaryMethod_index)) {
+		return fmt.Sprintf("KDEBoundaryMethod(%d)", i)
+	}
+	return _KDEBoundaryMethod_name[_KDEBoundaryMethod_index[i]:_KDEBoundaryMethod_index[i+1]]
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/kdekernel_string.go b/vendor/github.com/aclements/go-moremath/stats/kdekernel_string.go
new file mode 100644
index 0000000..20a23fd
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/kdekernel_string.go
@@ -0,0 +1,16 @@
+// generated by stringer -type=KDEKernel; DO NOT EDIT
+
+package stats
+
+import "fmt"
+
+const _KDEKernel_name = "GaussianKernelDeltaKernel"
+
+var _KDEKernel_index = [...]uint8{0, 14, 25}
+
+func (i KDEKernel) String() string {
+	if i < 0 || i+1 >= KDEKernel(len(_KDEKernel_index)) {
+		return fmt.Sprintf("KDEKernel(%d)", i)
+	}
+	return _KDEKernel_name[_KDEKernel_index[i]:_KDEKernel_index[i+1]]
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/linearhist.go b/vendor/github.com/aclements/go-moremath/stats/linearhist.go
new file mode 100644
index 0000000..c36335a
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/linearhist.go
@@ -0,0 +1,42 @@
+// Copyright 2015 The Go 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 stats
+
+// LinearHist is a Histogram with uniformly-sized bins.
+type LinearHist struct {
+	min, max, delta float64
+	low, high       uint
+	bins            []uint
+}
+
+// NewLinearHist returns an empty histogram with nbins uniformly-sized
+// bins spanning [min, max].
+func NewLinearHist(min, max float64, nbins int) *LinearHist {
+	delta := float64(nbins) / (max - min)
+	return &LinearHist{min, max, delta, 0, 0, make([]uint, nbins)}
+}
+
+func (h *LinearHist) bin(x float64) int {
+	return int(h.delta * (x - h.min))
+}
+
+func (h *LinearHist) Add(x float64) {
+	bin := h.bin(x)
+	if bin < 0 {
+		h.low++
+	} else if bin >= len(h.bins) {
+		h.high++
+	} else {
+		h.bins[bin]++
+	}
+}
+
+func (h *LinearHist) Counts() (uint, []uint, uint) {
+	return h.low, h.bins, h.high
+}
+
+func (h *LinearHist) BinToValue(bin float64) float64 {
+	return h.min + bin*h.delta
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/locationhypothesis_string.go b/vendor/github.com/aclements/go-moremath/stats/locationhypothesis_string.go
new file mode 100644
index 0000000..ab0f26c
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/locationhypothesis_string.go
@@ -0,0 +1,17 @@
+// generated by stringer -type LocationHypothesis; DO NOT EDIT
+
+package stats
+
+import "fmt"
+
+const _LocationHypothesis_name = "LocationLessLocationDiffersLocationGreater"
+
+var _LocationHypothesis_index = [...]uint8{0, 12, 27, 42}
+
+func (i LocationHypothesis) String() string {
+	i -= -1
+	if i < 0 || i+1 >= LocationHypothesis(len(_LocationHypothesis_index)) {
+		return fmt.Sprintf("LocationHypothesis(%d)", i+-1)
+	}
+	return _LocationHypothesis_name[_LocationHypothesis_index[i]:_LocationHypothesis_index[i+1]]
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/loghist.go b/vendor/github.com/aclements/go-moremath/stats/loghist.go
new file mode 100644
index 0000000..937f62a
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/loghist.go
@@ -0,0 +1,83 @@
+// Copyright 2015 The Go 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 stats
+
+import "math"
+
+// LogHist is a Histogram with logarithmically-spaced bins.
+type LogHist struct {
+	b         int
+	m         float64
+	mOverLogb float64
+	low, high uint
+	bins      []uint
+}
+
+// NewLogHist returns an empty logarithmic histogram with bins for
+// integral values of m * log_b(x) up to x = max.
+func NewLogHist(b int, m float64, max float64) *LogHist {
+	// TODO(austin) Minimum value as well?  If the samples are
+	// actually integral, having fractional bin boundaries can
+	// mess up smoothing.
+	mOverLogb := m / math.Log(float64(b))
+	nbins := int(math.Ceil(mOverLogb * math.Log(max)))
+	return &LogHist{b: b, m: m, mOverLogb: mOverLogb, low: 0, high: 0, bins: make([]uint, nbins)}
+}
+
+func (h *LogHist) bin(x float64) int {
+	return int(h.mOverLogb * math.Log(x))
+}
+
+func (h *LogHist) Add(x float64) {
+	bin := h.bin(x)
+	if bin < 0 {
+		h.low++
+	} else if bin >= len(h.bins) {
+		h.high++
+	} else {
+		h.bins[bin]++
+	}
+}
+
+func (h *LogHist) Counts() (uint, []uint, uint) {
+	return h.low, h.bins, h.high
+}
+
+func (h *LogHist) BinToValue(bin float64) float64 {
+	return math.Pow(float64(h.b), bin/h.m)
+}
+
+func (h *LogHist) At(x float64) float64 {
+	bin := h.bin(x)
+	if bin < 0 || bin >= len(h.bins) {
+		return 0
+	}
+	return float64(h.bins[bin])
+}
+
+func (h *LogHist) Bounds() (float64, float64) {
+	// XXX Plot will plot this on a linear axis.  Maybe this
+	// should be able to return the natural axis?
+	// Maybe then we could also give it the bins for the tics.
+	lowbin := 0
+	if h.low == 0 {
+		for bin, count := range h.bins {
+			if count > 0 {
+				lowbin = bin
+				break
+			}
+		}
+	}
+	highbin := len(h.bins)
+	if h.high == 0 {
+		for bin := range h.bins {
+			if h.bins[len(h.bins)-bin-1] > 0 {
+				highbin = len(h.bins) - bin
+				break
+			}
+		}
+	}
+	return h.BinToValue(float64(lowbin)), h.BinToValue(float64(highbin))
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/normaldist.go b/vendor/github.com/aclements/go-moremath/stats/normaldist.go
new file mode 100644
index 0000000..d00f96a
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/normaldist.go
@@ -0,0 +1,141 @@
+// Copyright 2015 The Go 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 stats
+
+import (
+	"math"
+	"math/rand"
+)
+
+// NormalDist is a normal (Gaussian) distribution with mean Mu and
+// standard deviation Sigma.
+type NormalDist struct {
+	Mu, Sigma float64
+}
+
+// StdNormal is the standard normal distribution (Mu = 0, Sigma = 1)
+var StdNormal = NormalDist{0, 1}
+
+// 1/sqrt(2 * pi)
+const invSqrt2Pi = 0.39894228040143267793994605993438186847585863116493465766592583
+
+func (n NormalDist) PDF(x float64) float64 {
+	z := x - n.Mu
+	return math.Exp(-z*z/(2*n.Sigma*n.Sigma)) * invSqrt2Pi / n.Sigma
+}
+
+func (n NormalDist) pdfEach(xs []float64) []float64 {
+	res := make([]float64, len(xs))
+	if n.Mu == 0 && n.Sigma == 1 {
+		// Standard normal fast path
+		for i, x := range xs {
+			res[i] = math.Exp(-x*x/2) * invSqrt2Pi
+		}
+	} else {
+		a := -1 / (2 * n.Sigma * n.Sigma)
+		b := invSqrt2Pi / n.Sigma
+		for i, x := range xs {
+			z := x - n.Mu
+			res[i] = math.Exp(z*z*a) * b
+		}
+	}
+	return res
+}
+
+func (n NormalDist) CDF(x float64) float64 {
+	return math.Erfc(-(x-n.Mu)/(n.Sigma*math.Sqrt2)) / 2
+}
+
+func (n NormalDist) cdfEach(xs []float64) []float64 {
+	res := make([]float64, len(xs))
+	a := 1 / (n.Sigma * math.Sqrt2)
+	for i, x := range xs {
+		res[i] = math.Erfc(-(x-n.Mu)*a) / 2
+	}
+	return res
+}
+
+func (n NormalDist) InvCDF(p float64) (x float64) {
+	// This is based on Peter John Acklam's inverse normal CDF
+	// algorithm: http://home.online.no/~pjacklam/notes/invnorm/
+	const (
+		a1 = -3.969683028665376e+01
+		a2 = 2.209460984245205e+02
+		a3 = -2.759285104469687e+02
+		a4 = 1.383577518672690e+02
+		a5 = -3.066479806614716e+01
+		a6 = 2.506628277459239e+00
+
+		b1 = -5.447609879822406e+01
+		b2 = 1.615858368580409e+02
+		b3 = -1.556989798598866e+02
+		b4 = 6.680131188771972e+01
+		b5 = -1.328068155288572e+01
+
+		c1 = -7.784894002430293e-03
+		c2 = -3.223964580411365e-01
+		c3 = -2.400758277161838e+00
+		c4 = -2.549732539343734e+00
+		c5 = 4.374664141464968e+00
+		c6 = 2.938163982698783e+00
+
+		d1 = 7.784695709041462e-03
+		d2 = 3.224671290700398e-01
+		d3 = 2.445134137142996e+00
+		d4 = 3.754408661907416e+00
+
+		plow  = 0.02425
+		phigh = 1 - plow
+	)
+
+	if p < 0 || p > 1 {
+		return nan
+	} else if p == 0 {
+		return -inf
+	} else if p == 1 {
+		return inf
+	}
+
+	if p < plow {
+		// Rational approximation for lower region.
+		q := math.Sqrt(-2 * math.Log(p))
+		x = (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q + c6) /
+			((((d1*q+d2)*q+d3)*q+d4)*q + 1)
+	} else if phigh < p {
+		// Rational approximation for upper region.
+		q := math.Sqrt(-2 * math.Log(1-p))
+		x = -(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q + c6) /
+			((((d1*q+d2)*q+d3)*q+d4)*q + 1)
+	} else {
+		// Rational approximation for central region.
+		q := p - 0.5
+		r := q * q
+		x = (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r + a6) * q /
+			(((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r + 1)
+	}
+
+	// Refine approximation.
+	e := 0.5*math.Erfc(-x/math.Sqrt2) - p
+	u := e * math.Sqrt(2*math.Pi) * math.Exp(x*x/2)
+	x = x - u/(1+x*u/2)
+
+	// Adjust from standard normal.
+	return x*n.Sigma + n.Mu
+}
+
+func (n NormalDist) Rand(r *rand.Rand) float64 {
+	var x float64
+	if r == nil {
+		x = rand.NormFloat64()
+	} else {
+		x = r.NormFloat64()
+	}
+	return x*n.Sigma + n.Mu
+}
+
+func (n NormalDist) Bounds() (float64, float64) {
+	const stddevs = 3
+	return n.Mu - stddevs*n.Sigma, n.Mu + stddevs*n.Sigma
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/package.go b/vendor/github.com/aclements/go-moremath/stats/package.go
new file mode 100644
index 0000000..644b399
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/package.go
@@ -0,0 +1,25 @@
+// Copyright 2015 The Go 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 stats implements several statistical distributions,
+// hypothesis tests, and functions for descriptive statistics.
+//
+// Currently stats is fairly small, but for what it does implement, it
+// focuses on high quality, fast implementations with good, idiomatic
+// Go APIs.
+package stats // import "github.com/aclements/go-moremath/stats"
+
+import (
+	"errors"
+	"math"
+)
+
+var inf = math.Inf(1)
+var nan = math.NaN()
+
+// TODO: Put all errors in the same place and maybe unify them.
+
+var (
+	ErrSamplesEqual = errors.New("all samples are equal")
+)
diff --git a/vendor/github.com/aclements/go-moremath/stats/sample.go b/vendor/github.com/aclements/go-moremath/stats/sample.go
new file mode 100644
index 0000000..0b5d23e
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/sample.go
@@ -0,0 +1,335 @@
+// Copyright 2015 The Go 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 stats
+
+import (
+	"math"
+	"sort"
+
+	"github.com/aclements/go-moremath/vec"
+)
+
+// Sample is a collection of possibly weighted data points.
+type Sample struct {
+	// Xs is the slice of sample values.
+	Xs []float64
+
+	// Weights[i] is the weight of sample Xs[i].  If Weights is
+	// nil, all Xs have weight 1.  Weights must have the same
+	// length of Xs and all values must be non-negative.
+	Weights []float64
+
+	// Sorted indicates that Xs is sorted in ascending order.
+	Sorted bool
+}
+
+// Bounds returns the minimum and maximum values of xs.
+func Bounds(xs []float64) (min float64, max float64) {
+	if len(xs) == 0 {
+		return math.NaN(), math.NaN()
+	}
+	min, max = xs[0], xs[0]
+	for _, x := range xs {
+		if x < min {
+			min = x
+		}
+		if x > max {
+			max = x
+		}
+	}
+	return
+}
+
+// Bounds returns the minimum and maximum values of the Sample.
+//
+// If the Sample is weighted, this ignores samples with zero weight.
+//
+// This is constant time if s.Sorted and there are no zero-weighted
+// values.
+func (s Sample) Bounds() (min float64, max float64) {
+	if len(s.Xs) == 0 || (!s.Sorted && s.Weights == nil) {
+		return Bounds(s.Xs)
+	}
+
+	if s.Sorted {
+		if s.Weights == nil {
+			return s.Xs[0], s.Xs[len(s.Xs)-1]
+		}
+		min, max = math.NaN(), math.NaN()
+		for i, w := range s.Weights {
+			if w != 0 {
+				min = s.Xs[i]
+				break
+			}
+		}
+		if math.IsNaN(min) {
+			return
+		}
+		for i := range s.Weights {
+			if s.Weights[len(s.Weights)-i-1] != 0 {
+				max = s.Xs[len(s.Weights)-i-1]
+				break
+			}
+		}
+	} else {
+		min, max = math.Inf(1), math.Inf(-1)
+		for i, x := range s.Xs {
+			w := s.Weights[i]
+			if x < min && w != 0 {
+				min = x
+			}
+			if x > max && w != 0 {
+				max = x
+			}
+		}
+		if math.IsInf(min, 0) {
+			min, max = math.NaN(), math.NaN()
+		}
+	}
+	return
+}
+
+// Sum returns the (possibly weighted) sum of the Sample.
+func (s Sample) Sum() float64 {
+	if s.Weights == nil {
+		return vec.Sum(s.Xs)
+	}
+	sum := 0.0
+	for i, x := range s.Xs {
+		sum += x * s.Weights[i]
+	}
+	return sum
+}
+
+// Weight returns the total weight of the Sasmple.
+func (s Sample) Weight() float64 {
+	if s.Weights == nil {
+		return float64(len(s.Xs))
+	}
+	return vec.Sum(s.Weights)
+}
+
+// Mean returns the arithmetic mean of xs.
+func Mean(xs []float64) float64 {
+	if len(xs) == 0 {
+		return math.NaN()
+	}
+	m := 0.0
+	for i, x := range xs {
+		m += (x - m) / float64(i+1)
+	}
+	return m
+}
+
+// Mean returns the arithmetic mean of the Sample.
+func (s Sample) Mean() float64 {
+	if len(s.Xs) == 0 || s.Weights == nil {
+		return Mean(s.Xs)
+	}
+
+	m, wsum := 0.0, 0.0
+	for i, x := range s.Xs {
+		// Use weighted incremental mean:
+		//   m_i = (1 - w_i/wsum_i) * m_(i-1) + (w_i/wsum_i) * x_i
+		//       = m_(i-1) + (x_i - m_(i-1)) * (w_i/wsum_i)
+		w := s.Weights[i]
+		wsum += w
+		m += (x - m) * w / wsum
+	}
+	return m
+}
+
+// GeoMean returns the geometric mean of xs. xs must be positive.
+func GeoMean(xs []float64) float64 {
+	if len(xs) == 0 {
+		return math.NaN()
+	}
+	m := 0.0
+	for i, x := range xs {
+		if x <= 0 {
+			return math.NaN()
+		}
+		lx := math.Log(x)
+		m += (lx - m) / float64(i+1)
+	}
+	return math.Exp(m)
+}
+
+// GeoMean returns the geometric mean of the Sample. All samples
+// values must be positive.
+func (s Sample) GeoMean() float64 {
+	if len(s.Xs) == 0 || s.Weights == nil {
+		return GeoMean(s.Xs)
+	}
+
+	m, wsum := 0.0, 0.0
+	for i, x := range s.Xs {
+		w := s.Weights[i]
+		wsum += w
+		lx := math.Log(x)
+		m += (lx - m) * w / wsum
+	}
+	return math.Exp(m)
+}
+
+// Variance returns the sample variance of xs.
+func Variance(xs []float64) float64 {
+	if len(xs) == 0 {
+		return math.NaN()
+	} else if len(xs) <= 1 {
+		return 0
+	}
+
+	// Based on Wikipedia's presentation of Welford 1962
+	// (http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm).
+	// This is more numerically stable than the standard two-pass
+	// formula and not prone to massive cancellation.
+	mean, M2 := 0.0, 0.0
+	for n, x := range xs {
+		delta := x - mean
+		mean += delta / float64(n+1)
+		M2 += delta * (x - mean)
+	}
+	return M2 / float64(len(xs)-1)
+}
+
+func (s Sample) Variance() float64 {
+	if len(s.Xs) == 0 || s.Weights == nil {
+		return Variance(s.Xs)
+	}
+	// TODO(austin)
+	panic("Weighted Variance not implemented")
+}
+
+// StdDev returns the sample standard deviation of xs.
+func StdDev(xs []float64) float64 {
+	return math.Sqrt(Variance(xs))
+}
+
+// StdDev returns the sample standard deviation of the Sample.
+func (s Sample) StdDev() float64 {
+	if len(s.Xs) == 0 || s.Weights == nil {
+		return StdDev(s.Xs)
+	}
+	// TODO(austin)
+	panic("Weighted StdDev not implemented")
+}
+
+// Quantile returns the sample value X at which q*weight of the sample
+// is <= X. This uses interpolation method R8 from Hyndman and Fan
+// (1996).
+//
+// q will be capped to the range [0, 1]. If len(xs) == 0 or all
+// weights are 0, returns NaN.
+//
+// Quantile(0.5) is the median. Quantile(0.25) and Quantile(0.75) are
+// the first and third quartiles, respectively. Quantile(P/100) is the
+// P'th percentile.
+//
+// This is constant time if s.Sorted and s.Weights == nil.
+func (s Sample) Quantile(q float64) float64 {
+	if len(s.Xs) == 0 {
+		return math.NaN()
+	} else if q <= 0 {
+		min, _ := s.Bounds()
+		return min
+	} else if q >= 1 {
+		_, max := s.Bounds()
+		return max
+	}
+
+	if !s.Sorted {
+		// TODO(austin) Use select algorithm instead
+		s = *s.Copy().Sort()
+	}
+
+	if s.Weights == nil {
+		N := float64(len(s.Xs))
+		//n := q * (N + 1) // R6
+		n := 1/3.0 + q*(N+1/3.0) // R8
+		kf, frac := math.Modf(n)
+		k := int(kf)
+		if k <= 0 {
+			return s.Xs[0]
+		} else if k >= len(s.Xs) {
+			return s.Xs[len(s.Xs)-1]
+		}
+		return s.Xs[k-1] + frac*(s.Xs[k]-s.Xs[k-1])
+	} else {
+		// TODO(austin): Implement interpolation
+
+		target := s.Weight() * q
+
+		// TODO(austin) If we had cumulative weights, we could
+		// do this in log time.
+		for i, weight := range s.Weights {
+			target -= weight
+			if target < 0 {
+				return s.Xs[i]
+			}
+		}
+		return s.Xs[len(s.Xs)-1]
+	}
+}
+
+// IQR returns the interquartile range of the Sample.
+//
+// This is constant time if s.Sorted and s.Weights == nil.
+func (s Sample) IQR() float64 {
+	if !s.Sorted {
+		s = *s.Copy().Sort()
+	}
+	return s.Quantile(0.75) - s.Quantile(0.25)
+}
+
+type sampleSorter struct {
+	xs      []float64
+	weights []float64
+}
+
+func (p *sampleSorter) Len() int {
+	return len(p.xs)
+}
+
+func (p *sampleSorter) Less(i, j int) bool {
+	return p.xs[i] < p.xs[j]
+}
+
+func (p *sampleSorter) Swap(i, j int) {
+	p.xs[i], p.xs[j] = p.xs[j], p.xs[i]
+	p.weights[i], p.weights[j] = p.weights[j], p.weights[i]
+}
+
+// Sort sorts the samples in place in s and returns s.
+//
+// A sorted sample improves the performance of some algorithms.
+func (s *Sample) Sort() *Sample {
+	if s.Sorted || sort.Float64sAreSorted(s.Xs) {
+		// All set
+	} else if s.Weights == nil {
+		sort.Float64s(s.Xs)
+	} else {
+		sort.Sort(&sampleSorter{s.Xs, s.Weights})
+	}
+	s.Sorted = true
+	return s
+}
+
+// Copy returns a copy of the Sample.
+//
+// The returned Sample shares no data with the original, so they can
+// be modified (for example, sorted) independently.
+func (s Sample) Copy() *Sample {
+	xs := make([]float64, len(s.Xs))
+	copy(xs, s.Xs)
+
+	weights := []float64(nil)
+	if s.Weights != nil {
+		weights = make([]float64, len(s.Weights))
+		copy(weights, s.Weights)
+	}
+
+	return &Sample{xs, weights, s.Sorted}
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/stream.go b/vendor/github.com/aclements/go-moremath/stats/stream.go
new file mode 100644
index 0000000..0deb904
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/stream.go
@@ -0,0 +1,100 @@
+// Copyright 2015 The Go 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 stats
+
+import (
+	"fmt"
+	"math"
+)
+
+// TODO(austin) Unify more with Sample interface
+
+// StreamStats tracks basic statistics for a stream of data in O(1)
+// space.
+//
+// StreamStats should be initialized to its zero value.
+type StreamStats struct {
+	Count           uint
+	Total, Min, Max float64
+
+	// Numerically stable online mean
+	mean          float64
+	meanOfSquares float64
+
+	// Online variance
+	vM2 float64
+}
+
+// Add updates s's statistics with sample value x.
+func (s *StreamStats) Add(x float64) {
+	s.Total += x
+	if s.Count == 0 {
+		s.Min, s.Max = x, x
+	} else {
+		if x < s.Min {
+			s.Min = x
+		}
+		if x > s.Max {
+			s.Max = x
+		}
+	}
+	s.Count++
+
+	// Update online mean, mean of squares, and variance.  Online
+	// variance based on Wikipedia's presentation ("Algorithms for
+	// calculating variance") of Knuth's formulation of Welford
+	// 1962.
+	delta := x - s.mean
+	s.mean += delta / float64(s.Count)
+	s.meanOfSquares += (x*x - s.meanOfSquares) / float64(s.Count)
+	s.vM2 += delta * (x - s.mean)
+}
+
+func (s *StreamStats) Weight() float64 {
+	return float64(s.Count)
+}
+
+func (s *StreamStats) Mean() float64 {
+	return s.mean
+}
+
+func (s *StreamStats) Variance() float64 {
+	return s.vM2 / float64(s.Count-1)
+}
+
+func (s *StreamStats) StdDev() float64 {
+	return math.Sqrt(s.Variance())
+}
+
+func (s *StreamStats) RMS() float64 {
+	return math.Sqrt(s.meanOfSquares)
+}
+
+// Combine updates s's statistics as if all samples added to o were
+// added to s.
+func (s *StreamStats) Combine(o *StreamStats) {
+	count := s.Count + o.Count
+
+	// Compute combined online variance statistics
+	delta := o.mean - s.mean
+	mean := s.mean + delta*float64(o.Count)/float64(count)
+	vM2 := s.vM2 + o.vM2 + delta*delta*float64(s.Count)*float64(o.Count)/float64(count)
+
+	s.Count = count
+	s.Total += o.Total
+	if o.Min < s.Min {
+		s.Min = o.Min
+	}
+	if o.Max > s.Max {
+		s.Max = o.Max
+	}
+	s.mean = mean
+	s.meanOfSquares += (o.meanOfSquares - s.meanOfSquares) * float64(o.Count) / float64(count)
+	s.vM2 = vM2
+}
+
+func (s *StreamStats) String() string {
+	return fmt.Sprintf("Count=%d Total=%g Min=%g Mean=%g RMS=%g Max=%g StdDev=%g", s.Count, s.Total, s.Min, s.Mean(), s.RMS(), s.Max, s.StdDev())
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/tdist.go b/vendor/github.com/aclements/go-moremath/stats/tdist.go
new file mode 100644
index 0000000..29bbb1a
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/tdist.go
@@ -0,0 +1,42 @@
+// Copyright 2015 The Go 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 stats
+
+import (
+	"math"
+
+	"github.com/aclements/go-moremath/mathx"
+)
+
+// A TDist is a Student's t-distribution with V degrees of freedom.
+type TDist struct {
+	V float64
+}
+
+func lgamma(x float64) float64 {
+	y, _ := math.Lgamma(x)
+	return y
+}
+
+func (t TDist) PDF(x float64) float64 {
+	return math.Exp(lgamma((t.V+1)/2)-lgamma(t.V/2)) /
+		math.Sqrt(t.V*math.Pi) * math.Pow(1+(x*x)/t.V, -(t.V+1)/2)
+}
+
+func (t TDist) CDF(x float64) float64 {
+	if x == 0 {
+		return 0.5
+	} else if x > 0 {
+		return 1 - 0.5*mathx.BetaInc(t.V/(t.V+x*x), t.V/2, 0.5)
+	} else if x < 0 {
+		return 1 - t.CDF(-x)
+	} else {
+		return math.NaN()
+	}
+}
+
+func (t TDist) Bounds() (float64, float64) {
+	return -4, 4
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/ttest.go b/vendor/github.com/aclements/go-moremath/stats/ttest.go
new file mode 100644
index 0000000..8742298
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/ttest.go
@@ -0,0 +1,147 @@
+// Copyright 2015 The Go 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 stats
+
+import (
+	"errors"
+	"math"
+)
+
+// A TTestResult is the result of a t-test.
+type TTestResult struct {
+	// N1 and N2 are the sizes of the input samples. For a
+	// one-sample t-test, N2 is 0.
+	N1, N2 int
+
+	// T is the value of the t-statistic for this t-test.
+	T float64
+
+	// DoF is the degrees of freedom for this t-test.
+	DoF float64
+
+	// AltHypothesis specifies the alternative hypothesis tested
+	// by this test against the null hypothesis that there is no
+	// difference in the means of the samples.
+	AltHypothesis LocationHypothesis
+
+	// P is p-value for this t-test for the given null hypothesis.
+	P float64
+}
+
+func newTTestResult(n1, n2 int, t, dof float64, alt LocationHypothesis) *TTestResult {
+	dist := TDist{dof}
+	var p float64
+	switch alt {
+	case LocationDiffers:
+		p = 2 * (1 - dist.CDF(math.Abs(t)))
+	case LocationLess:
+		p = dist.CDF(t)
+	case LocationGreater:
+		p = 1 - dist.CDF(t)
+	}
+	return &TTestResult{N1: n1, N2: n2, T: t, DoF: dof, AltHypothesis: alt, P: p}
+}
+
+// A TTestSample is a sample that can be used for a one or two sample
+// t-test.
+type TTestSample interface {
+	Weight() float64
+	Mean() float64
+	Variance() float64
+}
+
+var (
+	ErrSampleSize        = errors.New("sample is too small")
+	ErrZeroVariance      = errors.New("sample has zero variance")
+	ErrMismatchedSamples = errors.New("samples have different lengths")
+)
+
+// TwoSampleTTest performs a two-sample (unpaired) Student's t-test on
+// samples x1 and x2. This is a test of the null hypothesis that x1
+// and x2 are drawn from populations with equal means. It assumes x1
+// and x2 are independent samples, that the distributions have equal
+// variance, and that the populations are normally distributed.
+func TwoSampleTTest(x1, x2 TTestSample, alt LocationHypothesis) (*TTestResult, error) {
+	n1, n2 := x1.Weight(), x2.Weight()
+	if n1 == 0 || n2 == 0 {
+		return nil, ErrSampleSize
+	}
+	v1, v2 := x1.Variance(), x2.Variance()
+	if v1 == 0 && v2 == 0 {
+		return nil, ErrZeroVariance
+	}
+
+	dof := n1 + n2 - 2
+	v12 := ((n1-1)*v1 + (n2-1)*v2) / dof
+	t := (x1.Mean() - x2.Mean()) / math.Sqrt(v12*(1/n1+1/n2))
+	return newTTestResult(int(n1), int(n2), t, dof, alt), nil
+}
+
+// TwoSampleWelchTTest performs a two-sample (unpaired) Welch's t-test
+// on samples x1 and x2. This is like TwoSampleTTest, but does not
+// assume the distributions have equal variance.
+func TwoSampleWelchTTest(x1, x2 TTestSample, alt LocationHypothesis) (*TTestResult, error) {
+	n1, n2 := x1.Weight(), x2.Weight()
+	if n1 <= 1 || n2 <= 1 {
+		// TODO: Can we still do this with n == 1?
+		return nil, ErrSampleSize
+	}
+	v1, v2 := x1.Variance(), x2.Variance()
+	if v1 == 0 && v2 == 0 {
+		return nil, ErrZeroVariance
+	}
+
+	dof := math.Pow(v1/n1+v2/n2, 2) /
+		(math.Pow(v1/n1, 2)/(n1-1) + math.Pow(v2/n2, 2)/(n2-1))
+	s := math.Sqrt(v1/n1 + v2/n2)
+	t := (x1.Mean() - x2.Mean()) / s
+	return newTTestResult(int(n1), int(n2), t, dof, alt), nil
+}
+
+// PairedTTest performs a two-sample paired t-test on samples x1 and
+// x2. If μ0 is non-zero, this tests if the average of the difference
+// is significantly different from μ0. If x1 and x2 are identical,
+// this returns nil.
+func PairedTTest(x1, x2 []float64, μ0 float64, alt LocationHypothesis) (*TTestResult, error) {
+	if len(x1) != len(x2) {
+		return nil, ErrMismatchedSamples
+	}
+	if len(x1) <= 1 {
+		// TODO: Can we still do this with n == 1?
+		return nil, ErrSampleSize
+	}
+
+	dof := float64(len(x1) - 1)
+
+	diff := make([]float64, len(x1))
+	for i := range x1 {
+		diff[i] = x1[i] - x2[i]
+	}
+	sd := StdDev(diff)
+	if sd == 0 {
+		// TODO: Can we still do the test?
+		return nil, ErrZeroVariance
+	}
+	t := (Mean(diff) - μ0) * math.Sqrt(float64(len(x1))) / sd
+	return newTTestResult(len(x1), len(x2), t, dof, alt), nil
+}
+
+// OneSampleTTest performs a one-sample t-test on sample x. This tests
+// the null hypothesis that the population mean is equal to μ0. This
+// assumes the distribution of the population of sample means is
+// normal.
+func OneSampleTTest(x TTestSample, μ0 float64, alt LocationHypothesis) (*TTestResult, error) {
+	n, v := x.Weight(), x.Variance()
+	if n == 0 {
+		return nil, ErrSampleSize
+	}
+	if v == 0 {
+		// TODO: Can we still do the test?
+		return nil, ErrZeroVariance
+	}
+	dof := n - 1
+	t := (x.Mean() - μ0) * math.Sqrt(n) / math.Sqrt(v)
+	return newTTestResult(int(n), 0, t, dof, alt), nil
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/udist.go b/vendor/github.com/aclements/go-moremath/stats/udist.go
new file mode 100644
index 0000000..06c34ad
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/udist.go
@@ -0,0 +1,389 @@
+// Copyright 2015 The Go 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 stats
+
+import (
+	"math"
+
+	"github.com/aclements/go-moremath/mathx"
+)
+
+// A UDist is the discrete probability distribution of the
+// Mann-Whitney U statistic for a pair of samples of sizes N1 and N2.
+//
+// The details of computing this distribution with no ties can be
+// found in Mann, Henry B.; Whitney, Donald R. (1947). "On a Test of
+// Whether one of Two Random Variables is Stochastically Larger than
+// the Other". Annals of Mathematical Statistics 18 (1): 50–60.
+// Computing this distribution in the presence of ties is described in
+// Klotz, J. H. (1966). "The Wilcoxon, Ties, and the Computer".
+// Journal of the American Statistical Association 61 (315): 772-787
+// and Cheung, Ying Kuen; Klotz, Jerome H. (1997). "The Mann Whitney
+// Wilcoxon Distribution Using Linked Lists". Statistica Sinica 7:
+// 805-813 (the former paper contains details that are glossed over in
+// the latter paper but has mathematical typesetting issues, so it's
+// easiest to get the context from the former paper and the details
+// from the latter).
+type UDist struct {
+	N1, N2 int
+
+	// T is the count of the number of ties at each rank in the
+	// input distributions. T may be nil, in which case it is
+	// assumed there are no ties (which is equivalent to an M+N
+	// slice of 1s). It must be the case that Sum(T) == M+N.
+	T []int
+}
+
+// hasTies returns true if d has any tied samples.
+func (d UDist) hasTies() bool {
+	for _, t := range d.T {
+		if t > 1 {
+			return true
+		}
+	}
+	return false
+}
+
+// p returns the p_{d.N1,d.N2} function defined by Mann, Whitney 1947
+// for values of U from 0 up to and including the U argument.
+//
+// This algorithm runs in Θ(N1*N2*U) = O(N1²N2²) time and is quite
+// fast for small values of N1 and N2. However, it does not handle ties.
+func (d UDist) p(U int) []float64 {
+	// This is a dynamic programming implementation of the
+	// recursive recurrence definition given by Mann and Whitney:
+	//
+	//   p_{n,m}(U) = (n * p_{n-1,m}(U-m) + m * p_{n,m-1}(U)) / (n+m)
+	//   p_{n,m}(U) = 0                           if U < 0
+	//   p_{0,m}(U) = p{n,0}(U) = 1 / nCr(m+n, n) if U = 0
+	//                          = 0               if U > 0
+	//
+	// (Note that there is a typo in the original paper. The first
+	// recursive application of p should be for U-m, not U-M.)
+	//
+	// Since p{n,m} only depends on p{n-1,m} and p{n,m-1}, we only
+	// need to store one "plane" of the three dimensional space at
+	// a time.
+	//
+	// Furthermore, p_{n,m} = p_{m,n}, so we only construct values
+	// for n <= m and obtain the rest through symmetry.
+	//
+	// We organize the computed values of p as followed:
+	//
+	//       n →   N
+	//     m *
+	//     ↓ * *
+	//       * * *
+	//       * * * *
+	//       * * * *
+	//     M * * * *
+	//
+	// where each * is a slice indexed by U. The code below
+	// computes these left-to-right, top-to-bottom, so it only
+	// stores one row of this matrix at a time. Furthermore,
+	// computing an element in a given U slice only depends on the
+	// same and smaller values of U, so we can overwrite the U
+	// slice we're computing in place as long as we start with the
+	// largest value of U. Finally, even though the recurrence
+	// depends on (n,m) above the diagonal and we use symmetry to
+	// mirror those across the diagonal to (m,n), the mirrored
+	// indexes are always available in the current row, so this
+	// mirroring does not interfere with our ability to recycle
+	// state.
+
+	N, M := d.N1, d.N2
+	if N > M {
+		N, M = M, N
+	}
+
+	memo := make([][]float64, N+1)
+	for n := range memo {
+		memo[n] = make([]float64, U+1)
+	}
+
+	for m := 0; m <= M; m++ {
+		// Compute p_{0,m}. This is zero except for U=0.
+		memo[0][0] = 1
+
+		// Compute the remainder of this row.
+		nlim := N
+		if m < nlim {
+			nlim = m
+		}
+		for n := 1; n <= nlim; n++ {
+			lp := memo[n-1] // p_{n-1,m}
+			var rp []float64
+			if n <= m-1 {
+				rp = memo[n] // p_{n,m-1}
+			} else {
+				rp = memo[m-1] // p{m-1,n} and m==n
+			}
+
+			// For a given n,m, U is at most n*m.
+			//
+			// TODO: Actually, it's at most ⌈n*m/2⌉, but
+			// then we need to use more complex symmetries
+			// in the inner loop below.
+			ulim := n * m
+			if U < ulim {
+				ulim = U
+			}
+
+			out := memo[n] // p_{n,m}
+			nplusm := float64(n + m)
+			for U1 := ulim; U1 >= 0; U1-- {
+				l := 0.0
+				if U1-m >= 0 {
+					l = float64(n) * lp[U1-m]
+				}
+				r := float64(m) * rp[U1]
+				out[U1] = (l + r) / nplusm
+			}
+		}
+	}
+	return memo[N]
+}
+
+type ukey struct {
+	n1   int // size of first sample
+	twoU int // 2*U statistic for this permutation
+}
+
+// This computes the cumulative counts of the Mann-Whitney U
+// distribution in the presence of ties using the computation from
+// Cheung, Ying Kuen; Klotz, Jerome H. (1997). "The Mann Whitney
+// Wilcoxon Distribution Using Linked Lists". Statistica Sinica 7:
+// 805-813, with much guidance from appendix L of Klotz, A
+// Computational Approach to Statistics.
+//
+// makeUmemo constructs a table memo[K][ukey{n1, 2*U}], where K is the
+// number of ranks (up to len(t)), n1 is the size of the first sample
+// (up to the n1 argument), and U is the U statistic (up to the
+// argument twoU/2). The value of an entry in the memo table is the
+// number of permutations of a sample of size n1 in a ranking with tie
+// vector t[:K] having a U statistic <= U.
+func makeUmemo(twoU, n1 int, t []int) []map[ukey]float64 {
+	// Another candidate for a fast implementation is van de Wiel,
+	// "The split-up algorithm: a fast symbolic method for
+	// computing p-values of distribution-free statistics". This
+	// is what's used by R's coin package. It's a comparatively
+	// recent publication, so it's presumably faster (or perhaps
+	// just more general) than previous techniques, but I can't
+	// get my hands on the paper.
+	//
+	// TODO: ~40% of this function's time is spent in mapassign on
+	// the assignment lines in the two loops and another ~20% in
+	// map access and iteration. Improving map behavior or
+	// replacing the maps altogether with some other constant-time
+	// structure could double performance.
+	//
+	// TODO: The worst case for this function is when there are
+	// few ties. Yet the best case overall is when there are *no*
+	// ties. Can we get the best of both worlds? Use the fast
+	// algorithm for the most part when there are few ties and mix
+	// in the general algorithm just where we need it? That's
+	// certainly possible for sub-problems where t[:k] has no
+	// ties, but that doesn't help if t[0] has a tie but nothing
+	// else does. Is it possible to rearrange the ranks without
+	// messing up our computation of the U statistic for
+	// sub-problems?
+
+	K := len(t)
+
+	// Compute a coefficients. The a slice is indexed by k (a[0]
+	// is unused).
+	a := make([]int, K+1)
+	a[1] = t[0]
+	for k := 2; k <= K; k++ {
+		a[k] = a[k-1] + t[k-2] + t[k-1]
+	}
+
+	// Create the memo table for the counts function, A. The A
+	// slice is indexed by k (A[0] is unused).
+	//
+	// In "The Mann Whitney Distribution Using Linked Lists", they
+	// use linked lists (*gasp*) for this, but within each K it's
+	// really just a memoization table, so it's faster to use a
+	// map. The outer structure is a slice indexed by k because we
+	// need to find all memo entries with certain values of k.
+	//
+	// TODO: The n1 and twoU values in the ukeys follow strict
+	// patterns. For each K value, the n1 values are every integer
+	// between two bounds. For each (K, n1) value, the twoU values
+	// are every integer multiple of a certain base between two
+	// bounds. It might be worth turning these into directly
+	// indexible slices.
+	A := make([]map[ukey]float64, K+1)
+	A[K] = map[ukey]float64{ukey{n1: n1, twoU: twoU}: 0}
+
+	// Compute memo table (k, n1, twoU) triples from high K values
+	// to low K values. This drives the recurrence relation
+	// downward to figure out all of the needed argument triples.
+	//
+	// TODO: Is it possible to generate this table bottom-up? If
+	// so, this could be a pure dynamic programming algorithm and
+	// we could discard the K dimension. We could at least store
+	// the inputs in a more compact representation that replaces
+	// the twoU dimension with an interval and a step size (as
+	// suggested by Cheung, Klotz, not that they make it at all
+	// clear *why* they're suggesting this).
+	tsum := sumint(t) // always ∑ t[0:k]
+	for k := K - 1; k >= 2; k-- {
+		tsum -= t[k]
+		A[k] = make(map[ukey]float64)
+
+		// Construct A[k] from A[k+1].
+		for A_kplus1 := range A[k+1] {
+			rkLow := maxint(0, A_kplus1.n1-tsum)
+			rkHigh := minint(A_kplus1.n1, t[k])
+			for rk := rkLow; rk <= rkHigh; rk++ {
+				twoU_k := A_kplus1.twoU - rk*(a[k+1]-2*A_kplus1.n1+rk)
+				n1_k := A_kplus1.n1 - rk
+				if twoUmin(n1_k, t[:k], a) <= twoU_k && twoU_k <= twoUmax(n1_k, t[:k], a) {
+					key := ukey{n1: n1_k, twoU: twoU_k}
+					A[k][key] = 0
+				}
+			}
+		}
+	}
+
+	// Fill counts in memo table from low K values to high K
+	// values. This unwinds the recurrence relation.
+
+	// Start with K==2 base case.
+	//
+	// TODO: Later computations depend on these, but these don't
+	// depend on anything (including each other), so if K==2, we
+	// can skip the memo table altogether.
+	if K < 2 {
+		panic("K < 2")
+	}
+	N_2 := t[0] + t[1]
+	for A_2i := range A[2] {
+		Asum := 0.0
+		r2Low := maxint(0, A_2i.n1-t[0])
+		r2High := (A_2i.twoU - A_2i.n1*(t[0]-A_2i.n1)) / N_2
+		for r2 := r2Low; r2 <= r2High; r2++ {
+			Asum += mathx.Choose(t[0], A_2i.n1-r2) *
+				mathx.Choose(t[1], r2)
+		}
+		A[2][A_2i] = Asum
+	}
+
+	// Derive counts for the rest of the memo table.
+	tsum = t[0] // always ∑ t[0:k-1]
+	for k := 3; k <= K; k++ {
+		tsum += t[k-2]
+
+		// Compute A[k] counts from A[k-1] counts.
+		for A_ki := range A[k] {
+			Asum := 0.0
+			rkLow := maxint(0, A_ki.n1-tsum)
+			rkHigh := minint(A_ki.n1, t[k-1])
+			for rk := rkLow; rk <= rkHigh; rk++ {
+				twoU_kminus1 := A_ki.twoU - rk*(a[k]-2*A_ki.n1+rk)
+				n1_kminus1 := A_ki.n1 - rk
+				x, ok := A[k-1][ukey{n1: n1_kminus1, twoU: twoU_kminus1}]
+				if !ok && twoUmax(n1_kminus1, t[:k-1], a) < twoU_kminus1 {
+					x = mathx.Choose(tsum, n1_kminus1)
+				}
+				Asum += x * mathx.Choose(t[k-1], rk)
+			}
+			A[k][A_ki] = Asum
+		}
+	}
+
+	return A
+}
+
+func twoUmin(n1 int, t, a []int) int {
+	K := len(t)
+	twoU := -n1 * n1
+	n1_k := n1
+	for k := 1; k <= K; k++ {
+		twoU_k := minint(n1_k, t[k-1])
+		twoU += twoU_k * a[k]
+		n1_k -= twoU_k
+	}
+	return twoU
+}
+
+func twoUmax(n1 int, t, a []int) int {
+	K := len(t)
+	twoU := -n1 * n1
+	n1_k := n1
+	for k := K; k > 0; k-- {
+		twoU_k := minint(n1_k, t[k-1])
+		twoU += twoU_k * a[k]
+		n1_k -= twoU_k
+	}
+	return twoU
+}
+
+func (d UDist) PMF(U float64) float64 {
+	if U < 0 || U >= 0.5+float64(d.N1*d.N2) {
+		return 0
+	}
+
+	if d.hasTies() {
+		// makeUmemo computes the CDF directly. Take its
+		// difference to get the PMF.
+		p1, ok1 := makeUmemo(int(2*U)-1, d.N1, d.T)[len(d.T)][ukey{d.N1, int(2*U) - 1}]
+		p2, ok2 := makeUmemo(int(2*U), d.N1, d.T)[len(d.T)][ukey{d.N1, int(2 * U)}]
+		if !ok1 || !ok2 {
+			panic("makeUmemo did not return expected memoization table")
+		}
+		return (p2 - p1) / mathx.Choose(d.N1+d.N2, d.N1)
+	}
+
+	// There are no ties. Use the fast algorithm. U must be integral.
+	Ui := int(math.Floor(U))
+	// TODO: Use symmetry to minimize U
+	return d.p(Ui)[Ui]
+}
+
+func (d UDist) CDF(U float64) float64 {
+	if U < 0 {
+		return 0
+	} else if U >= float64(d.N1*d.N2) {
+		return 1
+	}
+
+	if d.hasTies() {
+		// TODO: Minimize U?
+		p, ok := makeUmemo(int(2*U), d.N1, d.T)[len(d.T)][ukey{d.N1, int(2 * U)}]
+		if !ok {
+			panic("makeUmemo did not return expected memoization table")
+		}
+		return p / mathx.Choose(d.N1+d.N2, d.N1)
+	}
+
+	// There are no ties. Use the fast algorithm. U must be integral.
+	Ui := int(math.Floor(U))
+	// The distribution is symmetric around U = m * n / 2. Sum up
+	// whichever tail is smaller.
+	flip := Ui >= (d.N1*d.N2+1)/2
+	if flip {
+		Ui = d.N1*d.N2 - Ui - 1
+	}
+	pdfs := d.p(Ui)
+	p := 0.0
+	for _, pdf := range pdfs[:Ui+1] {
+		p += pdf
+	}
+	if flip {
+		p = 1 - p
+	}
+	return p
+}
+
+func (d UDist) Step() float64 {
+	return 0.5
+}
+
+func (d UDist) Bounds() (float64, float64) {
+	// TODO: More precise bounds when there are ties.
+	return 0, float64(d.N1 * d.N2)
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/utest.go b/vendor/github.com/aclements/go-moremath/stats/utest.go
new file mode 100644
index 0000000..c31e54c
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/utest.go
@@ -0,0 +1,276 @@
+// Copyright 2015 The Go 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 stats
+
+import (
+	"math"
+	"sort"
+
+	"github.com/aclements/go-moremath/mathx"
+)
+
+// A LocationHypothesis specifies the alternative hypothesis of a
+// location test such as a t-test or a Mann-Whitney U-test. The
+// default (zero) value is to test against the alternative hypothesis
+// that they differ.
+type LocationHypothesis int
+
+//go:generate stringer -type LocationHypothesis
+
+const (
+	// LocationLess specifies the alternative hypothesis that the
+	// location of the first sample is less than the second. This
+	// is a one-tailed test.
+	LocationLess LocationHypothesis = -1
+
+	// LocationDiffers specifies the alternative hypothesis that
+	// the locations of the two samples are not equal. This is a
+	// two-tailed test.
+	LocationDiffers LocationHypothesis = 0
+
+	// LocationGreater specifies the alternative hypothesis that
+	// the location of the first sample is greater than the
+	// second. This is a one-tailed test.
+	LocationGreater LocationHypothesis = 1
+)
+
+// A MannWhitneyUTestResult is the result of a Mann-Whitney U-test.
+type MannWhitneyUTestResult struct {
+	// N1 and N2 are the sizes of the input samples.
+	N1, N2 int
+
+	// U is the value of the Mann-Whitney U statistic for this
+	// test, generalized by counting ties as 0.5.
+	//
+	// Given the Cartesian product of the two samples, this is the
+	// number of pairs in which the value from the first sample is
+	// greater than the value of the second, plus 0.5 times the
+	// number of pairs where the values from the two samples are
+	// equal. Hence, U is always an integer multiple of 0.5 (it is
+	// a whole integer if there are no ties) in the range [0, N1*N2].
+	//
+	// U statistics always come in pairs, depending on which
+	// sample is "first". The mirror U for the other sample can be
+	// calculated as N1*N2 - U.
+	//
+	// There are many equivalent statistics with slightly
+	// different definitions. The Wilcoxon (1945) W statistic
+	// (generalized for ties) is U + (N1(N1+1))/2. It is also
+	// common to use 2U to eliminate the half steps and Smid
+	// (1956) uses N1*N2 - 2U to additionally center the
+	// distribution.
+	U float64
+
+	// AltHypothesis specifies the alternative hypothesis tested
+	// by this test against the null hypothesis that there is no
+	// difference in the locations of the samples.
+	AltHypothesis LocationHypothesis
+
+	// P is the p-value of the Mann-Whitney test for the given
+	// null hypothesis.
+	P float64
+}
+
+// MannWhitneyExactLimit gives the largest sample size for which the
+// exact U distribution will be used for the Mann-Whitney U-test.
+//
+// Using the exact distribution is necessary for small sample sizes
+// because the distribution is highly irregular. However, computing
+// the distribution for large sample sizes is both computationally
+// expensive and unnecessary because it quickly approaches a normal
+// approximation. Computing the distribution for two 50 value samples
+// takes a few milliseconds on a 2014 laptop.
+var MannWhitneyExactLimit = 50
+
+// MannWhitneyTiesExactLimit gives the largest sample size for which
+// the exact U distribution will be used for the Mann-Whitney U-test
+// in the presence of ties.
+//
+// Computing this distribution is more expensive than computing the
+// distribution without ties, so this is set lower. Computing this
+// distribution for two 25 value samples takes about ten milliseconds
+// on a 2014 laptop.
+var MannWhitneyTiesExactLimit = 25
+
+// MannWhitneyUTest performs a Mann-Whitney U-test [1,2] of the null
+// hypothesis that two samples come from the same population against
+// the alternative hypothesis that one sample tends to have larger or
+// smaller values than the other.
+//
+// This is similar to a t-test, but unlike the t-test, the
+// Mann-Whitney U-test is non-parametric (it does not assume a normal
+// distribution). It has very slightly lower efficiency than the
+// t-test on normal distributions.
+//
+// Computing the exact U distribution is expensive for large sample
+// sizes, so this uses a normal approximation for sample sizes larger
+// than MannWhitneyExactLimit if there are no ties or
+// MannWhitneyTiesExactLimit if there are ties. This normal
+// approximation uses both the tie correction and the continuity
+// correction.
+//
+// This can fail with ErrSampleSize if either sample is empty or
+// ErrSamplesEqual if all sample values are equal.
+//
+// This is also known as a Mann-Whitney-Wilcoxon test and is
+// equivalent to the Wilcoxon rank-sum test, though the Wilcoxon
+// rank-sum test differs in nomenclature.
+//
+// [1] Mann, Henry B.; Whitney, Donald R. (1947). "On a Test of
+// Whether one of Two Random Variables is Stochastically Larger than
+// the Other". Annals of Mathematical Statistics 18 (1): 50–60.
+//
+// [2] Klotz, J. H. (1966). "The Wilcoxon, Ties, and the Computer".
+// Journal of the American Statistical Association 61 (315): 772-787.
+func MannWhitneyUTest(x1, x2 []float64, alt LocationHypothesis) (*MannWhitneyUTestResult, error) {
+	n1, n2 := len(x1), len(x2)
+	if n1 == 0 || n2 == 0 {
+		return nil, ErrSampleSize
+	}
+
+	// Compute the U statistic and tie vector T.
+	x1 = append([]float64(nil), x1...)
+	x2 = append([]float64(nil), x2...)
+	sort.Float64s(x1)
+	sort.Float64s(x2)
+	merged, labels := labeledMerge(x1, x2)
+
+	R1 := 0.0
+	T, hasTies := []int{}, false
+	for i := 0; i < len(merged); {
+		rank1, nx1, v1 := i+1, 0, merged[i]
+		// Consume samples that tie this sample (including itself).
+		for ; i < len(merged) && merged[i] == v1; i++ {
+			if labels[i] == 1 {
+				nx1++
+			}
+		}
+		// Assign all tied samples the average rank of the
+		// samples, where merged[0] has rank 1.
+		if nx1 != 0 {
+			rank := float64(i+rank1) / 2
+			R1 += rank * float64(nx1)
+		}
+		T = append(T, i-rank1+1)
+		if i > rank1 {
+			hasTies = true
+		}
+	}
+	U1 := R1 - float64(n1*(n1+1))/2
+
+	// Compute the smaller of U1 and U2
+	U2 := float64(n1*n2) - U1
+	Usmall := math.Min(U1, U2)
+
+	var p float64
+	if !hasTies && n1 <= MannWhitneyExactLimit && n2 <= MannWhitneyExactLimit ||
+		hasTies && n1 <= MannWhitneyTiesExactLimit && n2 <= MannWhitneyTiesExactLimit {
+		// Use exact U distribution. U1 will be an integer.
+		if len(T) == 1 {
+			// All values are equal. Test is meaningless.
+			return nil, ErrSamplesEqual
+		}
+
+		dist := UDist{N1: n1, N2: n2, T: T}
+		switch alt {
+		case LocationDiffers:
+			if U1 == U2 {
+				// The distribution is symmetric about
+				// Usmall. Since the distribution is
+				// discrete, the CDF is discontinuous
+				// and if simply double CDF(Usmall),
+				// we'll double count the
+				// (non-infinitesimal) probability
+				// mass at Usmall. What we want is
+				// just the integral of the whole CDF,
+				// which is 1.
+				p = 1
+			} else {
+				p = dist.CDF(Usmall) * 2
+			}
+
+		case LocationLess:
+			p = dist.CDF(U1)
+
+		case LocationGreater:
+			p = 1 - dist.CDF(U1-1)
+		}
+	} else {
+		// Use normal approximation (with tie and continuity
+		// correction).
+		t := tieCorrection(T)
+		N := float64(n1 + n2)
+		μ_U := float64(n1*n2) / 2
+		σ_U := math.Sqrt(float64(n1*n2) * ((N + 1) - t/(N*(N-1))) / 12)
+		if σ_U == 0 {
+			return nil, ErrSamplesEqual
+		}
+		numer := U1 - μ_U
+		// Perform continuity correction.
+		switch alt {
+		case LocationDiffers:
+			numer -= mathx.Sign(numer) * 0.5
+		case LocationLess:
+			numer += 0.5
+		case LocationGreater:
+			numer -= 0.5
+		}
+		z := numer / σ_U
+		switch alt {
+		case LocationDiffers:
+			p = 2 * math.Min(StdNormal.CDF(z), 1-StdNormal.CDF(z))
+		case LocationLess:
+			p = StdNormal.CDF(z)
+		case LocationGreater:
+			p = 1 - StdNormal.CDF(z)
+		}
+	}
+
+	return &MannWhitneyUTestResult{N1: n1, N2: n2, U: U1,
+		AltHypothesis: alt, P: p}, nil
+}
+
+// labeledMerge merges sorted lists x1 and x2 into sorted list merged.
+// labels[i] is 1 or 2 depending on whether merged[i] is a value from
+// x1 or x2, respectively.
+func labeledMerge(x1, x2 []float64) (merged []float64, labels []byte) {
+	merged = make([]float64, len(x1)+len(x2))
+	labels = make([]byte, len(x1)+len(x2))
+
+	i, j, o := 0, 0, 0
+	for i < len(x1) && j < len(x2) {
+		if x1[i] < x2[j] {
+			merged[o] = x1[i]
+			labels[o] = 1
+			i++
+		} else {
+			merged[o] = x2[j]
+			labels[o] = 2
+			j++
+		}
+		o++
+	}
+	for ; i < len(x1); i++ {
+		merged[o] = x1[i]
+		labels[o] = 1
+		o++
+	}
+	for ; j < len(x2); j++ {
+		merged[o] = x2[j]
+		labels[o] = 2
+		o++
+	}
+	return
+}
+
+// tieCorrection computes the tie correction factor Σ_j (t_j³ - t_j)
+// where t_j is the number of ties in the j'th rank.
+func tieCorrection(ties []int) float64 {
+	t := 0
+	for _, tie := range ties {
+		t += tie*tie*tie - tie
+	}
+	return float64(t)
+}
diff --git a/vendor/github.com/aclements/go-moremath/vec/package.go b/vendor/github.com/aclements/go-moremath/vec/package.go
new file mode 100644
index 0000000..6bd2061
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/vec/package.go
@@ -0,0 +1,6 @@
+// Copyright 2015 The Go 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 vec provides functions for float64 vectors.
+package vec // import "github.com/aclements/go-moremath/vec"
diff --git a/vendor/github.com/aclements/go-moremath/vec/vec.go b/vendor/github.com/aclements/go-moremath/vec/vec.go
new file mode 100644
index 0000000..970d441
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/vec/vec.go
@@ -0,0 +1,76 @@
+// Copyright 2015 The Go 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 vec
+
+import "math"
+
+// Vectorize returns a function g(xs) that applies f to each x in xs.
+//
+// f may be evaluated in parallel and in any order.
+func Vectorize(f func(float64) float64) func(xs []float64) []float64 {
+	return func(xs []float64) []float64 {
+		return Map(f, xs)
+	}
+}
+
+// Map returns f(x) for each x in xs.
+//
+// f may be evaluated in parallel and in any order.
+func Map(f func(float64) float64, xs []float64) []float64 {
+	// TODO(austin) Parallelize
+	res := make([]float64, len(xs))
+	for i, x := range xs {
+		res[i] = f(x)
+	}
+	return res
+}
+
+// Linspace returns num values spaced evenly between lo and hi,
+// inclusive. If num is 1, this returns an array consisting of lo.
+func Linspace(lo, hi float64, num int) []float64 {
+	res := make([]float64, num)
+	if num == 1 {
+		res[0] = lo
+		return res
+	}
+	for i := 0; i < num; i++ {
+		res[i] = lo + float64(i)*(hi-lo)/float64(num-1)
+	}
+	return res
+}
+
+// Logspace returns num values spaced evenly on a logarithmic scale
+// between base**lo and base**hi, inclusive.
+func Logspace(lo, hi float64, num int, base float64) []float64 {
+	res := Linspace(lo, hi, num)
+	for i, x := range res {
+		res[i] = math.Pow(base, x)
+	}
+	return res
+}
+
+// Sum returns the sum of xs.
+func Sum(xs []float64) float64 {
+	sum := 0.0
+	for _, x := range xs {
+		sum += x
+	}
+	return sum
+}
+
+// Concat returns the concatenation of its arguments. It does not
+// modify its inputs.
+func Concat(xss ...[]float64) []float64 {
+	total := 0
+	for _, xs := range xss {
+		total += len(xs)
+	}
+	out := make([]float64, total)
+	pos := 0
+	for _, xs := range xss {
+		pos += copy(out[pos:], xs)
+	}
+	return out
+}
diff --git a/vendor/vendor.json b/vendor/vendor.json
index fa0972c..12623ea 100644
--- a/vendor/vendor.json
+++ b/vendor/vendor.json
@@ -27,6 +27,54 @@
 			"revisionTime": "2017-02-02T22:37:47Z"
 		},
 		{
+			"checksumSHA1": "b95iQdItSVgUj7mm2lGPq6bAKuI=",
+			"path": "github.com/aclements/go-gg/generic",
+			"revision": "6dbb4e4fefb062aa608dfc5abc3279753d400000",
+			"revisionTime": "2016-09-21T14:28:07Z"
+		},
+		{
+			"checksumSHA1": "AnfMX9UQwGTIk1lOXvHC5neV/Ys=",
+			"path": "github.com/aclements/go-gg/generic/slice",
+			"revision": "6dbb4e4fefb062aa608dfc5abc3279753d400000",
+			"revisionTime": "2016-09-21T14:28:07Z"
+		},
+		{
+			"checksumSHA1": "dCOm7590kLFskNaN1qxnZqv3Uoc=",
+			"path": "github.com/aclements/go-gg/ggstat",
+			"revision": "6dbb4e4fefb062aa608dfc5abc3279753d400000",
+			"revisionTime": "2016-09-21T14:28:07Z"
+		},
+		{
+			"checksumSHA1": "iSkv+GYexi78UZka28H+2tdC5S0=",
+			"path": "github.com/aclements/go-gg/table",
+			"revision": "6dbb4e4fefb062aa608dfc5abc3279753d400000",
+			"revisionTime": "2016-09-21T14:28:07Z"
+		},
+		{
+			"checksumSHA1": "0R5tItJTud6pMwB/vVna0FdHsFg=",
+			"path": "github.com/aclements/go-moremath/fit",
+			"revision": "0ff62e0875ff52f29a54f430cefa15e4775c99b2",
+			"revisionTime": "2016-10-14T18:41:02Z"
+		},
+		{
+			"checksumSHA1": "eKoKQST+S7X0jcu8IEa/URMOsbY=",
+			"path": "github.com/aclements/go-moremath/mathx",
+			"revision": "0ff62e0875ff52f29a54f430cefa15e4775c99b2",
+			"revisionTime": "2016-10-14T18:41:02Z"
+		},
+		{
+			"checksumSHA1": "QJmKpM8L5l7aLi6zIRXBArf9mPk=",
+			"path": "github.com/aclements/go-moremath/stats",
+			"revision": "0ff62e0875ff52f29a54f430cefa15e4775c99b2",
+			"revisionTime": "2016-10-14T18:41:02Z"
+		},
+		{
+			"checksumSHA1": "DZGrSYl93Bo71dvmtmvdxwi6rKM=",
+			"path": "github.com/aclements/go-moremath/vec",
+			"revision": "0ff62e0875ff52f29a54f430cefa15e4775c99b2",
+			"revisionTime": "2016-10-14T18:41:02Z"
+		},
+		{
 			"checksumSHA1": "q3Bc7JpLWBqhZ4M7oreGo34RSkc=",
 			"path": "github.com/golang/protobuf/proto",
 			"revision": "8ee79997227bf9b34611aee7946ae64735e6fd93",