function [Sigma,B,X,z,e] = generate_data(dataset,B,n,k,noise_power)

    if(dataset == "all" || ...
       dataset == "aiyou_heavy_tail" || ...
       dataset == "chf_better" || ...
       dataset == "sample1" || ...
       dataset == "aiyou_total")
        d = 1;
    else
        d = k;
    end

    if nargin == 0
        k = 7;     % Dimension of source signals
        n = 20000;  % Number of data points
    end
    
    % B      = eye(k);
    %B = rand(k,k);
    % B      = normc(B);
    %B=eye(k);
    % B      = randn(k,k);
    % B      = B*diag(1./sqrt(diag(B'*B)));
    %B=eye(k);       
    U      = randn(k,k);
    Sigma  = noise_power*(1/k)*(U*U');
    % Sigma  = zeros(k,k);
    % rho   = 0.67;
    % Sigma = rho*(10*eye(k) - B*B');

    % Bernoulli
    p  = 1/2; 
    % z1          = (binornd(1,p,n,k)-p)/sqrt(p*(1-p));
    z1 = (binornd(1,p,n,d)-p);
    z_bernoulli = z1;

    % Uniform 
    z_uniform = unifrnd(-sqrt(3),sqrt(3),n,d);

    % Poisson
    z_poisson = poissrnd(1,n,d);
    % z_poisson = (z_poisson-ones(n,1)*mean(z_poisson))./std(z_poisson);
    z_poisson = (z_poisson-ones(n,d)*mean(z_poisson,"all"));

    % Exponential
    z_exp     = exprnd(1,n,d);
    z_exp     = (z_exp-ones(n,1)*mean(z_exp))./std(z_exp);
    % z_exp     = (z_exp-ones(n,d)*mean(z_exp,"all"));

    % Gamma
    z_gamma     = gamrnd(4,2,n,d);
    % z_gamma     = (z_gamma-ones(n,1)*mean(z_gamma))./std(z_gamma);
    z_gamma     = (z_gamma-ones(n,d)*mean(z_gamma,"all"));

    % lognormal
    z_lognormal  = random(makedist('Lognormal','mu',0,'sigma',1),n,d);
    z_lognormal  = (z_lognormal-ones(n,d)*mean(z_lognormal,"all"));

    % weibull
    z_weibull    = wblrnd(1,1,n,d);
    z_weibull    = (z_weibull-ones(n,d)*mean(z_weibull,"all"));

    % Logistic
    z_logistic = random(makedist('Logistic','mu',0,'sigma',1),n,d);
    z_logistic = (z_logistic-ones(n,d)*mean(z_logistic,"all"));
    % z_logistic = (z_logistic-ones(n,1)*mean(z_logistic))./std(z_logistic);
    % z_logistic = (z_logistic-ones(n,1)*mean(z_logistic))./std(z_logistic);

    % Zero-Kurtosis-bernoulli
    p                           = 1/2 + sqrt(1/12); 
    % z_zero_kurtosis_bernoulli = (binornd(1,p,n,k)-p)/sqrt(p*(1-p));
    z_zero_kurtosis_bernoulli   = binornd(1,p,n,d)-p;
    % z_zero_kurtosis_bernoulli = binornd(1,p,n,k)-p;
    % p                         = 1/2 - sqrt(1/12); 
    % z2                        = (binornd(1,p,n,1)-p)/sqrt(p*(1-p));
    % z_zero_kurtosis_bernoulli = cat(2, z1, z2);

    % Zero-Kurtosis-poisson
    lambda                      = 1/2; 
    multiplier                  = binornd(1,1/2,n,d);
    z_zero_kurtosis_poisson     = multiplier.*(sqrt(poissrnd(lambda,n,d))) + ...
                                  (1-multiplier).*(-1*sqrt(poissrnd(lambda,n,d)));
    z_zero_kurtosis_poisson     = z_zero_kurtosis_poisson - ...
                                  ones(n,d)*mean(z_zero_kurtosis_poisson,"all");

    % Zero-Kurtosis Uniform
    u1 = unifrnd(-1,1,n,d);
    a = sqrt(5 + sqrt(24));
    u2 = unifrnd(-a,a,n,d);
    p = binornd(1,1/2,n,1);R=unifrnd(n,1);
    z_zero_kurtosis_uniform = (R<p).*u1 + (R>=p).*u2;

%     % Cauchy 
%     z_cauchy = cauchyrnd(0,1,n,d);
%     z_cauchy  = (z_cauchy-ones(n,d)*mean(z_cauchy,"all"));

    % Normal
    z_normal = randn(n,d);

    % t-distribution with parameter=3
    z_t_3     = trnd(3, n, d);
    z_t_3     = (z_t_3-ones(n,1)*mean(z_t_3))./std(z_t_3);

    % t-distribution with parameter=5
    z_t_5     = trnd(5, n, d);
    z_t_5     = (z_t_5-ones(n,1)*mean(z_t_5))./std(z_t_5);

    % z_exp + normal   
    z_exp1       = exprnd(1,n,d);
    z_exp1       = (z_exp1-ones(n,1)*mean(z_exp1))./std(z_exp1);
    z_normal1    = randn(n,d);
    z_exp_normal = z_exp1 + z_normal1;

    %z_zero_kurtosis_uniform = 0.5*u1 + 0.5*u2;

    % z  = cat(2, z_zero_kurtosis, z_bernoulli, z_uniform, ...
    %        z_poisson, z_exp, z_logistic, z_gamma);
    
    % z  = cat(2, z_bernoulli, z_zero_kurtosis);
    
    % z_exp  = exprnd(1,n,d/4);
    % z_exp1 = (z_exp-ones(n,1)*mean(z_exp));
    % z_exp  = exprnd(10,n,d/4);
    % z_exp2 = (z_exp-ones(n,1)*mean(z_exp));
    % 
    % z_poisson1 = poissrnd(1,n,d/4);
    % z_poisson1 = (z_poisson1-ones(n,1)*mean(z_poisson1));
    % z_poisson2 = poissrnd(10,n,d/4);
    % z_poisson2 = (z_poisson2-ones(n,1)*mean(z_poisson2));
    
    % z = evrnd(0,.5,n,d);
    % z = (z-ones(n,1)*mean(z));
    % z = cat(2, z_exp1, z_exp2);
    % z = cat(2, z_poisson1, z_poisson2, z_exp1, z_exp2);
   
    % p     = rand(n,1);
    % z_mix = diag(p<.5)*(2*randn(n,d))+diag(p>=.5)*(randn(n,d));
    % z     = z_mix;
    % z     = cat(2, z_logistic, z_zero_kurtosis);

    if(dataset == "bernoulli")
        fprintf("Running on %s data\n", "bernoulli");
        z = z_bernoulli;
    elseif(dataset == "uniform")
         fprintf("Running on %s data\n", "uniform");
         z = z_uniform;
    elseif(dataset == "poisson")
        fprintf("Running on %s data\n", "poisson");
        z = z_poisson; 
    elseif(dataset == "exp")
        fprintf("Running on %s data\n", "exp");
        z = z_exp;
    elseif(dataset == "gamma")
        fprintf("Running on %s data\n", "gamma");
        z = z_gamma;
    elseif(dataset == "lognormal")
        fprintf("Running on %s data\n", "lognormal");
        z = z_lognormal;
    elseif(dataset ==  "weibull")
        fprintf("Running on %s data\n", "weibull");
        z = z_weibull;
    elseif(dataset == "logistic")
        fprintf("Running on %s data\n", "logistic");
        z = z_logistic;
    elseif(dataset == "zero_kurtosis_bernoulli")
        fprintf("Running on %s data\n", "zero_kurtosis_bernoulli");
        z = z_zero_kurtosis_bernoulli;
    elseif(dataset == "zero_kurtosis_poisson")
        fprintf("Running on %s data\n", "zero_kurtosis_poisson");
        z = z_zero_kurtosis_poisson;
    elseif(dataset == "zero_kurtosis_uniform")
        fprintf("Running on %s data\n", "zero_kurtosis_uniform");
        z = z_zero_kurtosis_uniform;
    elseif(dataset == "cauchy")
        fprintf("Running on %s data\n", "z_cauchy");
        z = z_cauchy;
    elseif(dataset == "aiyou_heavy_tail")
        fprintf("Running on %s data\n", "aiyou_heavy_tail");
        z = cat(2,z_cauchy,z_uniform);
    elseif(dataset == "chf_better")
        fprintf("Running on %s data\n", "chf_better");
        z = cat(2,z_zero_kurtosis_bernoulli,z_poisson);
    elseif(dataset == "sample1")
        fprintf("Running on %s data\n", "sample1");
        z = cat(2,z_exp,z_gamma);
    elseif(dataset == "aiyou_total")
        fprintf("Running on %s data\n", "aiyou_total");
        z  = cat(2, z_normal, z_exp, z_t_3, z_lognormal, ...
                 z_t_5, z_logistic, z_weibull, ...
                 z_exp_normal);
    else
        fprintf("Running on %s data\n", "all");
        z  = cat(2, z_bernoulli, z_uniform, z_poisson, z_exp, ...
            z_logistic, z_gamma, z_lognormal, ...
            z_weibull, z_zero_kurtosis_bernoulli, ...
            z_zero_kurtosis_poisson, z_zero_kurtosis_uniform);
    end
    
    % z  = heavy_tail(5,n,d);
    % z  = (z-ones(n,1)*mean(z))./std(z);

    e   = mvnrnd(zeros(k,1), Sigma, n);
    X   = z*B' + e;
end