Compare commits
12 Commits
a4fdbc1e2c
...
308eecba98
Author | SHA1 | Date | |
---|---|---|---|
308eecba98 | |||
b1a58f9b74 | |||
015d33adca | |||
fa714f91ae | |||
0bcd9f9b62 | |||
66cbb4132f | |||
a4389e605f | |||
0d31a4f1ab | |||
f1834341a9 | |||
ccad14b318 | |||
56ab018469 | |||
79dfcf79db |
22
README.md
22
README.md
|
@ -2,6 +2,8 @@
|
||||||
|
|
||||||
squiggle.c is a self-contained C99 library that provides functions for simple Monte Carlo estimation, based on [Squiggle](https://www.squiggle-language.com/).
|
squiggle.c is a self-contained C99 library that provides functions for simple Monte Carlo estimation, based on [Squiggle](https://www.squiggle-language.com/).
|
||||||
|
|
||||||
|
![](./core.png)
|
||||||
|
|
||||||
## Why C?
|
## Why C?
|
||||||
|
|
||||||
- Because it is fast
|
- Because it is fast
|
||||||
|
@ -12,7 +14,7 @@ squiggle.c is a self-contained C99 library that provides functions for simple Mo
|
||||||
- Because if you can implement something in C, you can implement it anywhere else
|
- Because if you can implement something in C, you can implement it anywhere else
|
||||||
- Because it can be made faster if need be
|
- Because it can be made faster if need be
|
||||||
- e.g., with a multi-threading library like OpenMP,
|
- e.g., with a multi-threading library like OpenMP,
|
||||||
- or by implementing faster but more complex algorithms
|
- o by implementing faster but more complex algorithms
|
||||||
- or more simply, by inlining the sampling functions (adding an `inline` directive before their function declaration)
|
- or more simply, by inlining the sampling functions (adding an `inline` directive before their function declaration)
|
||||||
- **Because there are few abstractions between it and machine code** (C => assembly => machine code with gcc, or C => machine code, with tcc), leading to fewer errors beyond the programmer's control.
|
- **Because there are few abstractions between it and machine code** (C => assembly => machine code with gcc, or C => machine code, with tcc), leading to fewer errors beyond the programmer's control.
|
||||||
|
|
||||||
|
@ -34,7 +36,7 @@ You can follow some example usage in the examples/ folder
|
||||||
|
|
||||||
### squiggle.c is short
|
### squiggle.c is short
|
||||||
|
|
||||||
[squiggle.c](squiggle.c) is less than 500 lines of C. The reader could just read it and grasp its contents.
|
[squiggle.c](squiggle.c) is less than 600 lines of C, with a core of <250 lines. The reader could just read it and grasp its contents.
|
||||||
|
|
||||||
### Core strategy
|
### Core strategy
|
||||||
|
|
||||||
|
@ -341,4 +343,18 @@ It emits one warning about something I already took care of, so by default I've
|
||||||
- [x] Have some more complicated & realistic example
|
- [x] Have some more complicated & realistic example
|
||||||
- [x] Add summarization functions: 90% ci (or all c.i.?)
|
- [x] Add summarization functions: 90% ci (or all c.i.?)
|
||||||
- [x] Link to the examples in the examples section.
|
- [x] Link to the examples in the examples section.
|
||||||
- [x] Add a few functions for doing simple algebra on normals, and lognormals?
|
- [x] Add a few functions for doing simple algebra on normals, and lognormals
|
||||||
|
- [x] Add prototypes
|
||||||
|
- [x] Use named structs
|
||||||
|
- [x] Add to header file
|
||||||
|
- [x] Provide example algebra
|
||||||
|
- [x] Add conversion between 90% ci and parameters.
|
||||||
|
- [x] Use that conversion in conjuction with small algebra.
|
||||||
|
- [x] Consider ergonomics of using ci instead of c_i
|
||||||
|
- [x] use named struct instead
|
||||||
|
- [x] demonstrate and document feeding a struct directly to a function; my_function((struct c_i){.low = 1, .high = 2});
|
||||||
|
- [ ] Consider desirability of defining shortcuts for those functions. Adds a level of magic, though.
|
||||||
|
- [ ] Test results
|
||||||
|
- [x] Move to own file? Or signpost in file? => signposted in file.
|
||||||
|
- [ ] Disambiguate sample_laplace--successes vs failures || successes vs total trials as two distinct and differently named functions
|
||||||
|
- [ ] Write twitter thread.
|
||||||
|
|
|
@ -1,7 +1,7 @@
|
||||||
#include <stdint.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <stdio.h>
|
|
||||||
#include "../../squiggle.h"
|
#include "../../squiggle.h"
|
||||||
|
#include <stdint.h>
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
|
||||||
// Estimate functions
|
// Estimate functions
|
||||||
double sample_0(uint64_t* seed)
|
double sample_0(uint64_t* seed)
|
||||||
|
@ -24,10 +24,11 @@ double sample_many(uint64_t* seed)
|
||||||
return sample_to(2, 10, seed);
|
return sample_to(2, 10, seed);
|
||||||
}
|
}
|
||||||
|
|
||||||
int main(){
|
int main()
|
||||||
|
{
|
||||||
// set randomness seed
|
// set randomness seed
|
||||||
uint64_t* seed = malloc(sizeof(uint64_t));
|
uint64_t* seed = malloc(sizeof(uint64_t));
|
||||||
*seed = 1000; // xorshift can't start with 0
|
*seed = 1000; // xorshift can't start with 0
|
||||||
|
|
||||||
double p_a = 0.8;
|
double p_a = 0.8;
|
||||||
double p_b = 0.5;
|
double p_b = 0.5;
|
||||||
|
@ -37,18 +38,18 @@ int main(){
|
||||||
double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
|
double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
|
||||||
double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
|
double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
|
||||||
|
|
||||||
int n_samples = 1000000;
|
int n_samples = 1000000;
|
||||||
double* result_many = (double *) malloc(n_samples * sizeof(double));
|
double* result_many = (double*)malloc(n_samples * sizeof(double));
|
||||||
for(int i=0; i<n_samples; i++){
|
for (int i = 0; i < n_samples; i++) {
|
||||||
result_many[i] = sample_mixture(samplers, weights, n_dists, seed);
|
result_many[i] = sample_mixture(samplers, weights, n_dists, seed);
|
||||||
}
|
}
|
||||||
printf("Mean: %f\n", array_mean(result_many, n_samples));
|
printf("Mean: %f\n", array_mean(result_many, n_samples));
|
||||||
|
|
||||||
// printf("result_many: [");
|
|
||||||
// for(int i=0; i<100; i++){
|
|
||||||
// printf("%.2f, ", result_many[i]);
|
|
||||||
// }
|
|
||||||
// printf("]\n");
|
|
||||||
|
|
||||||
free(seed);
|
// printf("result_many: [");
|
||||||
|
// for(int i=0; i<100; i++){
|
||||||
|
// printf("%.2f, ", result_many[i]);
|
||||||
|
// }
|
||||||
|
// printf("]\n");
|
||||||
|
|
||||||
|
free(seed);
|
||||||
}
|
}
|
||||||
|
|
Binary file not shown.
|
@ -60,9 +60,9 @@ int main()
|
||||||
}
|
}
|
||||||
printf("... ]\n");
|
printf("... ]\n");
|
||||||
|
|
||||||
struct c_i c_i_90 = get_90_confidence_interval(mixture, seed);
|
ci ci_90 = get_90_confidence_interval(mixture, seed);
|
||||||
printf("mean: %f\n", array_mean(mixture_result, n));
|
printf("mean: %f\n", array_mean(mixture_result, n));
|
||||||
printf("90%% confidence interval: [%f, %f]\n", c_i_90.low, c_i_90.high);
|
printf("90%% confidence interval: [%f, %f]\n", ci_90.low, ci_90.high);
|
||||||
|
|
||||||
free(seed);
|
free(seed);
|
||||||
}
|
}
|
||||||
|
|
Binary file not shown.
|
@ -46,8 +46,8 @@ int main()
|
||||||
|
|
||||||
// Before a first nuclear collapse
|
// Before a first nuclear collapse
|
||||||
printf("## Before the first nuclear collapse\n");
|
printf("## Before the first nuclear collapse\n");
|
||||||
struct c_i c_i_90_2023 = get_90_confidence_interval(yearly_probability_nuclear_collapse_2023, seed);
|
ci ci_90_2023 = get_90_confidence_interval(yearly_probability_nuclear_collapse_2023, seed);
|
||||||
printf("90%% confidence interval: [%f, %f]\n", c_i_90_2023.low, c_i_90_2023.high);
|
printf("90%% confidence interval: [%f, %f]\n", ci_90_2023.low, ci_90_2023.high);
|
||||||
|
|
||||||
double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * num_samples);
|
double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * num_samples);
|
||||||
for (int i = 0; i < num_samples; i++) {
|
for (int i = 0; i < num_samples; i++) {
|
||||||
|
@ -57,8 +57,8 @@ int main()
|
||||||
|
|
||||||
// After the first nuclear collapse
|
// After the first nuclear collapse
|
||||||
printf("\n## After the first nuclear collapse\n");
|
printf("\n## After the first nuclear collapse\n");
|
||||||
struct c_i c_i_90_2070 = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_example, seed);
|
ci ci_90_2070 = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_example, seed);
|
||||||
printf("90%% confidence interval: [%f, %f]\n", c_i_90_2070.low, c_i_90_2070.high);
|
printf("90%% confidence interval: [%f, %f]\n", ci_90_2070.low, ci_90_2070.high);
|
||||||
|
|
||||||
double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * num_samples);
|
double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * num_samples);
|
||||||
for (int i = 0; i < num_samples; i++) {
|
for (int i = 0; i < num_samples; i++) {
|
||||||
|
@ -68,8 +68,8 @@ int main()
|
||||||
|
|
||||||
// After the first nuclear collapse (antiinductive)
|
// After the first nuclear collapse (antiinductive)
|
||||||
printf("\n## After the first nuclear collapse (antiinductive)\n");
|
printf("\n## After the first nuclear collapse (antiinductive)\n");
|
||||||
struct c_i c_i_90_antiinductive = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_antiinductive, seed);
|
ci ci_90_antiinductive = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_antiinductive, seed);
|
||||||
printf("90%% confidence interval: [%f, %f]\n", c_i_90_antiinductive.low, c_i_90_antiinductive.high);
|
printf("90%% confidence interval: [%f, %f]\n", ci_90_antiinductive.low, ci_90_antiinductive.high);
|
||||||
|
|
||||||
double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * num_samples);
|
double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * num_samples);
|
||||||
for (int i = 0; i < num_samples; i++) {
|
for (int i = 0; i < num_samples; i++) {
|
||||||
|
|
BIN
examples/11_algebra/example
Executable file
BIN
examples/11_algebra/example
Executable file
Binary file not shown.
26
examples/11_algebra/example.c
Normal file
26
examples/11_algebra/example.c
Normal file
|
@ -0,0 +1,26 @@
|
||||||
|
#include "../../squiggle.h"
|
||||||
|
#include <math.h>
|
||||||
|
#include <stdint.h>
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
|
||||||
|
int main()
|
||||||
|
{
|
||||||
|
// set randomness seed
|
||||||
|
uint64_t* seed = malloc(sizeof(uint64_t));
|
||||||
|
*seed = 1000; // xorshift can't start with 0
|
||||||
|
|
||||||
|
normal_params n1 = { .mean = 1.0, .std = 3.0 };
|
||||||
|
normal_params n2 = { .mean = 2.0, .std = 4.0 };
|
||||||
|
normal_params sn = algebra_sum_normals(n1, n2);
|
||||||
|
printf("The sum of Normal(%f, %f) and Normal(%f, %f) is Normal(%f, %f)\n",
|
||||||
|
n1.mean, n1.std, n2.mean, n2.std, sn.mean, sn.std);
|
||||||
|
|
||||||
|
lognormal_params ln1 = { .logmean = 1.0, .logstd = 3.0 };
|
||||||
|
lognormal_params ln2 = { .logmean = 2.0, .logstd = 4.0 };
|
||||||
|
lognormal_params sln = algebra_product_lognormals(ln1, ln2);
|
||||||
|
printf("The product of Lognormal(%f, %f) and Lognormal(%f, %f) is Lognormal(%f, %f)\n",
|
||||||
|
ln1.logmean, ln1.logstd, ln2.logmean, ln2.logstd, sln.logmean, sln.logstd);
|
||||||
|
|
||||||
|
free(seed);
|
||||||
|
}
|
53
examples/11_algebra/makefile
Normal file
53
examples/11_algebra/makefile
Normal file
|
@ -0,0 +1,53 @@
|
||||||
|
# Interface:
|
||||||
|
# make
|
||||||
|
# make build
|
||||||
|
# make format
|
||||||
|
# make run
|
||||||
|
|
||||||
|
# Compiler
|
||||||
|
CC=gcc
|
||||||
|
# CC=tcc # <= faster compilation
|
||||||
|
|
||||||
|
# Main file
|
||||||
|
SRC=example.c ../../squiggle.c
|
||||||
|
OUTPUT=example
|
||||||
|
|
||||||
|
## Dependencies
|
||||||
|
MATH=-lm
|
||||||
|
|
||||||
|
## Flags
|
||||||
|
DEBUG= #'-g'
|
||||||
|
STANDARD=-std=c99
|
||||||
|
WARNINGS=-Wall
|
||||||
|
OPTIMIZED=-O3 #-Ofast
|
||||||
|
# OPENMP=-fopenmp
|
||||||
|
|
||||||
|
## Formatter
|
||||||
|
STYLE_BLUEPRINT=webkit
|
||||||
|
FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT)
|
||||||
|
|
||||||
|
## make build
|
||||||
|
build: $(SRC)
|
||||||
|
$(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(MATH) -o $(OUTPUT)
|
||||||
|
|
||||||
|
format: $(SRC)
|
||||||
|
$(FORMATTER) $(SRC)
|
||||||
|
|
||||||
|
run: $(SRC) $(OUTPUT)
|
||||||
|
OMP_NUM_THREADS=1 ./$(OUTPUT) && echo
|
||||||
|
|
||||||
|
time-linux:
|
||||||
|
@echo "Requires /bin/time, found on GNU/Linux systems" && echo
|
||||||
|
|
||||||
|
@echo "Running 100x and taking avg time $(OUTPUT)"
|
||||||
|
@t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 1 thread: |" | sed 's|$$|ms|' && echo
|
||||||
|
|
||||||
|
## Profiling
|
||||||
|
|
||||||
|
profile-linux:
|
||||||
|
echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar"
|
||||||
|
echo "Must be run as sudo"
|
||||||
|
$(CC) $(SRC) $(MATH) -o $(OUTPUT)
|
||||||
|
sudo perf record ./$(OUTPUT)
|
||||||
|
sudo perf report
|
||||||
|
rm perf.data
|
BIN
examples/12_algebra_and_conversion/example
Executable file
BIN
examples/12_algebra_and_conversion/example
Executable file
Binary file not shown.
33
examples/12_algebra_and_conversion/example.c
Normal file
33
examples/12_algebra_and_conversion/example.c
Normal file
|
@ -0,0 +1,33 @@
|
||||||
|
#include "../../squiggle.h"
|
||||||
|
#include <math.h>
|
||||||
|
#include <stdint.h>
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
|
||||||
|
int main()
|
||||||
|
{
|
||||||
|
// set randomness seed
|
||||||
|
uint64_t* seed = malloc(sizeof(uint64_t));
|
||||||
|
*seed = 1000; // xorshift can't start with 0
|
||||||
|
|
||||||
|
// Convert to 90% confidence interval form and back
|
||||||
|
lognormal_params ln1 = { .logmean = 1.0, .logstd = 3.0 };
|
||||||
|
ci ln1_ci = convert_lognormal_params_to_ci(ln1);
|
||||||
|
printf("The 90%% confidence interval of Lognormal(%f, %f) is [%f, %f]\n",
|
||||||
|
ln1.logmean, ln1.logstd,
|
||||||
|
ln1_ci.low, ln1_ci.high);
|
||||||
|
lognormal_params ln1_ci_paramas = convert_ci_to_lognormal_params(ln1_ci);
|
||||||
|
printf("The lognormal which has 90%% confidence interval [%f, %f] is Lognormal(%f, %f)\n",
|
||||||
|
ln1_ci.low, ln1_ci.high,
|
||||||
|
ln1.logmean, ln1.logstd);
|
||||||
|
|
||||||
|
lognormal_params ln2 = convert_ci_to_lognormal_params((ci){.low = 1, .high = 10});
|
||||||
|
lognormal_params ln3 = convert_ci_to_lognormal_params((ci){.low = 5, .high = 50});
|
||||||
|
|
||||||
|
lognormal_params sln = algebra_product_lognormals(ln2, ln3);
|
||||||
|
ci sln_ci = convert_lognormal_params_to_ci(sln);
|
||||||
|
|
||||||
|
printf("Result of some lognormal products: to(%f, %f)\n", sln_ci.low, sln_ci.high);
|
||||||
|
|
||||||
|
free(seed);
|
||||||
|
}
|
53
examples/12_algebra_and_conversion/makefile
Normal file
53
examples/12_algebra_and_conversion/makefile
Normal file
|
@ -0,0 +1,53 @@
|
||||||
|
# Interface:
|
||||||
|
# make
|
||||||
|
# make build
|
||||||
|
# make format
|
||||||
|
# make run
|
||||||
|
|
||||||
|
# Compiler
|
||||||
|
CC=gcc
|
||||||
|
# CC=tcc # <= faster compilation
|
||||||
|
|
||||||
|
# Main file
|
||||||
|
SRC=example.c ../../squiggle.c
|
||||||
|
OUTPUT=example
|
||||||
|
|
||||||
|
## Dependencies
|
||||||
|
MATH=-lm
|
||||||
|
|
||||||
|
## Flags
|
||||||
|
DEBUG= #'-g'
|
||||||
|
STANDARD=-std=c99
|
||||||
|
WARNINGS=-Wall
|
||||||
|
OPTIMIZED=-O3 #-Ofast
|
||||||
|
# OPENMP=-fopenmp
|
||||||
|
|
||||||
|
## Formatter
|
||||||
|
STYLE_BLUEPRINT=webkit
|
||||||
|
FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT)
|
||||||
|
|
||||||
|
## make build
|
||||||
|
build: $(SRC)
|
||||||
|
$(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(MATH) -o $(OUTPUT)
|
||||||
|
|
||||||
|
format: $(SRC)
|
||||||
|
$(FORMATTER) $(SRC)
|
||||||
|
|
||||||
|
run: $(SRC) $(OUTPUT)
|
||||||
|
OMP_NUM_THREADS=1 ./$(OUTPUT) && echo
|
||||||
|
|
||||||
|
time-linux:
|
||||||
|
@echo "Requires /bin/time, found on GNU/Linux systems" && echo
|
||||||
|
|
||||||
|
@echo "Running 100x and taking avg time $(OUTPUT)"
|
||||||
|
@t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 1 thread: |" | sed 's|$$|ms|' && echo
|
||||||
|
|
||||||
|
## Profiling
|
||||||
|
|
||||||
|
profile-linux:
|
||||||
|
echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar"
|
||||||
|
echo "Must be run as sudo"
|
||||||
|
$(CC) $(SRC) $(MATH) -o $(OUTPUT)
|
||||||
|
sudo perf record ./$(OUTPUT)
|
||||||
|
sudo perf report
|
||||||
|
rm perf.data
|
BIN
examples/13_ergonomic_algebra/example
Executable file
BIN
examples/13_ergonomic_algebra/example
Executable file
Binary file not shown.
26
examples/13_ergonomic_algebra/example.c
Normal file
26
examples/13_ergonomic_algebra/example.c
Normal file
|
@ -0,0 +1,26 @@
|
||||||
|
#include "../../squiggle.h"
|
||||||
|
#include <math.h>
|
||||||
|
#include <stdint.h>
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
|
||||||
|
#define ln lognormal_params
|
||||||
|
#define to(...) convert_ci_to_lognormal_params((ci) __VA_ARGS__)
|
||||||
|
#define from(...) convert_lognormal_params_to_ci((ln) __VA_ARGS__)
|
||||||
|
#define times(a,b) algebra_product_lognormals(a,b)
|
||||||
|
|
||||||
|
int main()
|
||||||
|
{
|
||||||
|
// set randomness seed
|
||||||
|
uint64_t* seed = malloc(sizeof(uint64_t));
|
||||||
|
*seed = 1000; // xorshift can't start with 0
|
||||||
|
|
||||||
|
ln a = to({.low = 1, .high = 10});
|
||||||
|
ln b = to({.low = 5, .high = 500});
|
||||||
|
ln c = times(a, b);
|
||||||
|
|
||||||
|
printf("Result: to(%f, %f)\n", from(c).low, from(c).high);
|
||||||
|
printf("One sample from it is: %f\n", sample_lognormal(c.logmean, c.logstd, seed));
|
||||||
|
|
||||||
|
free(seed);
|
||||||
|
}
|
53
examples/13_ergonomic_algebra/makefile
Normal file
53
examples/13_ergonomic_algebra/makefile
Normal file
|
@ -0,0 +1,53 @@
|
||||||
|
# Interface:
|
||||||
|
# make
|
||||||
|
# make build
|
||||||
|
# make format
|
||||||
|
# make run
|
||||||
|
|
||||||
|
# Compiler
|
||||||
|
CC=gcc
|
||||||
|
# CC=tcc # <= faster compilation
|
||||||
|
|
||||||
|
# Main file
|
||||||
|
SRC=example.c ../../squiggle.c
|
||||||
|
OUTPUT=example
|
||||||
|
|
||||||
|
## Dependencies
|
||||||
|
MATH=-lm
|
||||||
|
|
||||||
|
## Flags
|
||||||
|
DEBUG= #'-g'
|
||||||
|
STANDARD=-std=c99
|
||||||
|
WARNINGS=-Wall
|
||||||
|
OPTIMIZED=-O3 #-Ofast
|
||||||
|
# OPENMP=-fopenmp
|
||||||
|
|
||||||
|
## Formatter
|
||||||
|
STYLE_BLUEPRINT=webkit
|
||||||
|
FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT)
|
||||||
|
|
||||||
|
## make build
|
||||||
|
build: $(SRC)
|
||||||
|
$(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(MATH) -o $(OUTPUT)
|
||||||
|
|
||||||
|
format: $(SRC)
|
||||||
|
$(FORMATTER) $(SRC)
|
||||||
|
|
||||||
|
run: $(SRC) $(OUTPUT)
|
||||||
|
OMP_NUM_THREADS=1 ./$(OUTPUT) && echo
|
||||||
|
|
||||||
|
time-linux:
|
||||||
|
@echo "Requires /bin/time, found on GNU/Linux systems" && echo
|
||||||
|
|
||||||
|
@echo "Running 100x and taking avg time $(OUTPUT)"
|
||||||
|
@t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 1 thread: |" | sed 's|$$|ms|' && echo
|
||||||
|
|
||||||
|
## Profiling
|
||||||
|
|
||||||
|
profile-linux:
|
||||||
|
echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar"
|
||||||
|
echo "Must be run as sudo"
|
||||||
|
$(CC) $(SRC) $(MATH) -o $(OUTPUT)
|
||||||
|
sudo perf record ./$(OUTPUT)
|
||||||
|
sudo perf report
|
||||||
|
rm perf.data
|
BIN
examples/14_twitter_thread_example/example
Executable file
BIN
examples/14_twitter_thread_example/example
Executable file
Binary file not shown.
43
examples/14_twitter_thread_example/example.c
Normal file
43
examples/14_twitter_thread_example/example.c
Normal file
|
@ -0,0 +1,43 @@
|
||||||
|
#include "../../squiggle.h"
|
||||||
|
#include <stdint.h>
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
|
||||||
|
double sample_0(uint64_t* seed){
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
double sample_1(uint64_t* seed){
|
||||||
|
return 1;
|
||||||
|
}
|
||||||
|
|
||||||
|
double sample_normal_mean_1_std_2(uint64_t* seed){
|
||||||
|
return sample_normal(1, 2, seed);
|
||||||
|
}
|
||||||
|
|
||||||
|
double sample_1_to_3(uint64_t* seed){
|
||||||
|
return sample_to(1, 3, seed);
|
||||||
|
}
|
||||||
|
|
||||||
|
int main()
|
||||||
|
{
|
||||||
|
// set randomness seed
|
||||||
|
uint64_t* seed = malloc(sizeof(uint64_t));
|
||||||
|
*seed = 1000; // xorshift can't start with 0
|
||||||
|
|
||||||
|
int n_dists = 4;
|
||||||
|
double weights[] = { 1, 2, 3, 4 };
|
||||||
|
double (*samplers[])(uint64_t*) = {
|
||||||
|
sample_0,
|
||||||
|
sample_1,
|
||||||
|
sample_normal_mean_1_std_2,
|
||||||
|
sample_1_to_3
|
||||||
|
};
|
||||||
|
|
||||||
|
int n_samples = 10;
|
||||||
|
for (int i = 0; i < n_samples; i++) {
|
||||||
|
printf("Sample #%d: %f\n", i, sample_mixture(samplers, weights, n_dists, seed));
|
||||||
|
}
|
||||||
|
|
||||||
|
free(seed);
|
||||||
|
}
|
53
examples/14_twitter_thread_example/makefile
Normal file
53
examples/14_twitter_thread_example/makefile
Normal file
|
@ -0,0 +1,53 @@
|
||||||
|
# Interface:
|
||||||
|
# make
|
||||||
|
# make build
|
||||||
|
# make format
|
||||||
|
# make run
|
||||||
|
|
||||||
|
# Compiler
|
||||||
|
CC=gcc
|
||||||
|
# CC=tcc # <= faster compilation
|
||||||
|
|
||||||
|
# Main file
|
||||||
|
SRC=example.c ../../squiggle.c
|
||||||
|
OUTPUT=example
|
||||||
|
|
||||||
|
## Dependencies
|
||||||
|
MATH=-lm
|
||||||
|
|
||||||
|
## Flags
|
||||||
|
DEBUG= #'-g'
|
||||||
|
STANDARD=-std=c99
|
||||||
|
WARNINGS=-Wall
|
||||||
|
OPTIMIZED=-O3 #-Ofast
|
||||||
|
# OPENMP=-fopenmp
|
||||||
|
|
||||||
|
## Formatter
|
||||||
|
STYLE_BLUEPRINT=webkit
|
||||||
|
FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT)
|
||||||
|
|
||||||
|
## make build
|
||||||
|
build: $(SRC)
|
||||||
|
$(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(MATH) -o $(OUTPUT)
|
||||||
|
|
||||||
|
format: $(SRC)
|
||||||
|
$(FORMATTER) $(SRC)
|
||||||
|
|
||||||
|
run: $(SRC) $(OUTPUT)
|
||||||
|
./$(OUTPUT) && echo
|
||||||
|
|
||||||
|
time-linux:
|
||||||
|
@echo "Requires /bin/time, found on GNU/Linux systems" && echo
|
||||||
|
|
||||||
|
@echo "Running 100x and taking avg time $(OUTPUT)"
|
||||||
|
@t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do ./$(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 1 thread: |" | sed 's|$$|ms|' && echo
|
||||||
|
|
||||||
|
## Profiling
|
||||||
|
|
||||||
|
profile-linux:
|
||||||
|
echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar"
|
||||||
|
echo "Must be run as sudo"
|
||||||
|
$(CC) $(SRC) $(MATH) -o $(OUTPUT)
|
||||||
|
sudo perf record ./$(OUTPUT)
|
||||||
|
sudo perf report
|
||||||
|
rm perf.data
|
27
scratchpad/core.c
Normal file
27
scratchpad/core.c
Normal file
|
@ -0,0 +1,27 @@
|
||||||
|
|
||||||
|
uint64_t xorshift64(uint64_t* seed)
|
||||||
|
{
|
||||||
|
// Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs"
|
||||||
|
// <https://en.wikipedia.org/wiki/Xorshift>
|
||||||
|
uint64_t x = *seed;
|
||||||
|
x ^= x << 13;
|
||||||
|
x ^= x >> 7;
|
||||||
|
x ^= x << 17;
|
||||||
|
return *seed = x;
|
||||||
|
}
|
||||||
|
|
||||||
|
double sample_unit_uniform(uint64_t* seed)
|
||||||
|
{
|
||||||
|
// samples uniform from [0,1] interval.
|
||||||
|
return ((double)xorshift64(seed)) / ((double)UINT64_MAX);
|
||||||
|
}
|
||||||
|
|
||||||
|
double sample_unit_normal(uint64_t* seed)
|
||||||
|
{
|
||||||
|
// // See: <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>
|
||||||
|
double u1 = sample_unit_uniform(seed);
|
||||||
|
double u2 = sample_unit_uniform(seed);
|
||||||
|
double z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2);
|
||||||
|
return z;
|
||||||
|
}
|
||||||
|
|
56
scratchpad/makefile
Normal file
56
scratchpad/makefile
Normal file
|
@ -0,0 +1,56 @@
|
||||||
|
# Interface:
|
||||||
|
# make
|
||||||
|
# make build
|
||||||
|
# make format
|
||||||
|
# make run
|
||||||
|
|
||||||
|
# Compiler
|
||||||
|
CC=gcc
|
||||||
|
# CC=tcc # <= faster compilation
|
||||||
|
|
||||||
|
# Main file
|
||||||
|
SRC=scratchpad.c ../squiggle.c
|
||||||
|
OUTPUT=scratchpad
|
||||||
|
|
||||||
|
## Dependencies
|
||||||
|
MATH=-lm
|
||||||
|
|
||||||
|
## Flags
|
||||||
|
DEBUG= #'-g'
|
||||||
|
STANDARD=-std=c99
|
||||||
|
WARNINGS=-Wall
|
||||||
|
OPTIMIZED=-O3 #-Ofast
|
||||||
|
# OPENMP=-fopenmp
|
||||||
|
|
||||||
|
## Formatter
|
||||||
|
STYLE_BLUEPRINT=webkit
|
||||||
|
FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT)
|
||||||
|
|
||||||
|
## make build
|
||||||
|
build: $(SRC)
|
||||||
|
$(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(MATH) -o $(OUTPUT)
|
||||||
|
|
||||||
|
format: $(SRC)
|
||||||
|
$(FORMATTER) $(SRC)
|
||||||
|
|
||||||
|
run: $(SRC) $(OUTPUT)
|
||||||
|
./$(OUTPUT)
|
||||||
|
|
||||||
|
verify: $(SRC) $(OUTPUT)
|
||||||
|
./$(OUTPUT) | grep "NOT passed" -A 2 --group-separator='' || true
|
||||||
|
|
||||||
|
time-linux:
|
||||||
|
@echo "Requires /bin/time, found on GNU/Linux systems" && echo
|
||||||
|
|
||||||
|
@echo "Running 100x and taking avg time $(OUTPUT)"
|
||||||
|
@t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 1 thread: |" | sed 's|$$|ms|' && echo
|
||||||
|
|
||||||
|
## Profiling
|
||||||
|
|
||||||
|
profile-linux:
|
||||||
|
echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar"
|
||||||
|
echo "Must be run as sudo"
|
||||||
|
$(CC) $(SRC) $(MATH) -o $(OUTPUT)
|
||||||
|
sudo perf record ./$(OUTPUT)
|
||||||
|
sudo perf report
|
||||||
|
rm perf.data
|
BIN
scratchpad/scratchpad
Executable file
BIN
scratchpad/scratchpad
Executable file
Binary file not shown.
20
scratchpad/scratchpad.c
Normal file
20
scratchpad/scratchpad.c
Normal file
|
@ -0,0 +1,20 @@
|
||||||
|
#include "../squiggle.h"
|
||||||
|
#include <math.h>
|
||||||
|
#include <stdint.h>
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
|
||||||
|
int main()
|
||||||
|
{
|
||||||
|
// set randomness seed
|
||||||
|
uint64_t* seed = malloc(sizeof(uint64_t));
|
||||||
|
*seed = 1000; // xorshift can't start with a seed of 0
|
||||||
|
|
||||||
|
for (int i = 0; i < 100; i++) {
|
||||||
|
double draw = sample_unit_uniform(seed);
|
||||||
|
printf("%f\n", draw);
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
free(seed);
|
||||||
|
}
|
139
squiggle.c
139
squiggle.c
|
@ -7,20 +7,28 @@
|
||||||
#include <sys/types.h>
|
#include <sys/types.h>
|
||||||
#include <time.h>
|
#include <time.h>
|
||||||
|
|
||||||
|
// Some error niceties; these won't be used until later
|
||||||
#define MAX_ERROR_LENGTH 500
|
#define MAX_ERROR_LENGTH 500
|
||||||
#define EXIT_ON_ERROR 0
|
#define EXIT_ON_ERROR 0
|
||||||
#define PROCESS_ERROR(error_msg) process_error(error_msg, EXIT_ON_ERROR, __FILE__, __LINE__)
|
#define PROCESS_ERROR(error_msg) process_error(error_msg, EXIT_ON_ERROR, __FILE__, __LINE__)
|
||||||
|
|
||||||
const double PI = 3.14159265358979323846; // M_PI in gcc gnu99
|
#define PI 3.14159265358979323846 // M_PI in gcc gnu99
|
||||||
|
#define NORMAL90CONFIDENCE 1.6448536269514722
|
||||||
|
|
||||||
|
// # Key functionality
|
||||||
|
// Define the minimum number of functions needed to do simple estimation
|
||||||
|
// Starts here, ends until the end of the mixture function
|
||||||
|
|
||||||
// Pseudo Random number generator
|
// Pseudo Random number generator
|
||||||
uint64_t xorshift32(uint32_t* seed)
|
uint64_t xorshift32(uint32_t* seed)
|
||||||
{
|
{
|
||||||
// Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs"
|
// Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs"
|
||||||
// See <https://stackoverflow.com/questions/53886131/how-does-xorshift64-works>
|
// See:
|
||||||
// https://en.wikipedia.org/wiki/Xorshift
|
// <https://en.wikipedia.org/wiki/Xorshift>
|
||||||
// Also some drama: <https://www.pcg-random.org/posts/on-vignas-pcg-critique.html>, <https://prng.di.unimi.it/>
|
// <https://stackoverflow.com/questions/53886131/how-does-xorshift32-works>,
|
||||||
// for floats
|
// Also some drama:
|
||||||
|
// <https://www.pcg-random.org/posts/on-vignas-pcg-critique.html>,
|
||||||
|
// <https://prng.di.unimi.it/>
|
||||||
uint64_t x = *seed;
|
uint64_t x = *seed;
|
||||||
x ^= x << 13;
|
x ^= x << 13;
|
||||||
x ^= x >> 17;
|
x ^= x >> 17;
|
||||||
|
@ -30,7 +38,7 @@ uint64_t xorshift32(uint32_t* seed)
|
||||||
|
|
||||||
uint64_t xorshift64(uint64_t* seed)
|
uint64_t xorshift64(uint64_t* seed)
|
||||||
{
|
{
|
||||||
// same as above, but for generating doubles
|
// same as above, but for generating doubles instead of floats
|
||||||
uint64_t x = *seed;
|
uint64_t x = *seed;
|
||||||
x ^= x << 13;
|
x ^= x << 13;
|
||||||
x ^= x >> 7;
|
x ^= x >> 7;
|
||||||
|
@ -48,7 +56,7 @@ double sample_unit_uniform(uint64_t* seed)
|
||||||
|
|
||||||
double sample_unit_normal(uint64_t* seed)
|
double sample_unit_normal(uint64_t* seed)
|
||||||
{
|
{
|
||||||
// See: <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>
|
// // See: <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>
|
||||||
double u1 = sample_unit_uniform(seed);
|
double u1 = sample_unit_uniform(seed);
|
||||||
double u2 = sample_unit_uniform(seed);
|
double u2 = sample_unit_uniform(seed);
|
||||||
double z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2);
|
double z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2);
|
||||||
|
@ -71,17 +79,40 @@ double sample_lognormal(double logmean, double logstd, uint64_t* seed)
|
||||||
return exp(sample_normal(logmean, logstd, seed));
|
return exp(sample_normal(logmean, logstd, seed));
|
||||||
}
|
}
|
||||||
|
|
||||||
|
inline double sample_normal_from_90_confidence_interval(double low, double high, uint64_t* seed)
|
||||||
|
{
|
||||||
|
// Explanation of key idea:
|
||||||
|
// 1. We know that the 90% confidence interval of the unit normal is
|
||||||
|
// [-1.6448536269514722, 1.6448536269514722]
|
||||||
|
// see e.g.: https://stackoverflow.com/questions/20626994/how-to-calculate-the-inverse-of-the-normal-cumulative-distribution-function-in-p
|
||||||
|
// 2. So if we take a unit normal and multiply it by
|
||||||
|
// L / 1.6448536269514722, its new 90% confidence interval will be
|
||||||
|
// [-L, L], i.e., length 2 * L
|
||||||
|
// 3. Instead, if we want to get a confidence interval of length L,
|
||||||
|
// we should multiply the unit normal by
|
||||||
|
// L / (2 * 1.6448536269514722)
|
||||||
|
// Meaning that its standard deviation should be multiplied by that amount
|
||||||
|
// see: https://en.wikipedia.org/wiki/Normal_distribution?lang=en#Operations_on_a_single_normal_variable
|
||||||
|
// 4. So we have learnt that Normal(0, L / (2 * 1.6448536269514722))
|
||||||
|
// has a 90% confidence interval of length L
|
||||||
|
// 5. If we want a 90% confidence interval from high to low,
|
||||||
|
// we can set mean = (high + low)/2; the midpoint, and L = high-low,
|
||||||
|
// Normal([high + low]/2, [high - low]/(2 * 1.6448536269514722))
|
||||||
|
double mean = (high + low) / 2.0;
|
||||||
|
double std = (high - low) / (2.0 * NORMAL90CONFIDENCE);
|
||||||
|
return sample_normal(mean, std, seed);
|
||||||
|
}
|
||||||
|
|
||||||
double sample_to(double low, double high, uint64_t* seed)
|
double sample_to(double low, double high, uint64_t* seed)
|
||||||
{
|
{
|
||||||
// Given a (positive) 90% confidence interval,
|
// Given a (positive) 90% confidence interval,
|
||||||
// returns a sample from a lognormal
|
// returns a sample from a lognorma with a matching 90% c.i.
|
||||||
// with a matching 90% c.i.
|
// Key idea: If we want a lognormal with 90% confidence interval [a, b]
|
||||||
const double NORMAL95CONFIDENCE = 1.6448536269514722;
|
// we need but get a normal with 90% confidence interval [log(a), log(b)].
|
||||||
|
// Then see code for sample_normal_from_90_confidence_interval
|
||||||
double loglow = logf(low);
|
double loglow = logf(low);
|
||||||
double loghigh = logf(high);
|
double loghigh = logf(high);
|
||||||
double logmean = (loglow + loghigh) / 2;
|
return exp(sample_normal_from_90_confidence_interval(loglow, loghigh, seed));
|
||||||
double logstd = (loghigh - loglow) / (2.0 * NORMAL95CONFIDENCE);
|
|
||||||
return sample_lognormal(logmean, logstd, seed);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
double sample_gamma(double alpha, uint64_t* seed)
|
double sample_gamma(double alpha, uint64_t* seed)
|
||||||
|
@ -129,14 +160,16 @@ double sample_gamma(double alpha, uint64_t* seed)
|
||||||
|
|
||||||
double sample_beta(double a, double b, uint64_t* seed)
|
double sample_beta(double a, double b, uint64_t* seed)
|
||||||
{
|
{
|
||||||
|
// See: https://en.wikipedia.org/wiki/Gamma_distribution#Related_distributions
|
||||||
double gamma_a = sample_gamma(a, seed);
|
double gamma_a = sample_gamma(a, seed);
|
||||||
double gamma_b = sample_gamma(b, seed);
|
double gamma_b = sample_gamma(b, seed);
|
||||||
return gamma_a / (gamma_a + gamma_b);
|
return gamma_a / (gamma_a + gamma_b);
|
||||||
}
|
}
|
||||||
|
|
||||||
double sample_laplace(double successes, double failures, uint64_t* seed){
|
double sample_laplace(double successes, double failures, uint64_t* seed)
|
||||||
// see <https://wikiless.esmailelbob.xyz/wiki/Beta_distribution?lang=en#Rule_of_succession>
|
{
|
||||||
return sample_beta(successes + 1, failures + 1, seed);
|
// see <https://en.wikipedia.org/wiki/Beta_distribution?lang=en#Rule_of_succession>
|
||||||
|
return sample_beta(successes + 1, failures + 1, seed);
|
||||||
}
|
}
|
||||||
|
|
||||||
// Array helpers
|
// Array helpers
|
||||||
|
@ -177,8 +210,7 @@ double array_std(double* array, int length)
|
||||||
// Mixture function
|
// Mixture function
|
||||||
double sample_mixture(double (*samplers[])(uint64_t*), double* weights, int n_dists, uint64_t* seed)
|
double sample_mixture(double (*samplers[])(uint64_t*), double* weights, int n_dists, uint64_t* seed)
|
||||||
{
|
{
|
||||||
// You can see a simpler version of this function in the git history
|
// Sample from samples with frequency proportional to their weights.
|
||||||
// or in C-02-better-algorithm-one-thread/
|
|
||||||
double sum_weights = array_sum(weights, n_dists);
|
double sum_weights = array_sum(weights, n_dists);
|
||||||
double* cumsummed_normalized_weights = (double*)malloc(n_dists * sizeof(double));
|
double* cumsummed_normalized_weights = (double*)malloc(n_dists * sizeof(double));
|
||||||
cumsummed_normalized_weights[0] = weights[0] / sum_weights;
|
cumsummed_normalized_weights[0] = weights[0] / sum_weights;
|
||||||
|
@ -203,7 +235,11 @@ double sample_mixture(double (*samplers[])(uint64_t*), double* weights, int n_di
|
||||||
return result;
|
return result;
|
||||||
}
|
}
|
||||||
|
|
||||||
// Sample from an arbitrary cdf
|
// # More cool stuff
|
||||||
|
// This is no longer necessary to do basic estimation,
|
||||||
|
// but is still cool
|
||||||
|
|
||||||
|
// ## Sample from an arbitrary cdf
|
||||||
struct box {
|
struct box {
|
||||||
int empty;
|
int empty;
|
||||||
double content;
|
double content;
|
||||||
|
@ -399,10 +435,10 @@ double sampler_danger(struct box cdf(double), uint64_t* seed)
|
||||||
|
|
||||||
// Get confidence intervals, given a sampler
|
// Get confidence intervals, given a sampler
|
||||||
|
|
||||||
struct c_i {
|
typedef struct ci_t {
|
||||||
float low;
|
float low;
|
||||||
float high;
|
float high;
|
||||||
};
|
} ci;
|
||||||
int compare_doubles(const void* p, const void* q)
|
int compare_doubles(const void* p, const void* q)
|
||||||
{
|
{
|
||||||
// https://wikiless.esmailelbob.xyz/wiki/Qsort?lang=en
|
// https://wikiless.esmailelbob.xyz/wiki/Qsort?lang=en
|
||||||
|
@ -418,7 +454,7 @@ int compare_doubles(const void* p, const void* q)
|
||||||
|
|
||||||
return 0;
|
return 0;
|
||||||
}
|
}
|
||||||
struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed)
|
ci get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed)
|
||||||
{
|
{
|
||||||
int n = 100 * 1000;
|
int n = 100 * 1000;
|
||||||
double* samples_array = malloc(n * sizeof(double));
|
double* samples_array = malloc(n * sizeof(double));
|
||||||
|
@ -427,7 +463,7 @@ struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* se
|
||||||
}
|
}
|
||||||
qsort(samples_array, n, sizeof(double), compare_doubles);
|
qsort(samples_array, n, sizeof(double), compare_doubles);
|
||||||
|
|
||||||
struct c_i result = {
|
ci result = {
|
||||||
.low = samples_array[5000],
|
.low = samples_array[5000],
|
||||||
.high = samples_array[94999],
|
.high = samples_array[94999],
|
||||||
};
|
};
|
||||||
|
@ -436,47 +472,54 @@ struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* se
|
||||||
return result;
|
return result;
|
||||||
}
|
}
|
||||||
|
|
||||||
// Do algebra over lognormals and normals
|
// # Small algebra manipulations
|
||||||
struct normal_parameters {
|
|
||||||
|
// here I discover named structs,
|
||||||
|
// which mean that I don't have to be typing
|
||||||
|
// struct blah all the time.
|
||||||
|
typedef struct normal_params_t {
|
||||||
double mean;
|
double mean;
|
||||||
double std;
|
double std;
|
||||||
};
|
} normal_params;
|
||||||
|
|
||||||
struct lognormal_parameters {
|
normal_params algebra_sum_normals(normal_params a, normal_params b)
|
||||||
double logmean;
|
|
||||||
double logstd;
|
|
||||||
};
|
|
||||||
|
|
||||||
struct normal_parameters algebra_sum_normals(struct normal_parameters a, struct normal_parameters b)
|
|
||||||
{
|
{
|
||||||
struct normal_parameters result = {
|
normal_params result = {
|
||||||
.mean = a.mean + b.mean,
|
.mean = a.mean + b.mean,
|
||||||
.std = sqrt((a.std * a.std) + (b.std * b.std)),
|
.std = sqrt((a.std * a.std) + (b.std * b.std)),
|
||||||
};
|
};
|
||||||
return result;
|
return result;
|
||||||
}
|
}
|
||||||
struct normal_parameters algebra_shift_normal(struct normal_parameters a, double shift)
|
|
||||||
{
|
|
||||||
struct normal_parameters result = {
|
|
||||||
.mean = a.mean + shift,
|
|
||||||
.std = a.std,
|
|
||||||
};
|
|
||||||
return result;
|
|
||||||
}
|
|
||||||
|
|
||||||
struct lognormal_parameters algebra_product_lognormals(struct lognormal_parameters a, struct lognormal_parameters b)
|
typedef struct lognormal_params_t {
|
||||||
|
double logmean;
|
||||||
|
double logstd;
|
||||||
|
} lognormal_params;
|
||||||
|
|
||||||
|
lognormal_params algebra_product_lognormals(lognormal_params a, lognormal_params b)
|
||||||
{
|
{
|
||||||
struct lognormal_parameters result = {
|
lognormal_params result = {
|
||||||
.logmean = a.logmean + b.logmean,
|
.logmean = a.logmean + b.logmean,
|
||||||
.logstd = sqrt((a.logstd * a.logstd) + (b.logstd * b.logstd)),
|
.logstd = sqrt((a.logstd * a.logstd) + (b.logstd * b.logstd)),
|
||||||
};
|
};
|
||||||
return result;
|
return result;
|
||||||
}
|
}
|
||||||
struct lognormal_parameters algebra_scale_lognormal(struct lognormal_parameters a, double k)
|
|
||||||
|
lognormal_params convert_ci_to_lognormal_params(ci x)
|
||||||
{
|
{
|
||||||
struct lognormal_parameters result = {
|
double loghigh = logf(x.high);
|
||||||
.logmean = a.logmean + k,
|
double loglow = logf(x.low);
|
||||||
.logstd = a.logstd,
|
double logmean = (loghigh + loglow) / 2.0;
|
||||||
};
|
double logstd = (loghigh - loglow) / (2.0 * NORMAL90CONFIDENCE);
|
||||||
|
lognormal_params result = { .logmean = logmean, .logstd = logstd };
|
||||||
|
return result;
|
||||||
|
}
|
||||||
|
|
||||||
|
ci convert_lognormal_params_to_ci(lognormal_params y)
|
||||||
|
{
|
||||||
|
double h = y.logstd * NORMAL90CONFIDENCE;
|
||||||
|
double loghigh = y.logmean + h;
|
||||||
|
double loglow = y.logmean - h;
|
||||||
|
ci result = { .low = exp(loglow), .high = exp(loghigh) };
|
||||||
return result;
|
return result;
|
||||||
}
|
}
|
||||||
|
|
26
squiggle.h
26
squiggle.h
|
@ -52,10 +52,30 @@ struct box sampler_cdf_double(double cdf(double), uint64_t* seed);
|
||||||
struct box sampler_cdf_box(struct box cdf(double), uint64_t* seed);
|
struct box sampler_cdf_box(struct box cdf(double), uint64_t* seed);
|
||||||
|
|
||||||
// Get 90% confidence interval
|
// Get 90% confidence interval
|
||||||
struct c_i {
|
typedef struct ci_t {
|
||||||
float low;
|
float low;
|
||||||
float high;
|
float high;
|
||||||
};
|
} ci;
|
||||||
struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed);
|
ci get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed);
|
||||||
|
|
||||||
|
// small algebra manipulations
|
||||||
|
|
||||||
|
typedef struct normal_params_t {
|
||||||
|
double mean;
|
||||||
|
double std;
|
||||||
|
} normal_params;
|
||||||
|
normal_params algebra_sum_normals(normal_params a, normal_params b);
|
||||||
|
|
||||||
|
typedef struct lognormal_params_t {
|
||||||
|
double logmean;
|
||||||
|
double logstd;
|
||||||
|
} lognormal_params;
|
||||||
|
lognormal_params algebra_product_lognormals(lognormal_params a, lognormal_params b);
|
||||||
|
|
||||||
|
|
||||||
|
lognormal_params convert_ci_to_lognormal_params(ci x);
|
||||||
|
|
||||||
|
|
||||||
|
ci convert_lognormal_params_to_ci(lognormal_params y);
|
||||||
|
|
||||||
#endif
|
#endif
|
||||||
|
|
Loading…
Reference in New Issue
Block a user