Path: news.mathworks.com!not-for-mail
From: "Dmitrey Yershov" <pierrevanstulov@mail.ru>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Generate random numbers with fixed sum and different constraints
Date: Mon, 19 Nov 2012 10:58:07 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 141
Message-ID: <k8d3bv$a2s$1@newscl01ah.mathworks.com>
References: <k80995$9eo$1@newscl01ah.mathworks.com> <k87s11$h6u$1@newscl01ah.mathworks.com> <k88l7g$7tm$1@newscl01ah.mathworks.com> <k88qag$nm1$1@newscl01ah.mathworks.com> <k89ehg$rm6$1@newscl01ah.mathworks.com> <k8b4e9$ja7$1@newscl01ah.mathworks.com> <k8badi$8rd$1@newscl01ah.mathworks.com> <k8ck6h$lhp$1@newscl01ah.mathworks.com>
Reply-To: "Dmitrey Yershov" <pierrevanstulov@mail.ru>
NNTP-Posting-Host: www-06-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1353322687 10332 172.30.248.38 (19 Nov 2012 10:58:07 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Mon, 19 Nov 2012 10:58:07 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 3853991
Xref: news.mathworks.com comp.soft-sys.matlab:783057

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8ck6h$lhp$1@newscl01ah.mathworks.com>...
> Sorry, I miss the definition of the variable 'tol' in my previous posted code; here we go again:
> 
> % lower + upper bounds
> lo = [1 0 2 1 0];
> up = [5 6 4 4 2];
> % target sum
> s = 15;
> % number of samples
> n = 1000;
> 
> %% Engine
> m = length(lo);
> slo = sum(lo);
> sup = sum(up);
> 
> % A particular solution, x0
> x0 = interp1([slo; sup], [lo; up], s);
> if any(isnan(x0))
>     error('no solution exists')
> end
> 
> % feasible set: S = { x=x0+A*y : A*y <= b }
> x0 = x0(:);
> B = null(ones(size(lo)));
> A = [-B; B];
> b = [x0-lo(:); up(:)-x0];
> 
> % FEX: http://www.mathworks.com/matlabcentral/fileexchange/30892
> tol = 1e-6;
> [V,nr,nre] = lcon2vert(A,b,[],[],tol);
> 
> % Split S in simplex
> T = delaunayn(V);
> P = V(T,:);
> P = reshape(P, [size(T) m-1]);
> P = permute(P, [3 2 1]);
> np = size(P,3);
> 
> % Compute the volume of simplexes
> Q = bsxfun(@minus, P, P(:,end,:));
> Q(:,end,:) = [];
> vol = arrayfun(@(k) det(Q(:,:,k)), 1:np);
> vol = abs(vol);
> 
> % generate n samples (i) of v with probability ~ vol
> vol = vol / sum(vol);
> c = cumsum(vol); c(end)=1;
> [~, i] = histc(rand(1,n),[0 c]);
> 
> % Random rarycentric coordinates
> % FEX: http://www.mathworks.com/matlabcentral/fileexchange/9700
> w = randfixedsum(m,n,1,0,1);
> 
> % Put together
> V = P(:,:,i);
> y = zeros(m-1, n);
> for k=1:n
>     y(:,k) = V(:,:,k)*w(:,k);
> end
> % final result, m x n array of random vector
> x = bsxfun(@plus, x0, B*y);
> 
> % Bruno

Thank you very much, Bruno! I think that your procedure is the card! BUT now  when I try to generate points with the following input parameters:

lo = [0 0 0];
up = [2 1 2];
s = 2;
n = 1000;

it gives me the error:

Error using qhullmx
QH6154 qhull precision error: initial facet 1 is coplanar with the interior point
ERRONEOUS FACET:

While executing:  | qhull d Qt Qbb Qc
Options selected for Qhull 2010.1 2010/01/14:
  run-id 1283078028  delaunay  Qtriangulate  Qbbound-last  Qcoplanar-keep
  _pre-merge  _zero-centrum  Qinterior-keep  Pgood  _max-width 2.7
  Error-roundoff 3.8e-15  _one-merge 2.7e-14  Visible-distance 7.6e-15
  U-coplanar-distance 7.6e-15  Width-outside 1.5e-14  _wide-facet 4.5e-14

The input to qhull appears to be less than 3 dimensional, or a
computation has overflowed.

Qhull could not construct a clearly convex simplex from points:

The center point is coplanar with a facet, or a vertex is coplanar
with a neighboring facet.  The maximum round off error for
computing distances is 3.8e-15.  The center point, facets and distances
to the center point are as follows:


facet p1 p3 p0 distance=    0
facet p2 p3 p0 distance= -1.1e-16
facet p2 p1 p0 distance= -2.2e-16
facet p2 p1 p3 distance= -1.1e-16

These points either have a maximum or minimum x-coordinate, or
they maximize the determinant for k coordinates.  Trial points
are first selected from points that maximize a coordinate.

The min and max coordinates for each dimension are:
  0:   -0.8392    0.8928  difference= 1.732
  1:    -1.239     1.493  difference= 2.732
  2:         0     2.732  difference= 2.732

If the input should be full dimensional, you have several options that
may determine an initial simplex:
  - use 'QJ'  to joggle the input and make it full dimensional
  - use 'QbB' to scale the points to the unit cube
  - use 'QR0' to randomly rotate the input for different maximum points
  - use 'Qs'  to search all points for the initial simplex
  - use 'En'  to specify a maximum roundoff error less than 3.8e-15.
  - trace execution with 'T3' to see the determinant for each point.

If the input is lower dimensional:
  - use 'QJ' to joggle the input and make it full dimensional
  - use 'Qbk:0Bk:0' to delete coordinate k from the input.  You should
    pick the coordinate with the least range.  The hull will have the
    correct topology.
  - determine the flat containing the points, rotate the points
    into a coordinate plane, and delete the other coordinates.
  - add one or more points to make the input full dimensional.


This is a Delaunay triangulation and the input is co-circular or co-spherical:
  - use 'Qz' to add a point "at infinity" (i.e., above the paraboloid)
  - or use 'QJ' to joggle the input and avoid co-circular data


Error in delaunayn (line 101)
t = qhullmx(x', 'd ', opt);

Error in Untitled (line 32)
T = delaunayn(V);

I will be obliged to you if you can fixe it.