Попытка определить степень полинома для полиномиальной регрессии

Я пытаюсь предсказать вес новорожденного, используя полиномиальную регрессионную модель. Сначала мне нужно понять, какая степень полинома лучше всего подходит для моих данных. Для этого я разделил свой набор данных на обучающую выборку (70%) и выборку перекрестной проверки (30%) и затем строю график ошибки по степени полинома. Я запускал свой скрипт 4 раза, случайным образом выбирая данные, но каждый раз получал столь различные кривые, как вы можете видеть

Я не знаю, почему это происходит и что я делаю неправильно.

Вы можете увидеть мой код ниже:

главный скрипт

%========== Начало - объявление констант ==========% 
x_training_percent = 0.7;
cv_set_percent = 0.3;
%========== Конец - объявление констант ==========% 

load('data/dataset.m');
[m, n] = size(X);% m: количество примеров, n: количество признаков.

%========== Начало - получение обучающих и CV наборов ==========% 
training_set_size = round(m * x_training_percent);
cv_set_size = round(m * cv_set_percent);
test_set_size = round(m * test_set_percent);

random_order_idx = randperm(m);

indexes = random_order_idx(1:training_set_size);
random_order_idx(1:training_set_size) = []; % Удаление использованных индексов
x_training_o = X(indexes, :); % Оригинальный обучающий набор X
y_training = y(indexes, :);   % Оригинальный обучающий набор y

indexes = random_order_idx(1:cv_set_size);
random_order_idx(1:cv_set_size) = []; % Удаление использованных индексов
x_cv_o = X(indexes, :); % Оригинальный набор перекрестной проверки X
y_cv = y(indexes, :);   % Оригинальный набор перекрестной проверки y

%========== Конец - получение обучающих и CV наборов ==========% 

max_p = 10; % максимальная степень полинома

cv_error = zeros(max_p, 1);
training_error = zeros(max_p, 1);

for p = 1:max_p

  % Обработка обучающей выборки
  x_training = x_training_o;
  x_training = polyFeatures(x_training, p); % Добавление полиномиальных членов от 1 до p
  [x_training, mu, sigma] = featureNormalize(x_training);
  x_training = [ones(training_set_size, 1) x_training];

  % Обработка выборки перекрестной проверки
  x_cv = x_cv_o;
  x_cv = polyFeatures(x_cv, p); % Добавление полиномиальных членов от 1 до p
  x_cv = bsxfun(@minus, x_cv, mu);
  x_cv = bsxfun(@rdivide, x_cv, sigma);
  x_cv = [ones(cv_set_size, 1) x_cv];

  %========== Начало - Обучение ==========% 
  lambda = 0;
  theta = trainLinearReg(x_training, y_training, lambda);
  %========== Конец - Обучение ==========% 

  %========== Начало - Вычисление ошибок предсказания с полиномиальной степенью p ==========% 
    predictions = x_training * theta; % Предсказания с обучающей выборкой
    training_error(p, :) = (1 / (2 * training_set_size)) * sum((predictions - y_training) .^ 2);

    cv_predictions = x_cv * theta; % Предсказания с выборкой перекрестной проверки
    cv_error(p, :) = (1 / (2 * cv_set_size)) * sum((cv_predictions - y_cv) .^ 2);
  %========== Конец - Вычисление ошибок предсказания ==========% 

end

plot(1:p, training_error, 1:p, cv_error);

title(sprintf('Кривая обучения для линейной регрессии с lambda = %f', lambda));
legend('Обучение', 'Перекрестная проверка')
xlabel('степень полинома')
ylabel('Ошибка')

polyFeatures

function [X_poly] = polyFeatures(X, p)
%POLYFEATURES Отображает X (вектор 1D) в p-ю степень
%   [X_poly] = POLYFEATURES(X, p) принимает матрицу данных X (размер m x 1) и
%   отображает каждый пример в его полиномиальные признаки, где
%   X_poly(i, :) = [X(i) X(i).^2 X(i).^3 ...  X(i).^p];
%

X_poly = X; % Для p = 1

for i = 2:p
    X_poly = [X_poly (X .^ i)];
end

end

featureNormalize

function [X_norm, mu, sigma] = featureNormalize(X)
%FEATURENORMALIZE Нормализует признаки в X 
%   FEATURENORMALIZE(X) возвращает нормализованную версию X, где
%   среднее значение каждого признака равно 0, а стандартное отклонение
%   равно 1. Это часто является хорошим этапом предварительной обработки при
%   работе с алгоритмами обучения.

mu = mean(X);
X_norm = bsxfun(@minus, X, mu);

sigma = std(X_norm);
X_norm = bsxfun(@rdivide, X_norm, sigma);

end

trainingLinearReg

function [theta] = trainLinearReg(X, y, lambda)
%TRAINLINEARREG Обучает линейную регрессию, учитывая набор данных (X, y) и
%параметр регуляризации lambda
%   [theta] = TRAINLINEARREG (X, y, lambda) обучает линейную регрессию с использованием
%   набора данных (X, y) и параметра регуляризации lambda. Возвращает
%   обученные параметры theta.
%

% Инициализация Theta
initial_theta = zeros(size(X, 2), 1); 

% Создание "сокращенной записи" для функции стоимости, которую нужно минимизировать
costFunction = @(t) linearRegCostFunction(X, y, t, lambda);

% Теперь costFunction - это функция, которая принимает только один аргумент
options = optimset('MaxIter', 400, 'GradObj', 'on');

% Минимизация с использованием fmincg
theta = fmincg(costFunction, initial_theta, options);

end

linearRegCostFunction

function [J, grad] = linearRegCostFunction(X, y, theta, lambda)
%LINEARREGCOSTFUNCTION Вычисляет стоимость и градиент для регуляризированной линейной 
%регрессии с несколькими переменными
%   [J, grad] = LINEARREGCOSTFUNCTION(X, y, theta, lambda) вычисляет 
%   стоимость использования theta в качестве параметра для линейной регрессии для 
%   подгонки точек данных в X и y. Возвращает стоимость в J и градиент в grad

% Инициализация некоторых полезных значений
m = length(y); % количество обучающих примеров

% Необходимо правильно вернуть следующие переменные 
J = 0;
grad = zeros(size(theta));

hx = X * theta; % Предсказание

J = (1 / (2 * m)) * sum((hx - y) .^ 2) + (lambda / (2 * m)) * sum([ 0; theta(2:end, :) ] .^ 2);

grad = (1 / m) * (hx - y)' * X + (lambda / m) * [ 0; theta(2:end, :) ]';

grad = grad(:);

end

Может кто-нибудь помочь мне?

Вы, похоже, реализуете многие части этого процесса. Если вы переключитесь на установленные пакеты, вы сможете получить надежные результаты. Затем сравните надежные результаты с вашими, чтобы увидеть, есть ли у вас ошибка в реализации.