diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c index 9ed24e9..120814c 100644 --- a/scratchpad/scratchpad.c +++ b/scratchpad/scratchpad.c @@ -3,6 +3,7 @@ #include #include // FLT_MAX, FLT_MIN #include // INT_MAX +#include #define VERBOSE 1 // to do: reuse more informative printing from build-your-own-lisp @@ -42,6 +43,12 @@ float cdf_squared_0_1(float x){ } } +float cdf_normal_0_1(float x){ + float mean = 0; + float std = 1; + return 0.5 * ( 1 + erf((x-mean)/(std * sqrt(2)) )); +} + // Inverse cdf struct box inverse_cdf(float cdf(float), float p){ // given a cdf: [-Inf, Inf] => [0,1] @@ -124,8 +131,37 @@ struct box inverse_cdf(float cdf(float), float p){ } +// Get random number between 0 and 1 +uint32_t xorshift32 +(uint32_t* seed) +{ + // Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" + // See + // https://en.wikipedia.org/wiki/Xorshift + // Also some drama: , + + uint32_t x = *seed; + x ^= x << 13; + x ^= x >> 17; + x ^= x << 5; + return *seed = x; +} + +// Distribution & sampling functions + +float rand_0_to_1(uint32_t* seed){ + return ((float) xorshift32(seed)) / ((float) UINT32_MAX); +} + // sampler based on inverse cdf -// to-do: integrals +struct box sampler(float cdf(float), uint32_t* seed){ + struct box result; + float p = rand_0_to_1(seed); + result = inverse_cdf(cdf, p); + return result; +} + +// to-do: integrals => beta distribution // main with an example int main(){ @@ -147,6 +183,10 @@ int main(){ }else{ printf("Inverse of the cdf at %f is: %f\n", 0.5, result2.content); } + + // set randomness seed + uint32_t* seed = malloc(sizeof(uint32_t)); + *seed = 1000; // xorshift can't start with 0 return 0; }