Generate gaussian

This is an example Go program for generate-gaussian.

package main

import (
        "fmt"
        "math"
        "math/rand"
        "os"

        "gonum.org/v1/plot"
        "gonum.org/v1/plot/plotter"
        "gonum.org/v1/plot/plotutil"
        "gonum.org/v1/plot/vg"
)

func gaussian(mu, sigma float64) (p, q float64) {
        var u, v, rsq float64
        for {
                u = 2*rand.Float64() - 1
                v = 2*rand.Float64() - 1
                rsq = u*u + v*v
                if rsq < 1 {
                        break
                }
        }
        r := math.Sqrt(-2 * math.Log(rsq) / rsq)
        p = mu + u*r*sigma
        q = mu + v*r*sigma
        return p, q
}

func main() {
        max := 1000
        gauss := make([]float64, max)
        mu := 0.0
        sigma := 1.0
        for i := 0; i < max/2; i++ {
                gauss[2*i], gauss[2*i+1] = gaussian(mu, sigma)
        }
        bins := 20
        binsize := sigma * 8 / float64(bins)
        mubin := float64(bins) / 2
        graph := make([]int, bins)
        for i := 0; i < max; i++ {
                bin := int(math.Round(gauss[i]/binsize + mubin))
                if bin >= 0 && bin < bins {
                        graph[bin]++
                } else {
                        fmt.Printf("Discarding %g bin %d\n", gauss[i], bin)
                }
        }
        points := make(plotter.XYs, bins)
        for i := 0; i < bins; i++ {
                f := float64(i)
                start := f*binsize - mubin*binsize
                end := start + binsize
                points[i].X = (start + end) / 2
                points[i].Y = float64(graph[i])
                fmt.Printf("%d [%.2f-%.2f] %d\n", i, start, end, graph[i])
        }
        file := "gaussian.png"
        p := plot.New()
        p.Title.Text = "Gaussian-distributed random points"
        p.X.Label.Text = "X"
        p.Y.Label.Text = "Y"
        err := plotutil.AddLinePoints(p, "first", points)
        if err != nil {
                fmt.Fprintf(os.Stderr, "Error plotting points: %s\n", err)
        }
        err = p.Save(4*vg.Inch, 4*vg.Inch, file)
        if err != nil {
                fmt.Fprintf(os.Stderr, "Error saving to %s: %s\n", file, err)
                os.Exit(1)
        }
}
func randomPoints(n int) plotter.XYs {
        pts := make(plotter.XYs, n)
        for i := range pts {
                if i == 0 {
                        pts[i].X = rand.Float64()
                } else {
                        pts[i].X = pts[i-1].X + rand.Float64()
                }
                pts[i].Y = pts[i].X + 10*rand.Float64()
        }
        return pts
}

(download)

The output of the example looks like this:

0 [-4.00--3.60] 0
1 [-3.60--3.20] 0
2 [-3.20--2.80] 0
3 [-2.80--2.40] 3
4 [-2.40--2.00] 9
5 [-2.00--1.60] 23
6 [-1.60--1.20] 43
7 [-1.20--0.80] 75
8 [-0.80--0.40] 128
9 [-0.40-0.00] 132
10 [0.00-0.40] 173
11 [0.40-0.80] 130
12 [0.80-1.20] 126
13 [1.20-1.60] 81
14 [1.60-2.00] 46
15 [2.00-2.40] 20
16 [2.40-2.80] 8
17 [2.80-3.20] 3
18 [3.20-3.60] 0
19 [3.60-4.00] 0