clc
clear
x = [0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30];
y = [5.13, 3.20, 2.46, 1.74, 1.18, 1.03, 0.574, 0.424, 0.293, 0.197, 0.155];
linear_coeffs = polyfit(x, y, 1);
y_linear = polyval(linear_coeffs, x);
x_nonzero = x(2:end); 
y_nonzero = y(2:end); 
log_x = log(x_nonzero);
log_y = log(y_nonzero);
power_coeffs = polyfit(log_x, log_y, 1);
a_power = exp(power_coeffs(2));
b_power = power_coeffs(1);
y_power = a_power * x.^b_power;
log_y_full = log(y);
exp_coeffs = polyfit(x, log_y_full, 1);
a_exp = exp(exp_coeffs(2));
b_exp = exp_coeffs(1);
y_exp = a_exp * exp(b_exp * x);
residuals_linear = sum((y - y_linear).^2);
residuals_power = sum((y_nonzero - a_power * x_nonzero.^b_power).^2);
residuals_exp = sum((y - a_exp * exp(b_exp * x)).^2);

[~, best_model] = min([residuals_linear, residuals_power, residuals_exp]);

x_pred = 16.7;
switch best_model
    case 1 
        y_pred = polyval(linear_coeffs, x_pred);
        best_fit_y = y_linear;
        func_desc = sprintf('Linear: y = %.3f*x + %.3f', linear_coeffs(1), linear_coeffs(2));
    case 2 
        y_pred = a_power * x_pred^b_power;
        best_fit_y = y_power;
        func_desc = sprintf('Power-law: y = %.3f*x^%.3f', a_power, b_power);
    case 3 
        y_pred = a_exp * exp(b_exp * x_pred);
        best_fit_y = y_exp;
        func_desc = sprintf('Exponential: y = %.3f*exp(%.3f*x)', a_exp, b_exp);
end
figure;
hold on;
plot(x, y, 'o', 'DisplayName', 'Original Data');
plot(x, best_fit_y, '-', 'DisplayName', 'Best Fit');
plot(x_pred, y_pred, 'rx', 'DisplayName', sprintf('y at x=%.1f', x_pred));
xlabel('x');
ylabel('y');
title('Data and Best-Fitting Function by kaliyah johnson');
legend('Location', 'northeast');
grid on;
disp(['Best-fitting model: ', func_desc]);
disp(['Predicted y at x = ', num2str(x_pred), ' is y = ', num2str(y_pred)]);