Compare commits

...

6 Commits

8 changed files with 222 additions and 203 deletions

1
.gitignore vendored Normal file
View File

@ -0,0 +1 @@
./fermi

BIN
fermi

Binary file not shown.

23
main/error.go Normal file
View File

@ -0,0 +1,23 @@
package main
import (
"fmt"
)
type FermiError struct {
Code string
Msg string
}
func (e FermiError) Error() string {
return e.Msg
}
func PrintError(e error) {
if fermi_error, ok := e.(FermiError); ok {
fmt.Printf("Error: %s: %s\n", fermi_error.Code, fermi_error.Msg)
fmt.Printf("Type \"help\" or \"h\" to see grammar and command flags\n")
} else {
fmt.Printf("Error: %s\n", e.Error())
}
}

View File

@ -2,14 +2,11 @@ package main
import (
"bufio"
"errors"
"flag"
"fmt"
"git.nunosempere.com/NunoSempere/fermi/pretty"
"git.nunosempere.com/NunoSempere/fermi/sample"
"math"
"os"
"sort"
"strings"
)
@ -46,11 +43,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 {
@ -86,7 +83,7 @@ const HELP_MSG = "1. Grammar:\n" +
" + 1 10\n" +
" * beta 1 10\n" +
" 1 10 (multiplication taken as default operation)\n" +
" =: x\n" +
" =: x\n" +
" .\n" +
" 1 100\n" +
" + x\n" +
@ -96,6 +93,9 @@ const HELP_MSG = "1. Grammar:\n" +
" 1 10\n" +
" + beta 1 100\n" +
" )\n" +
" =. y\n" +
" mx x 1 y 2.33\n" +
" + mx x 30% y 70%\n" +
" exit\n" +
"\n" +
"2. Command flags:\n" +
@ -105,101 +105,19 @@ const HELP_MSG = "1. Grammar:\n" +
" Specifies a file with a model to run\n" +
" -n int\n" +
" Specifies the number of samples to draw when using samples (default 100000)\n" +
" -h Shows help message\n"
" -h Shows help message"
const NORMAL90CONFIDENCE = 1.6448536269514727
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++ {
@ -210,16 +128,16 @@ func operateDistsAsSamples(dist1 Dist, dist2 Dist, op string) (Dist, error) {
if ys[0] != 0 {
zs[i] = xs[i] / ys[i]
} else {
fmt.Println("Error: When dividing as samples, division by zero")
return nil, errors.New("Division by zero")
// fmt.Println("Error: When dividing as samples, division by zero")
return nil, FermiError{Code: "Division by zero", Msg: "When operating on samples, division by zero"}
}
case "+":
zs[i] = xs[i] + ys[i]
case "-":
zs[i] = xs[i] - ys[i]
default:
fmt.Println("Error: Operation not recognized")
return nil, errors.New("Operation not recognized")
// fmt.Println("Error: Operation not recognized")
return nil, FermiError{Code: "Unknown operation", Msg: "When operating on samples, operation not recognized"}
}
}
@ -244,15 +162,11 @@ func multiplyLogDists(l1 Lognormal, l2 Lognormal) Lognormal {
}
func multiplyBetaDists(beta1 Beta, beta2 Beta) Beta {
return Beta{a: beta1.a + beta2.a, b: beta1.b + beta2.b}
}
func multiplyLogDistAndScalar(l Lognormal, s Scalar) (Dist, error) {
if s == 0.0 {
return Scalar(0.0), nil
} else if s < 0.0 {
return operateDistsAsSamples(s, l, "+")
return operateDistsAsSamples(s, l, "*")
} else {
return multiplyLogDists(l, Lognormal{low: float64(s), high: float64(s)}), nil
}
@ -285,11 +199,6 @@ func multiplyDists(old_dist Dist, new_dist Dist) (Dist, error) {
return Scalar(float64(o) * float64(n)), nil
}
}
case Beta:
switch n := new_dist.(type) {
case Beta:
return multiplyBetaDists(o, n), nil
}
}
return operateDistsAsSamples(old_dist, new_dist, "*")
}
@ -303,14 +212,12 @@ func divideDists(old_dist Dist, new_dist Dist) (Dist, error) {
switch n := new_dist.(type) {
case Lognormal:
if n.high == 0 || n.low == 0 {
fmt.Println("Error: Can't divide by 0.0")
return nil, errors.New("Error: division by zero")
return nil, FermiError{Code: "Division by zero", Msg: "When operating two lognormals, one of the parameters is zero, which would result in division by zero"}
}
return multiplyLogDists(o, Lognormal{low: 1.0 / n.high, high: 1.0 / n.low}), nil
case Scalar:
if n == 0.0 {
fmt.Println("Error: Can't divide by 0.0")
return nil, errors.New("Error: division by zero scalar")
return nil, FermiError{Code: "Division by zero", Msg: "When operating a lognormal with a scalar, trying to divide but the scalar is zero"}
}
return multiplyLogDistAndScalar(o, Scalar(1.0/n))
}
@ -322,8 +229,7 @@ func divideDists(old_dist Dist, new_dist Dist) (Dist, error) {
return multiplyLogDistAndScalar(Lognormal{low: 1.0 / n.high, high: 1.0 / n.low}, o)
case Scalar:
if n == 0.0 {
fmt.Println("Error: Can't divide by 0.0")
return nil, errors.New("Error: division by zero scalar")
return nil, FermiError{Code: "Division by zero", Msg: "When operating two scalars, trying to divide but the divisor is is zero"}
}
return Scalar(float64(o) / float64(n)), nil
}
@ -344,7 +250,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, FermiError{Code: "Unknown operation", Msg: "When operating distributions, operation not recognized"}
}
}
@ -352,9 +258,17 @@ func operateDists(old_dist Dist, new_dist Dist, op string) (Dist, error) {
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")
switch {
case len(words) < 1:
return nil, FermiError{Code: "Not a mixture", Msg: "Input can't be a mixture, since it doesn't have enough words"}
case words[0] != "mx":
return nil, FermiError{Code: "Not a mixture", Msg: "Input can't be a mixture, since it is not preceded by the mx keyword"}
case len(words)%2 != 1:
return nil, FermiError{Code: "Not a mixture", Msg: "When parsing a mixture, input doesn't have equal number of variables and weights. \nMixture syntax: \nmx x 20% y 70% z 10%\ni.e.: mx var weight var2 weight2 ... var_n weight_n"}
case len(words) < 5:
return nil, FermiError{Code: "Not a mixture", Msg: "When parsing a mixture, not enough words. \nMixture syntax: \nmx x 20% y 70% z 10%\ni.e.: mx var weight var2 weight2 ... var_n weight_n"}
}
words = words[1:] // crop "mx" at the beginning
var fs []func(int, sample.State) float64
var weights []float64
@ -363,49 +277,49 @@ 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, FermiError{Code: "Not a mixture variable", Msg: "When parsing a mixture, 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, FermiError{Code: "Not a mixture weight", Msg: "When parsing a mixture, expected mixture weight 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"}
}
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)
}
func parseWordsIntoOpAndDist(words []string, vars map[string]Dist) (string, Dist, error) {
parseWordsErr := func(msg string) (string, Dist, error) {
return "", nil, FermiError{Code: "Malformed input", Msg: "When parsing words into operator and distribution: " + msg}
}
op := ""
var dist Dist
switch words[0] {
case "*", "/", "+", "-", "mx":
case "*", "/", "+", "-":
op = words[0]
words = words[1:]
default:
op = "*"
}
switch len(words) {
case 0:
switch {
case len(words) == 0:
return parseWordsErr("Operator must have operand; can't operate on nothing")
case 1:
case len(words) == 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
@ -414,9 +328,9 @@ func parseWordsIntoOpAndDist(words []string, vars map[string]Dist) (string, Dist
case err1 != nil && !var_word_exists:
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])
case len(words) == 2:
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")
@ -428,11 +342,11 @@ func parseWordsIntoOpAndDist(words []string, vars map[string]Dist) (string, Dist
return parseWordsErr("Trying to parse two floats as a lognormal, but the first number is larger than the second number")
}
dist = Lognormal{low: new_low, high: new_high}
case 3:
case len(words) == 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")
}
@ -440,18 +354,12 @@ func parseWordsIntoOpAndDist(words []string, vars map[string]Dist) (string, Dist
default:
return parseWordsErr("Input not understood or not implemented yet")
}
default:
switch op {
case "mx":
tmp, err := parseMixture(words, vars)
if err != nil {
return parseWordsErr("Error parsing a mixture: " + err.Error())
}
dist = tmp
op = "*"
default:
return parseWordsErr("Input not understood or not implemented yet")
case len(words) >= 4: // four or more words
if words[0] == "mx" {
dist, err := parseMixture(words, vars)
return op, dist, err
}
return parseWordsErr("Input not understood or not implemented yet")
}
return op, dist, nil
}
@ -472,17 +380,11 @@ replForLoop:
words := strings.Split(new_line_trimmed, " ")
switch {
case strings.TrimSpace(new_line_trimmed) == "": /* Empty line case */
/* Empty line */
case strings.TrimSpace(new_line_trimmed) == "":
continue replForLoop
/* Parenthesis */
case len(words) == 2 && (words[0] == "*" || words[0] == "+" || words[0] == "-" || words[0] == "/") && words[1] == "(":
new_stack := runRepl(Stack{old_dist: INIT_DIST, vars: stack.vars}, reader, echo_flag)
combined_dist, err := operateDists(stack.old_dist, new_stack.old_dist, words[0])
if err == nil {
stack.old_dist = combined_dist
}
case len(words) == 1 && words[0] == ")":
return stack
/* Special operations */
case words[0] == "exit" || words[0] == "e":
os.Exit(0)
@ -494,7 +396,8 @@ 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
@ -502,19 +405,39 @@ 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
/* Parenthesis */
case len(words) == 2 && (words[0] == "*" || words[0] == "+" || words[0] == "-" || words[0] == "/") && words[1] == "(":
new_stack := runRepl(Stack{old_dist: INIT_DIST, vars: stack.vars}, reader, echo_flag)
combined_dist, err := operateDists(stack.old_dist, new_stack.old_dist, words[0])
if err != nil {
PrintError(err)
} else {
stack.old_dist = combined_dist
}
case len(words) == 1 && words[0] == ")":
return stack
/* Bread and butter distribution operations */
default:
op, new_dist, err := parseWordsIntoOpAndDist(words, stack.vars)
if err != nil {
PrintError(err)
PrettyPrintDist(stack.old_dist)
continue replForLoop
}
combined_dist, err := operateDists(stack.old_dist, new_dist, op)
if err == nil {
stack.old_dist = combined_dist
if err != nil {
PrintError(err)
PrettyPrintDist(stack.old_dist)
continue replForLoop
}
stack.old_dist = combined_dist
}
prettyPrintDist(stack.old_dist)
PrettyPrintDist(stack.old_dist)
}
}

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)
}

View File

@ -1,5 +1,5 @@
build:
go build fermi.go
go build main/fermi.go main/pretty.go main/error.go
run:
go run fermi.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)

View File

@ -173,10 +173,6 @@ func multiplyLogDists(l1 Lognormal, l2 Lognormal) Lognormal {
}
func multiplyBetaDists(beta1 Beta, beta2 Beta) Beta {
return Beta{a: beta1.a + beta2.a, b: beta1.b + beta2.b}
}
func multiplyDists(old_dist Dist, new_dist Dist) (Dist, error) {
switch o := old_dist.(type) {
@ -205,16 +201,8 @@ func multiplyDists(old_dist Dist, new_dist Dist) (Dist, error) {
return operateDistsAsSamples(old_dist, new_dist, "*")
}
}
case Beta:
switch n := new_dist.(type) {
case Beta:
return multiplyBetaDists(o, n), nil
default:
return operateDistsAsSamples(old_dist, new_dist, "*")
}
default:
return operateDistsAsSamples(old_dist, new_dist, "*")
}
return operateDistsAsSamples(old_dist, new_dist, "*")
}
func divideDists(old_dist Dist, new_dist Dist) (Dist, error) {