NunoSempere/fermi
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

feng-shui: move pretty printers to own files, simply sampling

Browse Source
This commit is contained in:
NunoSempere 2024-12-28 23:25:35 +01:00
parent 0c03b7e6bf
commit 3e4df77541
7 changed files with 318 additions and 135 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
fermi
View File

Binary file not shown.

10
main/errors.go Normal file
View File

@ -0,0 +1,10 @@
package main
type MyError struct {
Code string
Msg string
}
func (me MyError) Error() string {
return me.Msg
}

124
fermi.go → main/fermi.go
View File

@ -5,11 +5,9 @@ import (
"errors"
"flag"
"fmt"
"git.nunosempere.com/NunoSempere/fermi/pretty"
"git.nunosempere.com/NunoSempere/fermi/sample"
"math"
"os"
"sort"
"strings"
)
@ -46,11 +44,11 @@ func (p Scalar) Sampler(i int, r sample.State) float64 {
}
func (ln Lognormal) Sampler(i int, r sample.State) float64 {
return sample.Sample_to(ln.low, ln.high, r)
return sample.To(ln.low, ln.high, r)
}
func (beta Beta) Sampler(i int, r sample.State) float64 {
return sample.Sample_beta(beta.a, beta.b, r)
return sample.Beta(beta.a, beta.b, r)
}
func (fs FilledSamples) Sampler(i int, r sample.State) float64 {
@ -112,94 +110,12 @@ const INIT_DIST Scalar = Scalar(1)
var N_SAMPLES = 100_000
/* Printers */
func prettyPrintDist(dist Dist) {
switch v := dist.(type) {
case Lognormal:
fmt.Printf("=> ")
pretty.PrettyPrint2Floats(v.low, v.high)
fmt.Println()
case Beta:
fmt.Printf("=> beta ")
pretty.PrettyPrint2Floats(v.a, v.b)
fmt.Println()
case Scalar:
fmt.Printf("=> scalar ")
w := float64(v)
pretty.PrettyPrintFloat(w)
fmt.Println()
case FilledSamples:
n := len(v.xs)
sorted_xs := make([]float64, n)
copy(sorted_xs, v.xs)
sort.Slice(sorted_xs, func(i, j int) bool {
return sorted_xs[i] < sorted_xs[j]
})
low := sorted_xs[int(math.Round(float64(n)*0.05))]
high := sorted_xs[int(math.Round(float64(n)*0.95))]
fmt.Printf("=> ")
pretty.PrettyPrint2Floats(low, high)
fmt.Printf(" (")
pretty.PrettyPrintInt(N_SAMPLES)
fmt.Printf(" samples)")
fmt.Println()
default:
fmt.Printf("%v\n", v)
}
}
func printAndReturnErr(err_msg string) error {
fmt.Println(err_msg)
fmt.Println("Type \"help\" (without quotes) to see a pseudogrammar and examples")
return errors.New(err_msg)
}
func prettyPrintStats(dist Dist) {
xs := sample.Sample_serially(dist.Sampler, N_SAMPLES)
n := len(xs)
mean := 0.0
for i := 0; i < n; i++ {
mean += xs[i]
}
mean /= float64(n)
fmt.Printf("Mean: %f\n", mean)
stdev := 0.0
for i := 0; i < n; i++ {
stdev += math.Pow(xs[i]-mean, 2)
}
stdev = math.Sqrt(stdev / float64(n))
fmt.Printf("Stdev: %f\n", stdev)
sorted_xs := make([]float64, n)
copy(sorted_xs, xs)
sort.Slice(sorted_xs, func(i, j int) bool {
return sorted_xs[i] < sorted_xs[j]
})
print_ci := func(ci float64, prefix string) {
x := sorted_xs[int(math.Round(float64(n)*ci))]
fmt.Printf("%s%f\n", prefix, x)
}
print_ci(0.01, "ci 1%: ")
print_ci(0.05, "ci 5%: ")
print_ci(0.10, "ci 10%: ")
print_ci(0.25, "ci 25%: ")
print_ci(0.50, "ci 50%: ")
print_ci(0.75, "ci 75%: ")
print_ci(0.90, "ci 90%: ")
print_ci(0.95, "ci 95%: ")
print_ci(0.99, "ci 99%: ")
}
/* Operations */
// Generic operations with samples
func operateDistsAsSamples(dist1 Dist, dist2 Dist, op string) (Dist, error) {
xs := sample.Sample_serially(dist1.Sampler, N_SAMPLES)
ys := sample.Sample_serially(dist2.Sampler, N_SAMPLES)
xs := sample.Serially(dist1.Sampler, N_SAMPLES)
ys := sample.Serially(dist2.Sampler, N_SAMPLES)
zs := make([]float64, N_SAMPLES)
for i := 0; i < N_SAMPLES; i++ {
@ -335,7 +251,7 @@ func operateDists(old_dist Dist, new_dist Dist, op string) (Dist, error) {
case "-":
return operateDistsAsSamples(old_dist, new_dist, "-")
default:
return nil, printAndReturnErr("Can't combine distributions in this way")
return nil, PrintAndReturnErr("Can't combine distributions in this way")
}
}
@ -344,7 +260,7 @@ func parseMixture(words []string, vars map[string]Dist) (Dist, error) {
// mx, mix, var weight var weight var weight ...
// Check syntax
if len(words)%2 != 0 {
return nil, printAndReturnErr("Not a mixture. \nMixture syntax: \nmx x 2.5 y 8 z 10\ni.e.: mx var weight var2 weight2 ... var_n weight_n")
return nil, PrintAndReturnErr("Not a mixture. \nMixture syntax: \nmx x 2.5 y 8 z 10\ni.e.: mx var weight var2 weight2 ... var_n weight_n")
}
var fs []func(int, sample.State) float64
@ -354,29 +270,29 @@ func parseMixture(words []string, vars map[string]Dist) (Dist, error) {
if i%2 == 0 {
dist, exists := vars[word]
if !exists {
return nil, printAndReturnErr("Expected mixture variable but didn't get a variable. \nMixture syntax: \nmx x 2.5 y 8 z 10\ni.e.: mx var weight var2 weight2 ... var_n weight_n")
return nil, PrintAndReturnErr("Expected mixture variable but didn't get a variable. \nMixture syntax: \nmx x 2.5 y 8 z 10\ni.e.: mx var weight var2 weight2 ... var_n weight_n")
}
f := dist.Sampler
fs = append(fs, f)
} else {
weight, err := pretty.ParseFloat(word)
weight, err := ParseFloat(word)
if err != nil {
return nil, printAndReturnErr("Expected mixture weight but didn't get a float. \nMixture syntax: \nmx x 2.5 y 8 z 10\ni.e.: mx var weight var2 weight2 ... var_n weight_n")
return nil, PrintAndReturnErr("Expected mixture weight but didn't get a float. \nMixture syntax: \nmx x 2.5 y 8 z 10\ni.e.: mx var weight var2 weight2 ... var_n weight_n")
}
weights = append(weights, weight)
}
}
// Sample from mixture
xs, err := sample.Sample_mixture_serially_from_samplers(fs, weights, N_SAMPLES)
xs, err := sample.Mixture_serially_from_samplers(fs, weights, N_SAMPLES)
if err != nil {
return nil, printAndReturnErr(err.Error())
return nil, PrintAndReturnErr(err.Error())
}
return FilledSamples{xs: xs}, nil
}
/* Parser and repl */
func parseWordsErr(err_msg string) (string, Dist, error) {
return "", nil, printAndReturnErr(err_msg)
return "", nil, PrintAndReturnErr(err_msg)
}
func parseWordsIntoOpAndDist(words []string, vars map[string]Dist) (string, Dist, error) {
@ -396,7 +312,7 @@ func parseWordsIntoOpAndDist(words []string, vars map[string]Dist) (string, Dist
return parseWordsErr("Operator must have operand; can't operate on nothing")
case 1:
var_word, var_word_exists := vars[words[0]]
single_float, err1 := pretty.ParseFloat(words[0]) // abstract this away to search for K/M/B/T/etc.
single_float, err1 := ParseFloat(words[0]) // abstract this away to search for K/M/B/T/etc.
switch {
case var_word_exists:
dist = var_word
@ -406,8 +322,8 @@ func parseWordsIntoOpAndDist(words []string, vars map[string]Dist) (string, Dist
return parseWordsErr("Trying to operate on a scalar, but scalar is neither a float nor an assigned variable")
}
case 2:
new_low, err1 := pretty.ParseFloat(words[0])
new_high, err2 := pretty.ParseFloat(words[1])
new_low, err1 := ParseFloat(words[0])
new_high, err2 := ParseFloat(words[1])
switch {
case err1 != nil || err2 != nil:
return parseWordsErr("Trying to operate by a distribution, but distribution is not specified as two floats")
@ -422,8 +338,8 @@ func parseWordsIntoOpAndDist(words []string, vars map[string]Dist) (string, Dist
case 3:
switch {
case words[0] == "beta" || words[0] == "b":
a, err1 := pretty.ParseFloat(words[1])
b, err2 := pretty.ParseFloat(words[2])
a, err1 := ParseFloat(words[1])
b, err2 := ParseFloat(words[2])
if err1 != nil || err2 != nil {
return parseWordsErr("Trying to specify a beta distribution? Try beta 1 2")
}
@ -485,7 +401,7 @@ replForLoop:
stack.old_dist = INIT_DIST
fmt.Println()
case words[0] == "stats" || words[0] == "s":
prettyPrintStats(stack.old_dist)
PrettyPrintStats(stack.old_dist)
/* Variable assignment */
case words[0] == "=:" && len(words) == 2:
stack.vars[words[1]] = stack.old_dist
@ -493,7 +409,7 @@ replForLoop:
case words[0] == "=." && len(words) == 2:
stack.vars[words[1]] = stack.old_dist
fmt.Printf("%s ", words[1])
prettyPrintDist(stack.old_dist)
PrettyPrintDist(stack.old_dist)
stack.old_dist = INIT_DIST
default:
op, new_dist, err := parseWordsIntoOpAndDist(words, stack.vars)
@ -505,7 +421,7 @@ replForLoop:
stack.old_dist = combined_dist
}
}
prettyPrintDist(stack.old_dist)
PrettyPrintDist(stack.old_dist)
}
}

86
pretty/pretty.go → main/pretty.go
View File

@ -1,9 +1,11 @@
package pretty
package main
import (
"errors"
"fmt"
"git.nunosempere.com/NunoSempere/fermi/sample"
"math"
"sort"
"strconv"
)
@ -87,3 +89,85 @@ func ParseFloat(word string) (float64, error) {
}
}
/* Printers */
func PrettyPrintDist(dist Dist) {
switch v := dist.(type) {
case Lognormal:
fmt.Printf("=> ")
PrettyPrint2Floats(v.low, v.high)
fmt.Println()
case Beta:
fmt.Printf("=> beta ")
PrettyPrint2Floats(v.a, v.b)
fmt.Println()
case Scalar:
fmt.Printf("=> scalar ")
w := float64(v)
PrettyPrintFloat(w)
fmt.Println()
case FilledSamples:
n := len(v.xs)
sorted_xs := make([]float64, n)
copy(sorted_xs, v.xs)
sort.Slice(sorted_xs, func(i, j int) bool {
return sorted_xs[i] < sorted_xs[j]
})
low := sorted_xs[int(math.Round(float64(n)*0.05))]
high := sorted_xs[int(math.Round(float64(n)*0.95))]
fmt.Printf("=> ")
PrettyPrint2Floats(low, high)
fmt.Printf(" (")
PrettyPrintInt(N_SAMPLES)
fmt.Printf(" samples)")
fmt.Println()
default:
fmt.Printf("%v\n", v)
}
}
func PrettyPrintStats(dist Dist) {
xs := sample.Serially(dist.Sampler, N_SAMPLES)
n := len(xs)
mean := 0.0
for i := 0; i < n; i++ {
mean += xs[i]
}
mean /= float64(n)
fmt.Printf("Mean: %f\n", mean)
stdev := 0.0
for i := 0; i < n; i++ {
stdev += math.Pow(xs[i]-mean, 2)
}
stdev = math.Sqrt(stdev / float64(n))
fmt.Printf("Stdev: %f\n", stdev)
sorted_xs := make([]float64, n)
copy(sorted_xs, xs)
sort.Slice(sorted_xs, func(i, j int) bool {
return sorted_xs[i] < sorted_xs[j]
})
print_ci := func(ci float64, prefix string) {
x := sorted_xs[int(math.Round(float64(n)*ci))]
fmt.Printf("%s%f\n", prefix, x)
}
print_ci(0.01, "ci 1%: ")
print_ci(0.05, "ci 5%: ")
print_ci(0.10, "ci 10%: ")
print_ci(0.25, "ci 25%: ")
print_ci(0.50, "ci 50%: ")
print_ci(0.75, "ci 75%: ")
print_ci(0.90, "ci 90%: ")
print_ci(0.95, "ci 95%: ")
print_ci(0.99, "ci 99%: ")
}
func PrintAndReturnErr(err_msg string) error {
fmt.Println(err_msg)
fmt.Println("Type \"help\" (without quotes) to see a pseudogrammar and examples")
return errors.New(err_msg)
}

2
makefile
View File

@ -1,5 +1,5 @@
build:
go build fermi.go
go build main/fermi.go main/pretty.go
run:
go run fermi.go

173
pretty.go Normal file
View File

@ -0,0 +1,173 @@
package main
import (
"errors"
"fmt"
"git.nunosempere.com/NunoSempere/fermi/sample"
"math"
"sort"
"strconv"
)
func PrettyPrintInt(n int) {
switch {
case math.Abs(float64(n)) >= 1_000_000_000_000:
fmt.Printf("%.2fT", float64(n)/1_000_000_000_000.0)
case math.Abs(float64(n)) >= 1_000_000_000:
fmt.Printf("%.2fB", float64(n)/1_000_000_000.0)
case math.Abs(float64(n)) >= 1_000_000:
fmt.Printf("%.2fM", float64(n)/1_000_000.0)
case math.Abs(float64(n)) >= 1_000:
fmt.Printf("%.2fK", float64(n)/1_000.0)
default:
fmt.Printf("%d", n)
}
}
func PrettyPrintFloat(f float64) {
switch {
case math.Abs(f) >= 1_000_000_000_000:
fmt.Printf("%.2fT", f/1_000_000_000_000)
case math.Abs(f) >= 1_000_000_000:
fmt.Printf("%.2fB", f/1_000_000_000)
case math.Abs(f) >= 1_000_000:
fmt.Printf("%.2fM", f/1_000_000)
case math.Abs(f) >= 1_000:
fmt.Printf("%.2fK", f/1_000)
case math.Abs(f) <= 0.0001:
fmt.Printf("%.6f", f)
case math.Abs(f) <= 0.001:
fmt.Printf("%.5f", f)
case math.Abs(f) <= 0.01:
fmt.Printf("%.4f", f)
case math.Abs(f) <= 0.1:
fmt.Printf("%.3f", f)
default:
fmt.Printf("%.2f", f)
}
}
func PrettyPrint2Floats(low float64, high float64) {
PrettyPrintFloat(low)
fmt.Printf(" ")
PrettyPrintFloat(high)
}
func multiplyOrPassThroughError(a float64, b float64, err error) (float64, error) {
if err != nil {
return b, err
} else {
return a * b, nil
}
}
func ParseFloat(word string) (float64, error) {
// l = len(word) // assuming no UTF stuff
switch len(word) {
case 0:
return 0, errors.New("String to be parsed into float must not be the empty string")
case 1:
return strconv.ParseFloat(word, 64)
}
n := len(word) - 1
f, err := strconv.ParseFloat(word[:n], 64)
switch word[n] {
case '%':
return multiplyOrPassThroughError(0.01, f, err)
case 'K':
return multiplyOrPassThroughError(1_000, f, err)
case 'M':
return multiplyOrPassThroughError(1_000_000, f, err)
case 'B':
return multiplyOrPassThroughError(1_000_000_000, f, err)
case 'T':
return multiplyOrPassThroughError(1_000_000_000_000, f, err)
default:
return strconv.ParseFloat(word, 64)
}
}
/* Printers */
func PrettyPrintDist(dist Dist) {
switch v := dist.(type) {
case Lognormal:
fmt.Printf("=> ")
PrettyPrint2Floats(v.low, v.high)
fmt.Println()
case Beta:
fmt.Printf("=> beta ")
PrettyPrint2Floats(v.a, v.b)
fmt.Println()
case Scalar:
fmt.Printf("=> scalar ")
w := float64(v)
PrettyPrintFloat(w)
fmt.Println()
case FilledSamples:
n := len(v.xs)
sorted_xs := make([]float64, n)
copy(sorted_xs, v.xs)
sort.Slice(sorted_xs, func(i, j int) bool {
return sorted_xs[i] < sorted_xs[j]
})
low := sorted_xs[int(math.Round(float64(n)*0.05))]
high := sorted_xs[int(math.Round(float64(n)*0.95))]
fmt.Printf("=> ")
PrettyPrint2Floats(low, high)
fmt.Printf(" (")
PrettyPrintInt(N_SAMPLES)
fmt.Printf(" samples)")
fmt.Println()
default:
fmt.Printf("%v\n", v)
}
}
func PrettyPrintStats(dist Dist) {
xs := sample.Sample_serially(dist.Sampler, N_SAMPLES)
n := len(xs)
mean := 0.0
for i := 0; i < n; i++ {
mean += xs[i]
}
mean /= float64(n)
fmt.Printf("Mean: %f\n", mean)
stdev := 0.0
for i := 0; i < n; i++ {
stdev += math.Pow(xs[i]-mean, 2)
}
stdev = math.Sqrt(stdev / float64(n))
fmt.Printf("Stdev: %f\n", stdev)
sorted_xs := make([]float64, n)
copy(sorted_xs, xs)
sort.Slice(sorted_xs, func(i, j int) bool {
return sorted_xs[i] < sorted_xs[j]
})
print_ci := func(ci float64, prefix string) {
x := sorted_xs[int(math.Round(float64(n)*ci))]
fmt.Printf("%s%f\n", prefix, x)
}
print_ci(0.01, "ci 1%: ")
print_ci(0.05, "ci 5%: ")
print_ci(0.10, "ci 10%: ")
print_ci(0.25, "ci 25%: ")
print_ci(0.50, "ci 50%: ")
print_ci(0.75, "ci 75%: ")
print_ci(0.90, "ci 90%: ")
print_ci(0.95, "ci 95%: ")
print_ci(0.99, "ci 99%: ")
}
func PrintAndReturnErr(err_msg string) error {
fmt.Println(err_msg)
fmt.Println("Type \"help\" (without quotes) to see a pseudogrammar and examples")
return errors.New(err_msg)
}

58
sample/sample.go
View File

@ -17,49 +17,49 @@ type func64i = func(int, State) float64
var global_state = rand.New(rand.NewPCG(uint64(1), uint64(2)))
func Sample_int(n int, r State) int {
func Int(n int, r State) int {
return r.IntN(n)
}
func Sample_unit_uniform(r State) float64 {
func Unit_uniform(r State) float64 {
return r.Float64()
}
func Sample_unit_normal(r State) float64 {
func Unit_normal(r State) float64 {
return r.NormFloat64()
}
func Sample_uniform(start float64, end float64, r State) float64 {
return Sample_unit_uniform(r)*(end-start) + start
func Uniform(start float64, end float64, r State) float64 {
return Unit_uniform(r)*(end-start) + start
}
func Sample_normal(mean float64, sigma float64, r State) float64 {
return mean + Sample_unit_normal(r)*sigma
func Normal(mean float64, sigma float64, r State) float64 {
return mean + Unit_normal(r)*sigma
}
func Sample_lognormal(logmean float64, logstd float64, r State) float64 {
return (math.Exp(Sample_normal(logmean, logstd, r)))
func Lognormal(logmean float64, logstd float64, r State) float64 {
return (math.Exp(Normal(logmean, logstd, r)))
}
func Sample_normal_from_90_ci(low float64, high float64, r State) float64 {
func Normal_from_90_ci(low float64, high float64, r State) float64 {
var normal90 float64 = 1.6448536269514727
var mean float64 = (high + low) / 2.0
var std float64 = (high - low) / (2.0 * normal90)
return Sample_normal(mean, std, r)
return Normal(mean, std, r)
}
func Sample_to(low float64, high float64, r State) float64 {
func To(low float64, high float64, r State) float64 {
// Given a (positive) 90% confidence interval,
// returns a sample from a lognorma with a matching 90% c.i.
// Key idea: If we want a lognormal with 90% confidence interval [a, b]
// we need but get a normal with 90% confidence interval [log(a), log(b)].
// Then see code for Sample_normal_from_90_ci
// Then see code for Normal_from_90_ci
var loglow float64 = math.Log(low)
var loghigh float64 = math.Log(high)
return math.Exp(Sample_normal_from_90_ci(loglow, loghigh, r))
return math.Exp(Normal_from_90_ci(loglow, loghigh, r))
}
func Sample_gamma(alpha float64, r State) float64 {
func Gamma(alpha float64, r State) float64 {
// a simple method for generating gamma variables, marsaglia and wan tsang, 2001
// https://dl.acm.org/doi/pdf/10.1145/358407.358414
@ -80,7 +80,7 @@ func Sample_gamma(alpha float64, r State) float64 {
InnerLoop:
for {
x = Sample_unit_normal(r)
x = Unit_normal(r)
v = 1.0 + c*x
if v > 0.0 {
break InnerLoop
@ -88,7 +88,7 @@ func Sample_gamma(alpha float64, r State) float64 {
}
v = v * v * v
u = Sample_unit_uniform(r)
u = Unit_uniform(r)
if u < 1.0-0.0331*(x*x*x*x) { // Condition 1
// the 0.0331 doesn't inspire much confidence
@ -105,17 +105,17 @@ func Sample_gamma(alpha float64, r State) float64 {
}
} else {
return Sample_gamma(1.0+alpha, r) * math.Pow(Sample_unit_uniform(r), 1.0/alpha)
return Gamma(1.0+alpha, r) * math.Pow(Unit_uniform(r), 1.0/alpha)
}
}
func Sample_beta(a float64, b float64, r State) float64 {
gamma_a := Sample_gamma(a, r)
gamma_b := Sample_gamma(b, r)
func Beta(a float64, b float64, r State) float64 {
gamma_a := Gamma(a, r)
gamma_b := Gamma(b, r)
return gamma_a / (gamma_a + gamma_b)
}
func Sample_mixture_once(fs []func64, weights []float64, r State) float64 {
func Mixture_once(fs []func64, weights []float64, r State) float64 {
// fmt.Println("weights initially: ", weights)
var sum_weights float64 = 0
@ -149,7 +149,7 @@ func Sample_mixture_once(fs []func64, weights []float64, r State) float64 {
}
func Sample_serially(f func64i, n_samples int) []float64 {
func Serially(f func64i, n_samples int) []float64 {
xs := make([]float64, n_samples)
// var global_state = rand.New(rand.NewPCG(uint64(1), uint64(2)))
for i := 0; i < n_samples; i++ {
@ -158,7 +158,7 @@ func Sample_serially(f func64i, n_samples int) []float64 {
return xs
}
func Sample_mixture_serially_from_samples(fs [][]float64, weights []float64, n_samples int) ([]float64, error) {
func Mixture_serially_from_samples(fs [][]float64, weights []float64, n_samples int) ([]float64, error) {
// Checks
if len(weights) != len(fs) {
@ -206,7 +206,7 @@ func Sample_mixture_serially_from_samples(fs [][]float64, weights []float64, n_s
return xs, nil
}
func Sample_mixture_serially_from_samplers(fs []func64i, weights []float64, n_samples int) ([]float64, error) {
func Mixture_serially_from_samplers(fs []func64i, weights []float64, n_samples int) ([]float64, error) {
// Checks
if len(weights) != len(fs) {
@ -247,7 +247,7 @@ func Sample_mixture_serially_from_samplers(fs []func64i, weights []float64, n_sa
return xs, nil
}
func Sample_parallel(f func64, n_samples int) []float64 {
func Parallel(f func64, n_samples int) []float64 {
var num_threads = 16
var xs = make([]float64, n_samples)
var wg sync.WaitGroup
@ -283,9 +283,9 @@ func main() {
Sample_many := func(r State) float64 { return Sample_to(2, 10, r) }
fs := [4](func64){Sample_0, Sample_1, Sample_few, Sample_many}
model := func(r State) float64 { return Sample_mixture(fs[0:], ws[0:], r) }
model := func(r State) float64 { return Mixture(fs[0:], ws[0:], r) }
n_samples := 1_000_000
xs := Sample_parallel(model, n_samples)
xs := Parallel(model, n_samples)
var avg float64 = 0
for _, x := range xs {
avg += x
@ -298,7 +298,7 @@ func main() {
var r = rand.New(rand.NewPCG(uint64(1), uint64(2)))
var avg float64 = 0
for i := 0; i < n_samples; i++ {
avg += Sample_mixture(fs[0:], ws[0:], r)
avg += Mixture(fs[0:], ws[0:], r)
}
avg = avg / float64(n_samples)
fmt.Printf("Average: %v\n", avg)