Compare commits
6 Commits
e93316446c
...
b645f6648a
Author | SHA1 | Date | |
---|---|---|---|
b645f6648a | |||
a36808d8e6 | |||
dcd0db36eb | |||
3e4df77541 | |||
0c03b7e6bf | |||
5e0020a112 |
1
.gitignore
vendored
Normal file
1
.gitignore
vendored
Normal file
|
@ -0,0 +1 @@
|
|||
./fermi
|
23
main/error.go
Normal file
23
main/error.go
Normal 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())
|
||||
}
|
||||
}
|
|
@ -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)
|
||||
}
|
||||
}
|
||||
|
|
@ -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
2
makefile
|
@ -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
|
||||
|
|
|
@ -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)
|
||||
|
|
|
@ -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) {
|
||||
|
|
Loading…
Reference in New Issue
Block a user