squiggle.c

commit 4deb2044d675927730e97093001022166f8b554e
parent c7d9f87ad7dada50a314ee6aa159935dded998b7
Author: NunoSempere <nuno.sempere@protonmail.com>
Date:   Sun, 16 Jul 2023 12:09:41 +0200

add more examples, tweak code a bit

Diffstat:
M
scratchpad/scratchpad
 | 
0

M
scratchpad/scratchpad.c
 | 
52
+++++++++++++++++++++++++++++++++++++---------------

2 files changed, 37 insertions(+), 15 deletions(-)
diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad
Binary files differ.
diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c
@@ -6,8 +6,6 @@
 #include <math.h>
 #define VERBOSE 1
 
-// to do: reuse more informative printing from build-your-own-lisp
-
 // Errors
 // https://mccue.dev/pages/7-27-22-c-errors
 // Huh, first time I've really needed a struct
@@ -17,6 +15,7 @@
 // options: 
 //   - exit
 //   - pass structs
+// to do: reuse more informative printing from build-your-own-lisp?
 struct box {
 	int empty;
 	float content; 
@@ -93,12 +92,10 @@ struct box inverse_cdf(float cdf(float), float p){
 		int convergence_condition = 0;
 		int count = 0;
 		while(!convergence_condition && (count < (INT_MAX/2) )){
-			if(VERBOSE){
-				printf("while loop\n");
-			}
 			float mid = (high + low)/2;
 			int mid_not_new = (mid == low) || (mid == high);
-			if(VERBOSE){
+			if(0){
+				printf("while loop\n");
 				printf("low: %f, high: %f\n", low, high); 
 				printf("mid: %f\n", mid);
 			}
@@ -161,32 +158,57 @@ struct box sampler(float cdf(float), uint32_t* seed){
 	return result;
 }
 
-// to-do: integrals => beta distribution
+// ~~[-] to-do: integrals => beta distribution~~
+// => instead, use this incomplete beta function implementation, based on continuous fractions:
+// <https://codeplea.com/incomplete-beta-function-c>
+// <https://github.com/codeplea/incbeta>
 
 // main with an example
 int main(){
 
 	// Uniform:
-	struct box result = inverse_cdf(cdf_uniform_0_1, 0.5);
-	if(result.empty){
-		printf("Inverse not calculated\n");
+	struct box result_1 = inverse_cdf(cdf_uniform_0_1, 0.5);
+	char* name_1 = "cdf_uniform_0_1";
+	if(result_1.empty){
+		printf("Inverse for %s not calculated\n", name_1);
 		exit(1);
 	}else{
-		printf("Inverse of the cdf at %f is: %f\n", 0.5, result.content);
+		printf("Inverse of %s at %f is: %f\n", name_1, 0.5, result_1.content);
 	}
 
 	// Squared cdf
-	struct box result2 = inverse_cdf(cdf_squared_0_1, 0.5);
-	if(result2.empty){
-		printf("Inverse not calculated\n");
+	struct box result_2 = inverse_cdf(cdf_squared_0_1, 0.5);
+	char* name_2 = "cdf_squared_0_1";
+	if(result_2.empty){
+		printf("Inverse for %s not calculated\n", name_2);
+		exit(1);
+	}else{
+		printf("Inverse of %s at %f is: %f\n", name_2, 0.5, result_2.content);
+	}
+
+	// Normal cdf
+	struct box result_3 = inverse_cdf(cdf_normal_0_1, 0.5);
+	char* name_3 = "cdf_normal_0_1";
+	if(result_3.empty){
+		printf("Inverse for %s not calculated\n", name_3);
 		exit(1);
 	}else{
-		printf("Inverse of the cdf at %f is: %f\n", 0.5, result2.content);
+		printf("Inverse of %s at %f is: %f\n", name_3, 0.5, result_3.content);
 	}
 
 	// set randomness seed
 	uint32_t* seed = malloc(sizeof(uint32_t));
 	*seed = 1000; // xorshift can't start with 0
 	
+	printf("\n\nGetting some samples from %s:\n", name_3);
+	int n=100;
+	for(int i=0;i<n; i++){
+		struct box sample = sampler(cdf_normal_0_1, seed);
+		if(sample.empty){
+			printf("Error in sampler function");
+		}else{
+			printf("%f\n", sample.content);
+		}
+	}
 	return 0;
 }