function [Sigma,B,X,z,e] = generate_data_belkin(n,k,noise_power,B)
    if nargin == 0
        k = 12;     % Dimension of source signals
        n = 10000;  % Number of data points
        noise_power = 0.67;
    end

    addpath(genpath('../util'));
    
    % U      = gallery('qmult',k);
    % Lambda = diag(rand(k,1)*2 + 1);
    % V      = gallery('qmult',k);
    % B      = U*Lambda*V'; % Mixing matrix
    % B      = eye(k);
           
    rho   = noise_power;
    % Sigma = rho*(10*eye(k) - B*B');
    V     = randn(k,k);
    Sigma = rho*(1/k)*(V*V');
    
    % Bernoulli
    p           = 1/2; 
    z1          = (binornd(1,p,n,3)-p)/sqrt(p*(1-p));
    z_bernoulli = z1;
    
    % Uniform 
    z_uniform = unifrnd(-sqrt(3),sqrt(3),n,3);
    
    % Laplace
    z_laplace = laprnd(n, 1, 0, 0.05);
    z_laplace = (z_laplace-ones(n,1)*mean(z_laplace))./std(z_laplace);

    % Exponential
    z_exp     = exprnd(5,n,3);
    z_exp     = (z_exp-ones(n,1)*mean(z_exp))./std(z_exp);

    % t-distribution (3 DOF)
    z_t_3     = trnd(3,n,1);
    z_t_3     = (z_t_3-ones(n,1)*mean(z_t_3))./std(z_t_3);
    
    % t-distribution (5 DOF)
    z_t_5     = trnd(5,n,1);
    z_t_5     = (z_t_5-ones(n,1)*mean(z_t_5))./std(z_t_5);

    % zero-kurtosis bernoulli
    p                           = 1/2 + sqrt(1/12); 
    z_zero_kurtosis_bernoulli   = binornd(1,p,n,3)-p;
    z_zero_kurtosis_bernoulli   = z_zero_kurtosis_bernoulli/sqrt(p*(1-p));
    
    % z  = cat(2, z_bernoulli, z_uniform, z_laplace, ...
    %             z_exp, z_t_3, z_t_5, z_zero_kurtosis_bernoulli);

    z  = cat(2, z_uniform, z_zero_kurtosis_bernoulli, z_exp);
    e  = mvnrnd(zeros(k,1), Sigma, n);
    X  = z*B' + e;
end