Probability distributions (MATLAB)

From LiteratePrograms

Contents 1 Discrete Distributions 1.1 Multinomial 1.1.1 Parameters 1.1.2 Selection 1.1.3 Empirical distribution 2 Continuous distributions 2.1 Bimodal distribution 2.1.1 Parameters 2.1.2 Underlying distributions 2.1.3 Selection of each distribution 2.1.4 Whole script

Discrete Distributions

Multinomial

Parameters

The only parameters are the probability to obtain some events? Here, 5 events with probabiility (pct):

6.6667   13.3333   20.0000   26.6667   33.3333

<<params_multin>>=
p = [1 2 3 4 5];
p = p/sum(p)
q = cumsum(p)/sum(p)
nb = 200;

Selection

<<sript_multinomial.m>>=
params_multin
rs     = rand(nb,1);
chosen = sum(repmat(rs,1,length(q))>repmat(q,nb,1),2)+1
empirical_multin

Empirical distribution

Note the use of the sparse command; the empirical distribution is:

8.5000   13.5000   21.5000   26.0000   30.5000

<<empirical_multin>>=
rep = sparse(chosen,1,1);
full(rep)'/sum(rep)*100

Continuous distributions

Bimodal distribution

The question of how to generate bimodal distributions is often asked on the MATLAB newgroup. Here is a simple way, given that this example is made using two Gaussians, but any other distributions can be used and put in the X variable.

 

Parameters

The parameters of such a distribution are:

• the probability to be in one distribution or in the other
• their means
• thier standard deviations (for Gaussians)

<<gb_parameters>>=
q = [.4 .6];
m = [0 50];
s = [5 5];
distrib = struct('mu', m, 'sigma', s, 'weight', q)
nb = 200;

Underlying distributions

Their generation is quite direct.

<<gb_gaussians>>=
%% Construction of two scaled and translated gaussians
X       = randn(nb,length(q)).*repmat(s,nb,1)+repmat(m,nb,1);

Selection of each distribution

This is clearly the most difficult point.

<<gb_selection>>=
%% Selection of one gaussian or the other
rsel    = rand(nb,1);
idx1    = (repmat(rsel,1,length(q))>repmat(cumsum(q),nb,1));
idx2    = (repmat(rsel,1,length(q))<repmat(cumsum(q),nb,1));
idx1(:,2:end) = idx1(:,1:end-1).*idx2(:,2:end);
idx1(:,1) = idx2(:,1);
X       = sum(X.*idx1,2);

Whole script

<<script_gbimodal.m>>=
gb_parameters
gb_gaussians
gb_selection
%% Plot result
[y,x]=ksdensity(X);figure;plot(x,y,'linewidth',2)

Download code