
function [X,y,w] = reg_least_sq_data(n,model,sigma_x,sigma_y)

% random inputs
X = sigma_x * randn(n,model.m);

% linear prediction
y = X * model.b;

% plus quadratic terms
for i=1:n
  xi = X(i,:)';
  y(i) = y(i) + xi' * model.A * xi;
end

% plus Gaussian measurement noise
y = y + sigma_y * randn(n,1);

