Найти наилучшее значение веса для гладких ограниченных наименьших квадратов с Python? - PullRequest
1 голос
/ 16 апреля 2020

Мне нужно решить задачу наименьших квадратов без каких-либо известных оценок параметра. Я налагаю ограничение на то, чтобы желаемое решение было гладким (параметры модели меняются медленно), поэтому я минимизирую разницу между соседними параметрами (традиционное средство, используемое для решения этой геологической проблемы).

Ограничения реализуются путем размещения уравнений ограничения в виде строк в уравнении исходных данных d = Gm. Вспомогательный параметр w выбирается методом проб и ошибок (в некоторых учебниках w называется множителем Лагранжа).

У меня есть следующее:

G = np.array([[1,0,1,0,0,6],
             [1,0,0,1,0,6.708],
             [1,0,0,0,1,8.485],
             [0,1,1,0,0,7.616],
             [0,1,0,1,0,7],
             [0,1,0,0,1,7.616]])

d = np.array([[2.323],
             [2.543],
             [2.857],
             [2.64],
             [2.529],
             [2.553]])

Теперь добавляем ограничение произвольной w-взвешенной гладкости (w = 0,01):

w = 0.01
G = np.array([[1,0,1,0,0,6],
                 [1,0,0,1,0,6.708],
                 [1,0,0,0,1,8.485],
                 [0,1,1,0,0,7.616],
                 [0,1,0,1,0,7],
                 [0,1,0,0,1,7.616],
                 [w,-w,0,0,0,0],
                 [0,w,-w,0,0,0],
                 [0,0,w,-w,0,0],
                 [0,0,0,w,-w,0],
                 [0,0,0,0,w,-w]])

d = np.array([[2.323],
             [2.543],
             [2.857],
             [2.64],
             [2.529],
             [2.553],
              [0],
              [0],
              [0],
              [0],
              [0]])

Однако, выбирая правильное значение w является ключевым шагом для ограничения хорошего решения для параметров модели.

Итак, мой вопрос: с Python есть ли способ, которым я могу l oop по многим расчетным решениям с различными значениями для w и выбрать значение, которое использовалось для достижения решения с наилучшим качеством

1 Ответ

0 голосов
/ 16 апреля 2020

В представленном решении я буду обозначать G_0 как G без дополнительного ограничения и аналогично d_0 равно d без дополнительных нулей. Я также предполагаю, что вы читаете G_0 и d_0 откуда-то, и я называю их известными.

import numpy as np

def create_W(n_rows, w):
    W = -np.diagflat(np.ones(n_rows), 1)
    np.fill_diagonal(W, 1)
    return W

def solution_quality_metric(m):
    # this need to be implemented to determine what you mean by "best"

n_rows = 5
d_w = np.zeros(n_rows)

# choose range for w values for example w_min = 0, w_max = 1, dw = 0.01

best_m = -np.inf
best_w = w_min

for w in np.arange(w_min, w_max, dw):
    W = create_W(n_rows, w)

    G = np.concatenate([G_0, W], axis=0)
    d = np.concatenate([d_0, d_w])
    m = np.lstsq(G, d)

    if solution_quality_metric(m) > best_m:
        best_m = solution_quality_metric(m)
        best_w = w 

Этот код, очевидно, не будет работать как есть, поскольку вы не указали, что вы подразумеваете под «решением с наилучшим качеством». Для этого вам нужно реализовать функцию solution_quality_metric

...