Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Generate random numbers with fixed sum and different constraints

Subject: Generate random numbers with fixed sum and different constraints

From: Dmitrey Yershov

Date: 14 Nov, 2012 14:19:17

Message: 1 of 26

Hello. I need to generate non-negative 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

http://www.mathworks.com/matlabcentral/fileexchange/9700

but in this alghorithm a<=xi<=b (a and b are the same for all xi). Any ideas?

Subject: Generate random numbers with fixed sum and different constraints

From: Greg Heath

Date: 17 Nov, 2012 11:22:09

Message: 2 of 26

"Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message <k80995$9eo$1@newscl01ah.mathworks.com>...
> Hello. I need to generate non-negative 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
>
> http://www.mathworks.com/matlabcentral/fileexchange/9700
>
> but in this alghorithm a<=xi<=b (a and b are the same for all xi). Any ideas?

Z = a + (b-a)*rand(m,n);

sumZ = repmat(sum(Z),m,1);

I'll let you figure out the rest.

Hope this helps

Greg

Subject: Generate random numbers with fixed sum and different constraints

From: Roger Stafford

Date: 17 Nov, 2012 18:32:16

Message: 3 of 26

"Greg Heath" <heath@alumni.brown.edu> wrote in message <k87s11$h6u$1@newscl01ah.mathworks.com>...
> "Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message <k80995$9eo$1@newscl01ah.mathworks.com>...
> > Hello. I need to generate non-negative 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
> >
> > http://www.mathworks.com/matlabcentral/fileexchange/9700
> >
> > but in this alghorithm a<=xi<=b (a and b are the same for all xi). Any ideas?
>
> Z = a + (b-a)*rand(m,n);
>
> sumZ = repmat(sum(Z),m,1);
>
> I'll let you figure out the rest.
>
> Hope this helps
>
> Greg
- - - - - - - - - -
  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.

  The space of points in this case having a sum of 1 within the permitted three-dimensional 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 area-wise uniform manner throughout the hexagon.

  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 n-dimensional rectangular box other than an n-dimensional 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.

Roger Stafford

Subject: Generate random numbers with fixed sum and different constraints

From: Bruno Luong

Date: 17 Nov, 2012 19:59:12

Message: 4 of 26

"Roger Stafford" wrote in message <k88l7g$7tm$1@newscl01ah.mathworks.com>...
> "Greg Heath" <heath@alumni.brown.edu> wrote in message <k87s11$h6u$1@newscl01ah.mathworks.com>...

>
> 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 n-dimensional rectangular box other than an n-dimensional 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.

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.

Bruno

Subject: Generate random numbers with fixed sum and different constraints

From: Roger Stafford

Date: 18 Nov, 2012 01:44:16

Message: 5 of 26

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k88qag$nm1$1@newscl01ah.mathworks.com>...
> 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.
>
> Bruno
- - - - - - - - - -
  Yes, that is a conceivable approach, Bruno. However it faces some formidable difficulties with large dimensionality. For an n-dimensional cube with n equal to, say, 51 the number of vertices in the n-1 dimensional polytope with a fixed sum set at or near the half-way 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 n-1 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
unequal bounding intervals.

Roger Stafford

Subject: Generate random numbers with fixed sum and different constraints

From: Greg Heath

Date: 18 Nov, 2012 02:05:23

Message: 6 of 26

"Roger Stafford" wrote in message <k88l7g$7tm$1@newscl01ah.mathworks.com>...
> "Greg Heath" <heath@alumni.brown.edu> wrote in message <k87s11$h6u$1@newscl01ah.mathworks.com>...
> > "Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message <k80995$9eo$1@newscl01ah.mathworks.com>...
> > > Hello. I need to generate non-negative 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
> > >
> > > http://www.mathworks.com/matlabcentral/fileexchange/9700
> > >
> > > but in this alghorithm a<=xi<=b (a and b are the same for all xi). Any ideas?
> >
> > Z = a + (b-a)*rand(m,n);
> >
> > sumZ = repmat(sum(Z),m,1);
> >
> > I'll let you figure out the rest.
> >
> > Hope this helps
> >
> > Greg
> - - - - - - - - - -
> 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.

WHOOPS! Silly mistake.

Please pardon my youthful exuberance.

Greg

P.S. The OP did not indicate that the number of points is specified. Was that an unintended omission?

Subject: Generate random numbers with fixed sum and different constraints

From: Dmitrey Yershov

Date: 18 Nov, 2012 17:04:09

Message: 7 of 26

"Roger Stafford" wrote in message <k89ehg$rm6$1@newscl01ah.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k88qag$nm1$1@newscl01ah.mathworks.com>...
> > 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.
> >
> > Bruno
> - - - - - - - - - -
> Yes, that is a conceivable approach, Bruno. However it faces some formidable difficulties with large dimensionality. For an n-dimensional cube with n equal to, say, 51 the number of vertices in the n-1 dimensional polytope with a fixed sum set at or near the half-way 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 n-1 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
> unequal bounding intervals.
>
> Roger Stafford

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).

Subject: Generate random numbers with fixed sum and different constraints

From: Bruno Luong

Date: 18 Nov, 2012 18:46:10

Message: 8 of 26

I have not fully tested, but it should go like this:

% lower + upper bounds
lo = [1 0 2 1 0];
up = [5 6 4 4 2];
% target sum
s = 15;
% number of samples
n = 1000;

m = length(lo);
slo = sum(lo);
sup = sum(up);

x0 = interp1([slo; sup], [lo; up], s);
if any(isnan(x0))
    error('no solution exists')
end

x0 = x0(:);
B = null(ones(size(lo)));
A = [-B; B];
b = [x0-lo(:); up(:)-x0];

% FEX: http://www.mathworks.com/matlabcentral/fileexchange/30892
% Thank you Matt
[V,nr,nre] = lcon2vert(A,b,[],[],tol);

T = delaunayn(V);
P = V(T,:);
P = reshape(P, [size(T) m-1]);
P = permute(P, [3 2 1]);
np = size(P,3);

% http://www.mathworks.com/matlabcentral/fileexchange/9700
% % Thank you Roger
Q = bsxfun(@minus, P, P(:,end,:));
Q(:,end,:) = [];
vol = arrayfun(@(k) det(Q(:,:,k)), 1:np);
vol = abs(vol);
w = randfixedsum(m,n,1,0,1);

% generate n samples of v with probability vol
vol = vol / sum(vol);
c = cumsum(vol); c(end)=1;
[~, i] = histc(rand(1,n),[0 c]);

V = P(:,:,i);
y = zeros(m-1, n);
for k=1:n
    y(:,k) = V(:,:,k)*w(:,k);
end
x = bsxfun(@plus, x0, B*y);

% Bruno

Subject: Generate random numbers with fixed sum and different constraints

From: Bruno Luong

Date: 19 Nov, 2012 06:39:13

Message: 9 of 26

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

Subject: Generate random numbers with fixed sum and different constraints

From: David Epstein

Date: 19 Nov, 2012 10:26:09

Message: 10 of 26

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8badi$8rd$1@newscl01ah.mathworks.com>...
> I have not fully tested, but it should go like this:
>
> % lower + upper bounds
> lo = [1 0 2 1 0];
> up = [5 6 4 4 2];
> % target sum
> s = 15;
> % number of samples
> n = 1000;
>
> m = length(lo);
> slo = sum(lo);
> sup = sum(up);
>
> x0 = interp1([slo; sup], [lo; up], s);
> if any(isnan(x0))
> error('no solution exists')
> end
>
> x0 = x0(:);
> B = null(ones(size(lo)));
> A = [-B; B];
> b = [x0-lo(:); up(:)-x0];
>
> % FEX: http://www.mathworks.com/matlabcentral/fileexchange/30892
> % Thank you Matt
> [V,nr,nre] = lcon2vert(A,b,[],[],tol);
>
> T = delaunayn(V);
> P = V(T,:);
> P = reshape(P, [size(T) m-1]);
> P = permute(P, [3 2 1]);
> np = size(P,3);
>
> % http://www.mathworks.com/matlabcentral/fileexchange/9700
> % % Thank you Roger
> Q = bsxfun(@minus, P, P(:,end,:));
> Q(:,end,:) = [];
> vol = arrayfun(@(k) det(Q(:,:,k)), 1:np);
> vol = abs(vol);
> w = randfixedsum(m,n,1,0,1);
>
> % generate n samples of v with probability vol
> vol = vol / sum(vol);
> c = cumsum(vol); c(end)=1;
> [~, i] = histc(rand(1,n),[0 c]);
>
> V = P(:,:,i);
> y = zeros(m-1, n);
> for k=1:n
> y(:,k) = V(:,:,k)*w(:,k);
> end
> x = bsxfun(@plus, x0, B*y);
>
> % Bruno

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:

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.
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.
3. Use Gram-Schmidt to construct linear isometry g that sends U into an n-dimensional euclidean space Z and f(V) into the subspace z_1=0. The input data for Gram-Schmidt has its first vector v a unit vector orthogonal to f(V), followed by the standard basis, orthogonal to the various (n-1)-dimensional faces of U. (Perhaps omit the basis vector nearest to v.)
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)).
5. Generate a uniform random sample in S and keep only those satisfying the appropriate 2*n inequalities.

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.

@Roger: what were your reasons for rejecting this approach in your randfixedsum package on FEX?

David

Subject: Generate random numbers with fixed sum and different constraints

From: Dmitrey Yershov

Date: 19 Nov, 2012 10:58:07

Message: 11 of 26

"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.

Subject: Generate random numbers with fixed sum and different constraints

From: Bruno Luong

Date: 19 Nov, 2012 12:04:13

Message: 12 of 26

"Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message
> 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.

Please try to replace the delaunayn() command with:

if m==3
    T = delaunay(V);
else
    T = delaunayn(V);
end

Obviously it doesn' work for 2-dimensional triangulation.

Bruno

Subject: Generate random numbers with fixed sum and different constraints

From: Dmitrey Yershov

Date: 19 Nov, 2012 12:29:08

Message: 13 of 26

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8d77s$lfd$1@newscl01ah.mathworks.com>...
> "Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message
> > 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.
>
> Please try to replace the delaunayn() command with:
>
> if m==3
> T = delaunay(V);
> else
> T = delaunayn(V);
> end
>
> Obviously it doesn' work for 2-dimensional triangulation.
>
> Bruno

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)!

P.S.: Special thanks to Roger Stafford and Matt J for basis procedures.

Subject: Generate random numbers with fixed sum and different constraints

From: Bruno Luong

Date: 19 Nov, 2012 12:39:13

Message: 14 of 26

"Dmitrey Yershov" <pierrevanstulov@mail.ru> wrote in message
>
> 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)!

You could further optimize the code by replacing the for-loop of V*w with James's MTIMESX

http://www.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-with-multi-dimensional-support

Bruno

Subject: Generate random numbers with fixed sum and different constraints

From: Bruno Luong

Date: 19 Nov, 2012 13:02:11

Message: 15 of 26

It seems to me this task
w = randfixedsum(m,n,1,0,1)

can be optimized as well. But let save this optimization topic for another day.

Bruno

Subject: Generate random numbers with fixed sum and different constraints

From: Roger Stafford

Date: 19 Nov, 2012 23:10:18

Message: 16 of 26

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8dakj$39u$1@newscl01ah.mathworks.com>...
> It seems to me this task
> w = randfixedsum(m,n,1,0,1)
> can be optimized as well. .....
- - - - - - - - -
  Yes, I agree, Bruno; 'randfixedsum' is rather cumbersome to use just for a single simplex. How about this instead:

 R = bsxfun(@power,rand(m-1,n),1./(m-1:-1:1)');
 X = cumprod([ones(1,n);R]).*[ones(m,n)-[R;zeros(1,n)]];

Roger Stafford

Subject: Generate random numbers with fixed sum and different constraints

From: Roger Stafford

Date: 20 Nov, 2012 01:54:09

Message: 17 of 26

"David Epstein" <David.Epstein.spam@remove.warwick.ac.uk> wrote in message <k8d1g1$4gr$1@newscl01ah.mathworks.com>...
> @Roger: what were your reasons for rejecting this approach in your randfixedsum package on FEX?
- - - - - - - - - -
  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.

  To get a feeling for this, consider an n-dimensional hypersphere of radius 1/2 enclosed in an n-dimensional cube with unit-length sides. The n-dimensional volume of the cube is 1 whereas that of the hypersphere for even n is (pi/4)^(n/2)/(n/2)! (See http://en.wikipedia.org/wiki/N-sphere.) For n = 50 this would be 1.53E-28, a small acceptance rate indeed.

Roger Stafford

Subject: Generate random numbers with fixed sum and different constraints

From: Bruno Luong

Date: 20 Nov, 2012 06:53:08

Message: 18 of 26

"Roger Stafford" wrote in message <k8ee8q$k6f$1@newscl01ah.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8dakj$39u$1@newscl01ah.mathworks.com>...
> > It seems to me this task
> > w = randfixedsum(m,n,1,0,1)
> > can be optimized as well. .....
> - - - - - - - - -
> Yes, I agree, Bruno; 'randfixedsum' is rather cumbersome to use just for a single simplex. How about this instead:
>
> R = bsxfun(@power,rand(m-1,n),1./(m-1:-1:1)');
> X = cumprod([ones(1,n);R]).*[ones(m,n)-[R;zeros(1,n)]];
>
> Roger Stafford

Really cool Roger! I would love if you could share the justification of this formula.

Bruno

Subject: Generate random numbers with fixed sum and different constraints

From: Bruno Luong

Date: 20 Nov, 2012 07:02:09

Message: 19 of 26

Just a small typo correction in a comment of my code. It should read:

% feasible set: S = { x=x0+B*y : A*y <= b }

% Bruno

Subject: Generate random numbers with fixed sum and different constraints

From: Bruno Luong

Date: 20 Nov, 2012 07:17:08

Message: 20 of 26

"Roger Stafford" wrote in message <k8ens1$kfr$1@newscl01ah.mathworks.com>...
> "David Epstein" <David.Epstein.spam@remove.warwick.ac.uk> wrote in message <k8d1g1$4gr$1@newscl01ah.mathworks.com>...
> > @Roger: what were your reasons for rejecting this approach in your randfixedsum package on FEX?
> - - - - - - - - - -
> 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.
>
> To get a feeling for this, consider an n-dimensional hypersphere of radius 1/2 enclosed in an n-dimensional cube with unit-length sides. The n-dimensional volume of the cube is 1 whereas that of the hypersphere for even n is (pi/4)^(n/2)/(n/2)! (See http://en.wikipedia.org/wiki/N-sphere.) For n = 50 this would be 1.53E-28, a small acceptance rate indeed.

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.

I have little experience in Delaunay in high dimensional space. But I could imagine it 's challenging too.

Bruno

Subject: Generate random numbers with fixed sum and different constraints

From: james bejon

Date: 20 Nov, 2012 08:16:14

Message: 21 of 26

"Roger Stafford" wrote in message <k8ee8q$k6f$1@newscl01ah.mathworks.com>...
> R = bsxfun(@power,rand(m-1,n),1./(m-1:-1:1)');
> X = cumprod([ones(1,n);R]).*[ones(m,n)-[R;zeros(1,n)]];

Very neat! (Though I have no idea how it works).

Subject: Generate random numbers with fixed sum and different constraints

From: Roger Stafford

Date: 20 Nov, 2012 19:54:09

Message: 22 of 26

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8f9ck$gng$1@newscl01ah.mathworks.com>...
> Really cool Roger! I would love if you could share the justification of this formula.

"james bejon" wrote in message <k8fe8e$2lo$1@newscl01ah.mathworks.com>...
> Very neat! (Though I have no idea how it works).
- - - - - - - -
  To Bruno, James, and anyone else interested:

  I recommended the two lines of code

 R = bsxfun(@power,rand(m-1,n),1./(m-1:-1:1)');
 X = cumprod([ones(1,n);R]).*[ones(m,n)-[R;zeros(1,n)]];

as a substitute for Bruno's

 w = randfixedsum(m,n,1,0,1)

(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 m-1 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:

 [1;0;0;0;...], [0;1;0;0;...], [0;0;1;0;...], ..., [0;0;0;...;0;0;1]

which we can call P1, P2, ..., Pm.

  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

 r1 = rand^(1/3); r2 = rand^(1/2); r3 = rand^(1);
 X = [ (1)*(1-r1) ; (1*r1)*(1-r2) ; (1*r1*r2)*(1-r3) ; (1*r1*r2*r3)*(1-0) ];

which with the definitions of the points is the same as

 X = (1-r1)*P1 + r1*((1-r2)*P2 + r2*((1-r3)*P3 + r3*P4));

  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

 (1-r3)*P3 + r3*P4

would distribute uniformly along the line between P3 and P4 using r3 = rand. The quantity

 (1-r2)*P2 + r2*((1-r3)*P3 + r3*P4)

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.

Roger Stafford

Subject: Generate random numbers with fixed sum and different constraints

From: Bruno Luong

Date: 21 Nov, 2012 07:38:17

Message: 23 of 26

Thanks Roger for the explanation.

Bruno

Subject: Generate random numbers with fixed sum and different constraints

From: james bejon

Date: 21 Nov, 2012 21:14:09

Message: 24 of 26

Very nice. I hadn't thought of it in those terms.

Subject: Generate random numbers with fixed sum and different constraints

From: JS Hong

Date: 17 Sep, 2013 05:52:09

Message: 25 of 26

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.

Subject: Generate random numbers with fixed sum and different constraints

From: JS Hong

Date: 17 Sep, 2013 06:28:06

Message: 26 of 26

Actually, I found the similar error with original parameters

lo = [1 0 2 1 0];
up = [5 6 4 4 2];

If s = any number between 2 and 6, the same error occurs.
2 and 6 is the max number whose both lower bounds are 0.

Can you tell me how to handle this error?
Thanks.

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us