%% Sampling from a probability distribution
% Toby Driscoll, December 2011

%%
% (Chebfun example stats/ResamplingRandomVariables.m)

%%
% A common problem in applications of random variables is to draw samples
% from a given distribution. It's easy to find computer codes for generating
% pseudorandom numbers that are distributed uniformly or normally, and
% these usually must be converted to simulate a different target
% distribution. The key steps are integration and function inversion, which
% Chebfun can do with great accuracy.

%%
%
LW = 'linewidth'; FS = 'fontsize'; MS = 'markersize';
splitting off

%% von Mises distribution
% The von Mises distribution is a periodic variant of the normal
% distribution. While the density is easily defined, it's otherwise not
% simple to work with analytically.
%
% We start with the density function, normalized to give total probability
% 1.
kappa = 1.5;
f = chebfun(@(x) exp(kappa*cos(x)),[-pi pi]);
density = f/sum(f);

%%
% Now we integrate to get the cumulative distribution function.
cdf = cumsum(density);
plot([density,cdf],LW,1.6), axis([-pi pi 0 1])
title('von Mises distribution',FS,12)
legend('density','distribution','Location','northwest')

%%
% Sampling from this distribution involves applying its inverse to
% uniformly sampled points. We could do this one-by-one using `roots`, but
% for a large number of points it is more efficient to find a chebfun
% for the inverse function with `inv`:
cdfinv = inv(cdf);
plot(cdfinv,LW,1.6)
title('Inverse of von Mises distribution',FS,12)

%%
% Now the resampling is easy. We compare the resulting histogram to the
% original von Mises density.
u = rand(1e4,1);                           % uniform
x = cdfinv(u);                             % von Mises
[count,bin] = hist(x,36);
count = count/sum(count*(bin(2)-bin(1)));  % renormalize, total area = 1
cla, bar(bin,count), hold on
plot(density,'r',LW,1.6), axis tight
title('Sampled points and the orignal density',FS,12)


%% Logit-normal distribution
% A more exotic and troublesome distribution is the logit-normal
% distribution. Its density and cdf are easy enough to define:
sig = 1.11;
f = @(x) exp( -(log(x./(1-x))).^2/(2*sig^2))./(x.*(1-x));
density = chebfun(f,[0 1]);
density = density/sum(density);
cdf = cumsum(density);
clf, plot([density,cdf],LW,1.6)
title('logit-normal distribution',FS,12)
legend('density','distribution','Location','northwest')

%%
% However, because $F'=f=0$ at the ends, the inverse function has 
% infinite slope at the ends, and a straightforward inversion will fail. To
% cope with this, we'll take some shortcuts. First, we'll use symmetry
% to restrict attention to $x> 1/2$. Second, we'll put Chebfun
% into splitting mode to help cope with the endpoint slope. Finally,
% we'll truncate the domain of the cdf slightly.
splitting on
cdfinv = inv( cdf{0.5,1-1e-3} );
clf, plot(cdfinv,LW,1.6)
title('Inverse of the logit-normal distribution',FS,12)

%%
% To apply the result for resampling, we have to reflect uniform values
% less than $1/2$ back into $[1/2,1]$, and reflect the results back.
u = rand(1e4,1);
flag = (u < 0.5);  u(flag) = 1-u(flag);
x = cdfinv( u );  x(flag) = 1-x(flag);
[count,bin] = hist(x,36);
count = count/sum(count*(bin(2)-bin(1)));  % renormalize, total area = 1
clf, bar(bin,count), hold on
plot(density,'r',LW,1.6), axis tight
title('Sampled points and the orignal density',FS,12)

%%
% We can see what our truncation of the original random variable costs us by
% looking at the domain of the inverse cdf:
cdfinv.ends.'
missing = 1 - ans(end)

%%
% Thus, a uniform variable that takes a value closer to 1 than this number
% won't be mapped accurately back to the logit-normal variable we want, unless we
% take further steps. Clearly, such events will be extremely rare.
