
Re: Generate random numbers with fixed sum and different constraints
Posted:
Nov 19, 2012 5:58 AM


"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 = [x0lo(:); up(:)x0]; > > % FEX: http://www.mathworks.com/matlabcentral/fileexchange/30892 > tol = 1e6; > [V,nr,nre] = lcon2vert(A,b,[],[],tol); > > % Split S in simplex > T = delaunayn(V); > P = V(T,:); > P = reshape(P, [size(T) m1]); > 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(m1, 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: runid 1283078028 delaunay Qtriangulate Qbboundlast Qcoplanarkeep _premerge _zerocentrum Qinteriorkeep Pgood _maxwidth 2.7 Errorroundoff 3.8e15 _onemerge 2.7e14 Visibledistance 7.6e15 Ucoplanardistance 7.6e15 Widthoutside 1.5e14 _widefacet 4.5e14
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.8e15. 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.1e16 facet p2 p1 p0 distance= 2.2e16 facet p2 p1 p3 distance= 1.1e16
These points either have a maximum or minimum xcoordinate, 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.8e15.  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 cocircular or cospherical:  use 'Qz' to add a point "at infinity" (i.e., above the paraboloid)  or use 'QJ' to joggle the input and avoid cocircular 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.

