Определение промежуточных данных

Определение промежуточных данных, Апроксимация, интерполяция и тому подобные шалости с табличными данными.

Пользователь

Сообщений: 732 Регистрация: 13.02.2019

25.10.2022 14:25:44

В данной теме попробую подобрать набор информации и ответов на вопросы:
Как найти промежуточные значения между заданными?
Как найти максимум функции заданной таблично?
и т.д.
~~Понятия Аппроксимация, Интерполяция и Экстраполяция гуглятся желающими самостоятельно.~~
Аппроксимация – под аппроксимационной кривой подразумевается некий полином (как правило, но не обязательно), график которого проходит наиболее близко к известным точкам (степень близости определяется по некоторому закону, как правило методу наименьших квадратов). При этом в известных значениях значения функции не обязательно совпадают с заданными значениями (в общем случае f(Xi) ≠ Yi).
Интерполяция – нахождение неизвестных промежуточных значений некоторой функции, по имеющемуся дискретному набору ее известных значений определенным способом. При этом в известных значениях значения функции совпадают с заданными значениями. Под интерполяционной кривой подразумевается некий полином (в нашем случае), график которого проходит через все известные точки.
Две особенности интерполяции:
- для получения полинома степени «n» требуется «n+1» заданная точка (например, полином первой степени f(x)=a·x+b требует две известные точки, или проще – линия строится по двум точкам, парабола по трём и т.д.);
- применение интерполяции методом «ближайшего соседа» в общем случае недопустимо.
Экстраполяция – особый тип аппроксимации, при котором функция аппроксимируется вне заданного интервала, а не между заданными значениями.

Часть 1. Аппроксимация с использованием встроенного функционала
Начнём с аппроксимации без использования макросов
Для начала - самый простой вариант:
1. Построить точечный график по имеющимся точкам;
2. На график добавить линию тренда с отображением уравнения на диаграмме;
3. Подобрать вид уравнения (степень полинома, вид уравнения...) который нравится/наиболее адекватно по субъективному восприятию отражает тенденцию изменения зависимости;
4. Скопировать уравнение линии тренда в ячейку и заменить "х" на "*А1^", где А1 - адрес ячейки в которой содержится значение аргумента, по которому требуется определить значение функции.
5. Пользоваться...
Если лень копировать уравнение с диаграммы, и хочется получить коэфф-ты полиномов, то можно сделать вот так:
Внимание! Не всегда коэфф-ты найденные нижеописанными способами будут соответствовать коэ-там на уравнении. Но об этом в 4-й части..
Аппроксимация полиномом

Код

Определение коэффициентов линейного уравнения
y=ax+b   
b=   1.397235    =ИНДЕКС(ЛИНЕЙН($B$2:$B$7;$A$2:$A$7;1);1;2)
a=   -0.027554    =ИНДЕКС(ЛИНЕЙН($B$2:$B$7;$A$2:$A$7;1);1;1)

Определение коэффициентов квадратичного уравнения
y=ax2+bx+c   
c=   1.662612    =ИНДЕКС(ЛИНЕЙН($B$2:$B$7;$A$2:$A$7^{1;2};);1;3)
b=   -0.071047    =ИНДЕКС(ЛИНЕЙН($B$2:$B$7;$A$2:$A$7^{1;2};);1;2)
a=   0.000849    =ИНДЕКС(ЛИНЕЙН($B$2:$B$7;$A$2:$A$7^{1;2};);1;1)

Определение коэффициентов кубического уравнения
y=ax3+bx²+cx+d   
d=   1.872900    =ИНДЕКС(ЛИНЕЙН($B$2:$B$7;$A$2:$A$7^{1;2;3};);1;4)
c=   -0.132192    =ИНДЕКС(ЛИНЕЙН($B$2:$B$7;$A$2:$A$7^{1;2;3};);1;3)
b=   0.004066    =ИНДЕКС(ЛИНЕЙН($B$2:$B$7;$A$2:$A$7^{1;2;3};);1;2)
a=   -0.000042    =ИНДЕКС(ЛИНЕЙН($B$2:$B$7;$A$2:$A$7^{1;2;3};);1;1)

Поиск коэффициентов для степеней аппроксимирующего полинома бОльших степений аналогичен.
Аппроксимация степенной функцией

Код
Коэф-ты уравнения y = а · Х ^ m m = -0.6320 =ИНДЕКС(ЛИНЕЙН(E3:E8;D3:D8);1) A = 1.1596 =ИНДЕКС(ЛИНЕЙН(E3:E8;D3:D8);2) a = 3.1888 =EXP(H4)

Аппроксимация логарифмической функцией

Код
Коэф-ты уравнения y = а · ln(x) + b a = -0.4676 =ИНДЕКС(ЛИНЕЙН(B3:B8;D3:D8);1) b = 2.0017 =ИНДЕКС(ЛИНЕЙН(B3:B8;D3:D8);2)

Использованием макросов
Преимуществом будет отсутствие необходимости использования ячеек листа. Макросы сохранённые в надстройку позволяют пользоваться ими без импорта в лист ну и т.д...
Недостатком - отсутствие визуализации решения.
Макрос Нумбер раз. Основной. Он собственно и находит все коэффициенты полинома.

Скрытый текст

Макрос нумбер 2. Необязательный. Подготовка исходных данных. Ввиду того что у меня данные могут браться из разных мест я их привожу к единому виду. Не стал выкидывать, ибо у меня 100% всё работает в такой связке, посему пускай будет.

Скрытый текст

Код

' Подпрограмма переопределения значений исходных данных
Function Подготовка_данных(Xvh As Variant, Yvh As Variant, _
                           Xvih() As Double, Yvih() As Double) As Variant    
    Dim NumYRows As Long, I As Long
    If TypeName(Xvh) = "Range" Then
    ' Если данные взяты с листа эксель
        NumXRows = Xvh.Count
        NumYRows = Yvh.Count        
        If NumXRows <> NumYRows Then MsgBox "Число Х и У не равно!"       
        ReDim Xvih(0 To NumXRows - 1)
        ReDim Yvih(0 To NumXRows - 1)
        If Xvh.Count = UBound(Xvh.Value2) Then
        ' Если данные X в столбцах
        For I = 1 To NumXRows
                Xvih(I - 1) = Xvh(I, 1).Value2
            Next I
        Else
        ' Если данные X в строках
            For I = 1 To NumXRows
                Xvih(I - 1) = Xvh(1, I).Value2
            Next I
        End If                
        If Yvh.Count = UBound(Yvh.Value2) Then
        ' Если данные Y в столбцах
            For I = 1 To NumYRows
                Yvih(I - 1) = Yvh(I, 1).Value2
            Next I
        Else
        ' Если данные Y в строках
            For I = 1 To NumYRows
                Yvih(I - 1) = Yvh(1, I).Value2
            Next I
        End If
    Else
    ' Если данные переданы из функции
        NumXRows = UBound(Xvh) + 1 : NumYRows = UBound(Yvh) + 1
        ReDim Xvih(0 To NumXRows - 1) : ReDim Yvih(0 To NumXRows - 1)

        For I = 0 To NumXRows - 1
            Xvih(I) = Xvh(I): Yvih(I) = Yvh(I)
        Next I
    End If
End Function

Макрос 3. Получение значения в указанной промежуточной точке. Большой вариант.

Скрытый текст

Код

' Аппроксимация для всего массива исходных данных
' В подпрограмму передаются все заданные точки,
' и аппроксимация ведётся по всем точкам!
' Xish   - массив данных по оси Х
' Yish   - массив данных по оси Y
' Xisk   - заданное значение Х, при котором требуется поиск соответствующего Y
' tipApr - признак аппроксимации.
'                                P - полиномом                           y = а · x^stepen + b · x^stepen-1 + ... + c
'                                S - степенной функцией            y = а · Х ^ m
'                                L - логарифмической функцией  y = а · ln(x)^stepen + b · ln(x)^stepen-1+ ... + c
'                                E - экспоненциальной ф-й         y = c *e ^(b * x)
' stepen - степень аппроксимирующего полинома или степень Х.

Public Function approx(Xish As Variant, Yish As Variant, _
                        ByVal Xisk As Double, _
                        Optional ByVal tipApr As String = "P", _
                        Optional ByVal stepen As Double = 1) As Variant
Dim i As Integer
Dim X() As Double
Dim Y() As Double
Dim Rez(0 To 8)
Dim cd() As Double
Подготовка_данных Xish, Yish, X, Y

Select Case tipApr
    Case "P" ' полиномиальная аппроксимация y = а1 · x^stepen + a2 · ln(x)^stepen-1 + ... + an
        stepen = Int(stepen)
        If stepen < 1 Then stepen = 1
        If stepen > 7 Then stepen = 7
        If (UBound(Y) - LBound(Y)) < stepen Then
            stepen = UBound(Y) - LBound(Y)
        End If
        Rez(0) = 0#
        Linia_trenda Y, X, stepen, cd
        For i = 1 To stepen + 1
            Rez(0) = Rez(0) + cd(i) * Xisk ^ (stepen - i + 1)
            Rez(i) = cd(i) ' Вывод коэффициентов a1...an
        Next i
        approx = Rez
               
    Case "S" ' Степенная аппроксимация y = а · Х ^ m
        For i = LBound(Y) To UBound(Y)
            X(i) = WorksheetFunction.Ln(X(i))
            Y(i) = WorksheetFunction.Ln(Y(i))
        Next i
        Linia_trenda Y, X, 1, cd
        Rez(0) = Exp(cd(2)) * Xisk ^ cd(1)
        Rez(1) = cd(2) ' Вывод "а"
        Rez(2) = cd(1) ' Вывод "m"
        approx = Rez
        
    Case "L" ' Логарифмическая аппроксимация 
             ' y = а1 · ln(x)^stepen + a2 · ln(x)^stepen-1+ ... + an
        For i = LBound(Y) To UBound(Y)
            X(i) = WorksheetFunction.Ln(X(i))
        Next i
        stepen = Int(stepen)
        If stepen < 1 Then stepen = 1
        If stepen > 7 Then stepen = 7
        If (UBound(Y) - LBound(Y)) < stepen Then
            stepen = UBound(Y) - LBound(Y)
        End If
        Rez(0) = 0
        Linia_trenda Y, X, stepen, cd
        For i = 1 To stepen + 1
            Rez(0) = Rez(0) + cd(i) * WorksheetFunction.Ln(Xisk) ^ (stepen - i + 1)
            Rez(i) = cd(i) ' Вывод коэффициентов a1...an
        Next i
        approx = Rez
        
    Case "E" ' Экспоненциальная аппроксимация y = c *e ^(b * x)
        For i = LBound(Y) To UBound(Y)
            Y(i) = WorksheetFunction.Ln(Y(i))
        Next i
        Linia_trenda Y, X, 1, cd
        Rez(0) = Exp(cd(2)) * Exp(cd(1) * Xisk)
        Rez(1) = Exp(cd(2)) ' Вывод "c"
        Rez(2) = cd(1)      ' Вывод "b"
        approx = Rez
               
    Case Else
        approx = "не верно задан тип аппроксимации"
        
End Select
End Function

Короткий вариант. Только полином

Скрытый текст

Как можно понять "большой вариант" не сложно дорабатывается под требуемые виды аппроксимирующих уравнений. Например у меня сейчас используются:
y = 1/(а · x^stepen + b · x^stepen-1 + ... + c)
y = a + b / x^stepen
y = 1/(а+b·x^stepen)
y = 1/(а+b·1/х)
y = а·x^2+b·1/х+c и т.д ...

За сим первую часть закрываю. Вторая часть будет о кусочной интерполяции ... чуть позже.

Изменено: tutochkin - 25.10.2022 16:54:21

Jack Famous

Пользователь

Сообщений: 12060 Регистрация: 07.11.2014

OS: Win 10 Корп. x64 | Excel 2016 x64: | Browser: Chrome

25.10.2022 15:43:43

Родственная тема: Прогноз. Как заполнить пропуски в данных по неизвестной закономерности при имеющихся соседних данных

Во всех делах очень полезно периодически ставить знак вопроса к тому, что вы с давних пор считали не требующим доказательств (Бертран Рассел) ►Благодарности сюда◄

tutochkin

Пользователь

Сообщений: 732 Регистрация: 13.02.2019

25.10.2022 16:32:04

Часть 2. Принятие существования недостоверности определения значений
На практике имеется следующее: есть некая точечно (не путать с «точно») заданная зависимость, и требуется определить значение между заданных точек (для аппроксимации или интерполяции). Например, есть заданные значения при значениях аргумента 0, 1, 2…9, а узнать надо при 0.5 и 8.5.

Как видно, в данном случае ни один из существующих полиномов линии тренда Excel не даёт интерполяционной кривой (проходящей через все точки), а наша цель:
- совпадение значений используемой функции значениям в реперных (заданных изначально) точках;
- совпадение или близость тенденции изменения параметров, т.е. в данном примере между т.1 и 2 функция должна иметь постоянно уменьшающееся значение, а не как для полинома 6-й степени: для значений менее 0,5 - уменьшение, а от 0,5 до 1 - значительное возрастание.
Одним из вариантов решения является применение разбиения известного количества заданных точек на несколько групп с малым количеством таким образом, чтобы:
- искомое значение Х содержалось в выбранном диапазоне;
- для выбранного количества точек можно было найти функцию, график которой будет проходить через все заданные точки.
Т.е. применение кусочной интерполяции.

Скрытый текст

Код

' Программа кусочного интерполирования точечно заданной функции.
' X - вектор исходных данных Х
' Y - вектор исходных данных Y
' Xisk - Х, при котором требуется отпределить Y
' toch - число точек по которым производится интерполяция после второй 
'        и до предпоследней точки.
' hvost - Указание количества точек используемых для начала и конца.
Public Function kus_interp(Xish As Variant, Yish As Variant, _
                           ByRef Xisk As Variant, _
                           Optional ByVal toch As Integer = 2, _
                           Optional ByVal hvost As Integer = 2) As Single
Dim i As Long
Dim j As Long
Dim k As Long
Dim xd() As Double
Dim yd() As Double
Dim cd() As Double
Dim X() As Double
Dim Y() As Double
Dim Xhvost() As Double
Dim Yhvost() As Double
Подготовка_данных Xish, Yish, X, Y
Imin = LBound(X)
Imax = UBound(X)
If toch > 8 Then toch = 8
If toch < 2 Then toch = 2
If (UBound(Y) - LBound(Y) + 1) < toch Then
    toch = UBound(Y) - LBound(Y) + 1
End If
If hvost > 8 Then hvost = 8
If hvost < 2 Then hvost = 2
If (UBound(Y) - LBound(Y) + 1) < hvost Then
    hvost = UBound(Y) - LBound(Y) + 1
End If
kus_interp = 0
ReDim xd(1 To toch) As Double
ReDim yd(1 To toch) As Double
ReDim Xhvost(1 To hvost) As Double
ReDim Yhvost(1 To hvost) As Double
If Xisk < X(Imin + 1) Then 
    For i = 1 To hvost
        Xhvost(i) = X(Imin - 1 + i)
        Yhvost(i) = Y(Imin - 1 + i)
    Next i
    Linia_trenda Yhvost, Xhvost, hvost - 1, cd
    For i = 1 To hvost
        kus_interp = kus_interp + cd(i) * Xisk ^ (hvost - i)
    Next i
Else
    If Xisk >= X(Imax - 1) Then 
        For i = 1 To hvost
            Xhvost(i) = X(Imax - hvost + i)
            Yhvost(i) = Y(Imax - hvost + i)
        Next i
        Linia_trenda Yhvost, Xhvost, hvost - 1, cd
        For i = 1 To hvost
            kus_interp = kus_interp + cd(i) * Xisk ^ (hvost - i)
        Next i
    Else 
        If toch = 2 Then
            For j = Imin To Imax - 1
             If Xisk < X(j + 1) Then
                    xd(1) = X(j): xd(2) = X(j + 1)
                    yd(1) = Y(j): yd(2) = Y(j + 1)
                    Linia_trenda yd, xd, 1, cd
                    For k = 1 To toch
                        kus_interp = kus_interp + cd(k) * Xisk ^ (toch - k)
                    Next k
                    Exit For
               Exit For
             End If
            Next j
        End If
        If toch > 2 Then
            For i = Imin + toch - 2 To Imax - 1
                If Xisk < X(i) Then
                    For k = 1 To toch
                        xd(k) = X(i + k - toch + 1)
                        yd(k) = Y(i + k - toch + 1)
                    Next k
                    Linia_trenda yd, xd, toch - 1, cd
                    For j = 1 To toch
                        kus_interp = kus_interp + cd(j) * Xisk ^ (toch - j)
                    Next j
                    Exit For
                End If
            Next i
        End If
     End If
  End If
End Function

Для нашего примера для определения значения функции:
- при Х=0,5 при линейной интерполяции для построения полинома используются 1-я и 2-я точки, а для интерполяции полиномом второй степени используются 1-я 2-я и 3-я точки;
- при Х=8,5 при линейной интерполяции для построения полинома используются 9-я и 10-я точки, а для интерполяции полиномом второй степени используются 8-я 9-я и 10-я точки.
Как видно из рисунка ниже, результат зависит от используемого метода интерполяции как внутри заданных значений, так и, в значительной мере, при экстраполяции.

Возможно использование и полиномов более высоких степеней, но как правило, достаточно полинома первого/второго порядка, т.к. зависимости далеко не всегда являются полиномными.
При поиске решения с применением полиномов 3-й и более степеней следует понимать, что в зависимости от выбранных точек результат будет разным. Т.е. если есть 4-ре точки (х1, х2, х3, х4), то результат, полученный по интерполяции по х1, х2, х3, будет отличаться от результата по х2,х3,х4.
Например, при выборе полинома второй степени и поиске значений при Х=1,5 данные, полученные по полиному, построенному точкам х=1 – 2 – 3, будут значительно отличаться от данных, полученных по полиному, построенному по точкам х=0 – 1 – 2. И так для большинства вариантов. В приведённом примере только на участке 4-5 есть совпадение полиномов.

Аналогичные проблемы будут и при использовании полиномов более высоких порядков. Например, для анализируемого случая строятся графики с использованием полиномов 7-го порядка (а·х7 + …), включающие первые и последние заданные точки. Первый интерполяционный график построен по точкам: х1-2-3-4-5-6-7-8 , второй – по х3 4 5 6 7 8 9 10. Наличие сильных расхождений в зоне пересечения (между точками х3-4-5-6-7-8 ) очевидно, но главную проблему представляет собой поведение 1-го инт.графика в зонах х1-2 и х7-8.

В условиях ограниченности исходных полиномы высоких порядков могут принести необоснованную погрешность расчёта. А для некоторых вариантов применение отличного от кусочного метода интерполяции просто невозможно.
Например, есть точки графика функции полинома 9-й степени. Провести интерполяцию данных штатными средствами Excel невозможно (ограничение полинома – 7-я степень). Кусочно-заданная функция потребует разбиение на несколько участков (более 3-х), и всё равно приведёт к значительным погрешностям даже в реперных точках. Либо у полученного графика будут аналогичные проблемы поиска данных в промежуточных точках:

Вывод – Не существует 100% достоверного математического способа определения промежуточных значений точечно заданной функции, за исключением частных вариантов. Любой из способов имеет некоторую погрешность (или допущение), и об этом следует помнить при расчётах.
В качестве ещё одного примера можно обратиться к графику в первом сообщении. На нём представлены 5 интерполяционных кривых, построенных по одним и тем же исходным точкам. Разница определения данных между точками в зависимости от вида интерполяции, очевидна.
Так же стоит отметить наличие ошибки расчётов величин с плавающей запятой (вещественных). Подробнее с действующим стандартом можно ознакомиться по адресу https://www.softelectro.ru/ieee754.html
Отдельным пунктом при кусочной интерполяции лежит понимание того, что имеется возможность управлять поведением функции при экстраполяции. Т.е. например в при кусочном интерполировании по всей функции использовать группы по 4-ре точки (интерполяция полиномом 3-й степени), а последний и первый участки, а так же зоны ДО и ЗА имеющимися данными, определять на основании полиномов первой степени (по парам первых и последних двух точек соответственно)....

Изменено: tutochkin - 25.10.2022 16:50:10

tutochkin

Пользователь

Сообщений: 732 Регистрация: 13.02.2019

27.10.2022 15:50:09

Не хочу копипастить данные от ZVI, МатросНаЗебре, MCH и других порядочных форумистов, посему немного вокруг и около

tutochkin

Пользователь

Сообщений: 732 Регистрация: 13.02.2019

14.12.2023 12:43:07

По мотивам задачки из раздела "работа"
Сразу скажу - на 100% я её не решил, а именно было требование:
1. получение массива промежуточных данных между точками графика максимально близко к тому что отрисовывает эксель. Понятие "максимально" размыто. Алгоритм по которому сглаживается соединяющая линия мне не известна. Однако пару методов я покажу ниже, при этом можно сказать"а вот совпадает то не совсем точно! И это будет правдой.
При этом "трассировки перемещения самолёта по небу" не оцифровывал никогда - задач таких нет и не будет. Чисто ради интереса занялся.
Итак, вот что у меня получилось:

Красная линия - исходная.
Синяя - кубический сплайн. В сети много решений, например тут
Зелёная - применение кусочной интерполяции (код выше в постах).
При этом надо для построения зелёной:
1. сделать две интерполяции, отдельно по Х и отдельно по У. введя доп.столбец (в ячейках А) нумерацию исходных данных. Да, можно было и переделать код макроса, но мне лень (я же не исполнитель решения задачи, а так... хобби). А так - реально все доп.столбцы можно убрать в код.
2. просчитать два столбца
3. Строить по данным, новым столбцам.

К неожиданным плюсам построения по кусочному методу могу отнести возможность применения разных законов для интерполяции по данным Х и У, что приводит к интересным результатам.

Упд. А вот это "На выходе таблица с назначаемым шагом по Х." для меня не совсем понятно. Шаг по Х, при том что исходные точки могут в эту сетку не влезть, и наличия трёх У для одного Х (при петлях) для меня не совсем понятно, и не рассматривал. Ну и за решение не брался.
Вот тут что в У писать при Х = 0,0038 ? Там 5 линий на данной вертикали

Изменено: tutochkin - 14.12.2023 13:11:13

Jack Famous

Пользователь

Сообщений: 12060 Регистрация: 07.11.2014

OS: Win 10 Корп. x64 | Excel 2016 x64: | Browser: Chrome

14.12.2023 13:26:08

tutochkin, приветствую!
Большое спасибо! Я думал насчёт интерполяции кусками, но сходу применить кубический сплайн (функция spline в файле от MCh) не вышло.

Цитата
tutochkin: Синяя - кубический сплайн. В сети много решений, например тут

инструмент по ссылке — странный. Сохранил, но пока не разобрался. Использовал spline.

Цитата
tutochkin: "На выходе таблица с назначаемым шагом по Х."

всё просто: пользователь задаёт шаг и между каждой парой X из исходной таблицы добавляются X с заданным шагом. Например, если X1 = 2, X2 = 6, Step = 1.5, то получится ряд: 2, 3.5, 5, 6 — в котором первый и последний X являются исходными.

Разберусь и покажу здесь, что у меня вышло с новыми данными от вас.

tutochkin

Пользователь

Сообщений: 732 Регистрация: 13.02.2019

14.12.2023 14:00:34

Цитата

Jack Famous написал:
всё просто: пользователь задаёт шаг и между каждой парой X из исходной таблицы добавляются X с заданным шагом. Например, если X1 = 2, X2 = 6, Step = 1.5, то получится ряд: 2, 3.5, 5, 6 — в котором первый и последний X являются исходными.

Ну задать Х то вполне можно. И даже в виде ряда. Но при этом как получить значение У ?
Повторюсь, вот тут возможны 5-ть вариантов...

В общем хотелось бы посмотреть как сие сделать. Без сарказма.

Прикрепленные файлы

2023-12-14_145310.jpg (214.63 КБ)

tutochkin Пользователь Сообщений: 732 Регистрация: 13.02.2019	#8 14.12.2023 21:06:34 Выдался свободный часик, налил крепкого чайку и поборол лень. Как то так в общем: Прикрепленные файлы Сплайн_через_кусочную_интерполяцию.xlsm (38.93 КБ)

Jack Famous Пользователь Сообщений: 12060 Регистрация: 07.11.2014 OS: Win 10 Корп. x64 \| Excel 2016 x64: \| Browser: Chrome	#9 15.12.2023 10:02:57 tutochkin, спасибо! Изучу вопрос. Во всех делах очень полезно периодически ставить знак вопроса к тому, что вы с давних пор считали не требующим доказательств (Бертран Рассел) ►Благодарности сюда◄

tutochkin Пользователь Сообщений: 732 Регистрация: 13.02.2019	#10 15.12.2023 10:36:55 Зачем надо не знаю, но залипательно...