| 1 | %X is generated following a gaussian mixture \sum pi_r N(meanX_k, covX_r) |
| 2 | %Y is generated then, with Y_i ~ \sum pi_r N(Beta_r.X_i, covY_r) |
| 3 | function[X,Y,Z] = generateIO(meanX, covX, covY, pi, beta, n) |
| 4 | |
| 5 | [p, ~, k] = size(covX); |
| 6 | [m, ~, ~] = size(covY); |
| 7 | if exist('octave_config_info') |
| 8 | %Octave statistics package doesn't have gmdistribution() |
| 9 | X = zeros(n, p); |
| 10 | Z = zeros(n); |
| 11 | cs = cumsum(pi); |
| 12 | for i=1:n |
| 13 | %TODO: vectorize ? http://stackoverflow.com/questions/2977497/weighted-random-numbers-in-matlab |
| 14 | tmpRand01 = rand(); |
| 15 | [~,Z(i)] = min(cs - tmpRand01 >= 0); |
| 16 | X(i,:) = mvnrnd(meanX(Z(i),:), covX(:,:,Z(i)), 1); |
| 17 | end |
| 18 | else |
| 19 | gmDistX = gmdistribution(meanX, covX, pi); |
| 20 | [X, Z] = random(gmDistX, n); |
| 21 | end |
| 22 | |
| 23 | Y = zeros(n, m); |
| 24 | BX = zeros(n,m,k); |
| 25 | for i=1:n |
| 26 | for r=1:k |
| 27 | %compute beta_r . X_i |
| 28 | BXir = zeros(1, m); |
| 29 | for mm=1:m |
| 30 | BXir(mm) = dot(X(i,:), beta(:,mm,r)); |
| 31 | end |
| 32 | %add pi(r) * N(beta_r . X_i, covY) to Y_i |
| 33 | Y(i,:) = Y(i,:) + pi(r) * mvnrnd(BXir, covY(:,:,r), 1); |
| 34 | end |
| 35 | end |
| 36 | |
| 37 | end |