IGPR.m

function [y_pred_mean, y_pred_var] = IGPR(x_train, y_train, x_test, theta)
    % Define the inverse Gaussian transform function
    inverse_gaussian_transform = @(y) exp(y);

    % Transform the output data
    y_transformed = log(y_train);

    % Define the Gaussian kernel function
    kernel = @(x1, x2, theta) exp(-theta*(x1 - x2)'*(x1 - x2));

    % Compute the kernel matrices
    K = kernel(x_train, x_train, theta);
    K_star = kernel(x_train, x_test, theta);
    K_star_star = kernel(x_test, x_test, theta);

    % Perform Cholesky decomposition
    L = chol(K + 1e-6*eye(size(K)), 'lower');

    % Compute the alpha vector
    alpha = L'\(L\y_transformed);

    % Compute the mean predictions
    y_pred_mean = zeros(size(x_test));
    for i = 1:size(x_test, 1)
        y_pred_mean(i) = inverse_gaussian_transform(K_star(:, i)' * alpha);
    end

    % Compute the predictive variance
    v = L\K_star;
    y_pred_var = zeros(size(x_test));
    for i = 1:size(x_test, 1)
        y_pred_var(i) = inverse_gaussian_transform(K_star_star(i, i) - v(:, i)' * v(:, i));
    end
end