% Example 7.2: Maximum entropy distribution % Section 7.2, Figures 7.2-7.3 % Boyd & Vandenberghe, "Convex Optimization" % Originally by Lieven Vandenberghe % Adapted for CVX by Michael Grant 4/11/06 % % We consider a probability distribution on 100 equidistant points in the % interval [-1,1]. We impose the following prior assumptions: % % -0.1 <= E(X) <= +0.1 % +0.5 <= E(X^2) <= +0.6 % -0.3 <= E(3*X^3-2*X) <= -0.2 % +0.3 <= Pr(X<0) <= 0.4 % % Along with the constraints sum(p) == 1, p >= 0, these constraints % describe a polyhedron of probability distrubtions. In the first figure, % the distribution that maximizes entropy is computed. In the second % figure, we compute upper and lower bounds for Prob(X<=a_i) for each % point -1 <= a_i <= +1 in the distribution, as well as the maximum % entropy CDF. % % Represent the polyhedron as follows: % A * p <= b % sum( p ) == 1 % p >= 0 % n = 100; a = linspace(-1,1,n); a2 = a .^ 2; a3 = 3 * ( a.^ 3 ) - 2 * a; ap = +( a < 0 ); A = [ a ; -a ; a2 ; -a2 ; a3 ; -a3 ; ap ; -ap ]; b = [ 0.1 ; 0.1 ;0.5 ; -0.5 ; -0.2 ; 0.3 ; 0.4 ; -0.3 ]; % % Compute the maximum entropy distribution % cvx_expert true cvx_begin variables pent(n) maximize( sum(entr(pent)) ) sum(pent) == 1; A * pent <= b; cvx_end % % Compute the bounds on Prob(X<=a_i), i=1,...,n % Ubnds = zeros(1,n); Lbnds = zeros(1,n); for t = 1 : n, cvx_begin quiet variable p( n ) minimize sum( p(1:t) ) p >= 0; sum( p ) == 1; A * p <= b; cvx_end Lbnds(t) = cvx_optval; cvx_begin quiet variable p( n ) maximize sum( p(1:t) ) p >= 0; sum( p ) == 1; A * p <= b; cvx_end Ubnds(t) = cvx_optval; disp( sprintf( '%g <= Prob(x<=%g) <= %g', Lbnds(t), a(t), Ubnds(t) ) ); end % % Generate the figures % figure( 1 ) stairs( a, pent ); xlabel( 'x' ); ylabel( 'PDF( x )' ); figure( 2 ) stairs( a, cumsum( pent ) ); grid on; hold on d = stairs(a, Lbnds,'r-'); set(d,'Color',[0 0.5 0]); d = stairs(a, Ubnds,'r-'); set(d,'Color',[0 0.5 0]); d = plot([-1,-1], [Lbnds(1), Ubnds(1)],'r-'); set(d,'Color',[0 0.5 0]); axis([-1.1 1.1 -0.1 1.1]); xlabel( 'x' ); ylabel( 'CDF( x )' ); hold off