/* nag_rand_gen_multinomial (g05tgc) Example Program. * * Copyright 2008, Numerical Algorithms Group. * * Mark 9, 2009. */ /* Pre-processor includes */ #include #include #include #include #include #define X(I, J) x[(order == Nag_ColMajor)?(J*pdx + I):(I*pdx + J)] int main(void) { /* Integer scalar and array declarations */ Integer exit_status = 0; Integer lr, x_size, i, j, lstate, pdx; Integer *state = 0, *x = 0; /* NAG structures */ NagError fail; Nag_ModeRNG mode; /* Double scalar and array declarations */ double p_max; double *r = 0; /* Set the distribution parameters */ Integer k = 4; Integer m = 6000; double p[] = { 0.08e0, 0.1e0, 0.8e0, 0.02e0 }; /* Set the sample size */ Integer n = 20; /* Return the results in column major order */ Nag_OrderType order = Nag_ColMajor; /* Choose the base generator */ Nag_BaseRNG genid = Nag_Basic; Integer subid = 0; /* Set the seed */ Integer seed[] = { 1762543 }; Integer lseed = 1; /* Initialise the error structure */ INIT_FAIL(fail); printf( "nag_rand_gen_multinomial (g05tgc) Example Program Results\n\n"); /* Get the length of the state array */ lstate = -1; nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n", fail.message); exit_status = 1; goto END; } pdx = (order == Nag_ColMajor)?n:k; x_size = (order == Nag_ColMajor)?pdx * k:pdx * n; /* Calculate the size of the reference vector */ p_max = 0.0; for (i = 1; i < k; i++) p_max = (p_max < p[i])?p[i]:p_max; lr = 30 + 20 * sqrt(m * p_max * (1 - p_max)); /* Allocate arrays */ if (!(r = NAG_ALLOC(lr, double)) || !(state = NAG_ALLOC(lstate, Integer)) || !(x = NAG_ALLOC(x_size, Integer))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Initialise the generator to a repeatable sequence */ nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Generate the variates, initialising the reference vector at the same time */ mode = Nag_InitializeAndGenerate; nag_rand_gen_multinomial(order, mode, n, m, k, p, r, lr, state, x, pdx, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_rand_gen_multinomial (g05tgc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Display the variates*/ for (i = 0; i < n; i++) { for (j = 0; j < k; j++) printf("%12ld", X(i, j)); printf("\n"); } END: if (r) NAG_FREE(r); if (state) NAG_FREE(state); if (x) NAG_FREE(x); return exit_status; }