Conference rusure::math

>    assume it does.)  You have the results of the last 'x' operations,
>    where 'x' is at least 100 but less than 5000.  (Does this even
>    matter, or is this statement a ringer?)

For many statisticians, there is a lower bound of sample size, and below that
bound you don't do statistics, because there is "insufficient data".  There
is also an upper bound where simplifying assumptions make most techniques
unnecessary.  The numbers above say that this is a problem designed for 
such statisticians: there is "enough data" but not "too much data."

>    Now I understand that statistically speaking, each of the 12
>    possible results should have occurred one-twelfth of the time;
>    or if X is the number of historical operations we have on file,
>    then each value's statistically accurate appearance could be
>    computed by (X * .0833).  Right?  

As .1 says, random is not necessarily uniform.  The common assumption in 
this type of problem is that your experiment produces random but not 
necessarily uniform results.

>    What else can be said about this at this point in time, assuming
>    we have X operation results on file?  The text makes reference to
>    "inferential statistics" a few times, but nowhere does it explain

This is sounding to me like a traditional class, so I will give the traditional 
approach.

First you compute some descriptive statistics, like this.  Define XA as the 
number of times result A appears, XB similarly and so forth.  Then define a
histogram of the X's.  Then compute FA = XA/X, FB = XB/X and so forth.  Note
that the F's are all bounded by zero and one, and must sum to one.  These F's
summarize the X data, and by the assumption of randomness, summarize it
completely.

These F's are also an approximation to the underlying probabilities of the
experiment.  Call PA the probability of getting result A in a given test, then
FA is an approximation to PA.

>    exactly what that is?  Is there a way to "accurately predict"
>    (given some "percentage of assurance") what the next result should
>    be?  It is clearly (?) more than just determining which result
>    has occurred the least to date?  The text makes reference to

Some statisticians, like me, would allow you to use this approximation to say 
the probability that the next result will be A is PA, approximately FA.  Note 
that this says the result occuring most to date is most likely.  Note also 
that to do this right, you need to analyze the approximation, which is too
complex for this reply.  Note finally that many statisticans consider this a
completely illegitimate statement.

To apply inferential statistics, you have to decide what inference (or 
hypothesis) you want to test, and then pick the appropriate statistic to test 
it.

A typical hypothesis to introduce here is that all the PAs are equal.  Then one
asks whether the FAs are such as to rule out that hypothesis, to some degree
of certainty.  The most common statistic to test this is the chi-square,
computed by the sum of 

	( FA - 0.0833 ) ^2 / FA

over A, B, and so forth.  You also need the "degrees of freedom", which for 
this type of problem is one less than the number of categories, or 11, and your
desired level of confidence, for example, 0.01.  You look in the back of the 
book for a table of chi-square, and compare the entry there to the number you
computed above.  If the number you computed is larger, than you can rule out
the hypothesis that the PAs are equal.

>    taking the standard deviations, then grouping, allowing us to
>    determine 67% or 98% sureness, but I don't understand that, nor
>    why/how that works.

This is good inferential statistics, but for another problem.  It does not 
apply here unless your course is being very sloppy.

    OK, here's what I did to shed additional light on this.  (Whereas
    in fact, it actually made the whole issue muddier for me...)

    I wrote a quickie program to generate 5000 variables, where the
    value of the random variable must be a whole number between 1
    and 12 inclusive.  (VAX BASIC RANDOMIZE and RND functions -- does
    that mean *really* random?  Does that mean "uniform"?)

    I then went thru random samplings of the data, taking anywhere
    from the last 200 sampling from any one point in time to about
    1200 previous samplings from any one point in time.  As .2
    suggested, I applied the formula there (FA - 0.0833)^2 / FA)
    for FA, FB, FC, etc.  What did that tell me?  The values
    computed for each of the 12 variables was not equal, really
    not even looking close.  Was I to be looking for something
    specific when doing this?  Did this actually tell me something
    that I'm just not seeing?

    By picking a point in time "X" and being able to look back
    "Y" samplings, but already knowing the result of event "X+1",
    I'm trying to determine what I would be looking for in the
    data that would have led me to believe that the actual result
    achieved in sample "X+1" was what I would have predicted based
    on that data.

    Did I miss something, is this more statistical than I can hope
    to understand, or am I mis-applying statistical rules?

    Thanks for any & all help
    Jon

>    I wrote a quickie program to generate 5000 variables, where the
>    value of the random variable must be a whole number between 1
>    and 12 inclusive.  (VAX BASIC RANDOMIZE and RND functions -- does
>    that mean *really* random?  Does that mean "uniform"?)

Assuming that you did something like N = FIX(1 + 12*RND), then that is nearly
uniform.  Technically the series you generated is pseudo-random, but last I
heard it would take a roomful of VAX-9000s and a few people-years of analysis
to detect the flaws in the current Random function in the RTL.  So don't worry
about it, and assume you have a uniform, random distribution.

>    I then went thru random samplings of the data, taking anywhere
>    from the last 200 sampling from any one point in time to about
>    1200 previous samplings from any one point in time.  As .2
>    suggested, I applied the formula there (FA - 0.0833)^2 / FA)
>    for FA, FB, FC, etc.  What did that tell me?  The values

Nothing yet.

>    computed for each of the 12 variables was not equal, really
>    not even looking close.  Was I to be looking for something
>    specific when doing this?  Did this actually tell me something
>    that I'm just not seeing?

You need to add the twelve values you get, then compare them to the chi-square 
value.  Call the sum of the twelve values CHI-SQUARE.  I predict that nine times
out of ten, CHI-SQUARE is less than 5.  This will confirm that RND is pretty
good at producing a uniform distribution.

Oops, I saw an error in the formula above.  It should be

	SUM of (FA - 0.0833...)^2 / 0.0833...

>    By picking a point in time "X" and being able to look back
>    "Y" samplings, but already knowing the result of event "X+1",
>    I'm trying to determine what I would be looking for in the
>    data that would have led me to believe that the actual result
>    achieved in sample "X+1" was what I would have predicted based
>    on that data.

Once you know that this is a uniform, random distribution, there is no more 
information to be obtained from any past sequence.  In principle, you could
do much more careful analyses of the sequences to detect and exploit the 
deviations of RND from a true random number, but then we are back to the
roomful of VAX 9000s.

        re .5,
        
>>		 If you plot the positions of N-tuples in N-space all the
>>    points fall in nice neat hyperplanes.
        
        I've always wondered whether for these N-tuples they
        intend
        
        	(x1,x2,...,xN), (x2,x3,...,xN+1), (x3,x4,...,xN+2),...
        
        or
        
        	(x1,...,xN),(xN+1,...x2N),(x2N+1,...,x3N),...
        
        Dan

    Re .6

    They mean the former (cyclic shift thru the sequence) but in
    practice it really doesn't matter...

    - Jim

RE: .6 (Dan)

    They don't really mean either one, although the first comes closer to
    the truth, the second comes closer to practical application.

    Any n-tuple generated by a linear congruential generator will occur on
    one of a limited number of hyperplanes.  The location of the hyperplane
    does not depend on the "phase" of the n-tuple relative to any
    particular, arbitrary starting point in the PRNGs cycle (for that
    matter, for n r.p to the PRNGs cycle length, any n-tuple may be
    considered to be at all phases).

    If you want to "see" the hyperplanes by some plotting technique, you
    can choose overlapping or non-overlapping n-tuples, or choose a random
    x1 by an appropriate independent means and use that an the next n-1
    values for the point.

    In considering how this flaw might influence actual applications, you
    would normally be worried about generating n values for one simulation
    "trial" and then generating n more for the next, and therefore the
    sequencial non-overlapping case would be what you would actually be
    worrying about.

    The spectral test operates algebraically on all the n-tuples (and
    therefore on overlapping ones) to characterize how much space there is
    between the planes.

					Topher

.5>    A bit of a nit.  Unless the BASIC RTL uses something other than the
>    standard VMS math RTL or the "current" RTL is different from the one we
>    have on our system, you have overstated the case by quite a bit.
>    The RTL uses a LCG, and it is pretty easy to detect the flaws if you
>    know how.  If you plot the positions of N-tuples in N-space all the
>    points fall in nice neat hyperplanes.  For small N, and properly chosen
>    parameters, the planes are very close together and so constitute a
>    minor deviation in the low-order bits.  For even moderate N, however
>    (say 10) the effect is unavoidably quite severe.

Yes, I knew about this.  In fact, I once plotted on a GT40 the PDP-11 random
number generator output using N=3.  It fell into about 5 planes.  The VMS 
RTL people were aware of this problem, and chose a much different multiplier,
69069.  I never checked, but I think it avoided this problem, in the sense that
the number of planes for N=3 became very large.

The nit here is that you don't get all the information: the floating point 
values have been scaled and truncated into the 12 bins.  If N=10, you will have
12^10 = 60 billion bins to fill.  

I believe there are ways of determining the seed from the data, but the only 
one I can think of is brute force: generate the entire 2^32 sequence of 
integers, and then pattern match your sequence against it.  I think that a 
sequence of about ten integers is likely to be unique.

RE: .9 (Wally)

    Yes, the multiplier they used (69069) is quite good for small N.  But
    one only checks LCG parameters for small N -- up to 6 or 8.  The
    behavior gets worse pretty quickly with increasing N.  As I remember
    (don't ask why I should remember this), the upper-bound for the number
    of hyperplanes for an LCG is:

			1/N
		(N! * M)

    with N being the dimensionality (number of elements in the n-tuple) and
    M begin the modulus used (which for the VMS PRNG is 2^31).  When N is
    3, the upper bound comes out to 2344 (non-hyper)planes, and the
    particular choice of multiplier does a good job of near-optimality.
    For N=10, however, the best that any multiplier could accomplish would
    be about 38 planes, and no one bothered to see how well the multiplier
    actually approached that limit.

    Yes, if you are talking about calculations after you have reduced the
    number to an integer between 0 and 12, then if the multiplier happens
    to be nearly optimum for this dimensionality, then the effect will be
    slight (since each "bin" will contain multiple 9-planes), but I was
    talking about weaknesses in the underlying routine (as I thought you
    were also).  (And, of course, we have virtually no assurance that this
    multiplier comes anywhere close to optimality for N=10).

>I believe there are ways of determining the seed from the data, but the only 
>one I can think of is brute force: generate the entire 2^32 sequence of 
>integers, and then pattern match your sequence against it.  I think that a 
>sequence of about ten integers is likely to be unique.

    So far I have limited myself to reasonable statistical flaws, i.e.,
    flaws that a reasonable statistical test might reveal without any
    special knowledge of the underlying algorithm.  You bring up algebraic
    "weakness", however.  People interested in cryptography have developed
    quite efficient techniques for determining a seed for an LCG given a
    few adjacent numbers derived from the high-order bits.  Theoretically,
    10 values in the range 0-11 should be enough or almost enough to
    determine the seed if the multiplier, modulus and additive factor are
    known.  In practice a few more might be needed.

    Basically, each new 0-11 value gives you lg(12) bits of information
    (where lg is log-base-2), or about 3.6 bits of information.  Nine
    values would give us 32.4 bits of information which is slightly more
    than enough (since we need 31 bits).  As it happens, though, a
    particular sequence of values may overdetermine parts of the answer,
    and therefore underdetermine others.  An additional value or two may
    be necessary to disambiguate.

    So 10 values are probably unique, but may not be.  Nine values may be
    unique.

					    Topher

!    Yes, the multiplier they used (69069) is quite good for small N. 

The current rtl RNG has a rather large failing. If one looks at the seg of
numbers generated by bit position, they cycle cooresponding to their 
bit position. bit 0 cycles from 0 1 0 1 0 1 0 ...
bit 1 cycle every four numbers.
bit 24 cycles 2**25 numbers.
bit n cycles every 2**n+1 numbers.

>    with N being the dimensionality (number of elements in the n-tuple) and
>    M begin the modulus used (which for the VMS PRNG is 2^31).  When N is
                                                         ^^^^
    2^32...

    The number 69069 was chosen by Marsaglia after an exhaustive
    search using the spectral test for k-tuples for k = 1 to 5.
    It falls down at 6 dimensions - I once wrote a little
    program that sent the normal to that hyperplane in dimension
    6 to the x axis by a a Householder reflection and plotted
    the x and y components of 6-tuples rotated by that transformation.

    A bunch of lines showed right up!

    It's no problem to fix this - just use another generator to
    shuffle the sequence returned by your LCG and you will avoid the
    hyperplane problem.

    - Jim

    I hate to break up the fun in this rathole, but back to one of the
    original questions: is there a statistical method by which one could
    say, for example, that based on the history of events to date, event
    'A' is more likely to occur next than event 'B'?
    
    Or putting a new twist on the question, if you were willing to assume
    (for whatever reason) that events 'A', 'B', 'C', and 'D' were NOT going
    to occur as the result of the next experiment, does that lend any help
    in determining which of the remaining 'E' thru 'L' events should occur
    as the result of the next experiment?
    
    I promise that after I get an answer here, you may return undisturbed
    to this rathole...   :-)
    
    Thanks
    Jon

    Re .14:
    
    I think the text you were reporting about in .0 was trying to make a
    point like this:
    
    Suppose you have a set of trials and the results of each, and you know
    there are twelve possible results.  Now the thing is you do NOT know
    how likely each result is.  Yes, they are random, but you do not know
    that each result has an equal probability.  Your goal is to try to
    figure out something about the probabilities.
    
    Well, one way you can do this is to assume, for the moment, that each
    event has an equal probability.  Then you look at the set of results
    you have.  To "look at" them, you analyze them with various statistical
    tests to see if the results are, considered all together, something you
    are likely to get if each event has equal probability.  If you got
    equal numbers for each event or close to equal numbers for each event
    (and "close" is told to you by the statistical tests you apply), then
    you say "Okay, the results are close to what I would expect most of the
    time if each event had equal probability.".  But if the results aren't
    close to what would happen most of the time, then you say "Hey, the
    results aren't close.  Now, this could be a fluke, since these are
    random trials, but this is so unlikely I suspect that the probabilities
    for the different results aren't equal.".
    
    There, you've done statistics.
    
    You could also work with different assumptions.  Instead of assuming
    the probabilities are equal, you could assume they are some other
    specific values, and then you would run statistical tests to see if the
    results are within what would happen most of the time with those
    probabilities.
    
    One point with statistics is that it can never really tell you that A
    is more likely to occur than B.  It does tell you that A has been
    occurring a lot more than you would expect if it were equally likely. 
    Now that is not a deductive proof, but it is good enough for human
    beings.  At that point, we tend to say "All right, something's fishy
    here, I bet A is actually more likely than B.".  That's not because it
    is more likely just on the next trial, but because it has been more
    likely all along and now you have enough information to suspect that.
    
    
    				-- edp

RE: .13 (Jim)

>    2^32...

    Right you are, Jim, I had misremembered.  Increase the best possible
    number of hyperplanes for 3-tuples by 2^(1/3) and the maximum number of
    hyperplanes for 10-typles by 2^(1/10) in my previous statements.


>    It's no problem to fix this - just use another generator to
>    shuffle the sequence returned by your LCG and you will avoid the
>    hyperplane problem.

    If you are going to use another generator, you might as well use a good
    one in the first place and not bother with the mixing.  Alternately,
    you can use the "self-shuffle" recommended in the second edition of
    Knuth vol 2 -- which is faster, doesn't need an "additional" generator,
    and seems to be just as good.

					Topher

Note 95 discusses a problem centered around a pseudo random digit generator 
that appears to be non-uniform, for reasons which are not yet known.

Lynn 

.14>    I hate to break up the fun in this rathole, but back to one of the
>    original questions: is there a statistical method by which one could
>    say, for example, that based on the history of events to date, event
>    'A' is more likely to occur next than event 'B'?

Unfortunately, this 'rathole' is providing the best answer to your question.  .1
.2 and .15 offer some other answers.

Perhaps a classification will make the situation clearer.  Assume you have a 
sequence of integers chosen from the set {1 ... 12 }, which shows no obvious
pattern.  If you allow me to use an imprecise word like 'random', then there
are four possibilities:

1. the sequence is uniform and random

2. the sequence is random but not uniform

3. the sequence is uniform but not random

4. the sequence is neither uniform nor random

Now lets look at the answer to your question for these four cases:

1.  In this case, analysis is a waste of time.  The probability of any integer
from the set coming up next is 1/12.  No amount of analysis and no statistical
technique will give you a better answer.

2.  In this case you can get enough sample data to estimate the probability 
that each integer will come up, called the probability distribution.  Once 
you have a good estimate of this distribution, you are done.  No amount of 
analysis and no statistical technique will give you a better answer.  There are
some cute techniques for putting bounds around your guesses for the next value,
but I'm not sure you are interested.

Note, while we are here, that the descriptive statistics of .2 apply to case 2 
above, and the inferential statistics of .2, .4 and .15 are a way of 
distinguishing between case 1 and case 2.

3. In this case you need to look for the deviations from randomness and 
exploit them.  There are basically two ways to do this: clear box (the good way)
and black box (the bad way).  The clear box method looks at the mechanism 
generating the sequence and searches for relationships among the elements.
The 'rathole' discussion here is just such a search for a particular mechanism.
If this is the mechanism of interest, then the 'rathole' provides the answer.

The black box method assumes that we cannot look directly at the mechanism, so
we analyze the sequence itself to see what clues it provides about deviations
from randomness.  A lot of stock-market and economic forecasting is done this 
way.  This is a well-studied but complex field.  If you are interested, check
the subject index of a large library under Statistics, Time Series Analysis.

4. This case introduces nothing new, it is just a combination of .2 and .3, so
techniques from both those cases apply.  In practice, case 4 is the one treated,
since it is much more common than case 3.

.14>    Or putting a new twist on the question, if you were willing to assume
>    (for whatever reason) that events 'A', 'B', 'C', and 'D' were NOT going
>    to occur as the result of the next experiment, does that lend any help
>    in determining which of the remaining 'E' thru 'L' events should occur
>    as the result of the next experiment?

All depends.  If you have no rational reason for that assumption, then it does 
not help.  If you have a reason, then it must be based on experience (black box)
or mechanism (clear box).  There are various statistical techniques for 
converting these rational reasons into probability distributions on E thru L
events.

RE: .20 (Jon)

>     * I understand when you say that "N" is the sample size, but is that
>     the same as "n" (lower case vs upper case), or is "n" supposed to
>     represent something else?

    They're the same.  Its a typo.  You don't think I actually took the
    time to proofread all that, do you? :-)

>     For example (and I haven't done the math yet, just pulling numbers out
>     of the air), say I get values of 1.23 and 2.34 for "A" and 1.88 and
>     2.99 for "B" after doing the math shown above.  What did those values
>     tell me about the statistical likelihood of "A" occuring as the result
>     of the next experiment over "B".  (If anything)?

    THESE particular values wouldn't tell you anything about the
    statistical likelihood of anything.  They would simply tell you that
    either I gave you the wrong formula, or you did your math wrong. :-)

    All bounds should come out between 0 and 1.  The assumptions made here
    mean that we can meaningfully talk about there being some real, specific
    probability that, for example, a particular trial will have outcome A.
    We don't know what that probability is, but, we can, on the basis of
    past behavior, estimate what that probability is.  The simplest thing
    we can do is make a point estimate and just say that "the probability
    of any particular trial having outcome A is probably approximately fa".
    This is frequently done, but the problem is that that statement doesn't
    tell us how probably or how approximately.

    The alternative is to use a confidence interval.  If you use the
    formula I gave with P=.95 and you had 100 trials, of which 40 came out
    with outcome A, then your 95% confidence interval for Fa (the unknown
    actual probability for outcome A) is roughly (.3, .5).  Ignoring the
    denials of the tradionalists that such a statement actually means
    anything at all, this means that you can be 95% sure that Fa is
    actually between .3 and .5.

    The same applies to the other outcomes.  You have to be a little
    careful, however, in interpretting all of the numbers at once.  There
    is some degree of interaction between them.  For example, not all of
    the outcomes "actual" values can be near their respective lower bounds
    or the total probability of any of the outcomes occuring would be less
    than 1.  Also, with 12 95% confidence intervals, there is close to a
    50% chance that at least one of the actual probabilities will be
    outside its calculated confidence interval.

    Nevertheless, these 12 confidence intervals will give you a good
    feeling for what the relative likelihoods of the various outcomes
    are.  I suggest "plotting" the results.  The x axis should range
    between 0 and 1 (although you might change that upper limit if your
    largest upper limit is relatively small).  The y axis should be your
    twelve outcomes (make two plots, one in which they are arranged A, B,
    C, D etc; and one in which the one with the smallest point estimate is
    at the bottom, the next largest above it, and so on).  For each outcome
    draw a horizontal bar from the lower bound to the upper bound with a
    dot at the center point (the point estimate).

					Topher

    Topher,

    Thanks much for your clarifications.  Now I've setup a couple of Q&D
    (quick & dirty) programs to do the math based on my test file, going
    back 100, 200, 300, 400, & 500 experiments, seeing what values I get
    each time for each "A" thru "L" variable.  (I picked those numbers
    arbitrarily; figured less than 100 experiments would not be enough to
    base a good guess on, and more than 500 might tend to be too many?)

    I must have done the math right, because I almost always get lower and
    upper bound values between 0 and 1.  (Although I occasionally get a
    negative value?)  These values didn't appear to tell me much (maybe I'm
    still not sure what to look for?), so I took your suggestion and
    graphed the values, connecting the midpoints (midpoint = halfway
    between upper & lower bounds) of Fa thru Fl.  Nice graph, but is it
    supposed to tell me anything significant?  Just wondering if I'm
    looking for anything in particular; ie:

                           Z(.9)          Z(.95)           Z(.99)
                           -----          ------           ------
    Fa              0.0923  0.1797     0.0855  0.1865    0.0722  0.1998
    Fb              0.0539  0.1301     0.0482  0.1358    0.0369  0.1471
    Fc              0.0404  0.1116     0.0352  0.1168    0.0248  0.1272
       etc.

    If this is telling me something that I'm missing, please let me
    know what it is.

    Thanks
    Jon

RE: .22 (Jon)

    > ... and more than 500 might tend to be too many?)

    I wouldn't worry about "too many" too much.  What is too many depends
    on context.  The more samples you take the more sensitive your tests
    become.  If you take "too many" samples your tests become majorly
    influenced by effects that are so small that you want to ignore them. 
    Your computer simulation, however, is a pretty direct model of the
    process you are interested in.  The only extraneous factors (e.g.,
    slight deviations from ideal behavior of your pseudorandom number
    generator) are so small as to require a great many samples before they
    would interfere.  And there is less of an issue with descriptive
    statistics than with hypothesis testing in any case.

    You could unquestionably increase that to 50,000 samples (which would
    give you roughly one more decimal place of accuracy) and could probably
    go to 5,000,000 samples (roughly two more decimal places of accuracy).

    > (Although I occasionally get a negative value?)

    You probably get the negative values with the P=.99, N=100 case.  The
    formula I gave is an approximation.  The usual rule of thumb given for
    this (and related) approximations is that the Nx (the number of samples
    with outcome "x") should be at least 5 for it to be "reasonably
    accurate".  The approximation, however, is weaker in the "tails", and
    .99 is far enough out in the tail to really require more data.  Just
    don't take the P=.99/N=100 case too seriously.

    > These values don't tell me much.

    That is as it should be, since what they are intended to tell you is
    something that you already know: that the probability on each trial for
    the outcome to be A is approximately 1/12 (= .0833...), to be B is
    1/12, to be C is 1/12, etc.  You know that because you set up your
    program for that to be the case.  In "real life" you wouldn't know
    that, and the point is to determine if one outcome is more likely than
    another and by how much.

    Here are some observations you can make, based, in part, on your
    knowledge that the *real* underlying probability is 1/12.

	1) As the number of samples gets larger the confidence intervals
	   get narrower.

	2) As P increases from .9 to .95 to .99 the confidence intervals get
	   wider.

	3) The .9 confidence intervals contain the real probability value
           about 90% of the time.  The .95 confidence intervals contain the
           real probability value about 95% of the time and the .99
           confidence intervals contain the real probability value 99% of
           the time. (These are not independent observations, though, when
           the 1 time out of 100 event of the 99% confidence interval being
           "missed" occurs the .95 and .9 confidence intervals will
           obviously also be missed).

    Here are some experiments to try which will make the data "more
    interesting".

    EXPERIMENT 1:  Instead of generating a whole number between 1 and 12
    generate a number between 1 and 13 instead.  If the number is 13 change
    it to 1.  What this means is that outcomes 2 through 12 are still
    equally probable, but outcome 1 is twice as likely as any of them. Note
    the differences.

    EXPERIMENT 2:  Set up an array, P, indexed from 1 to 12.  Then do the
    following (in pseudocode -- sorry my BASIC is to rusty to be useful):

	Step 1:  FOR i FROM 1 TO 12 DO P[i] := urand(0.0, 1.0);
	    Where the call on "urand" returns a uniform random real value
	    between 0 and 1.

	Step 2: (Optional) sort_descending(P)
            This step will make the simulation run faster.  For this
            purpose, however, its probably not worth the trouble.  When
            this is finished, the largest of the 12 values generated in
            step 1 will be in P[1]; the next largest in P[2] and so on down
            to the smallest in P[12].

	Step 3: accum := 0;
		FOR i FROM 1 TO 12 DO
		    accum := accum + P[i];
		    P[i] := accum;
            Upon completion of this step P[1] will contain the value it had
            at the end of the previous step.  P[2] will have the sum of its
            old value and P[1]'s old value.  P[3] will have the sum of its
            old value and of P[2]'s and P[1]'s old value.  Etc. up to P[12]
            which will contain the total of all of the original values. 
            accum will also contain the total of all of the original
            values.

	Step 4: FOR i FROM 1 TO 12 DO P[i] := P[i]/accum;
            Upon completion of this step all the values will have been
            "normalized", with P[1] through P[11] less than 1 and P[12]
            equal to 1.  They will be in increasing order.

    What array P now contains is a representation for a "random"
    probability for each of the 12 outcomes.  Specifically, P[1] contains
    the probability that an outcome numbered 1 or less (i.e., outcome 1)
    will occur on any particular trial; P[2] contains the probability that
    an outcome numbered 2 or less (i.e., either outcome 1 or outcome 2)
    will occur on any particular trial; P[3] contains the probability that
    an outcome numbered 3 or less (i.e., outcome 1, outcome 2, or outcome
    3) will occur on any particular trial; up to P[12] containing the
    probability that an outcome numbered 12 or less (i.e., any of the
    outcomes) will occur on any particular trial (i.e., 1 since one of the
    12 outcomes must occur).  To use this array you do the following when
    you want to generate a random outcome.

	    X := urand(0.0, 1.0);
	    FOR i FROM 1 TO 12 DO
		IF (P[i] less-than-or-equal-to X)
		THEN
		    outcome := i;
		    LEAVE LOOP;

    The variable "outcome" will contain the number between 1 and 12
    representing the randomly chosen outcome.  Just conduct your simulated
    experiment as before but use this outcome rather than the all-equally-
    likely outcomes you used originally.  Now you *don't* know in advance
    what the values are meant to tell you and the simulation is more
    "realistic" in that regard.  (I wrote the above for clarity, there is
    obviously a lot of optimization that can be done in the initialization
    steps, but they are unlikely to make a noticeable difference in the
    overall simulation).

                                Topher

.23>    > (Although I occasionally get a negative value?)
>
>    You probably get the negative values with the P=.99, N=100 case.  The
>    formula I gave is an approximation.  The usual rule of thumb given for
>    this (and related) approximations is that the Nx (the number of samples
>    with outcome "x") should be at least 5 for it to be "reasonably
>    accurate".  The approximation, however, is weaker in the "tails", and
>    .99 is far enough out in the tail to really require more data.  Just
>    don't take the P=.99/N=100 case too seriously.

There is another way to look at the negative values, which may clarify why they
appear and what to do about them.

By definition, all the sample Fx's must be bounded by zero and one.  The 
Px's of the underlying distribution, by definition, are similarly bounded by
zero and one.  But the error bounds in .19 are based on an approximate
distribution function, the Z or normal distribution, which is not so bounded.

The result is that when Fx happens to be close to zero, and z(P) is too large,
the value computed for the lower confidence limit can become negative.

One solution is to quietly replace negative numbers with zero.

Another is to use a properly bounded distribution, in this case the beta 
distribution.  This means more computation, but MIPS are cheap these days.