Математическая задача: масштабируйте график так, чтобы он соответствовал другому

Question

У меня есть 2 таблицы значений, и я хочу масштабировать первую так, чтобы она максимально соответствовала 2-й. Оба имеют одинаковую длину. Если оба графика нарисованы в виде графиков на диаграмме, они должны быть как можно ближе друг к другу. Но я не хочу квадратичных, а простых линейных весов. Моя проблема в том, что я понятия не имею, как на самом деле вычислить лучший коэффициент масштабирования из-за функции Abs.

Какой-то псевдокод:

//given:
float[] table1= ...;
float[] table2= ...;

//wanted:
float factor= ???; // I have no idea how to compute this

float remainingDifference=0;
for(int i=0; i<length; i++)
{
    float scaledValue=table1[i] * factor;
    //Sum up the differences. I use the Abs function because negative differences are differences too.
    remainingDifference += Abs(scaledValue - table2[i]);
}

Я хочу вычислить коэффициент масштабирования, чтобы оставшееся различие было минимальным.

Answer 1

Простые линейные веса сложны, как вы сказали.

a_n = first sequence
b_n = second sequence
c = scaling factor

Ваша функция невязки (суммы от i = 1 до N, количество баллов):

SUM( |a_i - c*b_i| )

Взятие производной по c дает:

  d/dc SUM( |a_i - c*b_i| )
= SUM( b_i * (a_i - c*b_i)/|a_i - c*b_i| )

Установить на 0 и решить для c сложно. Я не думаю, что есть аналитический способ сделать это. Вы можете попробовать https://math.stackexchange.com/, чтобы увидеть, есть ли у них какие-нибудь яркие идеи.

Однако, если вы работаете с квадратичными весами, это становится значительно проще:

  d/dc SUM( (a_i - c*b_i)^2 )
= SUM( 2*(a_i - c*b_i)* -c )
= -2c * SUM( a_i - c*b_i ) = 0
=> SUM(a_i) - c*SUM(b_i) = 0
=> c = SUM(a_i) / SUM(b_i)

Я настоятельно рекомендую последний подход, если вы можете.

Answer 2

Я бы предложил попробовать какой-нибудь вариант с Ньютоном Рафсоном.

Создайте функцию Diff (k), которая смотрит на разницу в площади между вашими двумя графиками между фиксированными маркерами A и B.

математически, я думаю, это будет целое число (от x = A до B) {f (x) - k * g (x)} dx

в любом случае реально вы могли бы просто вычесть значения,

например, если вы находитесь в диапазоне от X = -10 до 10, и у вас есть точка данных для f (i) и g (i) для каждого целого числа i в [-10, 10] (то есть 21 точка данных)

тогда вы просто суммируете (i = -10 до 10) {f (i) - k * g (i)}

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

и чем больше разница, вы ожидаете, чем больше разрыв

Итак, это должна быть довольно плавная функция (если у вас много точек данных)

так что вы хотите минимизировать Diff (k)

, поэтому вы хотите узнать, является ли производная, т.е.

так что просто сделайте Ньютона Рафсона с этой новой функцией D '(k)

Начните с k = 1, и решение должно быть достаточно быстрым

это, вероятно, даст вам оптимальное время вычислений

если вы хотите что-то более простое, просто начните с некоторых k1 и k2, которые находятся по обе стороны от 0

так сказать Diff (1.5) = -3 и Diff (2.9) = 7

так что тогда вы бы выбрали k, скажем, 3/10 пути (10 = 7 - -3) между 1,5 и 2,9

и в зависимости от того, дает ли это положительное или отрицательное значение, используйте его в качестве нового k1 или k2, промойте и повторите

Answer 3

В случае, если кто-нибудь столкнется с этим в будущем, вот код (c ++). Хитрость заключается в том, чтобы сначала отсортировать выборки по коэффициенту масштабирования, который обеспечит наилучшее соответствие для 2 выборок каждая.Затем начните с обоих концов итерации с коэффициентом, который приводит к минимальному абсолютному отклонению (L1-норма).

Все, кроме сортировки, имеет линейное время выполнения => Время выполнения равно O (n * log n)

/*
 * Find x so that the sum over std::abs(pA[i]-pB[i]*x) from i=0 to (n-1) is minimal
 * Then return x
 */
float linearFit(const float* pA, const float* pB, int n)
{
    /*
    * Algebraic solution is not possible for the general case
    * => iterative algorithm
    */

    if (n < 0)
        throw "linearFit has invalid argument: expected n >= 0";
    if (n == 0)
        return 0;//If there is nothing to fit, any factor is a perfect fit (sum is always 0)
    if (n == 1)
        return pA[0] / pB[0];//return x so that pA[0] = pB[0]*x

    //If you don't like this , use a std::vector :P
    std::unique_ptr<float[]> targetValues_(new float[n]);
    std::unique_ptr<int[]> indices_(new int[n]);
    //Get proper pointers:
    float* targetValues = targetValues_.get();//The value for x that would cause pA[i] = pB[i]*x
    int*   indices      = indices_.get();     //Indices of useful (not nan and not infinity) target values
    //The code above guarantees n > 1, so it is safe to get these pointers:
    int m = 0;//Number of useful target values
    for (int i = 0; i < n; i++)
    {
        float a = pA[i];
        float b = pB[i];
        float targetValue = a / b;
        targetValues[i] = targetValue;
        if (std::isfinite(targetValue))
        {
            indices[m++] = i;
        }
    }
    if (m <= 0)
        return 0;
    if (m == 1)
        return targetValues[indices[0]];//If there is only one target value, then it has to be the best one.

    //sort the indices by target value
    std::sort(indices, indices + m, [&](int ia, int ib){
        return targetValues[ia] < targetValues[ib];
    });

    //Start from the extremes and meet at the optimal solution somewhere in the middle:
    int l = 0;
    int r = m - 1;

    // m >= 2 is guaranteed => l > r
    float penaltyFactorL = std::abs(pB[indices[l]]);
    float penaltyFactorR = std::abs(pB[indices[r]]);
    while (l < r)
    {
        if (l == r - 1 && penaltyFactorL == penaltyFactorR)
        {
            break;
        }
        if (penaltyFactorL < penaltyFactorR)
        {
            l++;
            if (l < r)
            {
                penaltyFactorL += std::abs(pB[indices[l]]);
            }
        }
        else
        {
            r--;
            if (l < r)
            {
                penaltyFactorR += std::abs(pB[indices[r]]);
            }
        }
    }

    //return the best target value
    if (l == r)
        return targetValues[indices[l]];
    else
        return (targetValues[indices[l]] + targetValues[indices[r]])*0.5;
}