function [x,y] = data(n, a, b, u)

% Generate feature points out to twice the standard deviation from the mean
x = sort(a + randn(n, 1) * 2 * b);

% Exact function values
y = f(x, a, b);

% Add zero-mean noise
y = y + rand(size(y)) .* u - u / 2;