expdesign.html
de096106547amaster

expdesign.html
View Options

			<!DOCTYPE HTML>
			<html>
			<head>
			<meta charset="UTF-8">
			<title>Section 7.5.2: Experiment design</title>
			<link rel="canonical" href="http://cvxr.com/cvx/examples/cvxbook/Ch07_statistical_estim/html/expdesign.html">
			<link rel="stylesheet" href="../../../examples.css" type="text/css">
			</head>
			<body>
			<div id="header">
			<h1>Section 7.5.2: Experiment design</h1>
			Jump to:
			<a href="#source">Source code</a>
			<a href="#output">Text output</a>

			<a href="#plots">Plots</a>
			<a href="../../../index.html">Library index</a>
			</div>
			<div id="content">
			<a id="source"></a>
			<pre class="codeinput">
			<span class="comment">% Boyd & Vandenberghe, "Convex Optimization"</span>
			<span class="comment">% Original version by Lieven Vandenberghe</span>
			<span class="comment">% Updated for CVX by Almir Mutapcic - Jan 2006</span>
			<span class="comment">% (a figure is generated)</span>
			<span class="comment">%</span>
			<span class="comment">% This is an example of D-optimal, A-optimal, and E-optimal</span>
			<span class="comment">% experiment designs.</span>

			<span class="comment">% problem data</span>
			m = 10;
			angles1 = linspace(3*pi/4,pi,m);
			angles2 = linspace(0,-pi/2,m);

			<span class="comment">% sensor positions</span>
			V = [3.0*[cos(angles1); sin(angles1)], <span class="keyword">...</span>
			1.5*[cos(angles2); sin(angles2)]];
			p = size(V,2);
			n = 2;
			noangles = 5000;

			<span class="comment">% D-optimal design</span>
			<span class="comment">%</span>
			<span class="comment">% maximize log det Vdiag(lambda)V'</span>
			<span class="comment">% subject to sum(lambda)=1, lambda >=0</span>
			<span class="comment">%</span>

			<span class="comment">% setup the problem and solve it</span>
			cvx_begin
			variable <span class="string">lambda(p)</span>
			maximize ( det_rootn( Vdiag(lambda)V' ) )
			subject <span class="string">to</span>
			sum(lambda) == 1;
			lambda >= 0;
			cvx_end
			lambdaD = lambda; <span class="comment">% save the solution for confidence ellipsoids</span>

			<span class="comment">% plot results</span>
			figure(1)
			<span class="comment">% draw ellipsoid v'Wv <= 2</span>
			W = inv(Vdiag(lambda)V');
			angles = linspace(0,2*pi,noangles);
			R = chol(W); <span class="comment">% W = R'*R</span>
			ellipsoid = sqrt(2)*(R\[cos(angles); sin(angles)]);
			d = plot(ellipsoid(1,:), ellipsoid(2,:), <span class="string">'--'</span>, 0,0,<span class="string">'+'</span>);
			set(d, <span class="string">'Color'</span>, [0 0.5 0]); set(d(2),<span class="string">'MarkerFaceColor'</span>,[0 0.5 0]);
			hold <span class="string">on</span>;

			dot=plot(V(1,:),V(2,:),<span class="string">'o'</span>);
			ind = find(lambda > 0.001);
			dots = plot(V(1,ind),V(2,ind),<span class="string">'o'</span>);
			set(dots,<span class="string">'MarkerFaceColor'</span>,<span class="string">'blue'</span>);

			<span class="comment">% print out nonzero lambda</span>
			disp(<span class="string">'Nonzero lambda values for D design:'</span>);
			<span class="keyword">for</span> i=1:length(ind)
			text(V(1,ind(i)),V(2,ind(i)), [<span class="string">'l'</span>,int2str(ind(i))]);
			disp([<span class="string">'lambda('</span>,int2str(ind(i)),<span class="string">') = '</span>, num2str(lambda(ind(i)))]);
			<span class="keyword">end</span>;

			<span class="comment">%axis([-4.5 4.5 -4.5 4.5])</span>
			axis([-5 5 -5 5])
			set(gca,<span class="string">'Xtick'</span>,[]);
			set(gca,<span class="string">'Ytick'</span>,[]);
			hold <span class="string">off</span>, axis <span class="string">off</span>
			<span class="comment">% print -deps Ddesign.eps</span>

			<span class="comment">% A-optimal design</span>
			<span class="comment">%</span>
			<span class="comment">% minimize Trace (sum_i lambdaivivi')^{-1}</span>
			<span class="comment">% subject to lambda >= 0, 1'*lambda = 1</span>
			<span class="comment">%</span>

			<span class="comment">% SDP formulation</span>
			e = eye(2,2);
			cvx_begin <span class="string">sdp</span>
			variables <span class="string">lambda(p)</span> <span class="string">u(n)</span>
			minimize ( sum(u) )
			subject <span class="string">to</span>
			<span class="keyword">for</span> k = 1:n
			[ Vdiag(lambda)V' e(:,k);
			e(k,:) u(k) ] >= 0;
			<span class="keyword">end</span>
			sum(lambda) == 1;
			lambda >= 0;
			cvx_end
			lambdaA = lambda; <span class="comment">% save the solution for confidence ellipsoids</span>

			<span class="comment">% plot results</span>
			figure(2)
			<span class="comment">% draw ellipsoid v'Wv <= mu</span>
			W = inv(Vdiag(lambda)V')^2;
			mu = diag(V'WV);
			mu = mean(mu(ind));
			angles = linspace(0,2*pi,noangles);
			R = chol(W); <span class="comment">% W = R'*R</span>
			ellipsoid = sqrt(mu)*(R\[cos(angles); sin(angles)]);
			d = plot(ellipsoid(1,:), ellipsoid(2,:), <span class="string">'--'</span>,0,0,<span class="string">'+'</span>);
			set(d, <span class="string">'Color'</span>, [0 0.5 0]);
			set(d(2), <span class="string">'MarkerFaceColor'</span>, [0 0.5 0]);
			hold <span class="string">on</span>

			dot = plot(V(1,:),V(2,:),<span class="string">'o'</span>);
			ind = find(lambda > 0.001);
			dots = plot(V(1,ind),V(2,ind),<span class="string">'o'</span>);
			set(dots,<span class="string">'MarkerFaceColor'</span>,<span class="string">'blue'</span>);

			disp(<span class="string">'Nonzero lambda values for A design:'</span>);
			<span class="keyword">for</span> i=1:length(ind)
			text(V(1,ind(i)),V(2,ind(i)), [<span class="string">'l'</span>,int2str(ind(i))]);
			disp([<span class="string">'lambda('</span>,int2str(ind(i)),<span class="string">') = '</span>, num2str(lambda(ind(i)))]);
			<span class="keyword">end</span>;
			<span class="comment">%axis([-4.5 4.5 -4.5 4.5])</span>
			axis([-5 5 -5 5])
			set(gca,<span class="string">'Xtick'</span>,[]);
			set(gca,<span class="string">'Ytick'</span>,[]);
			axis <span class="string">off</span>, hold <span class="string">off</span>
			<span class="comment">% print -deps Adesign.eps</span>

			<span class="comment">% E-optimal design</span>
			<span class="comment">%</span>
			<span class="comment">% minimize w</span>
			<span class="comment">% subject to sum_i lambda_ivivi' >= w*I</span>
			<span class="comment">% lambda >= 0, 1'*lambda = 1;</span>
			<span class="comment">%</span>

			cvx_begin <span class="string">sdp</span>
			variables <span class="string">t</span> <span class="string">lambda(p)</span>
			maximize ( t )
			subject <span class="string">to</span>
			Vdiag(lambda)V' >= t*eye(n,n);
			sum(lambda) == 1;
			lambda >= 0;
			cvx_end

			lambdaE = lambda; <span class="comment">% save the solution for confidence ellipsoids</span>

			figure(3)
			<span class="comment">% draw ellipsoid v'Wv <= mu</span>
			mu = diag(V'WV);
			mu = mean(mu(ind));
			angles = linspace(0,2*pi,noangles);
			R = chol(W); <span class="comment">% W = R'*R</span>
			ellipsoid = sqrt(mu)*(R\[cos(angles); sin(angles)]);
			d = plot(ellipsoid(1,:), ellipsoid(2,:), <span class="string">'--'</span>, 0, 0, <span class="string">'+'</span>);
			set(d, <span class="string">'Color'</span>, [0 0.5 0]);
			set(d(2), <span class="string">'MarkerFaceColor'</span>, [0 0.5 0]);
			hold <span class="string">on</span>

			dot = plot(V(1,:),V(2,:),<span class="string">'o'</span>);
			lambda = lambda(1:p);
			ind = find(lambda > 0.001);
			dots = plot(V(1,ind),V(2,ind),<span class="string">'o'</span>);
			set(dots,<span class="string">'MarkerFaceColor'</span>,<span class="string">'blue'</span>);

			disp(<span class="string">'Nonzero lambda values for E design:'</span>);
			<span class="keyword">for</span> i=1:length(ind)
			text(V(1,ind(i)),V(2,ind(i)), [<span class="string">'l'</span>,int2str(ind(i))]);
			disp([<span class="string">'lambda('</span>,int2str(ind(i)),<span class="string">') = '</span>, num2str(lambda(ind(i)))]);
			<span class="keyword">end</span>;
			<span class="comment">%axis([-4.5 4.5 -4.5 4.5])</span>
			axis([-5 5 -5 5])
			set(gca,<span class="string">'Xtick'</span>,[]);
			set(gca,<span class="string">'Ytick'</span>,[]);
			axis <span class="string">off</span>, hold <span class="string">off</span>
			<span class="comment">% print -deps Edesign.eps</span>


			<span class="comment">% confidence ellipsoids</span>
			eta = 6.2514; <span class="comment">% chi2inv(.9,3) value (command available in stat toolbox)</span>
			<span class="comment">% draw 90 percent confidence ellipsoid for D design</span>
			W = Vdiag(lambdaD)V';
			angles = linspace(0,2*pi,noangles);
			R = chol(W); <span class="comment">% W = R'*R</span>
			ellipsoid = sqrt(eta)*(R\[cos(angles); sin(angles)]);

			figure(4)
			plot(0,0,<span class="string">'ok'</span>,ellipsoid(1,:), ellipsoid(2,:), <span class="string">'-'</span>);
			text(ellipsoid(1,1100),ellipsoid(2,1100),<span class="string">'D'</span>);
			hold <span class="string">on</span>

			<span class="comment">% draw 90 percent confidence ellipsoid for A design</span>
			W = Vdiag(lambdaA)V';
			angles = linspace(0,2*pi,noangles);
			R = chol(W); <span class="comment">% W = R'*R</span>
			ellipsoid = sqrt(eta)*(R\[cos(angles); sin(angles)]);
			plot(0,0,<span class="string">'ok'</span>,ellipsoid(1,:), ellipsoid(2,:), <span class="string">'-'</span>);
			text(ellipsoid(1,1),ellipsoid(2,1),<span class="string">'A'</span>);

			<span class="comment">% draw 90 percent confidence ellipsoid for E design</span>
			W = Vdiag(lambdaE)V';
			angles = linspace(0,2*pi,noangles);
			R = chol(W); <span class="comment">% W = R'*R</span>
			ellipsoid = sqrt(eta)*(R\[cos(angles); sin(angles)]);
			d=plot(0,0,<span class="string">'ok'</span>,ellipsoid(1,:), ellipsoid(2,:), <span class="string">'-'</span>);
			set(d,<span class="string">'Color'</span>,[0 0.5 0]);
			text(ellipsoid(1,4000),ellipsoid(2,4000),<span class="string">'E'</span>);

			<span class="comment">% draw 90 percent confidence ellipsoid for uniform design</span>
			W_u = inv(V*V'/p);
			R = chol(W_u); <span class="comment">% W = R'*R</span>
			ellipsoid_u = sqrt(eta)*(R\[cos(angles); sin(angles)]);
			plot(ellipsoid_u(1,:), ellipsoid_u(2,:), <span class="string">'--'</span>);
			text(ellipsoid_u(1),ellipsoid_u(2),<span class="string">'U'</span>);
			set(gca,<span class="string">'Xtick'</span>,[]);
			set(gca,<span class="string">'Ytick'</span>,[]);
			axis <span class="string">off</span>
			<span class="comment">% print -deps confidence.eps</span>
			hold <span class="string">off</span>
			</pre>
			<a id="output"></a>
			<pre class="codeoutput">

			Calling SDPT3: 35 variables, 11 equality constraints
			------------------------------------------------------------

			num. of constraints = 11
			dim. of sdp var = 6, num. of sdp blk = 2
			dim. of linear var = 21
			dim. of free var = 1 *** convert ublk to lblk
			*******************************************************************
			SDPT3: Infeasible path-following algorithms
			*******************************************************************
			version predcorr gam expon scale_data
			HKM 1 0.000 1 0
			it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime
			-------------------------------------------------------------------
			0\|0.000\|0.000\|4.8e+02\|4.9e+01\|5.9e+03\|-2.875957e-10 0.000000e+00\| 0:0:00\| chol 1 1
			1\|0.875\|0.781\|6.0e+01\|1.1e+01\|9.0e+02\| 7.118271e+00 -2.554174e+01\| 0:0:00\| chol 1 1
			2\|0.943\|0.862\|3.5e+00\|1.5e+00\|1.0e+02\| 6.657847e+00 -3.837796e+01\| 0:0:00\| chol 1 1
			3\|0.930\|0.272\|2.4e-01\|1.1e+00\|4.3e+01\|-3.634210e-01 -3.600963e+01\| 0:0:00\| chol 1 1
			4\|1.000\|0.913\|3.3e-07\|9.4e-02\|4.6e+00\|-1.156221e+00 -5.373577e+00\| 0:0:00\| chol 1 1
			5\|0.975\|0.811\|2.8e-08\|1.8e-02\|1.3e+00\|-2.298010e+00 -3.488416e+00\| 0:0:00\| chol 1 1
			6\|0.796\|0.626\|6.6e-09\|6.7e-03\|6.3e-01\|-2.902415e+00 -3.493317e+00\| 0:0:00\| chol 1 1
			7\|0.990\|0.928\|1.2e-08\|4.8e-04\|3.6e-02\|-3.170678e+00 -3.203472e+00\| 0:0:00\| chol 1 1
			8\|0.985\|0.987\|1.0e-09\|6.1e-06\|4.7e-04\|-3.181819e+00 -3.182255e+00\| 0:0:00\| chol 1 1
			9\|0.989\|0.989\|2.2e-09\|3.5e-06\|1.4e-05\|-3.181979e+00 -3.181984e+00\| 0:0:00\| chol 1 1
			10\|1.000\|0.989\|2.0e-10\|1.0e-07\|4.2e-07\|-3.181980e+00 -3.181981e+00\| 0:0:00\| chol 1 1
			11\|1.000\|0.989\|2.4e-12\|3.2e-09\|1.3e-08\|-3.181981e+00 -3.181981e+00\| 0:0:00\|
			stop: max(relative gap, infeasibilities) < 1.49e-08
			-------------------------------------------------------------------
			number of iterations = 11
			primal objective value = -3.18198051e+00
			dual objective value = -3.18198052e+00
			gap := trace(XZ) = 1.28e-08
			relative gap = 1.74e-09
			actual relative gap = 6.93e-10
			rel. primal infeas = 2.41e-12
			rel. dual infeas = 3.15e-09
			norm(X), norm(y), norm(Z) = 1.7e+01, 4.4e+00, 7.9e+00
			norm(A), norm(b), norm(C) = 2.9e+01, 2.0e+00, 2.4e+00
			Total CPU time (secs) = 0.22
			CPU time per iteration = 0.02
			termination code = 0
			DIMACS: 2.4e-12 0.0e+00 3.8e-09 0.0e+00 6.9e-10 1.7e-09
			-------------------------------------------------------------------
			------------------------------------------------------------
			Status: Solved
			Optimal value (cvx_optval): +3.18198

			Nonzero lambda values for D design:
			lambda(1) = 0.5
			lambda(10) = 0.5

			Calling SDPT3: 32 variables, 11 equality constraints
			------------------------------------------------------------

			num. of constraints = 11
			dim. of sdp var = 6, num. of sdp blk = 2
			dim. of linear var = 20
			*******************************************************************
			SDPT3: Infeasible path-following algorithms
			*******************************************************************
			version predcorr gam expon scale_data
			HKM 1 0.000 1 0
			it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime
			-------------------------------------------------------------------
			0\|0.000\|0.000\|4.9e+02\|4.7e+01\|5.6e+03\| 2.000000e+01 0.000000e+00\| 0:0:00\| chol 1 1
			1\|0.926\|0.999\|3.6e+01\|1.2e-01\|4.2e+02\| 1.663772e+01 -2.905847e+01\| 0:0:00\| chol 1 1
			2\|0.950\|1.000\|1.8e+00\|1.0e-02\|4.1e+01\| 1.545943e+00 -2.335057e+01\| 0:0:00\| chol 1 1
			3\|1.000\|1.000\|3.0e-07\|1.0e-03\|6.4e+00\| 1.636822e+00 -4.755167e+00\| 0:0:00\| chol 1 1
			4\|0.917\|0.879\|1.2e-07\|2.1e-04\|8.7e-01\| 1.221243e+00 3.518119e-01\| 0:0:00\| chol 1 1
			5\|1.000\|1.000\|2.4e-08\|1.0e-05\|5.9e-01\| 1.101610e+00 5.104357e-01\| 0:0:00\| chol 1 1
			6\|0.973\|0.976\|2.4e-09\|1.2e-06\|9.8e-02\| 9.168349e-01 8.188437e-01\| 0:0:00\| chol 1 1
			7\|1.000\|1.000\|1.7e-09\|1.0e-07\|4.1e-02\| 8.744063e-01 8.336361e-01\| 0:0:00\| chol 1 1
			8\|0.971\|0.970\|8.9e-10\|1.3e-08\|1.5e-03\| 8.505205e-01 8.490159e-01\| 0:0:00\| chol 1 1
			9\|0.988\|0.986\|3.2e-10\|1.3e-09\|1.9e-05\| 8.495403e-01 8.495208e-01\| 0:0:00\| chol 1 1
			10\|0.999\|0.997\|2.4e-13\|6.9e-11\|3.3e-07\| 8.495281e-01 8.495278e-01\| 0:0:00\| chol 1 1
			11\|1.000\|1.000\|7.0e-12\|1.0e-12\|1.1e-08\| 8.495279e-01 8.495279e-01\| 0:0:00\|
			stop: max(relative gap, infeasibilities) < 1.49e-08
			-------------------------------------------------------------------
			number of iterations = 11
			primal objective value = 8.49527931e-01
			dual objective value = 8.49527921e-01
			gap := trace(XZ) = 1.08e-08
			relative gap = 3.99e-09
			actual relative gap = 3.99e-09
			rel. primal infeas = 7.02e-12
			rel. dual infeas = 1.00e-12
			norm(X), norm(y), norm(Z) = 8.1e+00, 1.7e+00, 2.8e+00
			norm(A), norm(b), norm(C) = 4.0e+01, 2.7e+00, 2.4e+00
			Total CPU time (secs) = 0.18
			CPU time per iteration = 0.02
			termination code = 0
			DIMACS: 9.6e-12 0.0e+00 1.2e-12 0.0e+00 4.0e-09 4.0e-09
			-------------------------------------------------------------------
			------------------------------------------------------------
			Status: Solved
			Optimal value (cvx_optval): +0.849528

			Nonzero lambda values for A design:
			lambda(1) = 0.2966
			lambda(10) = 0.37797
			lambda(20) = 0.32543

			Calling SDPT3: 24 variables, 4 equality constraints
			------------------------------------------------------------

			num. of constraints = 4
			dim. of sdp var = 2, num. of sdp blk = 1
			dim. of linear var = 20
			dim. of free var = 1 *** convert ublk to lblk
			*******************************************************************
			SDPT3: Infeasible path-following algorithms
			*******************************************************************
			version predcorr gam expon scale_data
			HKM 1 0.000 1 0
			it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime
			-------------------------------------------------------------------
			0\|0.000\|0.000\|4.8e+02\|4.7e+01\|5.2e+03\|-2.875957e-10 0.000000e+00\| 0:0:00\| chol 1 1
			1\|0.970\|0.720\|1.5e+01\|1.3e+01\|2.3e+02\| 5.239168e+00 -1.929069e+01\| 0:0:00\| chol 1 1
			2\|1.000\|0.979\|1.9e-06\|2.8e-01\|2.5e+01\| 3.638124e+00 -1.987258e+01\| 0:0:00\| chol 1 1
			3\|1.000\|0.493\|5.1e-07\|1.4e-01\|1.5e+01\|-2.120931e-01 -1.518800e+01\| 0:0:00\| chol 1 1
			4\|1.000\|0.880\|5.0e-06\|1.7e-02\|2.0e+00\|-6.785315e-01 -2.651558e+00\| 0:0:00\| chol 1 1
			5\|1.000\|0.093\|4.7e-08\|1.6e-02\|1.5e+00\|-1.131853e+00 -2.590894e+00\| 0:0:00\| chol 1 1
			6\|1.000\|0.631\|8.8e-08\|5.8e-03\|7.5e-01\|-1.431151e+00 -2.171634e+00\| 0:0:00\| chol 1 1
			7\|1.000\|0.583\|4.6e-08\|2.4e-03\|2.4e-01\|-1.724221e+00 -1.960989e+00\| 0:0:00\| chol 1 1
			8\|0.998\|0.959\|2.2e-08\|9.8e-05\|8.1e-03\|-1.797089e+00 -1.805145e+00\| 0:0:00\| chol 1 1
			9\|0.988\|0.988\|8.7e-10\|1.2e-06\|9.5e-05\|-1.799966e+00 -1.800060e+00\| 0:0:00\| chol 1 1
			10\|0.989\|0.989\|1.7e-10\|9.4e-07\|3.1e-06\|-1.800000e+00 -1.800001e+00\| 0:0:00\| chol 1 1
			11\|1.000\|0.989\|1.7e-11\|3.0e-08\|1.0e-07\|-1.800000e+00 -1.800000e+00\| 0:0:00\| chol 1 1
			12\|1.000\|0.989\|3.3e-12\|1.0e-09\|3.3e-09\|-1.800000e+00 -1.800000e+00\| 0:0:00\|
			stop: max(relative gap, infeasibilities) < 1.49e-08
			-------------------------------------------------------------------
			number of iterations = 12
			primal objective value = -1.80000000e+00
			dual objective value = -1.80000000e+00
			gap := trace(XZ) = 3.28e-09
			relative gap = 7.13e-10
			actual relative gap = 2.56e-10
			rel. primal infeas = 3.35e-12
			rel. dual infeas = 1.01e-09
			norm(X), norm(y), norm(Z) = 2.0e+00, 2.1e+00, 5.8e+00
			norm(A), norm(b), norm(C) = 2.9e+01, 2.0e+00, 2.4e+00
			Total CPU time (secs) = 0.20
			CPU time per iteration = 0.02
			termination code = 0
			DIMACS: 3.3e-12 0.0e+00 1.2e-09 0.0e+00 2.6e-10 7.1e-10
			-------------------------------------------------------------------
			------------------------------------------------------------
			Status: Solved
			Optimal value (cvx_optval): +1.8

			Nonzero lambda values for E design:
			lambda(10) = 0.2
			lambda(20) = 0.8
			</pre>
			<a id="plots"></a>
			<div id="plotoutput">
			<img src="expdesign__01.png" alt=""> <img src="expdesign__02.png" alt=""> <img src="expdesign__03.png" alt=""> <img src="expdesign__04.png" alt="">
			</div>
			</div>
			</body>
			</html>

R1252/simulation-code/cvx/examples/cvxbook/Ch07_statistical_estim/html/expdesign.htmlde096106547amaster

expdesign.htmlView Options

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

expdesign.html
View Options