Многофакторная обобщенная линейная модель / Хабр

Поставим задачу автоматического подбора весовых коэффициентов вместе с настройкой функции связи обобщенной линейной модели.

Регрессионная модель:

$y=f(w_1x_1+\ldots+w_dx_d)$

Неизвестными считаются не только весовые коэффициенты $w_1,\ldots,w_d$ , но и функция .

Метод

Как выбрать веса одновременно с функцией связи? Если перейти к одной координате $z=w_1x_1+\ldots+w_dx_d$ , потребуется провести через точки кривую, принимающую вблизи каждого очень близкие значения.

Для нахождения функции используем кубические сплайны класса . Известно, что (при соответствующем выборе краевых условий) кубические сплайны минимизируют функционал

$J(f)=\int_a^b (f''(x))^2\,dx$

Значение функционала на кубическом сплайне мы возьмем в качестве меры качества подбора весовых коэффициентов. Однако, поскольку при растяжении точек $\mathbf x_i$ в раз функционал уменьшается в раз, наша целевая функция будет такой:

$J_0(\mathbf w)=(z_{\max}-z_{\min})^3\int_{z_{\min}}^{z_{\max}} (f''(z))^2\,dz$

В этой формуле функция - это кубический сплайн, построенный по точкам .

Эксперименты показали, что для оптимизации целевой функции не подходят стандартные алгоритмы BFGS и Nelder-Mead. В связи с этим будем использовать генетический алгоритм, реализованный в библиотеке PyGAD.

Пример

Сгенерируем в 3-мерном пространстве набор случайных точек $\mathbf x_i$ . Зададим вектор и вычислим $y_i = f(\mathbf w \cdot \mathbf x_i)+\varepsilon$ , где $f(z) = \sqrt{z}$ , $\varepsilon$ - случайный шум. Посмотрим, сможет ли метод восстановить веса вместе с функцией .

Код

import numpy as np
from scipy.optimize import minimize
from scipy.interpolate import CubicSpline
import scipy.integrate as integrate
import matplotlib.pyplot as plt
import pygad
import math
import random


def J(z, y):
    ind = np.argsort(z)
    cs = CubicSpline(z[ind], y[ind])
    d2 = cs.derivative(nu=2)
    integrand = lambda x: d2(x) ** 2
    I = sum([integrate.quad(integrand, d2.x[i], d2.x[i + 1])[0] for i in range(len(d2.x) - 1)])
    return I * (d2.x[-1] - d2.x[0]) ** 3


def train_spline(x, y):
    N, d = x.shape
    def f(w):
        res = J(np.dot(x, w), y)
        #print(res)
        return res

    num_generations = 100
    num_parents_mating = 10

    sol_per_pop = 20
    num_genes = d

    # веса от 0 до 5 (для простоты)
    gene_space = [{'low': 0, 'high': 5} for _ in range(d)]

    def on_generation(ga_instance):
        print(f"Generation = {ga_instance.generations_completed}")
        print(f"Fitness    = {ga_instance.best_solution(pop_fitness=ga_instance.last_generation_fitness)[1]}")

    def fitness_func(ga_instance, solution, solution_idx):
        return -f(solution)

    ga_instance = pygad.GA(num_generations=num_generations,
                           num_parents_mating=num_parents_mating,
                           sol_per_pop=sol_per_pop,
                           num_genes=num_genes,
                           fitness_func=fitness_func,
                           on_generation=on_generation,
                           gene_space=gene_space)

    ga_instance.run()

    solution, solution_fitness, solution_idx = ga_instance.best_solution(ga_instance.last_generation_fitness)

    return solution


def show_spline(z, y):
    ind = np.argsort(z)
    cs = CubicSpline(z[ind], y[ind])
    z_grid = np.linspace(np.min(z), np.max(z), 1000)
    y_grid = cs(z_grid)
    vfunc = np.vectorize(test_f)
    y_test = vfunc(z_grid)
    plt.figure(figsize=(8, 5))
    plt.plot(z, y, 'o', label='Data Points')
    plt.plot(z_grid, y_grid, label='Cubic Spline Interpolation')
    plt.plot(z_grid, y_test, label='Test')
    plt.xlabel('z')
    plt.ylabel('y')
    plt.legend()
    plt.grid()
    plt.show()



test_f = lambda x: math.sqrt(x)
test_w = [1, 2, 3]
test_w /= np.linalg.norm(test_w)
N = 100
d = 3
x = np.zeros((N, d))
y = np.zeros(N)
np.random.seed(123)
random.seed(123)
for i in range(N):
    for j in range(d):
        x[i, j] = np.random.random()
    y[i] = test_f(np.dot(x[i, :], test_w)) + 0.01 * (np.random.random() - 0.5)

w = train_spline(x, y)
print(w / w[0])

show_spline(np.dot(x, w) / np.linalg.norm(w), y)