Skip to content
b-k edited this page Apr 3, 2013 · 4 revisions

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, 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__})

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

   \param n      The number of draws you want
   \param mean   The mean of the distribution
   \param sd     The standard deviation of the distribution.
   \param seed   The RNG seed


 */
apop_data *rnorm_base(rnorm_s);//but call via the macro.

// end rnorm.h


apop_data *rnorm_base(rnorm_s in){
    #define Send_error(E) {apop_data *d=apop_data_alloc(); d->error='E'; return d;}
    Apop_stopif(in.sd <0, Send_error(s), 0, "sd must be positive.");
    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);
    close_enough(apop_vector_mean(r01->vector), 0);
    close_enough(apop_vector_var(r01->vector), 1);

    apop_data *r3 = rnorm(.mean=3, .n=1e5);
    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);
    close_enough(apop_vector_mean(rwide->vector), 3);
    close_enough(apop_vector_var(rwide->vector), 9);
}
#endif

Clone this wiki locally