CC BY
20
Вычислительные технологии
Том 24, № 5, 2019
О мере вариабельности оценки параметров
>к
в статистике интервальных данных*
С. П. ШАРЫй
Институт вычислительных технологий СО РАН, Новосибирск, Россия Новосибирский государственный университет, Россия Контактный e-mail: shary@ict.nsc.ru
Представлена одна из возможных конструкций количественной меры вариабельности оценки параметров в задаче восстановления линейной зависимости по интервальным данным. Она показывает степень изменчивости и неоднозначности оценки, и необходимость ее введения диктуется тем обстоятельством, что в условиях интервальных данных ответ к задаче, как правило, не является единственным. Дается вывод новой меры вариабельности, и ее применение иллюстрируется на конкретных примерах. В заключение статьи обсуждаются содержательные мотивации и смысл новой величины с точки зрения других известных характеристик данных.
Ключевые слова: задача восстановления зависимостей, интервальная неопределенность данных, сильное согласование параметров и данных, мера вариабельности.
Библиографическая ссылка: Шарый С.П. О мере вариабельности оценки параметров в статистике интервальных данных // Вычислительные технологии. 2019. Т. 24, № 5. С. 90-108. DOI: 10.25743/ICT.2019.24.5.008.
Введение и постановка задачи
Цель работы — представить одну из возможных конструкций количественной меры вариабельности и чувствительности для оценки параметров функциональной зависимости в статистике интервальных данных. Это сравнительно молодая ветвь современной науки об обработке данных, которая не опирается на аппарат теории вероятностей, но широко использует методы интервального анализа [1-4].
Термином "вариабельность" мы называем степень изменчивости и неоднозначности оценки, и необходимость ее введения диктуется тем обстоятельством, что в задачах обработки интервальных данных ответ, как правило, неединствен. Обычно мы получаем целое множество различных оценок, одинаково пригодных в качестве ответов к задаче и согласующихся с ее данными. То, насколько мало или обширно это множество, как раз-таки и характеризуется термином "вариабельность".
Термин "чувствительность оценки" относится к ситуациям, когда находится одна оценка параметров (это так называемые задачи точечного оценивания), и нам необходимо определить, насколько сильно она будет меняться при вариациях (изменениях) входных данных задачи.
* Title translation and abstract in English can be found on page 108.
© ИВТ СО РАН, 2019.
В традиционной теоретико-вероятностной статистике оценки тех или иных параметров, как известно, сами являются случайными величинами, а мерой их чувствительности и вариабельности может служить дисперсия оценки (см., к примеру, [5]). Она выражает как меру рассеяния значений оценки, так и меру чувствительности этой оценки к изменениям входных данных задачи.
Что может служить аналогом этой величины в статистике интервальных данных, которая не оперирует вероятностными понятиями и в которой данные не несут вероятностных распределений?
В отношении вариабельности оценки параметров ответ на этот вопрос кажется достаточно очевидным: ею может стать любая величина, характеризующая размеры множества решений (называемого часто также "информационным множеством") задачи, если оно непусто (рис. 1). Можно даже просто брать интервальные оценки множества решений, получаемые теми или иными методами интервального анализа. Определенный недостаток этого варианта (хорошо заметный в сравнении с дисперсией) — излишняя детализация ответа, выдаваемого в виде бруса в Кга или интервального вектора какого-нибудь другого вида, большое количество информации, которую еще необходимо "переварить" и привести к компактной и выразительной форме. Другой недостаток — относительная сложность нахождения такой оценки.
Рис. 1. Мерой вариабельности может служить размер множества решений задачи
Требуется относительно несложная и эффективно вычислимая величина, выражаемая одним числом, которая давала бы общее агрегированное представление об интересующем нас предмете. Аналогично дисперсии она может служить "прикидочной" характеристикой качества оценок при практическом решении различных задач обработки интервальных данных.
Конкретная форма этой величины, конечно, должна зависеть как от задачи, так и от метода ее решения. Мы рассмотрим ниже возможную конструкцию меры вариабельности в задаче восстановления зависимостей, когда по данным измерений или наблюдений требуется построить функциональную зависимость заданного вида, которая наилучшим образом "приближает" эти данные или "согласуется" с ними. Более точно,
?
исследуется задача определения параметров Х\, х2, ... , хп линейной функции вида
Ь = ххах + ... + хпап (1)
по ряду измерений-наблюдений независимых переменных а1, а2,..., ап и соответствующих значений функции Ь. В качестве способа оценивания коэффициентов х1, х2,..., хп рассматривается так называемый метод максимума согласования, предложенный и развитый в работах [6-9] и других публикациях автора.
1. Формулировка основного результата
Исходные данные задачи восстановления зависимости (1) — это набор значений независимых (предикторных) переменных и зависимой (критериальной) переменной для функции (1), которые получены в результате т измерений или наблюдений:
О-Ш 0-12, ... a1n, Я-2Ъ «22, ... Q>2n,
&m1, &m2, . . . &mn,
Мы предполагаем, что эти данные неточны и имеют интервальную неопределенность, возникшую вследствие погрешностей измерений и т. п. (рис. 2). Данные (2) и другие интервальные величины всюду в тексте выделяются жирным математическим шрифтом согласно неформальному международному стандарту [10].
Ь1, Ь2,
Ьщ.
(2)
Рис. 2. Иллюстрация задачи восстановления линейной зависимости по интервальным данным
Первый индекс у величин из (2) означает номер измерения, а второй — номер независимой переменной. Подставляя формально данные (2) в равенство (1), получим для нахождения оценок (х1,х2,... , хп) параметров линейной зависимости (1) интервальную систему линейных алгебраических уравнений (ИСЛАУ)
аи^1 + 0,12X2 + ... + а 1пхп = 6Ь 0,21^1 + а,22Х2 + ... + а,2пХп = Ь2,
< . . .. . : (3)
+ «т.2^2 + ... + ^тп^п bm,
или, кратко,
Ах = b
с т х п-матрицей А = (а^) и т-вектором в правой части Ь.
Множества параметров, согласующихся с данными измерений (2) в том или ином смысле, образуют различные множества решений для системы уравнений (3), (4). Наиболее популярными являются объединенное множество решений
БиПг(А, Ь) = { х Е Ега | Ах = Ь для некоторых А Е А и Ь Е Ь },
которое соответствует так называемому слабому согласованию параметров зависимости (1) и данных (2), и допусковое множество решений
(А, Ь) = { х Е Rra | Ах Е b для любых матриц А Е А},
-tol\
которое соответствует так называемому сильному согласованию параметров зависимости (1) и данных (2) (см. подробности, к примеру, в [4]).
Далее будем считать, что решение рассматриваемой задачи восстановления зависимости (1) находится методом максимума согласования [6-9]. Напомним, что оценкой параметров зависимости (1) в нем берется точка максимума специального "распознающего функционала", т. е. некоторой функции, которая дает количественную "меру согласования" этой оценки с интервальными эмпирическими данными (2).
Метод максимума согласования имеет различные версии — "слабую" и "сильную", которые отличаются пониманием того, как интервальные данные задачи должны "согласовываться" с восстанавливаемой зависимостью. "Сильная" версия метода максимума согласования [7] имеет ряд теоретических и практических преимуществ перед обычной "слабой" версией. Это полиномиальная трудоемкость, конечная вариабельность оценки и ее робастность и др. По этой причине ниже рассмотрим сильную версию метода максимума согласования, которой соответствует допусковое множество решений ^^(А, Ь) интервальной системы уравнений (3), (4). Его распознающий функционал обычно обозначается символом "То1", и он имеет вид
Tol (х, А, Ь)
{
где
mm < rad bi —
1< i<m
mid bi — ^^ a
3 = 1
'ij xj
}
(5)
rad bi = -(bi — Ь,) и mid bi
2
+ ь,)
— соответственно радиусы и середины компонент вектора правой части Ь, а модуль интервала — это максимум модулей его точек (или концов). Таким образом, для решения задачи нахождения параметров линейной зависимости (1) по набору данных (2) нужно найти безусловный максимум по всем х из Rra —
Tol (х, A, b) ^ max,
и доставляющий этот максимум вектор х = argmax^Rn Tol (ж, А, Ь) принимается за искомую оценку параметров функции (1).
Если max Tol > 0, то множество решений Stoi(А, Ь), т.е. множество параметров, сильно согласующихся с данными, непусто, и х Е Stoi(А, Ь). Если max Tol < 0, то множество решений Etoi (А, Ь) пусто и наборов параметров, которые сильно согласуются
с данными (2), не существует, но х все равно обеспечивает наилучшее "согласование" построенной линейной функции с данными (2) (точнее, наименьшее несогласование).
В качестве величины, характеризующей чувствительность и вариабельность оценки вектора параметров х = (Xi, х2,... , хп) линейной зависимости (1), которая получена методом максимума согласования по данным (2), предлагаем
. || arg max Tol 112
IVE (A, b) = л/П max Tol • ( min cond2A ) •-R-. (6)
R" Vie^ 2 ) ||b||2 v 7
В формуле (6) переменные имеют следующий смысл. п — размерность вектора параметров восстанавливаемой функции. || • ||2 — евклидова норма векторов (2-норма), определяемая в Rra как
I 1/2
||X||2 = ( > >|2 1
/га \ 1/2
ОН ■
cond2(A) — спектральное число обусловленности матрицы А, определяемое как
cond2A = max( ) ,
т. е. как отношение максимального и минимального сингулярных чисел матрицы А. cond2(A) является распространением на прямоугольный случай хорошо известного в вычислительной линейной алгебре понятия числа обусловленности квадратной матрицы (см. подробности в [11, 12]). Точка b — это "наиболее представительная" точка интервального вектора b, определяемая как
b = 1 (| midb + radb| + |midb — radb|), (7)
где операции "mid" и "rad" применяются покомпонентно.
Несмотря на определенную формулу (7) для b, следует отметить, что введение этой точки в значительной степени является вопросом здравого смысла. Общий подход к определению b состоит в том, что это должна быть некоторая "наиболее представительная" точка из вектора правой части b и в особых ситуациях этот выбор может отличаться от формулы (7). Например, если мы решим, что ||midb||2 является "достаточно большим" по сравнению с нормой радиуса ||radb||2, то мы можем взять просто b = mid b, т. е. как среднюю точку интервального вектора b.
Строго говоря, равенство b = midb справедливо для интервальных векторов, в которых концы интервальных компонент имеют одинаковые знаки, но иногда мы можем решить, что равенство b = mid b должно выполняться и без этого условия. Кроме того, совершенно не воспрещен любой другой выбор точки b, если он имеет убедительное практическое обоснование.
Сам символ IVE образован как аббревиатура английской фразы Interval Variability of the Estimate. Ниже мы строго покажем, что величина IVE адекватно характеризует размеры непустого допускового множества решений. Но этот факт может быть также понят на основе аналитической интуиции и наглядных образных соображений.
Допусковое множество решений интервальной системы линейных алгебраических уравнений является, как известно, множеством нулевого уровня распознающего функционала Tol или, другими словами, пересечением подграфика этого функционала с координатной плоскостью Tol = 0 (рис. 3):
Etol (A, b) = {х е Rn | Tol (х, A,b) > 0 }.
max Tol
Рис. 3. Величина максимума распознающего функционала дает представление о размерах допускового множества решений
Рис. 4. Помимо максимума распознающего функционала на размеры множества решений влияет также "крутизна" графика функционала
Как следствие, величина максимума распознающего функционала может служить, при прочих равных условиях, мерилом того, насколько обширно или узко допусковое множество решений. Чем больше max Tol, тем больше размеры допускового множества решений, и наоборот. Другие факторы, обеспечивающие "прочие равные условия", — это наклон гиперплоскостей, из кусков которых составлен полиэдральный график функционала Tol (прямых в одномерном случае на рис. 3 и 4). Наклон гиперплоскостей определяется коэффициентами задающих их уравнений, которые являются концами интервалов данных (2). Величину этого наклона обобщенным образом характеризует число обусловленности матрицы данных. Наконец, множитель
— это масштабирующий коэффициент, с помощью которого обеспечивается соизмеримость окончательной величины решению и вектору правой части системы уравнений. Так и получается формула (6).
|| argmax Tol ||2 ||ж||2
Й2 N2
Помимо (6) в качестве меры относительной вариабельности решения определенный смысл может иметь величина
( • j Л max^Rn Tol
п\ min cond2A ---, (8)
Vaga ) ||5|2
которую мы никак специально не обозначаем.
2. Теория: обоснование меры вариабельности
Рассматривая общий случай, будем считать, что количество измерений т может не совпадать с числом п неизвестных параметров линейной функции (1), причем т > п. Иными словами, измерений (наблюдений) проводится не меньше, чем количество искомых параметров функции. Тогда интервальная система линейных алгебраических уравнений (3), (4) является либо квадратной, либо переопределенной, т.е. имеющей интервальную т х п-матрицу А с т > п.
2.1. Оценка возмущений решения прямоугольных линейных систем
Отправной точкой наших конструкций, обосновывающих выбор меры вариабельности именно в виде величин (6) и (8), будет известное неравенство, оценивающее возмущение Лх ненулевого решения системы линейных алгебраических уравнений Ах = Ь в зависимости от изменений ЛЬ в векторе правой части (см., к примеру, [11-14]):
l|Ax|o d ||ДЬ||2
и и < COnd2^ и,,, .
М2 ть
(9)
Обычно оно рассматривается для квадратных систем линейных алгебраических уравнений, когда т = п, но в случае евклидовой нормы векторов и спектрального числа обусловленности матриц это неравенство справедливо и для более общего случая, при т > п. Кратко напомним его вывод в этой ситуации.
Если
Ах = Ь и А(х + Лх) = Ь + ЛЬ,
то
АДх = ДЬ.
Далее,
||Д*||2 |Nb
||ДЬ||:
||Д^2 ||Ь»2 Nh ||ДЬ|2
||Д*||2
max
||ДЖ»2 ||АЖ|2 ||Х|2 ||АДЖ|2
||Аг||2
||Д^|2 ||Аг|
max
min
ЦАДхЬ |№ ЦАДхЬ
<
АЖ=0 ||АДж|2 |МЬ \Аж=° ||Д^12
(^min(^))-1 amax(A) = COnd2^
1
||Är||2 max —~—~—
в силу свойств сингулярных чисел матрицы [14]. Сравнение начала и конца этой выкладки приводит к неравенству (9), которое, как нетрудно понять, достижимо для некоторых х и Ах, или, что равносильно, для некоторых правых частей Ь и их возмущений АЬ. Естественно, что выписанные выкладки и полученная оценка имеют смысл лишь при ат1П(А) = 0.
2.2. Интервальные системы с точечной матрицей
Рассмотрим интервальную систему линейных алгебраических уравнений
Ах = b (10)
с точечной (неинтервальной) m х п-матрицей А и интервальным m-вектором b в правой части, причем m > п. Предположим также, что ее допусковое множество решений непусто, т. е.
БЫ(А, b) = { х е Rra | Ах е b } = 0
(при точечной матрице оно совпадает с объединенным множеством решений). Как можно быстро и несложно оценить размеры этого множества решений? Форма ответа на этот вопрос может быть различной, и мы будем строить его в виде оценки "типа радиуса" для множества решений. Более точно, мы получим оценку евклидовой нормы max ||х' — х||2 по всем х' е Etoi(А, b) для некоторой специальной фиксированной точки х е Etoi (А, b), в качестве которой берется
х = arg max Tol ( х, А, b)
хек™
— аргумент безусловного максимума распознающего функционала для системы (10), т. е. вектор оценки параметров зависимости (1) по соответствующим данным (2). Строго говоря, эта точка может определяться неединственным образом, но тогда х может быть любой из точек, на которых достигается рассматриваемый максимум.
Пусть х' — какая-то точка из допускового множества решений Stoi (А, b). Как оценить ||х' — х||2? Ясно, что х' их являются решениями систем линейных алгебраических уравнений с матрицей А и некоторыми правыми частями ' и соответственно из интервального вектора b. Если х = 0 и b = 0, то можно применить неравенство (9), рассмотрев возмущение решения х системы линейных алгебраических уравнений А х = . Тогда Ах = х' — х, Ab = Ь' — Ь, и мы имеем
< cond^-JM^,
||х||2 < ||Ь||2
откуда получается абсолютная оценка
||х' X12 < С<Э^2 А -||Х||2 ■ 11 Ь'и -Ь|12 . (11)
Шь
Точка X находится в результате максимизации распознающего функционала То1, для вычисления обусловленности соп^ А существуют хорошо разработанные стандартные процедуры, и потому для практической работы с неравенством (11) нужно как-то оценить значения | - | 2 и | | 2.
Наиболее просто мы поступим с ||Ь||2, приближенно взяв ||Ь||2 ~ ||Ь||2, т. е. как норму "наиболее представительной" точки Ь интервального вектора Ь, определенной в разд. 1:
\\b\\s
где b = 2(|midb + radb| + |midb — radb|).
Естественно, что при этом может допускаться некоторое огрубление, так что теперь вместо (11) будем писать
\\x' — x\\2 ^ cond2 A •
1|АЫ|;
\\ b \\ 2
'12)
Теперь необходимо определить приращение правой части ДЬ = b' — Ь. Его оценкой сверху является вектор радиусов rad b правой части, но она слишком груба. Чтобы получить более тонкую оценку ДЬ, наряду с системой (10) рассмотрим также систему линейных алгебраических уравнений
Ах = b, (13)
у которой правая часть b получена "сжатием" интервального вектора b, т. е. как
b := [ b + М, b — М ], (14)
где величина "сжатия" М определяется следующим образом:
М := max Tol (х, A, b) > 0.
Поскольку этот максимум М функционала Tol достигается при конкретном значении аргумента х, то
М = Tol (х, A, b)
in
г<т I
min -ч rad b,; —
Ki<m
mid bi — ^^ a-ijXj i=i
1
< min rad bi
Ki<m
Как следствие, b + M <b — M и интервальный вектор (14) — правильный, т.е. его концы не "перехлестывают" друг за друга.
Но из свойств распознающего функционала следует, что для интервальной системы линейных алгебраических уравнений (13) максимум распознающего функционала равен нулю, т. е.
max Tol (х, A, b) = 0.
хек"
В самом деле, величины rad bi входят слагаемыми во все выражения, стоящие в (5) под знаком минимума по г. Поэтому если одновременно увеличить или уменьшить все rad bi на одну и ту же величину при неизменных серединах mid bi компонент правой части, то и общее значение распознающего функционала увеличится или уменьшится на ту же самую величину. Иными словами, если взять константу С > 0 и интервальный т-вектор
е = ([—1,1],..., [—1,1])Т,
то для системы Ах = b + Се, в которой все компоненты правой части одновременно расширены на [—С, С], имеем
Tol (х, A, b + Се) = Tol (х, A, b) + С.
15)
2
Значит,
max Tol (x, A,b + Ce) = max Tol (x, A,b) + C (16)
(и это также верно для систем (3), (4) с произвольными матрицами А, а не только точечными). Аналогичным образом на допусковое множество решений и распознающий функционал действует и равномерное сужение вектора правой части. Если мы сузили все компоненты на одну и ту же величину М, то и максимум распознающего функционала новой интервальной системы уменьшится на М.
Итак, в силу свойств распознающего функционала допусковое множество решений Efoi(A,b) для системы (13) имеет пустую внутренность (такие множества часто называют "нетелесными"). Мы будем условно считать этот факт равносильным свойству "иметь нулевые размеры". Естественно, что это упрощающее допущение, так как в действительности допусковое множество решений при нулевом максимуме распознающего функционала может быть отнюдь не одноточечным. Тем не менее мы принимаем это упрощение, в пользу которого говорит также то обстоятельство, что ситуация с нулевым максимумом распознающего функционала неустойчива: соответствующее допус-ковое множество решений может стать пустым при сколь угодно малом возмущении данных.
Еще один факт, который касается вспомогательной системы (13) с суженной правой частью и следует из (15), (16), состоит в том, что определенная ранее точка X является аргументом максимума распознающего функционала также для (13):
X = arg max Tol (x, A, b).
xeR"
По этой причине точка b = AX лежит в b.
Из сказанного вытекает, что множество решений интервальной линейной системы A x = b получается из множества решений системы A x = b, которое имеет "прене-брежимые размеры" и для которого maxx€Rn Tol (x, A, b) = 0, с помощью расширения вектора b в каждой компоненте одновременно на [—М, М], где
М = max Tol (x,A,b).
xeR"
Как следствие, в неравенстве (12) мы можем взять ДЬ таким, что
|Д bi| = |b -6г|<М, i = 1, 2,... ,т, (17)
причем оценка сверху может достигаться, и потому
|| ДЬ|U = max Tol (x, A, b) (18)
x€Rn
в чебышевской векторной норме ||у||те := max1<i<m |г/j |.
Строго говоря, интервальный вектор b Э Ь может иметь ненулевые размеры, и тогда ||ДЬкак будто способна превзойти М. Но если мы хотим сделать неравенство (12) наиболее точным (достижимым), то его правая часть должна зануляться при max Tol = 0 и это предопределяет выбор (17) и (18).
Чтобы подставить найденное значение ||ДЬв (12), где фигурируют евклидовы нормы, воспользуемся неравенством эквивалентности векторных норм. Как известно (см. [11]), для любого вектора у
|ML <|Ы|2 < V^ |MU (19)
В целом получаем
\\х' — ж||2 ~ /п cond2A ■ llargmax Tol II ■--. (20
\\Ъ\\2
2.3. Общие интервальные системы
Рассмотрим, наконец, общие интервальные системы линейных алгебраических уравнений Ах = b с интервальной матрицей коэффициентов. Пусть, как и ранее,
х = argmax Tol (х, А, b)
хек™
— аргумент безусловного максимума распознающего функционала для системы Ах = b, т.е. вектор оценки параметров линейной зависимости (1) с данными (2). Пусть также х' — какая-то точка из непустого допускового множества решений Etoi(А, b). Как оценить \\х' — b\\2 ?
В силу известных свойств допускового множества решений (см., к примеру, [4]) его можно представить как
(А, b) = П { ж е R | Ах е b } = f| (A, b),
АеА АеА
т. е. как пересечение множеств решений отдельных ИСЛАУ, имеющих точечные матрицы А из А. Для каждой интервальной линейной системы Ах = b с А е А точки х' и х лежат во множестве решений Etoi(A,b), и мы имеем оценку (20). Следовательно, для допускового множества решений системы Ах = b, которое является пересечением множеств решений отдельных систем Ах = b с А е А, должно быть справедливо
II ' -II < • I Г Л Л II гр 1 II maxx€Rn To1 1
\\х — х\\2 ^ min < Jn cond2A • ||argmax Tol |L •--- > .
\\ W ~ АеА у II x€R" I'2 \\b\2 J
Правую часть в полученном неравенстве можно оценить приближенно, если внести минимум по A e А под фигурные скобки, оценив отдельные множители: число обусловленности cond2А, норму аргумента максимума \\ argmax^R" Tol \\ и максимум распознающего функционала max Tol . Полагая
I argmax Tol (х, A, b)y ~ const = У argmax Tol (х, А, b) У,
II хек™ 11 11 хек™ 11
будем иметь
min max Tol (х, A, b) ~ max Tol (х, A, b),
АеА xes.n xes.n
11 ж' — b| 2 ^ /п min cond2 А- У arg max Tol У- —a xeR—o-. (21)
11 112 ~ v АеА 2 11 хек™ 11 \\b\\ v y
Эта же оценка в силу неравенства эквивалентности норм (19) верна и для чебышевской нормы:
II ' -и < г ( • А л А У Till maxxeR™ To1 max \\x — ж Lo ^ v^ min cond2A • У argmax Tol IL •---,
x'estol(A,b)u ^ V АеА J И хек™ "2 \\b\\2
что оправдывает оценку 1УЕ (6). Если необходимо оценить относительные размеры до-пускового множества решений, выразив их в отношении к норме правой части системы уравнений, то, снова используя (19), получаем
шаХх',х"езш(А,ь) ||х' - х"||
оо ^ / \ шаХхекп То1
п ( шт cond2А V АеА 2
||Ь ||2
||х||~
Это дает оценку (8).
3. Два численных примера
В качестве конкретного примера применения меры вариабельности 1УЕ рассмотрим задачу восстановления линейной функции от двух переменных
Ь = х^! + х2а2
по данным трех измерений
№ а! 0,2
1 [98,100] [99,101] [190, 210]
2 [97, 99] [98,100] [200, 220]
3 [96, 98] [97, 99] [190, 210]
Для определения коэффициентов х! и х2 составляем интервальную 3 х 2-систему линейных алгебраических уравнений
([98,100] [99,101Д
[97, 99] [98,100] V [96, 98] [97, 99] )
0
^[190, 210]^ [200, 220] У[190, 210]у
(22)
Ее матрица имеет неполный ранг, поскольку содержит точечную матрицу ранга 1
98 99 98 99 98 99
(23)
Тем не менее интервальная матрица системы (22) не содержит линейно зависимых точечных столбцов, и потому согласно критерию И.А. Шарой [15, 16] допусковое множество решений — ограниченное. Оно изображено на рис. 5, построенном с помощью процедуры EqnTolR2 из пакета 1п^1п1псЯ2 [17]. Минимальное спектральное число обусловленности точечных матриц, содержащихся в матрице ИСЛАУ, равно 103.83 и достигается на матрице
100 99 97 100 96 99
В этом можно убедиться перебором всех "угловых" точечных матриц для интервальной матрицы системы (22), опираясь на квазивогнутость числа обусловленности.
Х\
Рис. 5. Допусковое множество решений системы уравнений (22) (узкий наклонный многоугольник) в сравнении с брусом х, построенным с помощью оценки ТУЕ
Численное решение задачи максимизации распознающего функционала Tol можно выполнить с помощью свободной программы tolsolvty, предназначенной для систем компьютерной математики Scilab, MATLAB, Octave.1 В этих системах встроенными процедурами удобно также считать число обусловленности матриц, которое требуется для наших расчетов. Используя установки точности, заданные в tolsolvty "по умолчанию", получаем max Tol = 1.9095, и он достигается в точке х = (5.1857 • 10-7, 2.0603)т. Соответственно,
IVE = V2 • 1.9095 • 103.83 • , = 1.6399.
V2002 + 2102 + 2002
Интервальная оболочка допускового множества решений для системы (22), т.е. его оптимальная внешняя интервальная оценка — это брус
(
[-0.9620, 3.0227]\ [-0.9320, 3.0127]/
который также находится процедурой EqnTolR2. Мы видим, что значение величины IVE неплохо согласуется с радиусом компонент оптимальной оценки множества решений, равной 1.9924 и 1.9724.
Так как в методе максимума согласования аргумент х = arg max^Rn Tol безусловного максимума распознающего функционала играет важнейшую роль и фактически относительно него строится наша оценка вариабельности IVE, то имеет смысл посмотреть на брус х с компонентами Х — IVE ,%i + IVE], i = 1, 2. Он тоже изображен на рис. 5, и существенная асимметрия его расположения относительно множества решений вызвана, конечно, специфическим положением центра — точки х, а также плохой
программа tolsolvty доступна на веб-сайте [1], раздел "Программное обеспечение".
обусловленностью задачи. С другими данными брус х оценивает допусковое множество решений существенно лучше.
Приведем пример противоположного свойства (в некотором смысле двойственный предыдущему примеру), когда линейная функция трех переменных
Ь = х1а1 + х2а2 + х2а3 восстанавливается по данным двух измерений, сведенным в таблицу.
(24)
№ а1 а2 з
1 [98,100] [97, 99] [96,98] [190, 210]
2 [99,101] [98,100] [97, 99] [200, 220]
Отметим, что в этих данных трехмерные брусы неопределенности замеров 1 и 2 существенно "налегают" друг на друга: их пересечением является брус с непустой внутренностью
( [99,100] \ [98, 99] [97, 98] \J200,210]/
размеры которого сравнимы с размерами исходных брусов данных. Для определения параметров зависимости приходим к интервальной системе линейных алгебраических уравнений
(хД ,г
11 [190,210]
/[98,100] [97, 99] [96, 98]' У [99,101] [98,100] [97, 99]
М 1 I ([190, 210] А
/ Х2 = V [200,220]/ 7 \хз/ 4 7
(25)
которая является неопределенной. Более того, она содержит матрицу неполного ранга 1, которая транспонирована к матрице (23).
Так как матрица этой системы является транспонированной матрицей системы (22), значение соп^А для нее то же самое.
Определение коэффициентов линейной функции (24) на основе понятия слабого согласования параметров и данных (см. разд. 1) представляется затруднительным (если вообще возможным), поскольку объединенное множество решений системы (25) является неограниченным. Тем не менее и в этих неблагоприятных условиях в силу критерия И.А. Шарой [15, 16] допусковое множество решений интервальной системы уравнений непусто и ограниченно (рис. 6, где изображена тонкая шестиугольная пластинка). Нахождение максимума распознающего функционала этой системы с помощью программы tolsolvty дает значение тахТо1 = 3.9698, который достигается в точке
Х = а^тахТо1 = (2.0603, 3 ■ 10-6, 2.1 ■ 10-6)т. Ее можно взять в качестве оценки коэффициентов. Тогда мера вариабельности
ТУЕ = л/2 ■ 3.9698 ■ 103.83 ■
| Х У 2
л/2002 + 2102
4.1413.
4
3
2
хз 1
0
-1
-2.
Х\
Рис. 6. Допусковое множество решений системы уравнений (25)
Интервальная оболочка допускового множества решений для системы (25), т.е. его оптимальная внешняя интервальная оценка — это брус
который тоже находится процедурой EqnTolR3. Радиусы компонент оптимальной внешней оценки допускового множества решений равны 3.0024, 3.0329, 3.0510, что также не сильно отличается от значения 1УЕ. Рассмотренный пример показывает работоспособность оценки даже в случае недоопределенных ИСЛАУ. Но строгое исследование и обоснование этого факта, как мы уже отмечали, еще ожидают своего продолжения.
4. Обсуждение
Мы уже показывали, что мера вариабельности 1УЕ имеет интуитивно ясный и даже наглядный смысл. Ее интерпретация, данная в разд. 1, частично объясняет также, почему 1УЕ работает для недоопределенных ИСЛАУ, когда измерений меньше, чем восстанавливаемых параметров. Но на меру вариабельности 1УЕ можно взглянуть еще с одной стороны.
В традиционной теоретико-вероятностной статистике большую роль играет явление коллинеарности данных (иногда называемое мультиколлинеарностью) — наличие линейной зависимости между входными (предикторными) переменными регрессионной модели. Обычно считается, что к переменных рассматриваемой модели являются коллинеарными, если представляющие их векторы лежат в линейном пространстве размерности меньше к [18], так что некоторые из этих векторов являются линейными
комбинациями других. На практике подобная "точная коллинеарность" данных встречается нечасто, но ее наличие совершенно не обязательно для существования реальной "проблемы коллинеарности", когда векторы данных "почти линейно зависимы". Тогда оценки параметров неустойчивы, что ведет к увеличенной статистической неопределенности — росту дисперсии оценок.
Согласно современным воззрениям коллинеарность данных наиболее адекватно описывается числом обусловленности матрицы из этих данных (см., к примеру, [18], гл. 3). В этом смысле введенная нами мера IVE является фактически модифицированной мерой коллинеарности данных, скорректированной с учетом относительного масштаба решения и правой части системы, а также конкретного значения максимума распознающего функционала. Минимум по всем cond2A для A G А берется в силу специфики сильного согласования параметров и данных, и он хорошо отражает регуляри-зующую роль допускового множества решений (см. подробное изложение этого вопроса в работе [19]).2
При такой трактовке предложенная в работе мера вариабельности IVE имеет смысл даже при отрицательном значении максимума распознающего функционала max Tol, когда допусковое множество решений пусто и параметров линейной функции (1), сильно согласующихся с данными, не существует. Абсолютное значение IVE все равно показывает, с точностью до некоторого масштабирующего коэффициента, меру коллинеарности данных (меру их обусловленности), а отрицательный знак сигнализирует о статусе решения задачи, т. е. о том, что найденный вектор параметров не является сильно согласующимся с данными, но лишь обеспечивает наилучшее возможное в условиях задачи приближение.
В целом предложенная в этой работе мера вариабельности IVE находится довольно просто, так как после решения задачи восстановления линейной зависимости с помощью метода максимума согласования мы уже знаем max Tol и arg max Tol. Наиболее сложным при вычислении IVE является нахождение минимума чисел обусловленности точечных матриц из заданной интервальной матрицы. В рассмотренных выше примерах мы использовали перебор угловых точечных матриц, что в общем случае, конечно же, непрактично и неэффективно.
Если интервальная матрица "достаточно узка" и не сильно отличается от точечной матрицы, то приближенно можно положить
min cond2(A) « cond2(midA). (26)
AeA
Но в общем случае этот рецепт может работать плохо, так как левая и правая части приближенного равенства (26) будут отличаться весьма сильно. Так, в примерах с системами (22) и (25) из разд. 3 обусловленность средней матрицы равна 2.38 • 104, и потому, воспользовавшись упрощенной формулой (26), мы изменим оценки размеров множества решений более чем в 200 раз. В общем случае для более точного оценивания min cond2(A) можно порекомендовать перебор некоторых наиболее хорошо обусловленных (по тем или иным внешним признакам) угловых матриц из А.
Для более серьезных задач можно организовать поиск минимума числа обусловленности каким-либо методом прямого поиска по множеству всех вершин интервальной
2Вообще говоря, согласно работе [19] одним из эффективных способов борьбы с явлением кол-
линеарности данных является их обынтерваливание и переход к сильному согласованию данных и
параметров. Этот прием равносилен развитой в [19] интервальной регуляризации систем линейных уравнений.
матрицы. Нахождение градиента целевой функции в этой задаче затруднено ее сложным характером, но вот вычисление значений целевой функции выполняется достаточно эффективно. Как следствие, наиболее подходящими инструментами решения этой задачи будут методы нулевого порядка, которые используют только значения целевой функции, т. е. числа обусловленности матриц.
Можно предложить, к примеру, следующий эвристический алгоритм оценивания min cond2(A). Задавшись некоторым натуральным числом N, организуем N шагов случайного блуждания по множеству угловых матриц интервальной матрицы А. Будем вычислять число обусловленности полученной случайной угловой матрицы Ä и ее противоположной относительно диагонали бруса А, следуя идее так называемых диагональных методов глобальной оптимизации [20]. Иными словами, мы вычисляем cond2(Ä) и еще число обусловленности матрицы, у которой элементами взяты противоположные к А концы интервалов А, что делается для усиления "репрезентативности" нашего случайного блуждания. Наименьшее из значений, полученных после таких N шагов, берется в качестве оценки для min cond2(A). Ясно, что конкретное значение числа шагов N нужно брать в зависимости от желаемой точности оценивания min cond2(A) и тех ресурсов, которые мы готовы на него потратить. Если нужно получить оценку для min cond2(A) как можно быстрее, то и N может быть невелико. При небольших N возможно также применение к результату какой-то корректирующей процедуры (например, умножение на поправочный коэффициент из открытого интервала ]0,1[).
Описанный выше простой рандомизированный алгоритм работает на удивление хорошо. Так, для рассмотренных в разд. 3 примеров интервальных систем (22) и (25) (матрицы которых совпадают с точностью до транспонирования) уже при N = 4, 5 на длительной серии прогонов (несколько десятков и более) вычисленная с помощью нашего алгоритма оценка для min cond2(A) не отличалась от минимального значения более чем в 2 раза, да и это случалось довольно редко. Таким образом, вместо вычисления обсловленности 26 = 64 штук "угловых матриц" мы получали практически приемлемую оценку не более чем за 10 вычислений.
Благодарности. Автор благодарен своим коллегам А.Н. Баженову, С.И. Жилину и А.А. Нелюбину за стимулирующие обсуждения и ценные замечания по статье.
Список литературы / References
[1] Интервальный анализ и его приложения. Тематический веб-сайт http://www.nsc.ru/ interval
Interval Analysis and its Applications. A thematic web-site http://www.nsc.ru/interval
[2] Ж!олен Л., Кифер М., Дидри О., Вальтер Э. Прикладной интервальный анализ. М.; Ижевск: Изд-во "РХД", 2007. 467 с.
Jaulin, L., Kieffer, M., Didrit, O., Walter, E. Applied interval analysis. London: Springer, 2001. 379 p.
[3] Moore, R.E., Kearfott, R.B., Cloud, M.J. Introduction to interval analysis. Philadelphia: SIAM, 2009. 223 p.
[4] Шарый С.П. Конечномерный интервальный анализ. Новосибирск: ИВТ СО РАН & XYZ, 2019. 630 с. http://www.nsc.ru/interval/?page=Library/InteBooks
Shary, S.P. Finite-Dimensional Interval Analysis. Novosibirsk: IVT SO RAN & XYZ, 2019. 630 p. http://www.nsc.ru/interval/?page=Library/InteBooks (In Russ.)
[5] Крамер Г. Математические методы статистики. М.: Мир, 1975. 648 с. Cramer, H. Mathematical methods of statistics. Princeton Univ. Press, 1999. 575 p.
[6] Ш!арый С.П. Разрешимость интервальных линейных уравнений и анализ данных с неопределенностями // Автоматика и телемеханика. 2012. № 2. С. 111-125.
Shary, S.P. Solvability of interval linear equations and data analysis under uncertainty // Automation and Remote Control. 2012. Vol. 73, No. 2. P. 310-322.
[7] ШШарый С.П. Сильная согласованность в задаче восстановления зависимостей при интервальной неопределенности данных // Вычисл. технологии. 2017. Т. 22, № 2. С. 150-172. Shary, S.P. Strong compatibility in data fitting problem under interval data uncertainty // Comput. Technologies. 2017. Vol. 22, № 2. P. 150-172. (In Russ.)
[8] ШШарый С.П. Метод максимума согласования для восстановления зависимостей по данным с интервальной неопределенностью // Изв. Акад. наук. Теория и системы управления. 2017. № 6. C. 3-19.
Shary, S.P. Maximum compatibility method for data fitting under interval uncertainty // J. of Comput. and Syst. Sci. Intern. 2017. Vol. 56, iss. 6. P. 897-913.
[9] ШШарый С.П., ШШарая И.А. Распознавание разрешимости интервальных уравнений и его приложения к анализу данных // Вычисл. технологии. 2013. Т. 18, № 3. С. 80-109. Shary, S.P., Sharaya, I.A. Recognizing solvability of interval equations and its application to data analysis // Comput. Technologies. 2013. Vol. 18, No. 3. P. 80-109. (In Russ.)
[10] Kearfott, R.B., Nakao, M., Neumaier, A., Rump, S., Shary, S.P., van Hentenryck, P. Standardized notation in interval analysis // Comput. Technologies. 2010. Vol. 15, No. 1. P. 7-13.
[11] Голуб Дж., ван Лоун Ч. Матричные вычисления. М.: Мир, 1999. 548 с.
Golub, G.H., Van Loan, Ch.F. Matrix Computations. Baltimore: The Johns Hopkins Univ. Press, 1996. 728 p.
[12] Уоткинс Д. Основы матричных вычислений. М.: Бином. Лаб. знаний, 2009. 664 с. Watkins, D.S. Fundamentals of matrix computations. New York: Wiley-Intersci., 2002. 618 p.
[13] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Бином. Лаб. знаний, 2003. 632 с.
Bakhvalov, N.S., Zhidkov, N.P., Kobelkov, G.M. Numerical methods. Moscow: Binom. Lab. Znaniy, 2003. 632 p. (In Russ.)
[14] Годунов С.К. Современные аспекты линейной алгебры. Новосибирск: Научная книга, 1997. 388 с.
Godunov, S.K. Modern aspects of linear algebra. Novosibirsk: Nauchnaya Kniga, 1997. 388 p. (In Russ.)
[15] ШШарая И.А. Ограничено ли допустимое множество решений интервальной системы? // Вычисл. технологии. 2004. Т. 9, № 3. С. 108-112.
Sharaya, I.A. Is tolerable solution set bounded? // Comput. Technologies. 2004. Vol. 9, No. 3. P. 108-112. (In Russ.)
[16] Sharaya, I.A. On unbounded tolerable solution sets // Reliable Computing. 2005. Vol. 11, iss. 5. P. 425-432.
[17] ШШарая И.А. Пакет IntLinIncR2 для визуализации множеств решений интервальных линейных систем с двумя неизвестными: Программное обеспечение. Новосибирск, 2014. Адрес доступа: http://www.nsc.ru/interval/sharaya
Sharaya, I.A. IntLinInc2D, a software package for visualization of solution sets for interval linear 2D systems. Novosibirsk, 2014. Available at: http://www.nsc.ru/interval/sharaya
[18] Belsley, D.A., Kuh, E., Welsch, R.E. Regerssion diagnostics. Hoboken, N.Y.: Wiley-Intersci., 1980. 292 p.
[19] Shary, S.P. Interval regularization for imprecise linear algebraic equations.
Available at: https://arxiv.org/abs/1810.01481 (Deposited in arXiv.org on 27 Sept. 2018, No. arXiv:1810.01481. 21 p.)
[20] Сергеев Я.Д., Квасов Д.Е. Диагональные методы глобальной оптимизации. М.: Физ-матлит, 2008. 352 с.
Sergeyev, Ya.D., Kvasov, D.E. Deterministic global optimization. An introduction to the diagonal approach. New York: Springer, 2017. 136 p.
Поступила в 'редакцию 9 апреля 2019 г., с доработки — 6 мая 2019 г.
On a variability measure for estimates of parameters in the statistics of interval data
Shary, Sergey P.
Institute of Computational Technologies SB RAS, Novosibirsk, 630090, Russia Novosibirsk State University, Novosibirsk, 630090, Russia E-mail: shary@ict.nsc.ru
The article presents a possible construction for quantitative variability measure of parameter estimates in the data fitting problem under interval uncertainty. It is designed to show the degree of variability and ambiguity of the estimate. The problem is important since the answer to the data fitting problem may be not unique under the uncertainty of interval data. Otherwise, the problem may have a whole set of such feasible answers, and an estimate of the size of this set can serve as a variability measure of the solution.
So-called maximum compatibility method is considered as a method for evaluating the parameters of the function. An estimate of the parameters is taken as the maximum point of a special "recognizing functional", i.e. a certain function that characterizes compatibility of the estimate with interval empirical data. Then the variability measure is simply computed as the product of several values found in the solution of the data fitting problem. These values are maximum of the recognizing functional, minimum of the condition number of the data matrix and the ratio of the norm of the parameter estimate to the norm of the vector of function values.
A derivation of the new variability measure is given, and its application is illustrated with specific examples. In conclusion, the article discusses the motivation and interpretation of the new variability measure from the point of view of other known data characteristics.
Keywords: data fitting problem, interval data uncertainty, strong compatibility of parameters and data, variability measure.
Cite: Shary, S.P. On a variability measure for estimates of parameters in the statistics of interval data // Computational Technologies. 2019. Vol. 24, No. 5. P. 90-108. (In Russ.) DOI: 10.25743/ICT.2019.24.5.008.
Acknowledgements. The author is grateful to his colleagues Alexander N. Bazhenov, Sergei I. Zhilin and Alexei A. Nelyubin for stimulating discussions and valuable comments on the article.
Received on April 9, 2019 Received in revised form on May 6, 2019
© ICT SB RAS, 2019
