• Генерация координат: Программа создает два вектора случайных координат: rx и ry, которые находятся в диапазоне от 0 до 1. Эти координаты представляют собой точки на плоскости.

• Проверка попадания: Точка считается находящейся внутри круга, если расстояние dist от начала координат (0, 0) меньше радиуса, равного 1. Это условие можно проверить с помощью следующей формулы: rx² + ry² < 1.

• Вычисление: Итоговая оценка числа π (pi_est) вычисляется как отношение количества точек, попавших внутрь круга, к общему количеству сгенерированных точек.

Что такое полная бета-функция?

Чтобы понять неполную версию, нужно вспомнить полную. Полная бета-функция $\(B(a, b)\)$ — это определенный интеграл от нуля до единицы, который зависит от двух параметров $\(a\)$ и $\(b\)$:

$$\(B(a, b) = \int_{0}^{1} t^{a-1} (1-t)^{b-1} \, dt\)$$

Что такое неполная бета-функция?

В неполной бета-функции верхний предел интеграла заменяется на переменную $\(x\)$ (где $\(0 \le x \le 1\))$. Это значит, что мы интегрируем функцию не до конца, а только на отрезке от $\(0\)$ до $\(x\)$.

Обозначается она как $\(B_x(a, b)\)$ и определяется следующим образом:

$$\(B_x(a, b) = \int_{0}^{x} t^{a-1} (1-t)^{b-1} \, dt\)$$

from tabulate import tabulate
import numpy as np

def generate_data(rank, num_points, noise_level):
    np.random.seed(42)
    x = np.linspace(-10, 10, num_points)
    y = np.linspace(-10, 10, num_points)
    if rank == 1:
        z = 3 * x + 5 + noise_level * np.random.randn(num_points)
        data = np.column_stack((x, y, z))
    elif rank == 2:
        z = 2 * x + 3 * y + 5 + noise_level * np.random.randn(num_points)
        data = np.column_stack((x, y, z))
    elif rank == 3:
        z = 2 * x**2 + 3 * y**2 + 5 + noise_level * np.random.randn(num_points)
        data = np.column_stack((x, y, z))
    return data

def GMT_trend2d(data, rank):
    import numpy as np
    from sklearn.linear_model import LinearRegression
    # scale factor for normally distributed data is 1.4826
    # https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.median_abs_deviation.html
    MAD_NORMALIZE = 1.4826
    # significance value
    sig_threshold = 0.51

    if rank not in [1,2,3]:
        raise Exception('Number of model parameters "rank" should be 1, 2, or 3')

    #see gmt_stat.c
    def gmtstat_f_q (chisq1, nu1, chisq2, nu2):
        import scipy.special as sc

        if chisq1 == 0.0:
            return 1
        if chisq2 == 0.0:
            return 0
        return sc.betainc(0.5*nu2, 0.5*nu1, chisq2/(chisq2+chisq1))

    if rank in [2,3]:
        x = data[:,0]
        x = np.interp(x, (x.min(), x.max()), (-1, +1))
    if rank == 3:
        y = data[:,1]
        y = np.interp(y, (y.min(), y.max()), (-1, +1))
    z = data[:,2]
    w = np.ones(z.shape)

    if rank == 1:
        xy = np.expand_dims(np.zeros(z.shape),1)
    elif rank == 2:
        xy = np.expand_dims(x,1)
    elif rank == 3:
        xy = np.stack([x,y]).transpose()

    # create linear regression object
    mlr = LinearRegression()

    chisqs = []
    coeffs = []
    while True:
        # fit linear regression
        mlr.fit(xy, z, sample_weight=w)

        r = np.abs(z - mlr.predict(xy))
        chisq = np.sum((r**2*w))/(z.size-3)
        chisqs.append(chisq)
        k = 1.5 * MAD_NORMALIZE * np.median(r)
        w = np.where(r <= k, 1, (2*k/r) - (k * k/(r**2)))
        sig = 1 if len(chisqs)==1 else gmtstat_f_q(chisqs[-1], z.size-3, chisqs[-2], z.size-3)
        # Go back to previous model only if previous chisq < current chisq
        if len(chisqs)==1 or chisqs[-2] > chisqs[-1]:
            coeffs = [mlr.intercept_, *mlr.coef_]

        #print ('chisq', chisq, 'significant', sig)
        if sig < sig_threshold:
            break

    # get the slope and intercept of the line best fit
    return (coeffs[:rank])

def calculate_mse(data, coeffs, rank):
    z_actual = data[:, 2]

    if rank == 1:
        z_predicted = coeffs[0]
    elif rank == 2:
        # Interpolate x the same way as in GMT_trend2d
        x = data[:, 0]
        x_interp = np.interp(x, (x.min(), x.max()), (-1, +1))
        z_predicted = coeffs[0] + coeffs[1] * x_interp
    elif rank == 3:
        # Interpolate x and y the same way as in GMT_trend2d
        x = data[:, 0]
        x_interp = np.interp(x, (x.min(), x.max()), (-1, +1))
        y = data[:, 1]
        y_interp = np.interp(y, (y.min(), y.max()), (-1, +1))
        z_predicted = coeffs[0] + coeffs[1] * x_interp + coeffs[2] * y_interp

    mse = np.mean((z_actual - z_predicted) ** 2)
    return mse

def test_mse(num_points = 100*1000, ranks = [1, 2, 3], noise_levels = [0, 1, 10, 50]):
    import warnings

    results = []
    # Suppress the specific warning
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", category=RuntimeWarning)
        for rank in ranks:
            for noise_level in noise_levels:
                data = generate_data(rank, num_points, noise_level)
                # round the output
                coeffs_gmt = [v.round(8) for v in GMT_trend2d(data, rank)]

                mse_gmt = np.round(calculate_mse(data, coeffs_gmt, rank), 0)

                results.append([rank, noise_level, mse_gmt])

    headers = ["Rank", "Noise Level", "GMT_trend2d, MSE"]
    print(tabulate(results, headers=headers))

test_mse()

using namespace np;
using namespace scipy;
using namespace sklearn;

auto generate_data(auto rank, auto num_points, auto noise_level) {
    random::seed(42);
    auto x = linspace(-10.0, 10.0, num_points);
    auto y = linspace(-10.0, 10.0, num_points);
    if (rank == 1) {
        auto z = 3 * x + 5 + noise_level * random::randn(num_points);
        return column_stack(x, y, z);
    }
    if (rank == 2) {
        auto z = 2 * x + 3 * y + 5 + noise_level * random::randn(num_points);
        return column_stack(x, y, z);
    }
    auto z = 2 * x * x + 3 * y * y + 5 + noise_level * random::randn(num_points);
    return column_stack(x, y, z);
}

auto GMT_trend2d(const Array<float_> &data, int rank) {
    float_ MAD_NORMALIZE = 1.4826;
    float_ sig_threshold = 0.51;

    if (rank != 1 && rank != 2 && rank != 3) {
        throw sklearn::RuntimeError("Number of model parameters \"rank\" should be 1, 2, or 3");
    }

    auto gmtstat_f_q = [](float_ chisq1, float_ nu1, float_ chisq2, float_ nu2) {
        if (chisq1 == 0.0) return 1.0;
        if (chisq2 == 0.0) return 0.0;
        return scipy::special::betainc(0.5 * nu2, 0.5 * nu1, chisq2 / (chisq2 + chisq1));
    };

    Array<float_> x;
    if (rank == 2 || rank == 3) {
        auto x_ = data[":,0"];
        x = interp(x_, Array<float_>{x_.min(), x_.max()}, Array<float_>{-1, +1});
    }
    Array<float_> y;
    if (rank == 3) {
        auto y_ = data[":,1"];
        y = interp(y_, Array<float_>{y_.min(), y_.max()}, Array<float_>{-1, +1});
    }
    auto z = data[":, 2"].copy();
    Array<float_> w = ones(z.shape()).copy();

    Array<float_> xy;
    if (rank == 1) {
        xy = expand_dims(zeros(z.shape()), 1);
    } else if (rank == 2) {
        xy = expand_dims(x, 1);
    } else if (rank == 3) {
        xy = stack(x, y).transpose();
    }

    auto mlr = linear_model::LinearRegression{};

    std::vector<float_> chisqs;
    Array<float_> coeffs;

    while (true) {
        mlr.fit(xy, z, w);

        auto r = abs_sub(z, mlr.predict(xy));
        auto chisq = sum_sq_weighted(r, w) / static_cast<float_>(z.size() - 3);
        chisqs.push_back(chisq);

        auto k = 1.5 * MAD_NORMALIZE * median(r);
        w = where_tukey(r, k);
        auto sig = (chisqs.size() == 1 ? 1 : gmtstat_f_q(chisqs[chisqs.size() - 1], static_cast<float_>(z.size() - 3), chisqs[chisqs.size() - 2], static_cast<float_>(z.size() - 3)));
        if (chisqs.size() == 1 or chisqs[chisqs.size() - 2] > chisqs[chisqs.size() - 1]) {
            coeffs = mlr.coeffs_();
        }

        if (sig < sig_threshold) {
            break;
        }
    }
    auto result = Array<float_>{};
    for (int i = 0; i < rank; ++i) {
        result = append(result, Array<float_>{coeffs.get(i)});
    }
    return result;

auto calculate_mse(const Array<float_> &data, const Array<float_> &coeffs, int rank) {
    auto z_actual = data[":,2"].copy();
    Array<float_> z_predicted;

    if (rank == 1) {
        z_predicted = coeffs.get(0) * ones(z_actual.shape());
    } else if (rank == 2) {
        // Interpolate x the same way as in GMT_trend2d
        auto x_ = data[":,0"];
        auto x = interp(x_, Array<float_>{x_.min(), x_.max()}, Array<float_>{-1, +1});
        z_predicted = coeffs.get(0) + coeffs.get(1) * x;
    } else if (rank == 3) {
        // Interpolate x and y the same way as in GMT_trend2d
        auto x_ = data[":,0"];
        auto x = interp(x_, Array<float_>{x_.min(), x_.max()}, Array<float_>{-1, +1});
        auto y_ = data[":,1"];
        auto y = interp(y_, Array<float_>{y_.min(), y_.max()}, Array<float_>{-1, +1});
        z_predicted = coeffs.get(0) + coeffs.get(1) * x + coeffs.get(2) * y;
    }
    using sklearn::metrics::mean_squared_error;
    using sklearn::metrics::MeanSquaredErrorParameters;
    MeanSquaredErrorParameters<Array<float_>> params{.y_true = z_actual, .y_pred = z_predicted};
    return mean_squared_error(params);
}

void test_mse(int num_points = 100 * 1000, const std::vector<int> &ranks = {1, 2, 3}, const std::vector<int> &noise_levels = {0, 1, 10, 50}) {
    std::vector<std::tuple<int, int, float_>> results;
    for (auto rank: ranks) {
        for (auto noise_level: noise_levels) {
            auto data = generate_data(rank, num_points, noise_level);
            auto coeffs_gmt = GMT_trend2d(data, rank);
            // round coefficients to 8 decimal places
            auto coeffs_rounded = coeffs_gmt.copy();
            for (std::size_t i = 0; i < coeffs_rounded.size(); ++i) {
                coeffs_rounded.set(i, std::round(coeffs_rounded.get(i) * 1e8) / 1e8);
            }
            auto mse_gmt = calculate_mse(data, coeffs_rounded, rank);
            // round MSE to zero decimal places
            mse_gmt = std::round(mse_gmt);
            results.emplace_back(rank, noise_level, mse_gmt);
        }
    }
    // print table
    std::cout << "Rank\tNoise Level\tGMT_trend2d, MSE\n";
    for (const auto &[rank, noise_level, mse]: results) {
        std::cout << rank << "\t" << noise_level << "\t" << mse << "\n";
    }
}