Bernoulli Trial problem

From Computational Statistics (CSE383M and CS395T)
Jump to navigation Jump to search

Team members: LoriL, Rcardenas, Dan.

This is the code our team developed for the analytical solution to the Bernoulli trials problem whose problem statement is located here. The runtime for N = 12 is rather large (~20-25 minutes), but our code computes the probabilities for all n which I have printed below.

function [probn] = Bernoulli(N,p,q)

num = 0; % counter for n
probn = zeros(1,N+1);
for pset = 0:2^N-1
    for qset = 0:2^(N-1)-1
        pbinary = dec2bin(pset,N); % make the binary set for p
        qbinary = dec2bin(qset,N-1); % make the binary set for q
        pones = length(find(pbinary == '1')); % number of 1's in set of p
        qones = length(find(qbinary == '1')); % number of 1's in set of q
        probPQset = p^(pones)*(1-p)^(N-pones)*q^(qones)*(1-q)^(N-1-qones); % probability for this set of p's and q's. 
        for i = 1:N-1
            if qbinary(i) == '1'
               pbinary(i+1) = pbinary(i);% if q = 1, repeat last trial
            end
        end
        for j = 1:N
            if pbinary(j) == '1'
                num = num + 1; % Total number of 1's.
            end
        end
        probn(num+1) = probn(num+1) + probPQset; 
        num = 0;
    end
   end
end

Bernoulli.png

For this bar chart, the x-axis are the n values plus 1. So the n value that has the greatest probability is 10 (not 11). The values for the probabilities are listed in the following table.

P(n): p=0.8, q=0.2, N=12
n values Probability of n
0 2.634e-6
1 3.0159e-5
2 1.9956e-4
3 9.6013e-4
4 0.0037
5 0.0115
6 0.0301
7 0.0668
8 0.1247
9 0.1924
10 0.2368
11 0.2153
12 0.1175