|
|
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 = [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.
|
|
