http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503
MATLAB Central Newsreader  Generate random numbers with fixed sum and different constraints
Feed for thread: Generate random numbers with fixed sum and different constraints
enus
©19942014 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Wed, 14 Nov 2012 14:19:17 +0000
Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891676
Dmitrey Yershov
Hello. I need to generate nonnegative rundom numbers sum of which is equal to 1. Each number xi is constrained: ai<=xi<=bi. How can I do this? Similar question was solved here<br>
<br>
<a href="http://www.mathworks.com/matlabcentral/fileexchange/9700">http://www.mathworks.com/matlabcentral/fileexchange/9700</a><br>
<br>
but in this alghorithm a<=xi<=b (a and b are the same for all xi). Any ideas?

Sat, 17 Nov 2012 11:22:09 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891879
Greg Heath
"Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message <k80995$9eo$1@newscl01ah.mathworks.com>...<br>
> Hello. I need to generate nonnegative rundom numbers sum of which is equal to 1. Each number xi is constrained: ai<=xi<=bi. How can I do this? Similar question was solved here<br>
> <br>
> <a href="http://www.mathworks.com/matlabcentral/fileexchange/9700">http://www.mathworks.com/matlabcentral/fileexchange/9700</a><br>
> <br>
> but in this alghorithm a<=xi<=b (a and b are the same for all xi). Any ideas?<br>
<br>
Z = a + (ba)*rand(m,n);<br>
<br>
sumZ = repmat(sum(Z),m,1);<br>
<br>
I'll let you figure out the rest.<br>
<br>
Hope this helps<br>
<br>
Greg

Sat, 17 Nov 2012 18:32:16 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891890
Roger Stafford
"Greg Heath" <heath@alumni.brown.edu> wrote in message <k87s11$h6u$1@newscl01ah.mathworks.com>...<br>
> "Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message <k80995$9eo$1@newscl01ah.mathworks.com>...<br>
> > Hello. I need to generate nonnegative rundom numbers sum of which is equal to 1. Each number xi is constrained: ai<=xi<=bi. How can I do this? Similar question was solved here<br>
> > <br>
> > <a href="http://www.mathworks.com/matlabcentral/fileexchange/9700">http://www.mathworks.com/matlabcentral/fileexchange/9700</a><br>
> > <br>
> > but in this alghorithm a<=xi<=b (a and b are the same for all xi). Any ideas?<br>
> <br>
> Z = a + (ba)*rand(m,n);<br>
> <br>
> sumZ = repmat(sum(Z),m,1);<br>
> <br>
> I'll let you figure out the rest.<br>
> <br>
> Hope this helps<br>
> <br>
> Greg<br>
         <br>
I'd like to know the answer to that too, Greg. Suppose your m = 3, n = 1, a = 0, and b = 2/3, and suppose z comes up randomly with z = [.6;.1;.1] as is possible. How is that point supposed to be projected onto a plane so as to have a sum of 1? A simple division by its sum(z) = .8 gives [.75;.125;.125] which exceeds the stated limit.<br>
<br>
The space of points in this case having a sum of 1 within the permitted threedimensional cube cuts it in half in a planar hexagon, and it is difficult to see how to project all points in the cube onto this hexagon in a simple manner using just the sum(z), never mind doing so in an areawise uniform manner throughout the hexagon.<br>
<br>
This is one of the reasons I went to the trouble of writing 'randfixedsum' which breaks such a hexagon into triangles and deals with each separately. However doing a similar thing with an ndimensional rectangular box other than an ndimensional cube as Dmitrey wishes to do is a project involving much more complicated simplex structure. At the moment I have no idea how to set about such a task.<br>
<br>
Roger Stafford

Sat, 17 Nov 2012 19:59:12 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891898
Bruno Luong
"Roger Stafford" wrote in message <k88l7g$7tm$1@newscl01ah.mathworks.com>...<br>
> "Greg Heath" <heath@alumni.brown.edu> wrote in message <k87s11$h6u$1@newscl01ah.mathworks.com>...<br>
<br>
> <br>
> This is one of the reasons I went to the trouble of writing 'randfixedsum' which breaks such a hexagon into triangles and deals with each separately. However doing a similar thing with an ndimensional rectangular box other than an ndimensional cube as Dmitrey wishes to do is a project involving much more complicated simplex structure. At the moment I have no idea how to set about such a task.<br>
<br>
Might be convert linear inequality into hull vertices. Use Delaunay to decompose the convex polytopes into union of simplexes then apply the barycenter coordinates to generate uniform distribution on simplexes. Put all that together, it should be able to generate the uniform distribution with required constraints.<br>
<br>
Bruno

Sun, 18 Nov 2012 01:44:16 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891906
Roger Stafford
"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k88qag$nm1$1@newscl01ah.mathworks.com>...<br>
> Might be convert linear inequality into hull vertices. Use Delaunay to decompose the convex polytopes into union of simplexes then apply the barycenter coordinates to generate uniform distribution on simplexes. Put all that together, it should be able to generate the uniform distribution with required constraints.<br>
> <br>
> Bruno<br>
         <br>
Yes, that is a conceivable approach, Bruno. However it faces some formidable difficulties with large dimensionality. For an ndimensional cube with n equal to, say, 51 the number of vertices in the n1 dimensional polytope with a fixed sum set at or near the halfway point would be 51!/25!^2 = 6,446,940,928,325,352. This would present quite a challenge to 'delaunayn'! Moreover, besides providing for the uniform distribution of each simplex, one must also choose different simplices in proportion to their n1 dimensional volumes, and certainly with the unequal bound values Dmitrey has requested there would be a vast number of different volumes among them to compute. For a practical program to be able to handle large values of n there should be some underlying symmetry principle that greatly simplifies such proceedings as these, and at the moment I can't think what that might be for <br>
unequal bounding intervals. <br>
<br>
Roger Stafford

Sun, 18 Nov 2012 02:05:23 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891907
Greg Heath
"Roger Stafford" wrote in message <k88l7g$7tm$1@newscl01ah.mathworks.com>...<br>
> "Greg Heath" <heath@alumni.brown.edu> wrote in message <k87s11$h6u$1@newscl01ah.mathworks.com>...<br>
> > "Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message <k80995$9eo$1@newscl01ah.mathworks.com>...<br>
> > > Hello. I need to generate nonnegative rundom numbers sum of which is equal to 1. Each number xi is constrained: ai<=xi<=bi. How can I do this? Similar question was solved here<br>
> > > <br>
> > > <a href="http://www.mathworks.com/matlabcentral/fileexchange/9700">http://www.mathworks.com/matlabcentral/fileexchange/9700</a><br>
> > > <br>
> > > but in this alghorithm a<=xi<=b (a and b are the same for all xi). Any ideas?<br>
> > <br>
> > Z = a + (ba)*rand(m,n);<br>
> > <br>
> > sumZ = repmat(sum(Z),m,1);<br>
> > <br>
> > I'll let you figure out the rest.<br>
> > <br>
> > Hope this helps<br>
> > <br>
> > Greg<br>
>          <br>
> I'd like to know the answer to that too, Greg. Suppose your m = 3, n = 1, a = 0, and b = 2/3, and suppose z comes up randomly with z = [.6;.1;.1] as is possible. How is that point supposed to be projected onto a plane so as to have a sum of 1? A simple division by its sum(z) = .8 gives [.75;.125;.125] which exceeds the stated limit.<br>
<br>
WHOOPS! Silly mistake.<br>
<br>
Please pardon my youthful exuberance.<br>
<br>
Greg<br>
<br>
P.S. The OP did not indicate that the number of points is specified. Was that an unintended omission?

Sun, 18 Nov 2012 17:04:09 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891929
Dmitrey Yershov
"Roger Stafford" wrote in message <k89ehg$rm6$1@newscl01ah.mathworks.com>...<br>
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k88qag$nm1$1@newscl01ah.mathworks.com>...<br>
> > Might be convert linear inequality into hull vertices. Use Delaunay to decompose the convex polytopes into union of simplexes then apply the barycenter coordinates to generate uniform distribution on simplexes. Put all that together, it should be able to generate the uniform distribution with required constraints.<br>
> > <br>
> > Bruno<br>
>          <br>
> Yes, that is a conceivable approach, Bruno. However it faces some formidable difficulties with large dimensionality. For an ndimensional cube with n equal to, say, 51 the number of vertices in the n1 dimensional polytope with a fixed sum set at or near the halfway point would be 51!/25!^2 = 6,446,940,928,325,352. This would present quite a challenge to 'delaunayn'! Moreover, besides providing for the uniform distribution of each simplex, one must also choose different simplices in proportion to their n1 dimensional volumes, and certainly with the unequal bound values Dmitrey has requested there would be a vast number of different volumes among them to compute. For a practical program to be able to handle large values of n there should be some underlying symmetry principle that greatly simplifies such proceedings as these, and at the moment I can't think what that might be <br>
for <br>
> unequal bounding intervals. <br>
> <br>
> Roger Stafford<br>
<br>
Thanks to everyone who reply! In fact, I'll be happy if I have method for n<=5 which is able to generate 1000 points in less than 2 sec (processor 2,4 Ghs).

Sun, 18 Nov 2012 18:46:10 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891933
Bruno Luong
I have not fully tested, but it should go like this:<br>
<br>
% lower + upper bounds<br>
lo = [1 0 2 1 0];<br>
up = [5 6 4 4 2];<br>
% target sum<br>
s = 15;<br>
% number of samples<br>
n = 1000;<br>
<br>
m = length(lo);<br>
slo = sum(lo);<br>
sup = sum(up);<br>
<br>
x0 = interp1([slo; sup], [lo; up], s);<br>
if any(isnan(x0))<br>
error('no solution exists')<br>
end<br>
<br>
x0 = x0(:);<br>
B = null(ones(size(lo)));<br>
A = [B; B];<br>
b = [x0lo(:); up(:)x0];<br>
<br>
% FEX: <a href="http://www.mathworks.com/matlabcentral/fileexchange/30892">http://www.mathworks.com/matlabcentral/fileexchange/30892</a><br>
% Thank you Matt<br>
[V,nr,nre] = lcon2vert(A,b,[],[],tol);<br>
<br>
T = delaunayn(V);<br>
P = V(T,:);<br>
P = reshape(P, [size(T) m1]);<br>
P = permute(P, [3 2 1]);<br>
np = size(P,3);<br>
<br>
% <a href="http://www.mathworks.com/matlabcentral/fileexchange/9700">http://www.mathworks.com/matlabcentral/fileexchange/9700</a><br>
% % Thank you Roger<br>
Q = bsxfun(@minus, P, P(:,end,:));<br>
Q(:,end,:) = [];<br>
vol = arrayfun(@(k) det(Q(:,:,k)), 1:np);<br>
vol = abs(vol);<br>
w = randfixedsum(m,n,1,0,1);<br>
<br>
% generate n samples of v with probability vol<br>
vol = vol / sum(vol);<br>
c = cumsum(vol); c(end)=1;<br>
[~, i] = histc(rand(1,n),[0 c]);<br>
<br>
V = P(:,:,i);<br>
y = zeros(m1, n);<br>
for k=1:n<br>
y(:,k) = V(:,:,k)*w(:,k);<br>
end<br>
x = bsxfun(@plus, x0, B*y);<br>
<br>
% Bruno

Mon, 19 Nov 2012 06:39:13 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891955
Bruno Luong
Sorry, I miss the definition of the variable 'tol' in my previous posted code; here we go again:<br>
<br>
% lower + upper bounds<br>
lo = [1 0 2 1 0];<br>
up = [5 6 4 4 2];<br>
% target sum<br>
s = 15;<br>
% number of samples<br>
n = 1000;<br>
<br>
%% Engine<br>
m = length(lo);<br>
slo = sum(lo);<br>
sup = sum(up);<br>
<br>
% A particular solution, x0<br>
x0 = interp1([slo; sup], [lo; up], s);<br>
if any(isnan(x0))<br>
error('no solution exists')<br>
end<br>
<br>
% feasible set: S = { x=x0+A*y : A*y <= b }<br>
x0 = x0(:);<br>
B = null(ones(size(lo)));<br>
A = [B; B];<br>
b = [x0lo(:); up(:)x0];<br>
<br>
% FEX: <a href="http://www.mathworks.com/matlabcentral/fileexchange/30892">http://www.mathworks.com/matlabcentral/fileexchange/30892</a><br>
tol = 1e6;<br>
[V,nr,nre] = lcon2vert(A,b,[],[],tol);<br>
<br>
% Split S in simplex<br>
T = delaunayn(V);<br>
P = V(T,:);<br>
P = reshape(P, [size(T) m1]);<br>
P = permute(P, [3 2 1]);<br>
np = size(P,3);<br>
<br>
% Compute the volume of simplexes<br>
Q = bsxfun(@minus, P, P(:,end,:));<br>
Q(:,end,:) = [];<br>
vol = arrayfun(@(k) det(Q(:,:,k)), 1:np);<br>
vol = abs(vol);<br>
<br>
% generate n samples (i) of v with probability ~ vol<br>
vol = vol / sum(vol);<br>
c = cumsum(vol); c(end)=1;<br>
[~, i] = histc(rand(1,n),[0 c]);<br>
<br>
% Random rarycentric coordinates<br>
% FEX: <a href="http://www.mathworks.com/matlabcentral/fileexchange/9700">http://www.mathworks.com/matlabcentral/fileexchange/9700</a><br>
w = randfixedsum(m,n,1,0,1);<br>
<br>
% Put together<br>
V = P(:,:,i);<br>
y = zeros(m1, n);<br>
for k=1:n<br>
y(:,k) = V(:,:,k)*w(:,k);<br>
end<br>
% final result, m x n array of random vector<br>
x = bsxfun(@plus, x0, B*y);<br>
<br>
% Bruno

Mon, 19 Nov 2012 10:26:09 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891964
David Epstein
"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8badi$8rd$1@newscl01ah.mathworks.com>...<br>
> I have not fully tested, but it should go like this:<br>
> <br>
> % lower + upper bounds<br>
> lo = [1 0 2 1 0];<br>
> up = [5 6 4 4 2];<br>
> % target sum<br>
> s = 15;<br>
> % number of samples<br>
> n = 1000;<br>
> <br>
> m = length(lo);<br>
> slo = sum(lo);<br>
> sup = sum(up);<br>
> <br>
> x0 = interp1([slo; sup], [lo; up], s);<br>
> if any(isnan(x0))<br>
> error('no solution exists')<br>
> end<br>
> <br>
> x0 = x0(:);<br>
> B = null(ones(size(lo)));<br>
> A = [B; B];<br>
> b = [x0lo(:); up(:)x0];<br>
> <br>
> % FEX: <a href="http://www.mathworks.com/matlabcentral/fileexchange/30892">http://www.mathworks.com/matlabcentral/fileexchange/30892</a><br>
> % Thank you Matt<br>
> [V,nr,nre] = lcon2vert(A,b,[],[],tol);<br>
> <br>
> T = delaunayn(V);<br>
> P = V(T,:);<br>
> P = reshape(P, [size(T) m1]);<br>
> P = permute(P, [3 2 1]);<br>
> np = size(P,3);<br>
> <br>
> % <a href="http://www.mathworks.com/matlabcentral/fileexchange/9700">http://www.mathworks.com/matlabcentral/fileexchange/9700</a><br>
> % % Thank you Roger<br>
> Q = bsxfun(@minus, P, P(:,end,:));<br>
> Q(:,end,:) = [];<br>
> vol = arrayfun(@(k) det(Q(:,:,k)), 1:np);<br>
> vol = abs(vol);<br>
> w = randfixedsum(m,n,1,0,1);<br>
> <br>
> % generate n samples of v with probability vol<br>
> vol = vol / sum(vol);<br>
> c = cumsum(vol); c(end)=1;<br>
> [~, i] = histc(rand(1,n),[0 c]);<br>
> <br>
> V = P(:,:,i);<br>
> y = zeros(m1, n);<br>
> for k=1:n<br>
> y(:,k) = V(:,:,k)*w(:,k);<br>
> end<br>
> x = bsxfun(@plus, x0, B*y);<br>
> <br>
> % Bruno<br>
<br>
I would be interested to know how this approach compares with straightforward rejection sampling, that is, not cutting up into simplices. If the dimension is n, the relevant subspace is defined using 2*n inequalities and one equality, so, in dimension 51 only 102 inequalities. I was thinking of the following procedure:<br>
<br>
1. Let f be the change of scale map from Dmitrey's box B to the unit cube U. f sends Dimitrey's probability measure to the standard probability measure.<br>
2. Let V be the subspace of B defined by one linear equality and 2n inequalities. Then f(V) in U is defined by one linear equality and 2n inequalities.<br>
3. Use GramSchmidt to construct linear isometry g that sends U into an ndimensional euclidean space Z and f(V) into the subspace z_1=0. The input data for GramSchmidt has its first vector v a unit vector orthogonal to f(V), followed by the standard basis, orthogonal to the various (n1)dimensional faces of U. (Perhaps omit the basis vector nearest to v.)<br>
4. Using max and min, find the smallest rectangle R in Z containing g(U), and hence a rectangle S, obtained from R by setting z_1=0. S contains g(f(V)).<br>
5. Generate a uniform random sample in S and keep only those satisfying the appropriate 2*n inequalities.<br>
<br>
I suppose that the problem with this method is that the volume of g(f(V)) may be small compared with the volume of S. But maybe the ratio could be improved by tricks, for example "adaptive sampling"whenever a random point p of S turns out not to lie in g(f(V)), then we can use the coordinates of p to try to chop off part of S while maintaining the rectangular shape relative to the coordinate axes.<br>
<br>
@Roger: what were your reasons for rejecting this approach in your randfixedsum package on FEX?<br>
<br>
David

Mon, 19 Nov 2012 10:58:07 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891966
Dmitrey Yershov
"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8ck6h$lhp$1@newscl01ah.mathworks.com>...<br>
> Sorry, I miss the definition of the variable 'tol' in my previous posted code; here we go again:<br>
> <br>
> % lower + upper bounds<br>
> lo = [1 0 2 1 0];<br>
> up = [5 6 4 4 2];<br>
> % target sum<br>
> s = 15;<br>
> % number of samples<br>
> n = 1000;<br>
> <br>
> %% Engine<br>
> m = length(lo);<br>
> slo = sum(lo);<br>
> sup = sum(up);<br>
> <br>
> % A particular solution, x0<br>
> x0 = interp1([slo; sup], [lo; up], s);<br>
> if any(isnan(x0))<br>
> error('no solution exists')<br>
> end<br>
> <br>
> % feasible set: S = { x=x0+A*y : A*y <= b }<br>
> x0 = x0(:);<br>
> B = null(ones(size(lo)));<br>
> A = [B; B];<br>
> b = [x0lo(:); up(:)x0];<br>
> <br>
> % FEX: <a href="http://www.mathworks.com/matlabcentral/fileexchange/30892">http://www.mathworks.com/matlabcentral/fileexchange/30892</a><br>
> tol = 1e6;<br>
> [V,nr,nre] = lcon2vert(A,b,[],[],tol);<br>
> <br>
> % Split S in simplex<br>
> T = delaunayn(V);<br>
> P = V(T,:);<br>
> P = reshape(P, [size(T) m1]);<br>
> P = permute(P, [3 2 1]);<br>
> np = size(P,3);<br>
> <br>
> % Compute the volume of simplexes<br>
> Q = bsxfun(@minus, P, P(:,end,:));<br>
> Q(:,end,:) = [];<br>
> vol = arrayfun(@(k) det(Q(:,:,k)), 1:np);<br>
> vol = abs(vol);<br>
> <br>
> % generate n samples (i) of v with probability ~ vol<br>
> vol = vol / sum(vol);<br>
> c = cumsum(vol); c(end)=1;<br>
> [~, i] = histc(rand(1,n),[0 c]);<br>
> <br>
> % Random rarycentric coordinates<br>
> % FEX: <a href="http://www.mathworks.com/matlabcentral/fileexchange/9700">http://www.mathworks.com/matlabcentral/fileexchange/9700</a><br>
> w = randfixedsum(m,n,1,0,1);<br>
> <br>
> % Put together<br>
> V = P(:,:,i);<br>
> y = zeros(m1, n);<br>
> for k=1:n<br>
> y(:,k) = V(:,:,k)*w(:,k);<br>
> end<br>
> % final result, m x n array of random vector<br>
> x = bsxfun(@plus, x0, B*y);<br>
> <br>
> % Bruno<br>
<br>
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:<br>
<br>
lo = [0 0 0];<br>
up = [2 1 2];<br>
s = 2;<br>
n = 1000;<br>
<br>
it gives me the error:<br>
<br>
Error using qhullmx<br>
QH6154 qhull precision error: initial facet 1 is coplanar with the interior point<br>
ERRONEOUS FACET:<br>
<br>
While executing:  qhull d Qt Qbb Qc<br>
Options selected for Qhull 2010.1 2010/01/14:<br>
runid 1283078028 delaunay Qtriangulate Qbboundlast Qcoplanarkeep<br>
_premerge _zerocentrum Qinteriorkeep Pgood _maxwidth 2.7<br>
Errorroundoff 3.8e15 _onemerge 2.7e14 Visibledistance 7.6e15<br>
Ucoplanardistance 7.6e15 Widthoutside 1.5e14 _widefacet 4.5e14<br>
<br>
The input to qhull appears to be less than 3 dimensional, or a<br>
computation has overflowed.<br>
<br>
Qhull could not construct a clearly convex simplex from points:<br>
<br>
The center point is coplanar with a facet, or a vertex is coplanar<br>
with a neighboring facet. The maximum round off error for<br>
computing distances is 3.8e15. The center point, facets and distances<br>
to the center point are as follows:<br>
<br>
<br>
facet p1 p3 p0 distance= 0<br>
facet p2 p3 p0 distance= 1.1e16<br>
facet p2 p1 p0 distance= 2.2e16<br>
facet p2 p1 p3 distance= 1.1e16<br>
<br>
These points either have a maximum or minimum xcoordinate, or<br>
they maximize the determinant for k coordinates. Trial points<br>
are first selected from points that maximize a coordinate.<br>
<br>
The min and max coordinates for each dimension are:<br>
0: 0.8392 0.8928 difference= 1.732<br>
1: 1.239 1.493 difference= 2.732<br>
2: 0 2.732 difference= 2.732<br>
<br>
If the input should be full dimensional, you have several options that<br>
may determine an initial simplex:<br>
 use 'QJ' to joggle the input and make it full dimensional<br>
 use 'QbB' to scale the points to the unit cube<br>
 use 'QR0' to randomly rotate the input for different maximum points<br>
 use 'Qs' to search all points for the initial simplex<br>
 use 'En' to specify a maximum roundoff error less than 3.8e15.<br>
 trace execution with 'T3' to see the determinant for each point.<br>
<br>
If the input is lower dimensional:<br>
 use 'QJ' to joggle the input and make it full dimensional<br>
 use 'Qbk:0Bk:0' to delete coordinate k from the input. You should<br>
pick the coordinate with the least range. The hull will have the<br>
correct topology.<br>
 determine the flat containing the points, rotate the points<br>
into a coordinate plane, and delete the other coordinates.<br>
 add one or more points to make the input full dimensional.<br>
<br>
<br>
This is a Delaunay triangulation and the input is cocircular or cospherical:<br>
 use 'Qz' to add a point "at infinity" (i.e., above the paraboloid)<br>
 or use 'QJ' to joggle the input and avoid cocircular data<br>
<br>
<br>
Error in delaunayn (line 101)<br>
t = qhullmx(x', 'd ', opt);<br>
<br>
Error in Untitled (line 32)<br>
T = delaunayn(V);<br>
<br>
I will be obliged to you if you can fixe it.

Mon, 19 Nov 2012 12:04:13 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891974
Bruno Luong
"Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message <br>
> Error in delaunayn (line 101)<br>
> t = qhullmx(x', 'd ', opt);<br>
> <br>
> Error in Untitled (line 32)<br>
> T = delaunayn(V);<br>
> <br>
> I will be obliged to you if you can fixe it.<br>
<br>
Please try to replace the delaunayn() command with:<br>
<br>
if m==3<br>
T = delaunay(V);<br>
else<br>
T = delaunayn(V);<br>
end<br>
<br>
Obviously it doesn' work for 2dimensional triangulation.<br>
<br>
Bruno

Mon, 19 Nov 2012 12:29:08 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891981
Dmitrey Yershov
"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8d77s$lfd$1@newscl01ah.mathworks.com>...<br>
> "Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message <br>
> > Error in delaunayn (line 101)<br>
> > t = qhullmx(x', 'd ', opt);<br>
> > <br>
> > Error in Untitled (line 32)<br>
> > T = delaunayn(V);<br>
> > <br>
> > I will be obliged to you if you can fixe it.<br>
> <br>
> Please try to replace the delaunayn() command with:<br>
> <br>
> if m==3<br>
> T = delaunay(V);<br>
> else<br>
> T = delaunayn(V);<br>
> end<br>
> <br>
> Obviously it doesn' work for 2dimensional triangulation.<br>
> <br>
> Bruno<br>
<br>
Thank you very VERY much, Bruno! Your procedure works fast enough to solve my problem and gives proper results (at least by sight for 3 dimensions)! <br>
<br>
P.S.: Special thanks to Roger Stafford and Matt J for basis procedures.

Mon, 19 Nov 2012 12:39:13 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891982
Bruno Luong
"Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message <br>
> <br>
> Thank you very VERY much, Bruno! Your procedure works fast enough to solve my problem and gives proper results (at least by sight for 3 dimensions)! <br>
<br>
You could further optimize the code by replacing the forloop of V*w with James's MTIMESX <br>
<br>
<a href="http://www.mathworks.com/matlabcentral/fileexchange/25977mtimesxfastmatrixmultiplywithmultidimensionalsupport">http://www.mathworks.com/matlabcentral/fileexchange/25977mtimesxfastmatrixmultiplywithmultidimensionalsupport</a><br>
<br>
Bruno

Mon, 19 Nov 2012 13:02:11 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#891985
Bruno Luong
It seems to me this task<br>
w = randfixedsum(m,n,1,0,1)<br>
<br>
can be optimized as well. But let save this optimization topic for another day.<br>
<br>
Bruno

Mon, 19 Nov 2012 23:10:18 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#892016
Roger Stafford
"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8dakj$39u$1@newscl01ah.mathworks.com>...<br>
> It seems to me this task<br>
> w = randfixedsum(m,n,1,0,1)<br>
> can be optimized as well. .....<br>
        <br>
Yes, I agree, Bruno; 'randfixedsum' is rather cumbersome to use just for a single simplex. How about this instead:<br>
<br>
R = bsxfun(@power,rand(m1,n),1./(m1:1:1)');<br>
X = cumprod([ones(1,n);R]).*[ones(m,n)[R;zeros(1,n)]];<br>
<br>
Roger Stafford

Tue, 20 Nov 2012 01:54:09 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#892023
Roger Stafford
"David Epstein" <David.Epstein.spam@remove.warwick.ac.uk> wrote in message <k8d1g1$4gr$1@newscl01ah.mathworks.com>...<br>
> @Roger: what were your reasons for rejecting this approach in your randfixedsum package on FEX?<br>
         <br>
I think you will find that as n increases the acceptance rate in such a procedure shrinks toward zero altogether too rapidly, thereby restricting one in practice to a rather small range for n.<br>
<br>
To get a feeling for this, consider an ndimensional hypersphere of radius 1/2 enclosed in an ndimensional cube with unitlength sides. The ndimensional volume of the cube is 1 whereas that of the hypersphere for even n is (pi/4)^(n/2)/(n/2)! (See <a href="http://en.wikipedia.org/wiki/Nsphere.">http://en.wikipedia.org/wiki/Nsphere.</a>) For n = 50 this would be 1.53E28, a small acceptance rate indeed.<br>
<br>
Roger Stafford

Tue, 20 Nov 2012 06:53:08 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#892042
Bruno Luong
"Roger Stafford" wrote in message <k8ee8q$k6f$1@newscl01ah.mathworks.com>...<br>
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8dakj$39u$1@newscl01ah.mathworks.com>...<br>
> > It seems to me this task<br>
> > w = randfixedsum(m,n,1,0,1)<br>
> > can be optimized as well. .....<br>
>         <br>
> Yes, I agree, Bruno; 'randfixedsum' is rather cumbersome to use just for a single simplex. How about this instead:<br>
> <br>
> R = bsxfun(@power,rand(m1,n),1./(m1:1:1)');<br>
> X = cumprod([ones(1,n);R]).*[ones(m,n)[R;zeros(1,n)]];<br>
> <br>
> Roger Stafford<br>
<br>
Really cool Roger! I would love if you could share the justification of this formula.<br>
<br>
Bruno

Tue, 20 Nov 2012 07:02:09 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#892043
Bruno Luong
Just a small typo correction in a comment of my code. It should read:<br>
<br>
% feasible set: S = { x=x0+B*y : A*y <= b }<br>
<br>
% Bruno

Tue, 20 Nov 2012 07:17:08 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#892046
Bruno Luong
"Roger Stafford" wrote in message <k8ens1$kfr$1@newscl01ah.mathworks.com>...<br>
> "David Epstein" <David.Epstein.spam@remove.warwick.ac.uk> wrote in message <k8d1g1$4gr$1@newscl01ah.mathworks.com>...<br>
> > @Roger: what were your reasons for rejecting this approach in your randfixedsum package on FEX?<br>
>          <br>
> I think you will find that as n increases the acceptance rate in such a procedure shrinks toward zero altogether too rapidly, thereby restricting one in practice to a rather small range for n.<br>
> <br>
> To get a feeling for this, consider an ndimensional hypersphere of radius 1/2 enclosed in an ndimensional cube with unitlength sides. The ndimensional volume of the cube is 1 whereas that of the hypersphere for even n is (pi/4)^(n/2)/(n/2)! (See <a href="http://en.wikipedia.org/wiki/Nsphere.">http://en.wikipedia.org/wiki/Nsphere.</a>) For n = 50 this would be 1.53E28, a small acceptance rate indeed.<br>
<br>
IMHO, the difficulty is in the same order than solving linear programing with the constraints A*y <= b, which have to find accurately the vertices without ambiguity. At least for the convex part.<br>
<br>
I have little experience in Delaunay in high dimensional space. But I could imagine it 's challenging too.<br>
<br>
Bruno

Tue, 20 Nov 2012 08:16:14 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#892049
james bejon
"Roger Stafford" wrote in message <k8ee8q$k6f$1@newscl01ah.mathworks.com>...<br>
> R = bsxfun(@power,rand(m1,n),1./(m1:1:1)');<br>
> X = cumprod([ones(1,n);R]).*[ones(m,n)[R;zeros(1,n)]];<br>
<br>
Very neat! (Though I have no idea how it works).

Tue, 20 Nov 2012 19:54:09 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#892079
Roger Stafford
"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8f9ck$gng$1@newscl01ah.mathworks.com>...<br>
> Really cool Roger! I would love if you could share the justification of this formula.<br>
<br>
"james bejon" wrote in message <k8fe8e$2lo$1@newscl01ah.mathworks.com>...<br>
> Very neat! (Though I have no idea how it works).<br>
       <br>
To Bruno, James, and anyone else interested:<br>
<br>
I recommended the two lines of code<br>
<br>
R = bsxfun(@power,rand(m1,n),1./(m1:1:1)');<br>
X = cumprod([ones(1,n);R]).*[ones(m,n)[R;zeros(1,n)]];<br>
<br>
as a substitute for Bruno's<br>
<br>
w = randfixedsum(m,n,1,0,1)<br>
<br>
(with X here instead of w) because 'randfixedsum' is designed for handling complicated polytopes which break down to into a great many individual simplices, whereas Bruno's call represents a single m1 dimensional simplex within R^m. (It's analogous to using a cannon to swat flies.) It is more efficient to compute the single simplex directly without having to bother with provisions for a decomposition into simplices. Its vertices are the m points:<br>
<br>
[1;0;0;0;...], [0;1;0;0;...], [0;0;1;0;...], ..., [0;0;0;...;0;0;1]<br>
<br>
which we can call P1, P2, ..., Pm.<br>
<br>
For example letting m = 4 we have a 3D tetrahedron simplex imbedded in R^4 space, and as can be seen the two code lines, in effect, carry out the computation<br>
<br>
r1 = rand^(1/3); r2 = rand^(1/2); r3 = rand^(1);<br>
X = [ (1)*(1r1) ; (1*r1)*(1r2) ; (1*r1*r2)*(1r3) ; (1*r1*r2*r3)*(10) ];<br>
<br>
which with the definitions of the points is the same as<br>
<br>
X = (1r1)*P1 + r1*((1r2)*P2 + r2*((1r3)*P3 + r3*P4));<br>
<br>
This distributes points within the tetrahedron in a statistically uniform manner. The product coefficients of the four points amount to barycentric coordinates with a sum of one. The expression<br>
<br>
(1r3)*P3 + r3*P4<br>
<br>
would distribute uniformly along the line between P3 and P4 using r3 = rand. The quantity<br>
<br>
(1r2)*P2 + r2*((1r3)*P3 + r3*P4)<br>
<br>
would then distribute uniformly within the triangle P2P3P4 with r2 = rand^(1/2) since the area of a triangle varies as the square of its height. Finally the full expression with r1 = rand^(1/3) gives a uniform distribution within the tetrahedron P1P2P3P4 since its volume is proportional to the cube of its height.<br>
<br>
Roger Stafford

Wed, 21 Nov 2012 07:38:17 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#892101
Bruno Luong
Thanks Roger for the explanation.<br>
<br>
Bruno

Wed, 21 Nov 2012 21:14:09 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#892142
james bejon
Very nice. I hadn't thought of it in those terms.

Tue, 17 Sep 2013 05:52:09 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#911350
JS Hong
Similar error happened to me with five variable. <br>
Could anybody help me to fix this?<br>
<br>
I copy your code and run it. <br>
But with following parameters, <br>
<br>
s = 183.8975;<br>
<br>
lo = [ 0 0 0 0 0];<br>
<br>
up = [420.0241 424.2856 648.9985 249.4882 113.4215];<br>
<br>
the same error occurs. <br>
<br>
Is this because s is small?<br>
<br>
<br>
<br>
<br>
"Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message <k8d3bv$a2s$1@newscl01ah.mathworks.com>...<br>
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8ck6h$lhp$1@newscl01ah.mathworks.com>...<br>
> > Sorry, I miss the definition of the variable 'tol' in my previous posted code; here we go again:<br>
> > <br>
> > % lower + upper bounds<br>
> > lo = [1 0 2 1 0];<br>
> > up = [5 6 4 4 2];<br>
> > % target sum<br>
> > s = 15;<br>
> > % number of samples<br>
> > n = 1000;<br>
> > <br>
> > %% Engine<br>
> > m = length(lo);<br>
> > slo = sum(lo);<br>
> > sup = sum(up);<br>
> > <br>
> > % A particular solution, x0<br>
> > x0 = interp1([slo; sup], [lo; up], s);<br>
> > if any(isnan(x0))<br>
> > error('no solution exists')<br>
> > end<br>
> > <br>
> > % feasible set: S = { x=x0+A*y : A*y <= b }<br>
> > x0 = x0(:);<br>
> > B = null(ones(size(lo)));<br>
> > A = [B; B];<br>
> > b = [x0lo(:); up(:)x0];<br>
> > <br>
> > % FEX: <a href="http://www.mathworks.com/matlabcentral/fileexchange/30892">http://www.mathworks.com/matlabcentral/fileexchange/30892</a><br>
> > tol = 1e6;<br>
> > [V,nr,nre] = lcon2vert(A,b,[],[],tol);<br>
> > <br>
> > % Split S in simplex<br>
> > T = delaunayn(V);<br>
> > P = V(T,:);<br>
> > P = reshape(P, [size(T) m1]);<br>
> > P = permute(P, [3 2 1]);<br>
> > np = size(P,3);<br>
> > <br>
> > % Compute the volume of simplexes<br>
> > Q = bsxfun(@minus, P, P(:,end,:));<br>
> > Q(:,end,:) = [];<br>
> > vol = arrayfun(@(k) det(Q(:,:,k)), 1:np);<br>
> > vol = abs(vol);<br>
> > <br>
> > % generate n samples (i) of v with probability ~ vol<br>
> > vol = vol / sum(vol);<br>
> > c = cumsum(vol); c(end)=1;<br>
> > [~, i] = histc(rand(1,n),[0 c]);<br>
> > <br>
> > % Random rarycentric coordinates<br>
> > % FEX: <a href="http://www.mathworks.com/matlabcentral/fileexchange/9700">http://www.mathworks.com/matlabcentral/fileexchange/9700</a><br>
> > w = randfixedsum(m,n,1,0,1);<br>
> > <br>
> > % Put together<br>
> > V = P(:,:,i);<br>
> > y = zeros(m1, n);<br>
> > for k=1:n<br>
> > y(:,k) = V(:,:,k)*w(:,k);<br>
> > end<br>
> > % final result, m x n array of random vector<br>
> > x = bsxfun(@plus, x0, B*y);<br>
> > <br>
> > % Bruno<br>
> <br>
> 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:<br>
> <br>
> lo = [0 0 0];<br>
> up = [2 1 2];<br>
> s = 2;<br>
> n = 1000;<br>
> <br>
> it gives me the error:<br>
> <br>
> Error using qhullmx<br>
> QH6154 qhull precision error: initial facet 1 is coplanar with the interior point<br>
> ERRONEOUS FACET:<br>
> <br>
> While executing:  qhull d Qt Qbb Qc<br>
> Options selected for Qhull 2010.1 2010/01/14:<br>
> runid 1283078028 delaunay Qtriangulate Qbboundlast Qcoplanarkeep<br>
> _premerge _zerocentrum Qinteriorkeep Pgood _maxwidth 2.7<br>
> Errorroundoff 3.8e15 _onemerge 2.7e14 Visibledistance 7.6e15<br>
> Ucoplanardistance 7.6e15 Widthoutside 1.5e14 _widefacet 4.5e14<br>
> <br>
> The input to qhull appears to be less than 3 dimensional, or a<br>
> computation has overflowed.<br>
> <br>
> Qhull could not construct a clearly convex simplex from points:<br>
> <br>
> The center point is coplanar with a facet, or a vertex is coplanar<br>
> with a neighboring facet. The maximum round off error for<br>
> computing distances is 3.8e15. The center point, facets and distances<br>
> to the center point are as follows:<br>
> <br>
> <br>
> facet p1 p3 p0 distance= 0<br>
> facet p2 p3 p0 distance= 1.1e16<br>
> facet p2 p1 p0 distance= 2.2e16<br>
> facet p2 p1 p3 distance= 1.1e16<br>
> <br>
> These points either have a maximum or minimum xcoordinate, or<br>
> they maximize the determinant for k coordinates. Trial points<br>
> are first selected from points that maximize a coordinate.<br>
> <br>
> The min and max coordinates for each dimension are:<br>
> 0: 0.8392 0.8928 difference= 1.732<br>
> 1: 1.239 1.493 difference= 2.732<br>
> 2: 0 2.732 difference= 2.732<br>
> <br>
> If the input should be full dimensional, you have several options that<br>
> may determine an initial simplex:<br>
>  use 'QJ' to joggle the input and make it full dimensional<br>
>  use 'QbB' to scale the points to the unit cube<br>
>  use 'QR0' to randomly rotate the input for different maximum points<br>
>  use 'Qs' to search all points for the initial simplex<br>
>  use 'En' to specify a maximum roundoff error less than 3.8e15.<br>
>  trace execution with 'T3' to see the determinant for each point.<br>
> <br>
> If the input is lower dimensional:<br>
>  use 'QJ' to joggle the input and make it full dimensional<br>
>  use 'Qbk:0Bk:0' to delete coordinate k from the input. You should<br>
> pick the coordinate with the least range. The hull will have the<br>
> correct topology.<br>
>  determine the flat containing the points, rotate the points<br>
> into a coordinate plane, and delete the other coordinates.<br>
>  add one or more points to make the input full dimensional.<br>
> <br>
> <br>
> This is a Delaunay triangulation and the input is cocircular or cospherical:<br>
>  use 'Qz' to add a point "at infinity" (i.e., above the paraboloid)<br>
>  or use 'QJ' to joggle the input and avoid cocircular data<br>
> <br>
> <br>
> Error in delaunayn (line 101)<br>
> t = qhullmx(x', 'd ', opt);<br>
> <br>
> Error in Untitled (line 32)<br>
> T = delaunayn(V);<br>
> <br>
> I will be obliged to you if you can fixe it.

Tue, 17 Sep 2013 06:28:06 +0000
Re: Generate random numbers with fixed sum and different constraints
http://www.mathworks.com/matlabcentral/newsreader/view_thread/324503#911351
JS Hong
Actually, I found the similar error with original parameters<br>
<br>
lo = [1 0 2 1 0];<br>
up = [5 6 4 4 2]; <br>
<br>
If s = any number between 2 and 6, the same error occurs. <br>
2 and 6 is the max number whose both lower bounds are 0. <br>
<br>
Can you tell me how to handle this error?<br>
Thanks.