Можно ли выполнить оценку максимального правдоподобия с ограничениями с помощью cvxpy, например, сделать это совместимым с DCP?

Вопрос или проблема

# Определите параметры распределения как переменные cp
# N = количество технологий, G = количество гауссовых распределений на технологию
means_NxG = cp.Variable((N, G))
stddevs_NxG = cp.Variable((N, G), nonneg=True)
scaling_factors_NxG = cp.Variable((N, G), nonneg=True)

# Вычислить W_NxP на основе параметров распределения
# P = количество точек в ценовом диапазоне
W_NxP = calculate_with_G_normal_distributions_W_NxP_cvxpy(
    price_domain_Px1=price_domain_Px1, 
    means_NxG=means_NxG, 
    stddevs_NxG=stddevs_NxG, 
    scaling_factors_NxG=scaling_factors_NxG
)

# Минимизировать сумму квадратов
objective_TxP = y_TxP - ACF_TxN @ W_NxP
objective_loss_1x1 = cp.norm(objective_TxP, 'fro')
objective = cp.Minimize(objective_loss_1x1)

prob = cp.Problem(
    objective=objective, 
    constraints=[некоторые ограничения, включая W_NxP]
)
prob.solve(verbose=True)

def calculate_with_G_normal_distributions_W_NxP_cvxpy(
        price_domain_Px1: Px1,
        means_NxG: cp.Variable, 
        stddevs_NxG: cp.Variable, 
        scaling_factors_NxG: cp.Variable
):
    N, G = means_NxG.shape  
    P = price_domain_Px1.shape[0]  # P = количество точек в ценовом диапазоне

    # Инициализируем W_NxP как список, который будет хранить результаты для каждой технологии
    W_NxP = []

    # Уплостите ценовой диапазон
    price_domain_Px_ = price_domain_Px1.flatten()

    # Построить W_NxP, используя выражения cvxpy
    for technology in range(N):  # Перебор технологий
        technology_density_1xP = 0  # Начните с нуля для каждой технологии
        for gaussian_parameter in range(G):  # Перебор гауссовых компонентов для каждой технологии
            mean = means_NxG[technology, gaussian_parameter]
            stddev = stddevs_NxG[technology, gaussian_parameter]
            scaling = scaling_factors_NxG[technology, gaussian_parameter]

            # Вычислить гауссовское pdf для каждой точки в ценовом диапазоне
            pdf_values_1xP = scaling * cp.exp(-0.5 * cp.square((price_domain_Px_ - mean) / stddev))

            # Суммируем значения pdf для каждого гаусса в этой технологии
            technology_density_1xP = technology_density_1xP + pdf_values_1xP

        # Измените размер результирующего вектора на 1xP и добавьте его в список
        W_NxP.append(cp.reshape(technology_density_1xP, (1, P)))

    # Соедините список векторов 1xP вертикально, чтобы создать окончательную матрицу NxP
    W_NxP = cp.vstack(W_NxP)
    return W_NxP

Ответ или решение

Для выполнения оценивания максимального правдоподобия с ограничениями, используя библиотеку CVXPY, полезно понимать принципы дисциплинированного выпуклого программирования (DCP). Ваша цель состоит в том, чтобы минимизировать сумму квадратов, отказавшись от использования градиентного спуска и обработав задачу в контексте CVXPY. Давайте разберёмся, как это сделать, соблюдая принципы DCP.

Понимание задачи

Сначала нужно выделить, что переменные means, stddevs и scaling_factors являются переменными, которые CVXPY должен интерпретировать как DCP-соответствующие. Необходимо гарантировать, что операции с ними сохраняют DCP-свойства. Ваша задача состоит в расчёте функции плотности распределения (PDF), которая подлежит оптимизации в рамках задачи минимизации.

Обновленный код

Вот пример, как вы можете изменить свой код, чтобы он был DCP-соответствующим, и исправить функцию calculate_with_G_normal_distributions_W_NxP_cvxpy.

import cvxpy as cp
import numpy as np

# Определите параметры вашей задачи
N = ...  # количество технологий
G = ...  # количество гауссовых распределений на технологию
P = ...  # количество точек в области цен
price_domain_Px1 = np.array([...])  # ваш массив цен
y_TxP = cp.Parameter((T, P))  # ожидаемые значения
ACF_TxN = cp.Parameter((T, N))  # некоторая матрица для перемножения

# Определите параметры распределения как переменные CVXPY
means_NxG = cp.Variable((N, G))
stddevs_NxG = cp.Variable((N, G), nonneg=True)
scaling_factors_NxG = cp.Variable((N, G), nonneg=True)

# Обновлённая функция для расчёта W_NxP
def calculate_with_G_normal_distributions_W_NxP_cvxpy(price_domain_Px1, means_NxG, stddevs_NxG, scaling_factors_NxG):
    N, G = means_NxG.shape  
    P = price_domain_Px1.shape[0]

    # Инициализация матрицы плотности
    W_NxP = []

    for technology in range(N):
        technology_density_1xP = cp.Constant(0)  # Начало с константы 0

        for gaussian_parameter in range(G):
            mean = means_NxG[technology, gaussian_parameter]
            stddev = stddevs_NxG[technology, gaussian_parameter]
            scaling = scaling_factors_NxG[technology, gaussian_parameter]

            # Расчёт pdf значений для гауссового распределения
            pdf_values_1xP = scaling * cp.exp(-0.5 * cp.square((price_domain_Px1 - mean) / stddev))

            technology_density_1xP += pdf_values_1xP  # Суммирование pdf значений

        W_NxP.append(technology_density_1xP)

    return cp.hstack(W_NxP)

# Используйте вычисленное W_NxP в оптимизационной задаче
W_NxP = calculate_with_G_normal_distributions_W_NxP_cvxpy(price_domain_Px1, means_NxG, stddevs_NxG, scaling_factors_NxG)

# Определите функцию потерь
objective_TxP = y_TxP - ACF_TxN @ W_NxP
objective_loss_1x1 = cp.norm(objective_TxP, 'fro')
objective = cp.Minimize(objective_loss_1x1)

# Определение проблемы
prob = cp.Problem(
    objective=objective,
    constraints=[W_NxP >= 0]  # Например, ваши ограничения
)

# Решение проблемы
prob.solve(verbose=True)

# Результаты
print("Оптимальные значения:", W_NxP.value)

Важные моменты

1. DCP-соответствие: Убедитесь, что операции с переменными производятся согласно правилам DCP. В данной реализации мы используем cp.Constant(0) для инициализации technology_density_1xP, что позволяет избежать проблем с преобразованием типов.

2. Непараметризованные переменные: Вы можете задать дополнительные ограничения для переменных, чтобы улучшить качество решения (например, на пороги).

3. Оптимизация результатов: Убедитесь, что используемые вами матрицы и данные корректны и соответствуют размерам. Возможно, вам потребуется провести предобработку данных.

Таким образом, используя преобразованные к DCP-совместимым функциям и правильные структуры, вы сможете эффективно решать задачи максимального правдоподобия с помощью CVXPY.