Homec4science
Diffusion EMPoWER (master) Browse

R1252/simulation-code/cvx/examples/cvxbook/Ch07_statistical_estim/html/probbounds.html
de096106547amaster

Last Change
Raw File

probbounds.html
View Options

<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Section 7.4.3: Probability bounds example with Voronoi diagram</title>
<link rel="canonical" href="http://cvxr.com/cvx/examples/cvxbook/Ch07_statistical_estim/html/probbounds.html">
<link rel="stylesheet" href="../../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Section 7.4.3: Probability bounds example with Voronoi diagram</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% Boyd &amp; Vandenberghe "Convex Optimization"</span>
<span class="comment">% Original version by Lieven Vandenberghe</span>
<span class="comment">% Adapted for CVX by Michael Grant, 2005/12/19</span>
<span class="comment">% Generates figures 7.5-7.8</span>
<span class="comment">% The constellation points. Feel free to change them, but they must</span>
<span class="comment">% produce a valid Voronoi diagram. Therefore, there must be three or</span>
<span class="comment">% more points, and no three can be collinear. To test your selected</span>
<span class="comment">% points, run VORONOI( Cs(:,1), Cs(:,2) ) and see if a complete diagram</span>
<span class="comment">% is drawn; if so, your points should work.</span>
Cs = [ <span class="keyword">...</span>
0, 0 ; <span class="keyword">...</span>
1.2, 2.4 ; <span class="keyword">...</span>
-3, +3 ; <span class="keyword">...</span>
-4, 0 ; <span class="keyword">...</span>
-1.6, -3.2 ; <span class="keyword">...</span>
1.84615384615385, -2.76923076923077 ; <span class="keyword">...</span>
2.35294117647059, 0.58823529411765 ];
Cmax = max(max(abs(Cs))) * 1.25;
<span class="comment">% Plot the constellation points and the Voronoi tesselation</span>
clf
Cx = Cs( :, 1 );
Cy = Cs( :, 2 );
m = length( Cx );
Cs = Cs';
[ Vx, Vy ] = voronoi( Cx, Cy );
plot( Vx, Vy, <span class="string">'b-'</span>, Cx, Cy, <span class="string">'o'</span> );
axis <span class="string">equal</span>
axis( Cmax * [ -1, 1, -1, 1 ] );
axis <span class="string">off</span>
hold <span class="string">on</span>
<span class="comment">% Draw unit circles around each constellation point</span>
noangles = 200;
angles = linspace( 0, 2 * pi, noangles );
crcpts = [ cos(angles) ; sin(angles) ];
<span class="keyword">for</span> i=1 : m,
text( Cx(i)+0.25, Cy(i)+0.25, [ <span class="string">'s'</span>, int2str(i) ] );
ellipse = [ cos(angles) ; sin(angles) ] + Cs(:,i) * ones(1,noangles);
plot( ellipse(1,:), ellipse(2,:), <span class="string">':'</span> );
<span class="keyword">end</span>;
<span class="comment">% print -deps chebbnds_example.eps</span>
<span class="comment">% Determine the polyhedrons for each Voronoi region by computing the</span>
<span class="comment">% Delaunay triangulation; that is, matrices A and b such that</span>
<span class="comment">% A * ( x - c ) &lt;= b</span>
<span class="comment">% where c is the constellation point. The faces of a polyhedron for a given</span>
<span class="comment">% point consist of the perpindicular bisectors of edges of the Delaunay</span>
<span class="comment">% triangles to which it belongs.</span>
m = size( Cs, 2 );
tri = delaunay( Cx, Cy );
ee = sparse( tri, tri( :, [ 3, 1, 2 ] ), 1, m, m );
ee = ee + ee';
<span class="keyword">for</span> k = 1 : m,
v2 = find( ee( :, k ) );
pk = Cs( :, v2 );
qk = Cs( :, k ) * ones( 1, length( v2 ) );
Ak = pk - qk;
bk = 0.5 * sum( Ak .* Ak, 1 );
As{k} = Ak';
bs{k} = bk';
<span class="keyword">end</span>
<span class="comment">% For each polyhedron, compute lower bounds on the probability of</span>
<span class="comment">% correct detection with sigma = 1. Check the results by plotting the</span>
<span class="comment">% ellipsoid x'*P*x + 2*q'*x + r = 1, which should inscribe the polyhedron.</span>
ints = 1 : m;
<span class="comment">% Uncomment to do only the first polyhedron, like the book does</span>
<span class="comment">% ints = 1;</span>
<span class="keyword">for</span> i = ints( : ).',
[ cd_cheb, P, q, r, X, lambda ] = cheb( As{i}, bs{i}, eye(2) );
ellipse = sqrt(1-r+q'*(P\q)) * P^(-1/2) * crcpts + <span class="keyword">...</span>
(-P\q + Cs(:,i)) * ones(1,noangles);
plot( ellipse(1,:), ellipse(2,:), <span class="string">'r-'</span> );
dots = plot( X(1,:)+Cx(i), X(2,:)+Cy(i), <span class="string">'ro'</span> );
set( dots, <span class="string">'MarkerFaceColor'</span>, <span class="string">'red'</span> );
set( dots, <span class="string">'MarkerSize'</span>, 4 );
<span class="keyword">end</span>
hold <span class="string">off</span>
<span class="comment">% print -deps chebbnds_example2.eps</span>
<span class="comment">% Compute Chebyshev lower bounds for Prob( As(i) * v &lt;= bs(i) )</span>
<span class="comment">% where v = N(Cs(i),sigma) for varying values of sigma</span>
nsigma = 500;
sigmas = linspace( 0.001, 6.0, nsigma )';
cd_cheb = zeros( nsigma, m );
fprintf( <span class="string">'Computing lower bounds'</span> );
<span class="comment">% Uncomment to match the book</span>
ints = 1 : m;
<span class="comment">% ints = 1 : 3;</span>
<span class="keyword">for</span> i = ints( : ).',
<span class="keyword">for</span> k = 1 : nsigma,
cd_cheb(k,i) = cheb( As{i}, bs{i}, sigmas(k) * eye(2) );
<span class="keyword">end</span>;
<span class="keyword">if</span> rem( k, 10 ) == 0,
fprintf( <span class="string">'.'</span> );
<span class="keyword">end</span>
<span class="keyword">end</span>;
fprintf( <span class="string">'done.\n'</span> );
figure(2)
mc = size( cd_cheb, 2 );
plot(sqrt(sigmas(:,ones(1,mc))), cd_cheb);
<span class="keyword">for</span> i = 1 : mc,
text( sqrt(sigmas(nsigma/4)), cd_cheb(nsigma/4,i), [<span class="string">'b'</span>,int2str(i)] );
<span class="keyword">end</span>;
xlabel(<span class="string">'x'</span>);
ylabel(<span class="string">'y'</span>);
axis( [ 0, 2.5, 0, 1 ] );
<span class="comment">% For the central set, compute Chebyshev lower bounds, Monte Carlo</span>
<span class="comment">% estimates, and Chernoff bounds.</span>
<span class="comment">% for central set, compute cheb lower bounds, mc estimates,</span>
<span class="comment">% and chernoff bounds</span>
<span class="comment">%</span>
nsigma = 50;
sigmas = linspace( 0.1, 0.5, nsigma );
cd1 = zeros( 1, nsigma ); <span class="comment">% lower bounds for prob( x in C1 )</span>
mc1 = zeros( 1, nsigma ); <span class="comment">% monte carlo estimates of prob( x in C1 )</span>
cher1 = zeros( m-1, nsigma ); <span class="comment">% chernoff upper bounds on Prob( x in Cj | s = s_1 )</span>
fprintf( <span class="string">'Computing lower bounds and Monte Carlo sims'</span> );
<span class="keyword">for</span> i = 1 : nsigma,
<span class="comment">% Compute the Chebyshev lower bound on Prob( As{1} * v &lt;= bs{1} )</span>
<span class="comment">% for v in N( 0, Sigma )</span>
Sigma = sigmas(i)^2 * eye(2);
cd1(i) = cheb( As{1}, bs{1}, Sigma );
mc1(i) = montecarlo( As{1}, bs{1}, Sigma, 10000 );
<span class="keyword">if</span> rem( i, 5 ) == 0,
fprintf( <span class="string">'.'</span> );
<span class="keyword">end</span>
<span class="keyword">end</span>
fprintf( <span class="string">'done.\nComputing upper bounds'</span> );
<span class="keyword">for</span> j = 2 : m,
A = As{j};
b = bs{j} - A * ( Cs(:,1) - Cs(:,j) );
<span class="comment">% Compute the Chernoff upper bound on</span>
<span class="comment">% Prob( As{j} * ( v + Cs{1} - Cs{j} ) &lt;= bs{j} )</span>
<span class="comment">% for v in N( 0, Sigma )</span>
<span class="keyword">for</span> i = 1 : nsigma,
cher1( j - 1, i ) = cher( A, b, sigmas(i)^2*eye(2) );
<span class="keyword">end</span>
fprintf( <span class="string">'.'</span> );
<span class="keyword">end</span>;
fprintf( <span class="string">'done.\n'</span> );
cher1 = max( 1 - sum( cher1 ), 0 );
figure(4)
plot( sigmas, cher1, <span class="string">'-'</span>, sigmas, mc1, <span class="string">'--'</span> );
axis( [ 0.2 0.5 0.9 1 ] );
xlabel( <span class="string">'x'</span> );
ylabel( <span class="string">'y'</span> );
<span class="comment">%print -deps chernoff_example.eps</span>
</pre>
<a id="output"></a>
<pre class="codeoutput">
Computing lower bounds.......done.
Computing lower bounds and Monte Carlo sims..........done.
Computing upper bounds......done.
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="probbounds__01.png" alt=""> <img src="probbounds__02.png" alt=""> <img src="probbounds__03.png" alt="">
</div>
</div>
</body>
</html>
c4science · Help