inital commit

.gitignore

@@ -0,0 +1,2 @@
+*.m~
+fig/

randfixedsum.m

@@ -0,0 +1,97 @@
+function [x,v] = randfixedsum(n,m,s,a,b)
+
+% [x,v] = randfixedsum(n,m,s,a,b)
+%
+%   This generates an n by m array x, each of whose m columns
+% contains n random values lying in the interval [a,b], but
+% subject to the condition that their sum be equal to s.  The
+% scalar value s must accordingly satisfy n*a <= s <= n*b.  The
+% distribution of values is uniform in the sense that it has the
+% conditional probability distribution of a uniform distribution
+% over the whole n-cube, given that the sum of the x's is s.
+%
+%   The scalar v, if requested, returns with the total
+% n-1 dimensional volume (content) of the subset satisfying
+% this condition.  Consequently if v, considered as a function
+% of s and divided by sqrt(n), is integrated with respect to s
+% from s = a to s = b, the result would necessarily be the
+% n-dimensional volume of the whole cube, namely (b-a)^n.
+%
+%   This algorithm does no "rejecting" on the sets of x's it
+% obtains.  It is designed to generate only those that satisfy all
+% the above conditions and to do so with a uniform distribution.
+% It accomplishes this by decomposing the space of all possible x
+% sets (columns) into n-1 dimensional simplexes.  (Line segments,
+% triangles, and tetrahedra, are one-, two-, and three-dimensional
+% examples of simplexes, respectively.)  It makes use of three
+% different sets of 'rand' variables, one to locate values
+% uniformly within each type of simplex, another to randomly
+% select representatives of each different type of simplex in
+% proportion to their volume, and a third to perform random
+% permutations to provide an even distribution of simplex choices
+% among like types.  For example, with n equal to 3 and s set at,
+% say, 40% of the way from a towards b, there will be 2 different
+% types of simplex, in this case triangles, each with its own
+% area, and 6 different versions of each from permutations, for
+% a total of 12 triangles, and these all fit together to form a
+% particular planar non-regular hexagon in 3 dimensions, with v
+% returned set equal to the hexagon's area.
+%
+% Roger Stafford - Jan. 19, 2006
+
+% Check the arguments.
+if (m~=round(m))|(n~=round(n))|(m<0)|(n<1)
+ error('n must be a whole number and m a non-negative integer.')
+elseif (s<n*a)|(s>n*b)|(a>=b)
+ error('Inequalities n*a <= s <= n*b and a < b must hold.')
+end
+
+% Rescale to a unit cube: 0 <= x(i) <= 1
+s = (s-n*a)/(b-a);
+
+% Construct the transition probability table, t.
+% t(i,j) will be utilized only in the region where j <= i + 1.
+k = max(min(floor(s),n-1),0); % Must have 0 <= k <= n-1
+s = max(min(s,k+1),k); % Must have k <= s <= k+1
+s1 = s - [k:-1:k-n+1]; % s1 & s2 will never be negative
+s2 = [k+n:-1:k+1] - s;
+w = zeros(n,n+1); w(1,2) = realmax; % Scale for full 'double' range
+t = zeros(n-1,n);
+tiny = 2^(-1074); % The smallest positive matlab 'double' no.
+for i = 2:n
+ tmp1 = w(i-1,2:i+1).*s1(1:i)/i;
+ tmp2 = w(i-1,1:i).*s2(n-i+1:n)/i;
+ w(i,2:i+1) = tmp1 + tmp2;
+ tmp3 = w(i,2:i+1) + tiny; % In case tmp1 & tmp2 are both 0,
+ tmp4 = (s2(n-i+1:n) > s1(1:i)); % then t is 0 on left & 1 on right
+ t(i-1,1:i) = (tmp2./tmp3).*tmp4 + (1-tmp1./tmp3).*(~tmp4);
+end
+
+% Derive the polytope volume v from the appropriate
+% element in the bottom row of w.
+v = n^(3/2)*(w(n,k+2)/realmax)*(b-a)^(n-1);
+
+% Now compute the matrix x.
+x = zeros(n,m);
+if m == 0, return, end % If m is zero, quit with x = []
+rt = rand(n-1,m); % For random selection of simplex type
+rs = rand(n-1,m); % For random location within a simplex
+s = repmat(s,1,m);
+j = repmat(k+1,1,m); % For indexing in the t table
+sm = zeros(1,m); pr = ones(1,m); % Start with sum zero & product 1
+for i = n-1:-1:1  % Work backwards in the t table
+ e = (rt(n-i,:)<=t(i,j)); % Use rt to choose a transition
+ sx = rs(n-i,:).^(1/i); % Use rs to compute next simplex coord.
+ sm = sm + (1-sx).*pr.*s/(i+1); % Update sum
+ pr = sx.*pr; % Update product
+ x(n-i,:) = sm + pr.*e; % Calculate x using simplex coords.
+ s = s - e; j = j - e; % Transition adjustment
+end
+x(n,:) = sm + pr.*s; % Compute the last x
+
+% Randomly permute the order in the columns of x and rescale.
+rp = rand(n,m); % Use rp to carry out a matrix 'randperm'
+[ig,p] = sort(rp); % The values placed in ig are ignored
+x = (b-a)*x(p+repmat([0:n:n*(m-1)],n,1))+a; % Permute & rescale x
+
+return

randomhermmat.m

@@ -0,0 +1,4 @@
+C = rand(3) + i*rand(3);
+H = (C + ctranspose(C))/2
+[V D]=eig(H)
+

stochmateig.m

@@ -0,0 +1,37 @@
+% http://www.mathworks.com/matlabcentral/fileexchange/9700-random-vectors-with-fixed-sum
+%function [] = stochmateig(N,dim,colsum,elmin,elmax)
+
+%% Initialize parameters (uncomment to run as script)
+
+N=100000;
+dim=3;
+colsum = 1;
+elmin=0;
+elmax=1;
+
+%%
+ss=zeros(dim,dim,N);
+eV=zeros(dim,N);
+
+for i=1:N
+    ss(:,:,i)=randfixedsum(dim,dim,colsum,elmin,elmax);
+    [V D]=eig(ss(:,:,i));
+    eV(:,i) = V(:,1);
+end
+
+% Settings
+%-------------------------------
+set(0,'defaultaxesfontsize',20);
+set(0,'defaulttextfontsize',20);
+scrsz = get(0,'ScreenSize');
+%-------------------------------
+
+% Plot data
+%-------------------------------
+%if opensavefigFlag
+%    h=figure('Visible','off','Position',[0 0 scrsz(3)/4.5 scrsz(4)/1.4]); set(h,'Color','w');
+%else
+    h=figure('Position',[0 0 scrsz(3)/4.5 scrsz(4)/1.4]); set(h,'Color','w');
+
+scatter3(eV(1,:)',eV(2,:)',eV(3,:)',5,'r');
+%plot3(eV(1,:)',eV(2,:)',eV(3,:)');