function [X] = generate_feature_matrix(n, p, Sigma, type, param)
   
   % inputs:
   % n: sample size
   % p: feature size
   % Sigma: feature covariance
   % type: type of distribution for entries of Z
   % param: parameter for distribution
   % bern: probability of success
   % t: degrees of freedom
   % output:
   % X = \Sigma^{1/2} * Z
   % scaling: each Z entry has mean 0 and variance 1
   
   if strcmp(type, 'gau')
       Z = randn(n, p);
   end
   
   if strcmp(type, 'bern')
       Z = (binornd(1, param, [n p]) - param) ...
           / sqrt(param * (1 - param));
   end
   
   if strcmp(type, 't')
       Z = trnd(param, [n p]) / sqrt(param / (param - 2));
   end
   X = Z * sqrtm(Sigma);
end