Re: a few questions about discrete rng

Re: a few questions about discrete rng

[Rodney Sparapani]

[-- Attachment #1: Type: TEXT/plain, Size: 3588 bytes --]

>
>Rodney Sparapani writes:
> > Brian:
> >  Apparently, I didn't understand the standards documentation with
> > respect to size_t.  But, it's too bad about gsl_vector.  I think
> > this decision should be reconsidered.  Also, I remember reading
> > that there have been some developments with respect to vectors in
> > C99.  Might this not also be a good reason to re-think the design?
>
>Nope ;-)
>

I was just reading the design part of the manual and it discusses size_t.  
However, I found it to be pretty vague.  Maybe you could elaborate a bit by 
adding your statement about it being the only portable way to index arrays 
across 32-bit, 64-bit, etc.  

> > It is alot easier to pass one parameter than 2 since you can easily
> > get them backwards.
>
>gsl's vectors use a stride, which adds unnecessary overhead to array
>based functions.

I don't understand this.  Would this tend to add overhead to linear algebra 
operations as well?  

>
> >  But, back to discrete with size_t and without gsl_vector.  I
> > thought that it would be better to streamline it a bit and make it
> > look more like the other rng's.  I'm attaching an example of that
> > which also contains a prototype for a multinomial randist.
> > Comments are welcome and please feel free to use this in GSL.
>
>Is there any algorithm better than O(n)?

Good question!  Since the discrete rng was already implemented in GSL, I didn't 
think to look for a faster method.  I checked Johnson which is the classic 
reference for multivariate distributions.  However, if there's material there, I 
missed it.  A more recent work, Random Number Generation and Monte Carlo Methods 
by JE Gentle (1998), seems to imply that there is:

"To generate a multinomial, a simple way is to work with the marginals; they are 
binomials.  The generation is done sequentially.  Each succeeding conditional 
marginal is binomial.  For efficiency, the first marginal considered would be 
the one with the largest probability.

Another interesting algorithm, due to Brown and Bromberg, is based on the fact 
that the conditional distribution of independent Poisson random variables, given 
their sum, is multinomial.  The use of this relationship requires construction 
of extensive tables.  Davis found this method to be slightly faster, once the 
setup operations are performed.  If multinomials are to be generated from 
distributions with different parameters, however, the simple method is more 
efficient."

MB Brown and J Bromberg, An efficient two-stage procedure for generating random 
variates from the multinomial distribution.  Am. Statistican(1984), 38, 216-9

CS Davis, The computer generation of multinomial random variates.  Computational 
Stat. and Data Analysis(1993), 16, 205-17.

I think Gentle's simple method would be equivalent to my simple method based on 
the discrete rng.  It appears that a method does exist that could beat it, but 
not by much and the overhead would probably make it slower if you needed 
simulations from multinomials with different probabilities.  Therefore, the 
simple methods appear to be a good compromise between efficiency and complexity. 
I'm attaching the latest which is now independent of discrete.c (due to 
borrowing and constructing a discrete rng).  Note that some rng features are 
still unimplemented, but I will add them.

Rodney Sparapani              Medical College of Wisconsin
Sr. Biostatistician           Patient Care & Outcomes Research (PCOR)
rsparapa@mcw.edu              http://www.mcw.edu/pcor
Was 'Name That Tune' rigged?  WWLD -- What Would Lombardi Do

[-- Attachment #2: multinomial.c --]
[-- Type: TEXT/x-sun-c-file, Size: 6445 bytes --]

#include <stdlib.h>
#include <gsl/gsl_vector_uint.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/*   Stack code borrowed from discrete.c                 */
/*** Begin Stack (this code is used just in this file) ***/

/* Stack code converted to use unsigned indices (i.e. s->i == 0 now
   means an empty stack, instead of -1), for consistency and to give a
   bigger allowable range. BJG */

typedef struct {
    size_t N;                      /* max number of elts on stack */
    size_t *v;                     /* array of values on the stack */
    size_t i;                      /* index of top of stack */
} gsl_stack_t;

static gsl_stack_t *
new_stack(size_t N) {
    gsl_stack_t *s;
    s = (gsl_stack_t *)malloc(sizeof(gsl_stack_t));
    s->N = N;
    s->i = 0;                  /* indicates stack is empty */
    s->v = (size_t *)malloc(sizeof(size_t)*N);
    return s;
}

static void
push_stack(gsl_stack_t *s, size_t v)
{
    if ((s->i) >= (s->N)) {
        fprintf(stderr,"Cannot push stack!\n");
        abort();                /* FIXME: fatal!! */
    }
    (s->v)[s->i] = v;
    s->i += 1;
}

static size_t pop_stack(gsl_stack_t *s)
{
    if ((s->i) == 0) {
        fprintf(stderr,"Cannot pop stack!\n");
        abort();               /* FIXME: fatal!! */
    }
    s->i -= 1;
    return ((s->v)[s->i]);
}

static inline size_t size_stack(const gsl_stack_t *s)
{
    return s->i;
}

static void free_stack(gsl_stack_t *s)
{
    free((char *)(s->v));
    free((char *)s);
}

/*** End Stack ***/


typedef struct {                /* struct for Walker algorithm */
    size_t K;
    size_t *A;
    double *F;
    gsl_rng *r;
} gsl_rng_discrete;

gsl_rng_discrete *
gsl_rng_discrete_alloc(const gsl_rng_type * T, size_t Kevents, const double *ProbArray)
{
    size_t k,b,s;
    gsl_rng_discrete *g;
    size_t nBigs, nSmalls;
    gsl_stack_t *Bigs;
    gsl_stack_t *Smalls;
    double *E;
    double pTotal = 0.0, mean, d;
    
    if (Kevents < 1) {
      /* Could probably treat Kevents=1 as a special case */

      GSL_ERROR_VAL ("number of events must be a positive integer", 
			GSL_EINVAL, 0);
    }

    /* Make sure elements of ProbArray[] are positive.
     * Won't enforce that sum is unity; instead will just normalize
     */

    for (k=0; k<Kevents; ++k) {
        if (ProbArray[k] < 0) {
	  GSL_ERROR_VAL ("probabilities must be non-negative",
			    GSL_EINVAL, 0) ;
        }
        pTotal += ProbArray[k];
    }

    /* Begin setting up the main "object" (just a struct, no steroids) */
    g = (gsl_rng_discrete *)malloc(sizeof(gsl_rng_discrete));
    g->r = gsl_rng_alloc (T);
    g->K = Kevents;
    g->F = (double *)malloc(sizeof(double)*Kevents);
    g->A = (size_t *)malloc(sizeof(size_t)*Kevents);

    E = (double *)malloc(sizeof(double)*Kevents);

    if (E==NULL) {
      GSL_ERROR_VAL ("Cannot allocate memory for randevent", ENOMEM, 0);
    }

    for (k=0; k<Kevents; ++k) {
        E[k] = ProbArray[k]/pTotal;
    }

    /* Now create the Bigs and the Smalls */
    mean = 1.0/Kevents;
    nSmalls=nBigs=0;
    for (k=0; k<Kevents; ++k) {
        if (E[k] < mean) ++nSmalls;
        else             ++nBigs;
    }
    Bigs   = new_stack(nBigs);
    Smalls = new_stack(nSmalls);
    for (k=0; k<Kevents; ++k) {
        if (E[k] < mean) {
            push_stack(Smalls,k);
        }
        else {
            push_stack(Bigs,k);
        }
    }
    /* Now work through the smalls */
    while (size_stack(Smalls) > 0) {
        s = pop_stack(Smalls);
        if (size_stack(Bigs) == 0) {
            /* Then we are on our last value */
            (g->A)[s]=s;
            (g->F)[s]=1.0;
            break;
        }
        b = pop_stack(Bigs);
        (g->A)[s]=b;
        (g->F)[s]=Kevents*E[s];
#if DEBUG
        fprintf(stderr,"s=%2d, A=%2d, F=%.4f\n",s,(g->A)[s],(g->F)[s]);
#endif        
        d = mean - E[s];
        E[s] += d;              /* now E[s] == mean */
        E[b] -= d;
        if (E[b] < mean) {
            push_stack(Smalls,b); /* no longer big, join ranks of the small */
        }
        else if (E[b] > mean) {
            push_stack(Bigs,b); /* still big, put it back where you found it */
        }
        else {
            /* E[b]==mean implies it is finished too */
            (g->A)[b]=b;
            (g->F)[b]=1.0;
        }
    }
    while (size_stack(Bigs) > 0) {
        b = pop_stack(Bigs);
        (g->A)[b]=b;
        (g->F)[b]=1.0;
    }
    /* Stacks have been emptied, and A and F have been filled */

    
#if 0
    /* if 1, then artificially set all F[k]'s to unity.  This will
     * give wrong answers, but you'll get them faster.  But, not
     * that much faster (I get maybe 20%); that's an upper bound
     * on what the optimal preprocessing would give.
     */
    for (k=0; k<Kevents; ++k) {
        (g->F)[k] = 1.0;
    }
#endif

#if KNUTH_CONVENTION
    /* For convenience, set F'[k]=(k+F[k])/K */
    /* This saves some arithmetic in gsl_rng_discrete(); I find that
     * it doesn't actually make much difference.
     */
    for (k=0; k<Kevents; ++k) {
        (g->F)[k] += k;
        (g->F)[k] /= Kevents;
    }
#endif    

    free_stack(Bigs);
    free_stack(Smalls);
    free((char *)E);

    return g;
}

size_t
gsl_rng_discrete_get(const gsl_rng_discrete *g)
{
    size_t c=0;
    double u,f;
    u = gsl_rng_uniform(g->r);
#if KNUTH_CONVENTION
    c = (u*(g->K));
#else
    u *= g->K;
    c = u;
    u -= c;
#endif
    f = (g->F)[c];
    /* fprintf(stderr,"c,f,u: %d %.4f %f\n",c,f,u); */
    if (f == 1.0) return c;

    if (u < f) {
        return c;
    }
    else {
        return (g->A)[c];
    }
}

void gsl_rng_discrete_free(gsl_rng_discrete *g)
{
    gsl_rng_free(g->r);
    free((char *)(g->r));
    free((char *)(g->A));
    free((char *)(g->F));
    free((char *)g);
}

gsl_vector_uint * gsl_ran_multinomial(const gsl_rng_discrete *d, const size_t n) 
{
  gsl_vector_uint * X = gsl_vector_uint_calloc(d->K);
  unsigned int i, j;

  for(i=0; i<n; i++) {
    j=gsl_rng_discrete_get(d);
    gsl_vector_uint_set(X, j, gsl_vector_uint_get(X, j)+1);
  }

  return X;
}

int main()
{
  double p[5]={0.05, 0.10, 0.15, 0.20, 0.5};
  
  gsl_rng_discrete * d = gsl_rng_discrete_alloc(gsl_rng_gfsr4, 5, p);

  gsl_vector_uint *x = gsl_ran_multinomial(d, 10000);
 
     FILE * f = fopen("test.txt", "w");
     gsl_vector_uint_fprintf (f, x, "%d");
     fclose (f);

     gsl_rng_discrete_free(d);

     return 0;
}

 

Re: a few questions about discrete rng

[Philip Kendall]

On Wed, Jul 10, 2002 at 02:38:55PM -0500, Rodney Sparapani wrote:
>
> I was just reading the design part of the manual and it discusses size_t.  
> However, I found it to be pretty vague.  Maybe you could elaborate a bit by 
> adding your statement about it being the only portable way to index arrays 
> across 32-bit, 64-bit, etc.  

On a 64-bit platform, arrays can have more than 2^32 entries, but an
unsigned int is often still only 32 bits long (this is certainly the
case for gcc/Alpha, Compaq cc/Alpha and Sun cc/Sparc).

Phil

-- 
  Philip Kendall <pak21@srcf.ucam.org>
  http://www.srcf.ucam.org/~pak21/

Re: a few questions about discrete rng

[Brian Gough]

Rodney Sparapani writes:
 > > > It is alot easier to pass one parameter than 2 since you can easily
 > > > get them backwards.
 > >
 > >gsl's vectors use a stride, which adds unnecessary overhead to array
 > >based functions.
 > 
 > I don't understand this.  Would this tend to add overhead to linear algebra 
 > operations as well?  

In linear algebra the stride is already used by BLAS.

Re: a few questions about discrete rng

[Brian Gough]

Rodney Sparapani writes:
 > Brian:
 >  Apparently, I didn't understand the standards documentation with
 > respect to size_t.  But, it's too bad about gsl_vector.  I think
 > this decision should be reconsidered.  Also, I remember reading
 > that there have been some developments with respect to vectors in
 > C99.  Might this not also be a good reason to re-think the design?

Nope ;-)

 > It is alot easier to pass one parameter than 2 since you can easily
 > get them backwards.

gsl's vectors use a stride, which adds unnecessary overhead to array
based functions.

 >  But, back to discrete with size_t and without gsl_vector.  I
 > thought that it would be better to streamline it a bit and make it
 > look more like the other rng's.  I'm attaching an example of that
 > which also contains a prototype for a multinomial randist.
 > Comments are welcome and please feel free to use this in GSL.

Is there any algorithm better than O(n)?

Re: a few questions about discrete rng

[Rodney Sparapani]

[-- Attachment #1: Type: TEXT/plain, Size: 995 bytes --]

Brian:

Apparently, I didn't understand the standards documentation with respect to 
size_t.  But, it's too bad about gsl_vector.  I think this decision should be 
reconsidered.  Also, I remember reading that there have been some developments 
with respect to vectors in C99.  Might this not also be a good reason to 
re-think the design?  It is alot easier to pass one parameter than 2 since you 
can easily get them backwards.

But, back to discrete with size_t and without gsl_vector.  I thought that it 
would be better to streamline it a bit and make it look more like the other 
rng's.  I'm attaching an example of that which also contains a prototype for a 
multinomial randist.  Comments are welcome and please feel free to use this in 
GSL.

Rodney Sparapani              Medical College of Wisconsin
Sr. Biostatistician           Patient Care & Outcomes Research (PCOR)
rsparapa@mcw.edu              http://www.mcw.edu/pcor
Was 'Name That Tune' rigged?  WWLD -- What Would Lombardi Do

[-- Attachment #2: multinomial.c --]
[-- Type: TEXT/x-sun-c-file, Size: 5034 bytes --]

#include <gsl/gsl_vector_uint.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "discrete.c"

typedef struct {                /* struct for Walker algorithm */
    size_t K;
    size_t *A;
    double *F;
    gsl_rng *r;
} gsl_rng_discrete;

gsl_rng_discrete *
gsl_rng_discrete_alloc(const gsl_rng_type * T, size_t Kevents, const double *ProbArray)
{
    size_t k,b,s;
    gsl_rng_discrete *g;
    size_t nBigs, nSmalls;
    gsl_stack_t *Bigs;
    gsl_stack_t *Smalls;
    double *E;
    double pTotal = 0.0, mean, d;
    
    if (Kevents < 1) {
      /* Could probably treat Kevents=1 as a special case */

      GSL_ERROR_VAL ("number of events must be a positive integer", 
			GSL_EINVAL, 0);
    }

    /* Make sure elements of ProbArray[] are positive.
     * Won't enforce that sum is unity; instead will just normalize
     */

    for (k=0; k<Kevents; ++k) {
        if (ProbArray[k] < 0) {
	  GSL_ERROR_VAL ("probabilities must be non-negative",
			    GSL_EINVAL, 0) ;
        }
        pTotal += ProbArray[k];
    }

    /* Begin setting up the main "object" (just a struct, no steroids) */
    g = (gsl_rng_discrete *)malloc(sizeof(gsl_rng_discrete));
    g->r = gsl_rng_alloc (T);
    g->K = Kevents;
    g->F = (double *)malloc(sizeof(double)*Kevents);
    g->A = (size_t *)malloc(sizeof(size_t)*Kevents);

    E = (double *)malloc(sizeof(double)*Kevents);

    if (E==NULL) {
      GSL_ERROR_VAL ("Cannot allocate memory for randevent", ENOMEM, 0);
    }

    for (k=0; k<Kevents; ++k) {
        E[k] = ProbArray[k]/pTotal;
    }

    /* Now create the Bigs and the Smalls */
    mean = 1.0/Kevents;
    nSmalls=nBigs=0;
    for (k=0; k<Kevents; ++k) {
        if (E[k] < mean) ++nSmalls;
        else             ++nBigs;
    }
    Bigs   = new_stack(nBigs);
    Smalls = new_stack(nSmalls);
    for (k=0; k<Kevents; ++k) {
        if (E[k] < mean) {
            push_stack(Smalls,k);
        }
        else {
            push_stack(Bigs,k);
        }
    }
    /* Now work through the smalls */
    while (size_stack(Smalls) > 0) {
        s = pop_stack(Smalls);
        if (size_stack(Bigs) == 0) {
            /* Then we are on our last value */
            (g->A)[s]=s;
            (g->F)[s]=1.0;
            break;
        }
        b = pop_stack(Bigs);
        (g->A)[s]=b;
        (g->F)[s]=Kevents*E[s];
#if DEBUG
        fprintf(stderr,"s=%2d, A=%2d, F=%.4f\n",s,(g->A)[s],(g->F)[s]);
#endif        
        d = mean - E[s];
        E[s] += d;              /* now E[s] == mean */
        E[b] -= d;
        if (E[b] < mean) {
            push_stack(Smalls,b); /* no longer big, join ranks of the small */
        }
        else if (E[b] > mean) {
            push_stack(Bigs,b); /* still big, put it back where you found it */
        }
        else {
            /* E[b]==mean implies it is finished too */
            (g->A)[b]=b;
            (g->F)[b]=1.0;
        }
    }
    while (size_stack(Bigs) > 0) {
        b = pop_stack(Bigs);
        (g->A)[b]=b;
        (g->F)[b]=1.0;
    }
    /* Stacks have been emptied, and A and F have been filled */

    
#if 0
    /* if 1, then artificially set all F[k]'s to unity.  This will
     * give wrong answers, but you'll get them faster.  But, not
     * that much faster (I get maybe 20%); that's an upper bound
     * on what the optimal preprocessing would give.
     */
    for (k=0; k<Kevents; ++k) {
        (g->F)[k] = 1.0;
    }
#endif

#if KNUTH_CONVENTION
    /* For convenience, set F'[k]=(k+F[k])/K */
    /* This saves some arithmetic in gsl_rng_discrete(); I find that
     * it doesn't actually make much difference.
     */
    for (k=0; k<Kevents; ++k) {
        (g->F)[k] += k;
        (g->F)[k] /= Kevents;
    }
#endif    

    free_stack(Bigs);
    free_stack(Smalls);
    free((char *)E);

    return g;
}

size_t
gsl_rng_discrete_get(const gsl_rng_discrete *g)
{
    size_t c=0;
    double u,f;
    u = gsl_rng_uniform(g->r);
#if KNUTH_CONVENTION
    c = (u*(g->K));
#else
    u *= g->K;
    c = u;
    u -= c;
#endif
    f = (g->F)[c];
    /* fprintf(stderr,"c,f,u: %d %.4f %f\n",c,f,u); */
    if (f == 1.0) return c;

    if (u < f) {
        return c;
    }
    else {
        return (g->A)[c];
    }
}

void gsl_rng_discrete_free(gsl_rng_discrete *g)
{
    gsl_rng_free(g->r);
    free((char *)(g->r));
    free((char *)(g->A));
    free((char *)(g->F));
    free((char *)g);
}

gsl_vector_uint * gsl_ran_multinomial(const gsl_rng_discrete *d, const size_t n) 
{
  gsl_vector_uint * X = gsl_vector_uint_calloc(d->K);
  unsigned int i, j;

  for(i=0; i<n; i++) {
    j=gsl_rng_discrete_get(d);
    gsl_vector_uint_set(X, j, gsl_vector_uint_get(X, j)+1);
  }

  return X;
}

int main()
{
  double p[5]={0.05, 0.10, 0.15, 0.20, 0.5};

  
  gsl_rng_discrete * d = gsl_rng_discrete_alloc(gsl_rng_gfsr4, 5, p);

  gsl_vector_uint *x = gsl_ran_multinomial(d, 10000);
 
     FILE * f = fopen("test.txt", "w");
     gsl_vector_uint_fprintf (f, x, "%d");
     fclose (f);

     gsl_rng_discrete_free(d);

     return 0;
}

 

Re: a few questions about discrete rng

[Brian Gough]

Rodney Sparapani writes:
 >  Thanks for your reply.  I don't know if I'm adding to this
 > discussion either :o) But, many other parts of GSL including
 > gsl_vector use size_t.  So, I can see where you got it from.
 > However, it's wide spread misuse is not good IMHO.  I suppose many
 > C programmers can "see through" typedefs and not be as easily
 > distracted as a bumbling statistician.  Something more
 > user-friendly would be nice though.  Also, after having learned a
 > complex API like gsl_vector, I think it best to propagate that
 > knowledge throughout.  Lastly, gsl_vector would have a certain
 > flexibility that pointers would not.  For example,
 > gsl_ran_discrete_t would not have to also store K.  Maybe not the
 > best example, but I think you can see what I mean.

Hi, 

size_t is the preferred way to represent the size of objects or arrays
in C.  Using "int" is not portable, e.g. to 64-bit platforms.  

We don't use gsl_vector for simple arrays.  One of the design
decisions was that it should be possible to use GSL in a minimal way,
without having to drag in too many header files if one only needs to
call a single function. See doc/gsl-design.texi for details.

Brian

Re: a few questions about discrete rng

[Rodney Sparapani]

James:

Thanks for your reply.  I don't know if I'm adding to this discussion either :o) 
But, many other parts of GSL including gsl_vector use size_t.  So, I can see 
where you got it from.  However, it's wide spread misuse is not good IMHO.  
I suppose many C programmers can "see through" typedefs and not be as easily
distracted as a bumbling statistician.  Something more user-friendly would be 
nice though.  Also, after having learned a complex API like gsl_vector, I think
it best to propagate that knowledge throughout.  Lastly, gsl_vector would have a
certain flexibility that pointers would not.  For example, gsl_ran_discrete_t
would not have to also store K.  Maybe not the best example, but I think you can
see what I mean.  

Rodney Sparapani              Medical College of Wisconsin
Sr. Biostatistician           Patient Care & Outcomes Research (PCOR)
rsparapa@mcw.edu              http://www.mcw.edu/pcor
Was 'Name That Tune' rigged?  WWLD -- What Would Lombardi Do

Re: a few questions about discrete rng

[James Theiler]

On Mon, 8 Jul 2002, Rodney Sparapani wrote:

] I'm working on a multinomial RNG for GSL.  I was looking at the code in
] randist/discrete.c and I'm perplexed.  For the number of categories, the type
] declaration is size_t.  This is the type returned by the sizeof operator
] according to the ANSI C standard for the number of bytes that a type has.  I
] don't think this is an appropriate definition.  I believe it should be an
] unsigned int, but not this one which has a special meaning and makes the source
] hard to read.  Also, the probability array is not of type gsl_vector.  Is there
] some reason for this?  If we all agree on these changes, then I will submit an
] updated discrete.c with my multinomial code.  Thanks.

I'll let Brian (or others) comment on size_t; I've only recently been
introduced to the convention that uses size_t for certain unsigned
integers, usually array sizes.  If I used it, I was copying somebody
else.

With respect to using gsl_vector versus an array, you ask: Is there
some reason for this?  The reason is that when I wrote it, gsl_vector
was not well established (at least not in my usage).  I don't have a
strong opinion.  You might argue that gsl_vector is more appropriately
limited to linear algebra operations.  In the early days of GSL, I
advocated an approach which required less dependencies on the various
types, but I've been won over to the single library approach, and in
that case the gsl_vector is a pretty useful type.  So, once again, I
realize I'm not actually contributing anything to the discussion, and
am deferring to Brian on whether this is an appropriate use of
gsl_vector.

regards,
jt

---------------------------------------------
James Theiler                     jt@lanl.gov
MS-D436, NIS-2, LANL        tel: 505/665-5682
Los Alamos, NM 87545        fax: 505/665-4414
----- Space and Remote Sensing Sciences -----

a few questions about discrete rng

[Rodney Sparapani]

I'm working on a multinomial RNG for GSL.  I was looking at the code in 
randist/discrete.c and I'm perplexed.  For the number of categories, the type 
declaration is size_t.  This is the type returned by the sizeof operator 
according to the ANSI C standard for the number of bytes that a type has.  I 
don't think this is an appropriate definition.  I believe it should be an 
unsigned int, but not this one which has a special meaning and makes the source 
hard to read.  Also, the probability array is not of type gsl_vector.  Is there 
some reason for this?  If we all agree on these changes, then I will submit an 
updated discrete.c with my multinomial code.  Thanks.

Rodney Sparapani              Medical College of Wisconsin
Sr. Biostatistician           Patient Care & Outcomes Research (PCOR)
rsparapa@mcw.edu              http://www.mcw.edu/pcor
Was 'Name That Tune' rigged?  WWLD -- What Would Lombardi Do