ex1.go

package main

import (
	"fmt"
	"image/color"
	"log"
	"math"
	"sort"

	"golang.org/x/exp/rand"
	"gonum.org/v1/gonum/stat"
	"gonum.org/v1/gonum/stat/distuv"
	"gonum.org/v1/plot"
	"gonum.org/v1/plot/plotter"
	"gonum.org/v1/plot/vg"
)

const mu float64 = 0.2
const sigma float64 = 0.1

func main() {
	drawBetaDistributionFunction()
	drawRandomHistogram()
}

func drawBetaDistributionFunction() {
	p, err := plot.New()
	if err != nil {
		panic(err)
	}
	p.Title.Text = "Beta distribution function"

	alpha := calcAlpha()
	beta := calcBeta(alpha)
	betaDistr := plotter.NewFunction(distuv.Beta{Alpha: alpha, Beta: beta}.Prob)

	p.Add(betaDistr)

	p.X.Min = 0
	p.X.Max = 1
	p.Y.Min = 0
	p.Y.Max = 5

	if err := p.Save(10*vg.Inch, 10*vg.Inch, "func.png"); err != nil {
		panic(err)
	}
}

func drawRandomHistogram() {
	p, err := plot.New()
	if err != nil {
		panic(err)
	}
	p.Title.Text = "Beta distribution and random sample histogram"

	alpha := calcAlpha()
	beta := calcBeta(alpha)

	src := rand.New(rand.NewSource(1))
	b := distuv.Beta{Alpha: alpha, Beta: beta, Src: src}

	v := make(plotter.Values, 10000)
	sum := 0.0
	for i := range v {
		r := b.Rand()
		v[i] = r
		sum = sum + r
	}
	sort.Float64s(v)

	h, err := plotter.NewHist(v, 16)
	if err != nil {
		panic(err)
	}
	h.Normalize(1)
	p.Add(h)

	labels, err := plotter.NewLabels(plotter.XYLabels{
		XYs: []plotter.XY{
			{X: 0.8, Y: 4.5},
			{X: 0.8, Y: 4.3},
			{X: 0.8, Y: 4.1},
			{X: 0.8, Y: 3.9},
			{X: 0.8, Y: 3.7},
			{X: 0.8, Y: 3.5},
		},
		Labels: []string{fmt.Sprintf("Sampled Mean: %f\n", sum/10000.0), fmt.Sprintf("True mean: %f\n", b.Mean()), fmt.Sprintf("Sampled Variance: %f\n", stat.Variance(v, nil)),
			fmt.Sprintf("True variance: %f\n", b.Variance()), fmt.Sprintf("Sampled Quantile: %f\n", stat.Quantile(0.95, stat.Empirical, v, nil)),
			fmt.Sprintf("True quantile: %f\n", b.Quantile(0.95))},
	},
	)
	if err != nil {
		log.Fatalf("could not creates labels plotter: %+v", err)
	}

	p.Add(labels)

	betaDistr := plotter.NewFunction(distuv.Beta{Alpha: alpha, Beta: beta}.Prob)
	betaDistr.Color = color.RGBA{R: 255, G: 100, A: 255}
	betaDistr.Width = vg.Points(3)
	p.Add(betaDistr)

	p.X.Min = 0
	p.X.Max = 1
	p.Y.Min = 0
	p.Y.Max = 5

	if err := p.Save(10*vg.Inch, 10*vg.Inch, "hist.png"); err != nil {
		panic(err)
	}
}

func calcAlpha() float64 {
	return mu * (((mu * (1 - mu)) / math.Pow(sigma, 2)) - 1)
}

func calcBeta(alpha float64) float64 {
	return alpha * (1 - mu) / mu
}