JJSierraM/squiggle.c
1
0
Fork 0
forked from personal/squiggle.c
Code Issues Pull Requests Packages Projects Releases Wiki Activity

add beta distribution samples!!

Browse Source
This commit is contained in:
NunoSempere 2023-07-16 14:38:12 +02:00
parent e94774ae3a
commit 607554f22b
2 changed files with 51 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
scratchpad/scratchpad
View File

Binary file not shown.

52
scratchpad/scratchpad.c
View File

@ -211,11 +211,11 @@ struct box incbeta(float a, float b, float x) {
if (x < 0.0 || x > 1.0){ if (x < 0.0 || x > 1.0){
if(EXIT_ON_ERROR){ if(EXIT_ON_ERROR){
printf("x out of bounds, in function incbeta, in %s (%d)", __FILE__, __LINE__); printf("x = %f, x out of bounds [0, 1], in function incbeta, in %s (%d)", __FILE__, __LINE__);
exit(1); exit(1);
}else{ }else{
char error_msg[200]; char error_msg[200];
snprintf(error_msg, 200, "x out of bounds, in function incbeta, in %s (%d)", __FILE__, __LINE__); snprintf(error_msg, 200, "x = %f, x out of bounds [0, 1], in function incbeta, in %s (%d)", x, __FILE__, __LINE__);
result.empty = 1; result.empty = 1;
result.error_msg = error_msg; result.error_msg = error_msg;
return result; return result;
@ -286,8 +286,16 @@ struct box incbeta(float a, float b, float x) {
} }
struct box cdf_beta(float x){ struct box cdf_beta(float x){
if(x < 0){
struct box result = { .empty = 0, .content = 0};
return result;
} else if(x > 1){
struct box result = { .empty = 0, .content = 1};
return result;
} else {
float successes = 1, failures = (2023-1945); float successes = 1, failures = (2023-1945);
return incbeta(successes, failures, x); return incbeta(successes, failures, x);
}
} }
float cdf_dangerous_beta(float x){ float cdf_dangerous_beta(float x){
@ -295,19 +303,31 @@ float cdf_dangerous_beta(float x){
// But it will propagate through the code // But it will propagate through the code
// So it doesn't feel like a great architectural choice; // So it doesn't feel like a great architectural choice;
// I prefer my choice of setting a variable which will determine whether to exit on failure or not. // I prefer my choice of setting a variable which will determine whether to exit on failure or not.
float successes = 1, failures = (2023-1945); // Ok, so the proper thing to do would be to refactor inverse_cdf
// but, I could also use a GOTO? <https://stackoverflow.com/questions/245742/examples-of-good-gotos-in-c-or-c>
// Ok, alternatives are:
// - Refactor inverse_cdf to take a box, take the small complexity + penalty. Add a helper
// - Duplicate the code, have a refactored inverse_cdf as well as a normal cdf
// - Do something hacky
// a. dangerous beta, which exits
// b. clever & hacky go-to statements
// i. They actually look fun to implement
// ii. But they would be hard for others to use.
if(x < 0){
return 0;
} else if(x > 1){
return 1;
} else {
float successes = 100, failures = 100;
struct box result = incbeta(successes, failures, x); struct box result = incbeta(successes, failures, x);
if(result.empty){ if(result.empty){
printf("%s", result.error_msg); printf("%s\n", result.error_msg);
exit(1); exit(1);
return 1;
}else{ }else{
return result.content; return result.content;
} }
} }
struct box dangerous_beta_sampler(uint32_t* seed)
// Think through what to do to feed the incbeta box into
{
return sampler(cdf_dangerous_beta, seed);
} }
int main() int main()
@ -374,5 +394,19 @@ int main()
float time_spent_2 = (float)(end_2 - begin_2) / CLOCKS_PER_SEC; float time_spent_2 = (float)(end_2 - begin_2) / CLOCKS_PER_SEC;
printf("Time spent: %f", time_spent_2); printf("Time spent: %f", time_spent_2);
// Get some beta samples
clock_t begin_3 = clock();
printf("\n\nGetting some samples from box sampler_dangerous_beta\n");
for (int i = 0; i < n; i++) {
struct box sample = sampler(cdf_dangerous_beta, seed);
if (sample.empty) {
printf("Error in sampler function");
} else {
printf("%f\n", sample.content);
}
}
clock_t end_3 = clock();
float time_spent_3 = (float)(end_3 - begin_3) / CLOCKS_PER_SEC;
printf("Time spent: %f", time_spent_3);
return 0; return 0;
} }