Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: Generate random numbers with fixed sum and different constraints
Replies: 25   Last Post: Sep 17, 2013 2:28 AM

 Messages: [ Previous | Next ]
 JS Hong Posts: 3 Registered: 9/6/10
Re: Generate random numbers with fixed sum and different constraints
Posted: Sep 17, 2013 1:52 AM

Similar error happened to me with five variable.
Could anybody help me to fix this?

I copy your code and run it.
But with following parameters,

s = 183.8975;

lo = [ 0 0 0 0 0];

up = [420.0241 424.2856 648.9985 249.4882 113.4215];

the same error occurs.

Is this because s is small?

"Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message <k8d3bv\$a2s\$1@newscl01ah.mathworks.com>...
> "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.

Date Subject Author
11/14/12 Dmitrey Yershov
11/17/12 Greg Heath
11/17/12 Roger Stafford
11/17/12 Bruno Luong
11/17/12 Roger Stafford
11/18/12 Dmitrey Yershov
11/18/12 Bruno Luong
11/19/12 Bruno Luong
11/19/12 Dmitrey Yershov
11/19/12 Bruno Luong
11/19/12 Dmitrey Yershov
11/19/12 Bruno Luong
11/19/12 Bruno Luong
11/19/12 Roger Stafford
11/20/12 Bruno Luong
11/20/12 Bruno Luong
11/20/12 Roger Stafford
11/21/12 Bruno Luong
11/21/12 james bejon
11/20/12 james bejon
9/17/13 JS Hong
9/17/13 JS Hong
11/19/12 David Epstein
11/19/12 Roger Stafford
11/20/12 Bruno Luong
11/17/12 Greg Heath