Calculating Approximate Distribution from Trials

Question

Update: Now answered below; if anyone has any better answers, then please let me know - I'd be most appreciative. :)

My question is the following, and here is the setup (fairly standard Polya's urn question - note, there is a very similar question, but no-one has answered it; however, the point of my question is very different):

I have a bucket with one white and one black ball in. I (randomly, uniformly) choose one of the balls from the bucket; I then replace that ball along with another of the same colour (choose a white, put that white back in with another white, so now three balls, two of which are white). I wish to determine the limit distribution.

I have showed that the proportion of white balls (or of black) is a martingale, so the martingale convergence theorem gives convergence, finite almost surely.

Here is where my question moves away from the standard question. I am looking to model this situation computationally. (I've written a quick Matlab code; I could include this, but I don't think it would really help.) Now, choosing the number of steps to be large (say, 10^7), I attain samples of this. My question is the following:

How can I use these samples to approximate the distribution?

Should I divide the unit interval into, say, 200 (or 201) subintervals, then test to see whether the answer falls in this interval, then plot an histogram of the number of times that it falls in this region? I'm really looking for some better way than this!

Any advice would be most appreciated! Thanks!

(Also, please don't transfer this to computer science or anything - it's a maths question, and the fact that I'm computationally finding the results is neither here nor there!) :)

Answer 1

After a bit more thought, I have realised that I can actually plot a pretty accurate histogram plot. I can't think of a better way of approximating the distribution curve.

Here is the Matlab code:

M = 10^3; % This is the number of times the situation will be run
N = 10^4; % This is the number of 'turns' per operation
P = zeros(M,1);
for m = 1:M
    W = 1; B = 1;
    for n = 1:N
        p = W/(W+B);    % this is the proportion that are WHITE
        r = rand(1);
        if r > p        % this simulates picking a BLACK ball
            W = W; B = B+1;     % increase BLACKS by 1
        else
            W = W+1; B = B;     % increase WHITES by 1
        end
    end
    P(m) = W/(W+B);
end
x = 0:0.05:1;
hist(P,x)

Unfortunately, the scale on the left hand side is number of entries in the given range. It would be better if I could divide it all by $M$, the total number of runs through - I don't know how I would implement this though, and it doesn't really matter.

If anyone has any ideas as to how I can improve on this, then I'd be most appreciative! :)