Probability distributions (MATLAB)

From LiteratePrograms
Jump to: navigation, search

Contents

[edit] Discrete Distributions

[edit] Multinomial

Density of the MATLAB-generated multinomial distribution

[edit] 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;

[edit] 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

[edit] 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

[edit] Continuous distributions

[edit] Bimodal distribution

Density of the MATLAB-generated 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.



 

[edit] 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;

[edit] 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);

[edit] 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);

[edit] 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
hijacker
hijacker
hijacker
hijacker