Rnorm

apop_draw makes a single draw from a distribution. This is a convenience function to make an arbitrary number of draws from a Normal distribution. The first part is what would go into a header file, including Doxygen-style documentation, then the function itself, then a test function, which demonstrates usage. If you define TESTING at compilation, then you can run this as a standalone program that tests the function.

#include <apop.h>

//////rnorm.h
typedef struct{
    size_t n;
    double mean, sd;
    size_t seed;
} rnorm_s;

#define rnorm(...) rnorm_base((rnorm_s){ __VA_ARGS__})

/** Convenience function to draw a list of random Normals, in honor of R's rnorm function.

   \param n      The number of draws you want. Default=1
   \param mean   The mean of the distribution. Default=0
   \param sd     The standard deviation of the distribution. Default=1
   \param seed   The RNG seed. Default: use apop_opts.rng_seed
   \return An apop_data set with the vector element filled with random draws.
   \exception If sd is negative, return a blank apop_data structure with error=='s'.
 */
apop_data *rnorm_base(rnorm_s);//but call via the macro.

// end rnorm.h


apop_data *rnorm_base(rnorm_s in){
    Apop_stopif(in.sd <0, {apop_data *d=apop_data_alloc(); d->error='s'; return d;}
                        , 0, "sd must be positive; you sent in %g.", in.sd);
    if (in.sd==0) in.sd=1;
    if (in.n==0) in.n=1;

    static gsl_rng *r;
    if (!r) r = apop_rng_alloc(in.seed ? in.seed : apop_opts.rng_seed++);
    apop_model *norm = apop_model_set_parameters(apop_normal, in.mean, in.sd);
    apop_data *out = apop_data_alloc(in.n);
    for (size_t i=0; i< in.n; i++)
        apop_draw(out->vector->data+i, r, norm);
    apop_model_free(norm);
    return out;
}


#ifdef TESTING
int main(){
    #define close_enough(a, b) assert(fabs(a-b) < 1e-2)
    apop_data *r = rnorm(5e5);
    assert(!r->error);
    close_enough(apop_vector_mean(r->vector), 0);
    close_enough(apop_vector_var(r->vector), 1);

    apop_data *r01 = rnorm(1e5, .seed=234);
    assert(!r01->error);
    close_enough(apop_vector_mean(r01->vector), 0);
    close_enough(apop_vector_var(r01->vector), 1);

    apop_data *r3 = rnorm(.mean=3, .n=1e5);
    assert(!r3->error);
    close_enough(apop_vector_mean(r3->vector), 3);
    close_enough(apop_vector_var(r3->vector), 1);

    apop_data *rwide = rnorm(.mean=3, .sd=3, .n=1e5);
    assert(!rwide->error);
    close_enough(apop_vector_mean(rwide->vector), 3);
    close_enough(apop_vector_var(rwide->vector), 9);

    apop_opts.log_file=fopen("/dev/null", "w"); //UNIX-specific to silence the warning.
    apop_data *bad = rnorm(.mean=3, .sd=-3, .n=1e5);
    assert(bad->error=='s');

}
#endif