ISSNG868-5886
НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2GG2, том 12, № 2, c. 3G-49
ОРИГИНАЛЬНЫЕ СТАТЬИ =
УДК519.245: 621.391.1
© А. Л. Буляница, В. Е. Курочкин
ИССЛЕДОВАНИЕ СВОЙСТВ И ПРОГРАММНО-АППАРАТНАЯ РЕАЛИЗАЦИЯ АЛГОРИТМА СТОХАСТИЧЕСКОЙ АППРОКСИМАЦИИ В МОДИФИКАЦИИ Я. З. ЦЫПКИНА
В работе рассматривается модификация Я.З. Цыпкина алгоритма стохастической аппроксимации Роббинса—Монро, позволяющая оценить постоянный сигнал в условиях аддитивной помехи с априорно неизвестным законом распределения, включая статистические выбросы. Проанализированы доказательства свойства несмещенности формируемой оценки и исследована эффективность оценивания. Рассмотрена экономичная аппаратная реализация алгоритма, а также его программная реализация, позволяющая, помимо оценивания информативного сигнала, осуществлять подбор параметров алгоритма и начальное приближение оценки, а также выявлять разладку в последовательности измерений.
ВВЕДЕНИЕ
В кинетических анализаторах различных субстанций, в том числе на основе хемо- и биосенсоров, выходной сигнал детектора характеризуется известной детерминированной формой с неизвестными случайными информативными параметрами, которые подлежат оцениванию. Форма сигнала определяется видом массопереноса, а также типом и, как правило, геометрией конкретного сенсора.
Существует широкий спектр сенсоров, работающих в условиях, когда скорости химических процессов ограничены скоростями диффузии реагентов. При отсутствии осложняющих взаимодействий (структурно-механические барьеры, расслоение и пробой фаз, накопление вещества на границе раздела) и учете объемных и поверхностных реакций первого порядка для оценки количества экстрагированного продукта используют решение уравнения диффузии в параболической форме для полубесконечных фаз [1-6]. При этом исчерпываемое вещество должно быть в избытке. Тогда зависимость количества экстрагируемого продукта q от времени примет вид линейного тренда первого порядка [2, 5] в специальной системе координат: q = а0 + Ь0^, где а0, Ь0 — информативные оцениваемые параметры.
В условиях, при которых извлекающее вещество практически не расходуется [2, 3, 5], временная зависимость q(t) также носит характер линейного тренда первого порядка
q = а0 + V- (1)
В зависимости от типа чувствительного эле-
мента информативным откликом может быть различная физическая величина — тепло, оптическая плотность и т.д. [7], которая, как правило, преобразуется в электрический сигнал. Например, величина электрического сигнала в прямых иммуносенсорах, основанных на взаимодействии антиген—антитело на границе жидкой и твердой фаз
[8], имеет вид (1). При этом Ь0 ~ 2С0>/В/ п , где С0 — концентрация искомых антител (антигенов). При использовании оптического хемосенсора величина оптической плотности Дd, измеренная в отраженном свете, также подчиняется закону (1) и, как показано в [9], информативный параметр имеет вид Ь0 ~ 2еКхС0^В/п . Здесь В — коэффициент диффузии продукта реакции, а £Кх — молярный коэффициент поглощения в фазе сенсора.
Таким образом, информативный отклик хемо-и биосенсора можно представить в общем виде как
q = а0 + Ь^а .
В системе координат ^“, q) последнее уравнение сводится к форме (1), т. е. является прямой (линейным трендом первого порядка).
С помощью величин а0, Ь0 можно однозначно определить искомую концентрацию С0. Однако эти величины зависят от способа измерения, по-разному влияют на оценку концентрации, и соответствующая методика оценивания указанных параметров нуждается в оптимизации, например по критерию минимальной ошибки оценивания концентрации [10].
Еще одним типом задач, в которых необходимо оценивание коэффициентов в зависимости типа
(1), является анализ калибровочных (градуировочных) зависимостей. В (1) коэффициент b0 однозначно связан с концентрацией анализируемой компоненты С0, так что калибровочная (градуировочная) характеристика прибора также может иметь схожий вид
bo =У + А(Со)5. (2)
Таким образом, и на этапе определения концентрации С0 по параметру b0 выходного сигнала чувствительного элемента также приходится решать аналогичную задачу. Кроме того, зависимость типа (2) может устанавливать связь между информативным параметром выходного сигнала (b0) и некоторыми параметрами, определяющими условия эксперимента (геометрия и размеры проточной аналитической системы, скорость транспортировки анализируемого вещества и т.д.).
Например, в [11, 12] рассмотрена проточная аналитическая система коаксиальной конструкции с амперометрическим чувствительным элементом. В условиях, когда анализируемые концентрации С0 не очень велики, некоторые из градуировочных характеристик системы имеют вид
b0 = aC0 и b0 = ^ < v >s . (3)
Здесь < v > — скорость транспортировки анализируемого вещества, показатель S в зависимости от диапазона скорости принимает значения 1/3, 1/2 или 1.
Зависимость информативного параметра b0 от скорости для детекторов типа отражающей стенки ("wall-jet") [13-15] также имеет вид (3) при 8 = 3/4.
Таким образом, на различных стадиях проведения химического или биологического экспресс-анализа необходимо решать сходные задачи и исследовать однотипные объекты независимо от того, рассматривается кинетическая кривая отклика чувствительного элемента либо калибровочная (градуировочная) характеристика прибора.
Основными задачами являются: анализ информативной составляющей и помех в сигналах типа линейного тренда, выбор алгоритмов оценивания информативных параметров сигналов типа линейного тренда нулевого и первого порядка, оптимизация алгоритмов по критериям как минимума ошибки оценивания информативного параметра, так и по критерию экономии вычислительных ресурсов алгоритма.
ИНТЕРВАЛЬНЫЕ И РЕКУРРЕНТНЫЕ АЛГОРИТМЫ ОБРАБОТКИ ИНФОРМАЦИИ
Как говорилось ранее, в сенсорных приборах химического и биологического (иммунного) экс-
пресс-анализа выходной сигнал с чувствительного элемента стремятся привести к одному из следующих стандартных типов:
а) сигнал постоянного уровня (линейный тренд нулевого порядка);
б) линейный тренд первого порядка;
в) линейный тренд первого порядка с насыщением;
г) сигнал типа "трапеция", т. е. совокупность линейных трендов первого и нулевого порядков.
Так как в большинстве случаев информативным является параметр положения линейного тренда первого порядка, однозначно связанный с концентрацией анализируемого вещества, то актуальной становится задача оценивания параметра положения (тангенса угла наклона) линейного тренда в присутствии помех различного происхождения.
В то же время переход к первому разностному сигналу сводит описанные выше сигналы четырех типов к совокупности сигналов постоянного уровня на фоне симметричной (и, следовательно, центрированной) помехи.
Решение поставленной задачи — оценивание постоянного сигнала с* по измерениям хь содержащим аддитивную помеху § с априорно неизвестным законом распределения, т. е.
х. = с * +§, / = 1,2,... — может осуществляться
различными методами. При этом следует учесть, что помеха симметрична (и центрирована) только в том случае, если был осуществлен переход к первой разности выходного сигнала чувствительного элемента сенсора.
Методы оценивания можно условно разделить на две группы: нерекуррентные (интервальные), к которым, в частности, относятся метод максимального правдоподобия, или метод наименьших квадратов, и рекуррентные (например, стохастической аппроксимации) [16]. Существует и рекуррентная модификация метода наименьших квадратов.
При использовании нерекуррентных методов для получения оценки сп+1 величины с* на (п+1)-м шаге оценивания необходимо произвести пересчет с использованием всех накопленных измерений х1,х2,...,хп,хп+1 и по крайней мере некоторых из предыдущих оценок с1,с2,...,сп, что предъявляет повышенные требования к объему оперативной памяти устройства оценивания и к его быстродействию. Однако в целом характерные времена химических и тем более биологических (иммунных) реакций достаточно велики (не менее
0.1-1 с, а для биологических реакций — до нескольких минут). Поэтому проблема быстродействия устройства оценивания достаточно легко решается с использованием практически любой современной элементной базы.
В рекуррентных методах оценивания оценка сп+1 пересчитывается из предыдущей оценки сп на основе новой информации, т. е. текущего измерения хп+1. Таким образом, сп+1 = /(сп, хп+1) или, что то же самое, сп+1 = / * (сп, Е,п+1) . Появившиеся в последние годы микроконтроллеры, обладающие высоким быстродействием и низкой стоимостью, с небольшим объемом оперативного запоминающего устройства (ОЗУ) идеально подходят для реализации на их основе рекуррентных оценщиков.
Алгоритм стохастической аппроксимации был предложен в работе [17] для решения задачи поиска единственного корня уравнения регрессии при выполнении следующих предположений: у является случайной величиной с функцией распределения ¥ (у|с), зависящей от некоторого параметра с. Если
существует функция регрессии М (с) = | у d¥ (у|с)
при некотором вещественном а, уравнение М(с) = а имеет единственный корень с = в , то алгоритм поиска в имеет вид
■'п+1
= Сп + ап^(Сп )•
(4)
«п <^
п=1
п=1
ного тренда первого порядка [24, 25].
Особый интерес представляет модификация Я.З. Цыпкина [23] в следующей форме:
в
Сп+1 = Сп-^(Сп - Хп+1І
п
(5)
Здесь функция V — релейная функция с зоной нечувствительности 2Д. В простейшем случае Д = 0, и (5) сводится к виду сигнатурного алгоритма, допускающему простую аппаратную реализацию
Сп+1 = Сп -— ®^§п(сп - Хп+1І
(6)
С = X,.
Здесь сп — оценка величины в на п-м шаге оценивания, у(сп) — наблюдаемое значение величины у при с = сп , функция У(сп) = а- у(сп). Формула (4) называется алгоритмом стохастической аппроксимации Роббинса—Монро [17]. При выполнении условий, налагаемых на числовую последовательность {ап }:
(т. е. сам числовой ряд расходится, а ряд из квадратов сходится), была доказана несмещенность и состоятельность получаемой оценки в.
Состоятельность формулируется как Уе > 0 : Р(|сп - в < е) ^ 1, п ^ +^ (Р — вероятность указанного в скобках события). Несмещенность определяется как М{сп}= в , где М{..} обозначает математическое ожидание. В ряде случаев несмещенность также трактуют как асимптотическое свойство, выполняемое при п ^ +го .
В дальнейшем процедура Роббинса—Монро (4) получила свое развитие во многих работах, из которых можно выделить лишь основные направления, в которых проходила модификация.
Во-первых, предлагались различные подходы к выбору последовательности ап [18, 19]. Во-вторых, различным образом выбиралась функция V [20-23]. В-третьих, были предложены модификации, позволяющие оценивать величину линей-
Благодаря алгоритмической простоте, алгоритмы (5) и (6) получили широкое распространение для решения задач оценивания постоянного сигнала в присутствии аддитивных помех с априорно неизвестным законом распределения, в том числе в случае помехи типа "выброс".
СТРУКТУРА И ОСНОВНЫЕ СВОЙСТВА АЛГОРИТМА СТОХАСТИЧЕСКОЙ АППРОКСИМАЦИИ В ФОРМЕ Я.З. ЦЫПКИНА
Алгоритм стохастической аппроксимации Роббинса—Монро хорошо известен еще с пятидесятых годов, но его интенсивное изучение и модификация в нашей стране приходятся на семидесятые годы. Главные успехи в этом направлении связаны с именем академика Я.З. Цыпкина [23, 26] и с представителями его школы А.А. Бедельбаевой [27], Т.П. Красулиной [28] и др. Результатом исследований явилась модификация алгоритма Роббинса—Монро в форме, эквивалентной (5):
п+1
= Сп -—/п ■Х¥(Сп - Хп+1) .
(7)
При этом
У( г) =
- 1, для г < -Д,
0, для |г| < Д,
1, для г > Д;
хп+1 = с* + §п+1 — измерение; §п+1 — помеха; сп, сп+1 — оценки величины с* на п-м и (п+1)-м шагах оценивания; в — параметр алгоритма; 2Д — величина зоны нечувствительности [27]. В случае Д = 0 получается сигнатурный алгоритм Я.З. Цыпкина (6), т. к. ¥(г) = sign(z).
Важнейшим, с нашей точки зрения, теоретическим результатом явилось доказательство свойства минимаксности оценки алгоритма (7) в услови-
С1 = Х1
п
«п =“
ях априорной неопределенности данных о помехе: наихудшим (по дисперсии ошибки) являлось симметричное распределение Лапласа с плотностью распределения вероятности (ПРВ)
/(х) = ^ехр 2а
(
х
и, исходя из этого факта,
осуществлялся подбор функции V, что обеспечивало минимум дисперсии ошибки [23, 26]. Из всех возможных распределения вероятностей помехи _Дх) выбирались такие, для которых выполнялось условие /(0) = 1/(2а) > 0, т. е. невырожденные распределения. Таким образом, уже при достижении указанного результата (свойства минимаксно -сти) речь шла не о полной априорной неопределенности, а о более узком классе — классе невырожденных распределений.
К основным асимптотическим свойствам оценки, формируемой с помощью алгоритма (7), относятся состоятельность, несмещенность и эффективность.
Исследуя алгоритм (7), А.А. Бедельбаевой [27] удалось сформулировать условия, при которых оценка сп является несмещенной и состоятельной. Эти условия накладывали существенные ограничения на выбор параметров алгоритма (7). Приведенная в [27] зависимость дисперсии ошибки оценивания от плотности (р) и функции (¥) распределения помехи, а также от параметра Д, позволила сформулировать эти условия в форме:
в > втт , Атт =[4Р(Д)]"
(8)
Была проанализирована относительная эффективность алгоритма (7) по сравнению с методом наименьших квадратов (МНК) для некоторых видов помех (нормальной, равномерной, лаплассов-ской, а также модели типа "выбросов" [27]).
Основные выводы [27] широко известны: эффективность оценивания с помощью МНК существенно ухудшается по мере отклонения закона распределения помехи от нормального, оценка МНК не является робастной (устойчивой к выбросам), что выражается в резком увеличении дисперсии ошибки оценивания при наличии выбросов, представленных в форме смеси Тьюки. Под смесью Тьюки понимается случайная величина £, формируемая в виде суммы двух слагаемых вида
I = (1 -е^ (0,ст2) + е N(0, /л2а2).
Здесь N — символ нормального (гауссовского) распределения, величина е определяет долю выбросов (е обычно берут в пределах от 0.01 до 0.10) и параметр л много больше единицы (например, Л = 3 или даже /л = 10) [27, 30]. При этом эффективность самого алгоритма (7) существенно зависит от выбора параметров алгоритма (в, Д).
Минимальная дисперсия ошибки оценивания по параметру в достигается при условии
в0р, = 2вт
(9)
Кроме того, помеха должна быть симметрична.
Следует отметить, что для доказательства асимптотических свойств оценки была использована информация о плотности распределения вероятности помехи на границе зоны нечувствительности.
Полученные результаты [23, 27] относительно важнейших асимптотических свойств оценки — несмещенности и состоятельности — базируются на нетривиальной информации о величине плотности распределения помехи на границе зоны нечувствительности (в общем случае даже не в нуле). Кроме того, требование симметричности помехи, хотя и принимается многими исследователями априорно (например, так называемые "аксиомы Маликова" [29]), также носит ограничительный характер.
Таким образом, доказательство несмещенности и состоятельности оценки сп было проведено в условиях существенно отличных от анонсированной априорной неопределенности данных о помехе, даже если ограничиваться только невырожденными помехами. По этой причине подобный метод доказательства не представляется полностью убедительным.
Оптимум по параметру Д также существенным образом зависит от закона распределения помехи. Следовательно, подбор параметров алгоритма (7) может являться способом регулирования эффективности оценивания. Однако при отсутствии априорной информации о помехе заранее выбрать оптимальные параметры алгоритма и даже просто гарантированно выполнить условие (8) несмещенности и состоятельности оценки достаточно сложно. Исходя из зависимости дисперсии ошибки оценивания [27] от параметров (в, Д), можно сделать вывод о том, что дисперсия очень слабо возрастает вблизи своего минимума. Таким образом, использование точной процедуры оценивания и подбора параметра в ~ вор не требуется.
Значительно позднее работы [27] был рассмотрен интересный частный случай треугольной (симпсоновской) помехи [31]. Этот случай выделен особо, т. к. он вполне может быть реализован во многих задачах химии, биохимии и биофизики, когда оценивается первая разность выходного сигнала чувствительного элемента. При этом помеха исходного информативного сигнала удовлетворяет равномерному закону распределения.
Это тот случай, когда наличие в измерительной схеме аналого-цифровых элементов обусловливает квантование исходного измерения по уровню. В том случае, если указанный эффект является
главным источником погрешности, случайная величина (помеха исходного измерения) распределена равномерно [32] в диапазоне ±5/2, где 5 — единица младшего разряда. Тогда разностная помеха л удовлетворяет треугольному распределению Симпсона с размахом 8/Д1 при Д (— величине кванта времени. Разностная помеха будет симметричной относительно нуля и центрированной как интегральная свертка исходной помехи.
Если исходные информативные (измеряемые) сигналы представляют собой линейные тренды с информативным параметром — величиной тангенса угла наклона С*, пропорциональным концентрации определяемой субстанции [2, 33], то задачу оценивания С* можно свести к задаче оценивания сигнала постоянного уровня путем перехода к первому разностному сигналу. В этом случае для оценивания С* можно использовать алгоритм стохастической аппроксимации Роббинса—Монро в модификации Я.З. Цыпкина (7).
Модификация этого алгоритма, разработанная Дупачем [24] для отслеживания линейного тренда, задачи оценивания С* не решает, поскольку отслеживает положение точки на тренде, определяемое, помимо С*, и параметром сдвига.
В [27] показано, что минимальная дисперсия ошибки оценивания С* — о2ге1 обеспечивается при ворх = 1/(2 Р(Д)) и определяется выражением
^(Д) = ¥(-Д)/(2р2 (Д)).
(10)
как
Д* = а^тт{аг2е1(Д)|Д > 0} .
(11)
Как показано в [27], (11) может иметь как тривиальное (Д = 0), так и нетривиальное (Д > 0) решение в зависимости от закона распределения разностной помехи л-
Докажем, что в случае треугольной (симпсо-новской) помехи и при условии оптимального выбора параметра алгоритма в, согласно (9), дисперсия ошибки оценивания не зависит от величины зоны нечувствительности алгоритма, т. е. что
сгг2е1 (Д) = соп81;(Д). Тем самым будет доказана невозможность воздействия на величину о2ге1 с помощью параметра алгоритма Д, и бессмысленность постановки задачи (11), и тем самым возможность произвольного выбора этого параметра.
Простая проверка позволяет убедиться в выполнении указанного свойства: поскольку помеха удовлетворяет распределению Симпсона, то для положительных Д выполняется
р(Д) = (С -Д)/ С2 и ¥(-Д) = (С -Д)2/(2С2).
(12)
Здесь С — полуразмах помехи.
Подстановкой в (10) получаем, что
стг2е1(Д) = С2 /4 , т. е. константа не зависит от величины Д. Для обеспечения условия р(Д) > 0 требуется только Д < С.
Докажем, что независимость дисперсии ошибки от величины зоны нечувствительности выполняется только в случае симпсоновской помехи.
Случайная помеха характеризуется плотностью распределения вероятностей р(х) и функцией распределения ¥(х). Условие (10) а2(х) = соп81;(х) после дифференцирования по переменной х сводится к утверждению
Г ¥ (-х) ^
дх
= 0.
(13)
Очевидно, что (10) имеет смысл только при р(Д) > 0.
В общем случае варьирование величины Д в алгоритме (7) позволяет в соответствии с (10) влиять на величину дисперсии ошибки оценивания информативного параметра С* — о2ге1. Поэтому имеет смысл задача оптимального поиска величины параметра алгоритма (7) — Д*, формулируемая
Полагаем х > 0 и обозначим у = ¥(-х). Тогда р(х) = -у', р'(х) = -у" и (13) сводится к обыкновенному дифференциальному уравнению второго порядка
(У У - 2 • у • у ' = 0.
Особое решение у'(х) = 0 нарушает условие р > 0. Следовательно, у'/ у = 2 • у"/ у' и у' = -С1 • у[у . Произвольная константа С1 отлична от нуля. (В противном случае, получается уже отвергнутое ранее тривиальное решение). Повторное интегрирование приводит к выражению
у(х) = (С2 - С • х)2.
Последнее выражение у(х) представляет собой ПРВ треугольной (симпсоновской) помехи при выполнении дополнительного условия нормировки на единицу, благодаря которому произвольные константы интегрирования С1, С2 будут связаны друг с другом.
При х > 0 р(х) = 2С1 (С2 - С1 х); из-за неотрицательности р(х) диапазон допустимых положительных значений х ограничен величиной С2 / С1 . Условие нормировки сводится к форме:
С2/ С1
|р(х)6х = 1/2, или С2 = 1/л/2 . Тогда полу-
0
размах помехи С есть С2 / С1 = 1/(л[2С1), а у(х) представляется в следующем виде
у(х) = (1Д/2 - С1 х)2 = С121^л/2С1) - х). Так как
С,2 = 1/(2С2), то у(х) = ——х—, что полностью 1 2С2
совпадает с выражением (12).
Видно, что никакая другая функция у(х) не удовлетворяет соответствующему дифференциальному уравнению. Следовательно, распределение Симпсона является единственным распределением, обладающим свойством независимости дисперсии оценки, формируемой алгоритмом (7) при оптимальном выборе параметра в, от параметра Д, определяющего величину зоны нечувствительности.
Основные выводы.
1. Алгоритм оценивания С* — параметра положения (тангенса угла наклона) линейного тренда первого порядка может быть реализован в два этапа: вычисление первой разности, что обеспечивает симметричность и центрированность помехи Л, и оценивание сигнала постоянного уровня или линейного тренда нулевого порядка согласно алгоритму (7).
2. В случае преобладания погрешности квантования измеряемого сигнала по уровню над другими погрешностями измерения распределение разностной помехи л удовлетворяет треугольному закону Симпсона.
3. В этих условиях, если величина параметра алгоритма (7) в выбирается оптимальным образом в соответствии с (9), дисперсия ошибки оценивания информативного параметра С* определяется исключительно интервалом квантования измеряемого сигнала по времени Д1 и разрядностью аналого-цифровых элементов 5. Тем самым применение алгоритма (7) с любым параметром Д (в том числе и Д= 0) обеспечивает равную точность (дисперсию) оценивания.
АНАЛИЗ АСИМПТОТИЧЕСКИХ СВОЙСТВ ОЦЕНКИ РЕКУРСИВНОГО АЛГОРИТМА Я.З. ЦЫПКИНА С ПОЗИЦИЙ ТЕОРИИ АВТОМАТИЧЕСКОГО УПРАВЛЕНИЯ
Другим способом доказательства указанных асимптотических свойств оценки — несмещенности и состоятельности, — свободным от ограничительных предположений, является преобразование алгоритма Я.З. Цыпкина к системе автоматического управления и анализ асимптотической устойчивости ее положения равновесия.
Введем ошибку оценивания истинной измеряемой величины с* как £п = сп - с * . Тогда свойство несмещенности оценки формулируется как £п ^ 0 при п ^ ^. Под состоятельностью следует понимать достижение указанного результата при любом законе распределения помех, т. е. независимо от помехи. Алгоритм (1) приводится к системе автоматического управления (САУ), изображенной на рис.1.
Представленная САУ обладает интересными свойствами. Она дискретна, обладает двумя контурами обратной связи, одна из которых отрицательна, имеет нелинейный элемент НЭ (неидеальное реле с зоной нечувствительности 2Д, которое реализует функцию ^), передаточная характеристика которого К1 содержит звено с переменными параметрами, реализующее умножение на масштаб в / п , с передаточной характеристикой К2, а также линию задержки на 1 такт (2Г1) с передаточной характеристикой К3.
Само понятие "зоны нечувствительности" нуждается в некотором пояснении. С одной стороны, можно ввести "алгоритмическую зону нечувствительности" — формальное задание величины параметра Д алгоритма (7). С другой стороны, существует и "физическая зона нечувствительности",
х[п+1]
К2[п]
-► с[п+1]
Кз
Рис. 1. Структурная схема САУ, реализующей алгоритм (7) (пояснения в тексте)
определяемая условиями эксперимента и обработки данных (например, разрядностью аналогоцифровых преобразователей, микропроцессорными ошибками округления и т.п.). То есть случай отсутствия зоны нечувствительности является некоторой вычислительной абстракцией.
Нетрудно видеть, что сформулированное выше требование несмещенности оценки коррелирует с формулировкой асимптотической устойчивости САУ, хотя эти утверждения не являются равносильными. В принципе реализуемая ситуация сходимости оценки сп не к с*, а к какой-либо другой величине с0 также свидетельствует об асимптотической устойчивости САУ, но свойство несмещенности оценки при этом отсутствует. Таким образом, требование асимптотической устойчивости САУ является необходимым условием несмещенности оценки сп.
Обзор различных методик оценки устойчивости нелинейных САУ и анализ устойчивости данной структуры представлены в работе [34].
Оценка асимптотической устойчивости нелинейной САУ проводится, исходя из анализа структуры ее линейной части с помощью критерия Попова или его модификаций. О наличии асимптотической устойчивости САУ можно судить по существованию решений неравенства (14), которое должно выполняться для всех ю:
Н (р) =
+Г +~ 1
= |ехр(-р1) V — [1(1 - (т + 1)70) -1(1 - т7))].
0 т=1т
Заметим, что слагаемое суммы с номером т отлично от нуля только при 1 е [тТ0, (т + 1)Т0 ]. Меняя порядок операций интегрирования и суммиро-
1 (т+1)Т0
вания, получим Н (р) = V— I ехр(-р1 )й1.
т
т=1 т тТ0
Слагаемые с номером т — табличные интегралы, равные соответственно
ехр(-тТ0) - ехр(-(т + 1)Т0)
1 -ехр(-рТ0) ( Т ч
-------------ехр(-ртТ0).
Первый множитель не зависит от индекса суммирования т. Таким образом,
Н (р) =
1 - ехр( -рТ0) ^ ехр(-трТ)
т=1
т
Яе[(1 + 7юИ)Ж(;ю)]+ 1/к > 0.
(14)
Вошедшая в формулу бесконечная сумма заменяется в соответствии с разложением в ряд Мак-
т=1
Здесь Ж — передаточная характеристика ли- лорена 1п(1 - х) = -^
нейной части САУ, к — тангенс угла наклона передаточной характеристики НЭ. В нашем случае к = 1/А .
Передаточная характеристика звена с переменными параметрами строится в несколько этапов.
а) Представление передаточной характеристики Н(і) во временной области
хт
т
Таким образом, окончательный вид передаточной характеристики звена с переменными параметрами есть
К2 (р) = - 1п(1 - ехр(-рТ0 ))
1 - ехр(-рТ())
1
к(і) = У — [ - (т + 1)Т) - 1(і - тТ0)].
т=т
Здесь 1(і) — функция Хевисайда (функция единичного скачка), определяемая как к(і) = 0 для і < 0; к(і) = 1 для і > 1 . Таким образом, передаточная характеристика принимает значения 1/т, когда і є [тТ0,(т + 1)Т0 ], т. е. на т-м шаге оценивания; Т0 — шаг дискретизации по времени.
б) Представление передаточной характеристики в частотной области (на основе преобразования
Лапласа) Н(р) =|й(і)ехр(-рі)^ при р = ]ю .
0
С точностью до численного множителя в
В формулу К2 (]ю) входит комплексная экспонента, заменяемая по формуле Эйлера как ехр(]'ю) = С08(ю) + 7 8Іп(ю), и комплексный логарифм — многозначная функция, определяемая как Ьп(г) = 1п|г| + У(а^ г + 2п1), І є 2 (|г|, аг§ г —
модуль и аргумент комплексного числа г). В рассматриваемом случае |г| = |28Іп(юА і /2), и
аг§ г = (юА і) / 2 или П - (юА і) / 2 в зависимости
от знака синуса.
Передаточная характеристика линии задержки есть К3 (р) = ехр(-рТ0).
Выражение для приведенной передаточной характеристики линейной части системы Ж(]ю) строится с учетом контура обратной связи. Преобразование САУ и выделение передаточной характеристики иллюстрируются рис. 2 [35].
“*0
НЭ
/
Ж(р)
Ж(р)
Рис. 2. Схема выделения передаточной характеристики Жл(р) линейной части САУ (пояснения в тексте)
Передаточная характеристика Ж(р) определяется выражением (15)
(15)
где Е(р) — изображение выходного сигнала е линейной части системы, ¥(р) — изображение входного сигнала / линейной части системы. В явном виде (15) представляется в форме
Е (р )=-КЩр1 • ¥ )■
Таким образом, можно построить Ж(/ю) — передаточную характеристику линейной части системы с учетом контура обратной связи. Выражение для Ж представлено в [34]. После его подстановки в критерий (14) при отсутствии входного воздействия получим условие
вс
4П
(1 - 2П) + в 1п|28Іп(пП)|И + А > 0, (16)
где О = (юТ0)/(2п)— нормированная фаза.
Подбор Н заведомо осуществим, если второе слагаемое (16) отлично от нуля. В противном случае требуется выполнение ограничения
А>вТ0(2 -1/Й)/4.
(17)
и требуемая величина Д = вТ0 / 5 определяется при О = 5/ 6 ;
г) неравенство (17) допускает интерпретацию: Т0 / Д < /(в) — аналог критерия устойчивости явной разностной схемы первого порядка, поскольку Т0 и Д можно интерпретировать как кванты времени и уровня соответственно. В то же время очевидно, что поскольку в алгоритме Я.З. Цыпкина (7) явно связаны известная сп и вычисляемая сп+1 оценки величины, то он и представляет собой разновидность явной разностной схемы первого порядка [37].
МОДИФИКАЦИЯ ПОДХОДА М. АОКИ ДЛЯ АНАЛИЗА СХОДИМОСТИ ОЦЕНКИ АЛГОРИТМА СТОХАСТИЧЕСКОЙ АППРОКСИМАЦИИ В ФОРМЕ Я.З. ЦЫПКИНА
В данном разделе рассматривается еще один метод доказательства сходимости сп к с*, который основывается на подходе М. Аоки [38].
Поскольку помеха £п случайна, то случайными величинами являются сп и еп. Обозначим плотность распределения вероятностей (ПРВ) и функцию распределения (ФР) случайной величины еп как /п (х) и ¥п (х) соответственно. Очевидно, что ПРВ сп будет просто смещена относительно ПРВ £п на величину с*.
Идея М. Аоки [38] допускает следующую интерпретацию: если р(/п, /п+к) ^ 0 при п и
любом к е N (р — аналог расстояния), то еп асимптотически стремится к нулю. Выбор функции р жестко не задан, но в [38] выбрана М-оценка Хьюбера р( /п
, /п+к ) = эир|/п (х)
/п+к
(х) при
х е Я . В работе [39] рассматривались и другие определения расстояния между ПРВ /1 (х) и /2 ( х) . Например, информационное расхождение Кульбака [40] или квадратичная оценка вида
Р = |(/и(х) - /2(х))2дх .
Исследование неравенства (16) позволяет сделать выводы [34-36]:
а) при отсутствии зоны нечувствительности данная САУ не имеет асимптотически устойчивого положения равновесия (в ^ 0, иначе все оценки с* будут повторять с1);
б) в случае ограничения производной передаточной характеристики НЭ (например, путем введения зоны нечувствительности) асимптотически устойчивое положение равновесия достижимо;
в) следствием дискретности САУ будет О < 1,
Докажем, что при п ^ те разность
/п+1 (х) - /п (х) ^ 0 при любом х е Я .
Обозначим плотность распределения вероятностей и функцию распределения помехи £ соответственно как /х) и ¥(х). Рассмотрим элементарное событие йА п+1, заключающееся в том, что £п+1 = х, или £п+1 е [х, х + dx]. Его бесконечно малая вероятность Р в соответствии со свойствами ПРВ есть с точностью до бесконечно малых более высокого порядка /п+1 (х)<йх .
Сложное событие йА п+1 представляется в форме:
е
йАп+1 = (V = -1П £п ~ х - в / п) и
и (т = 0 П £п = х) и и (Т = +1 П &п “ х +в / п), (18)
где П и и обозначают произведение и сумму событий соответственно.
Необходимо учесть: а) независимость £п от £п+к, что приведет к неизменности ПРВ и ФР помехи на любом шаге оценивания; б) альтернативность событий V = -1, V = 0 и V = +1 — первых множителей в слагаемых суммы (18); в) независимость событий — первых множителей от событий — вторых множителей каждого слагаемого (18). Последняя часть утверждения верна, поскольку случайная помеха не должна зависеть от текущей ошибки оценивания.
Вероятности элементарных событий, входящих в (18), выражаются как
Р(£п = х -в / п) = /п (х -в / п)йх ,
Р(£п ~ х + в / п) = /п (х + в / п)<Эх ,
Р(£п = х) = /п(x)dx.
В то же время вероятность событий — первых множителей представляет собой условную вероятность реализации первого события при условии реализации второго события. Тем самым вероятность Р+ события V = +1 вычисляется по следующей схеме:
Р+ = Р(еп - £п+1 > Д) = Р(£п+1 < £п - Д) =
= Р(£п+1 < х + в /п -Д) = ¥(х + в /п -Д). (19)
Величины Р0 = Р( V = 0) и Р- = Р^ = -1) вычисляются по схемам, аналогичным (19), с использованием известных соотношений для ФР. Таким образом, Р0 = ¥(х + Д) - ¥(х - Д) и Р- = 1 - ¥ (х -в / п + Д).
Выполнив почленное сокращение на dx всех слагаемых левой и правой частей равенства (18), можно получить рекуррентную формулу связи ПРВ последовательных случайных величин еп и
£п+1 [41]:
/п+1 (х) = /п (х - в /п)[1 - ¥(х - в /п + Д)] +
+ /п (х)[¥(х + Д) - ¥(х -Д)]+ (20)
+ /п (х + в / п)¥(х + в /п -Д);
/1 (х) = /(х).
В работе [39] был рассмотрен случай сигнатурного
оценщика (Д = 0) и получена рекуррентная формула на основе аналогичных рассуждений. Она имеет более простой вид по сравнению с (20):
/п+1( х) = /п (х -в / п)[1 - ¥ (х -в / п)] +
+ /п (х + в / п)¥ (х +в / п). (21)
Безусловно, (21) может быть также получено из общей формулы (20) подстановкой Д = 0.
Поскольку при п ^ те в / п ^ 0, из (20) и (21) следует, что при произвольном х е Я разность между /п+1( х) и /п (х) также будет бесконечно малой.
Очевидно, что доказанное утверждение является необходимым условием асимптотического стремления р к нулю. Можно показать, что оно будет и достаточным.
Как утверждается в [38], при п ^ да ПРВ /п(х) ^5(х), где 5(х) — дельта-функция Дирака.
Следовательно, при использовании алгоритма Я.З. Цыпкина (7) в условиях искажения постоянного сигнала с* аддитивной центрированной помехой £ выполнено условие £п ^ 0 несмещенности оценки сп относительно величины с*.
Допустив помимо центрированности помехи ее симметричность, можно на основе метода полной математической индукции по номеру шага оценивания п и с учетом (20) доказать симметричность любой ПРВ /п+1 (х), п > 1.
Рекуррентные формулы (20) и (21) позволяют промоделировать ПРВ /п ( х) на начальных шагах оценивания, задавая в качестве /1(х) ПРВ случайной помехи £. Динамика ПРВ ошибки оценивания на первых шагах применения алгоритма (7) иллюстрируется рис. 3-5.
Анализ рис. 3-5 позволяет утверждать, что даже при столь различных закона распределения случайной помехи (равномерном, релеевском и лаплассовском) за несколько начальных шагов оценивания ПРВ ошибки оценивания приобретает ярко выраженный одномодальный и центрированный характер. Исходя из выражений (20) и (21) можно утверждать, что на характер динамики /п ( х) , помимо ПРВ случайной помехи, существенным образом влияют параметры алгоритма (7) в и Д. Однако качественно характер динамики /п (х) будет одинаков как для любого закона распределения случайной помехи, так и для любого набора параметров алгоритма (7).
3 О
Рис. 3. Равномерный закон распределения
-3 О
Рис. 4. Закон распределения Лапласа
■3 О
Рис. 5. Закон распределения Релея
ПРАКТИЧЕСКАЯ РЕАЛИЗАЦИЯ РЕКУРСИВНОГО АЛГОРИТМА ОЦЕНИВАНИЯ ПАРАМЕТРОВ ЛИНЕЙНЫХ ТРЕНДОВ НУЛЕВОГО И ПЕРВОГО ПОРЯДКОВ С КОНТУРОМ ОПРЕДЕЛЕНИЯ РАЗЛАДКИ
Как говорилось ранее, в приборах иммунного и химического экспресс-анализа чаще всего удается достичь того, чтобы выходной сигнал х(1) представлял собой какую-либо совокупность линейных трендов нулевого (х(1) = с *) и первого (х(1) = а + с * 1) порядков. Достаточно распространенными, в частности, являются сигналы типа "линейный тренд с насыщением" или "трапеция". При этом информативным параметром, как правило, является величина сигнала тренда нулевого порядка и тангенс угла наклона линейного тренда первого порядка — с*. Оценивание с* осуществляется при наличии аддитивной помехи в условиях минимума априорной информации о законе ее распределения, ограниченности времени оценивания и объема вычислительных ресурсов. Последние ограничения связаны с использованием Р1С-процессоров, применение которых в значительной степени позволяет уменьшить габариты вычислительных устройств и приборов в целом.
Известно, что понижения порядка тренда с одновременным преобразованием помехи к симметричному виду как результат свертки помехи с самой собой можно достичь переходом к первому разностному сигналу. При этом информативные параметры непосредственно или с помощью перехода к первому разностному сигналу можно оценивать как величину тренда нулевого порядка.
Подобная задача может быть решена с помощью модификации Я.З. Цыпкина алгоритма стохастической аппроксимации. Достоинства алгоритма (7) — экономичность и робастность — хорошо известны. Вместе с тем выполнение требования к выбору параметра в и методика выбора Д должны базироваться на априорно отсутствующей нетривиальной информации либо о законе распределения помехи, либо в крайнем случае о величине плотности распределения вероятностей и функции распределения помехи на границе зоны нечувствительности.
Для практической реализации устройства измерения постоянного сигнала на фоне аддитивной помехи требуется решить несколько задач [42]:
а) задание параметров алгоритма — в и Д,
б) выбор начального приближения (оценки с1),
в) формулирование критерия остановки оценивания.
При этом алгоритм, представляющий рекурсивную процедуру оценивания, становится вычислительным элементом устройства измерения постоянного сигнала. Важным требованием является
использование процедур, требующих малых вычислительных затрат. В частности, решено было отказаться от алгоритмов оценивания, предполагающих упорядочивание, например вычисление медианы и центра сгиба, оценок типа Тьюки и т.д.
Подбор параметра в
Для оценивания величины вшп необходимо приближенно рассчитать величину плотности распределения помехи на границе зоны нечувствительности. В простейшем случае идеального реле при Д = 0 необходимо оценить плотность распределения помехи в нуле — р(0).
Эта информация может быть получена с помощью введения начального (холостого) участка из N измерений, на котором можно оценить размах помехи [43] и среднюю плотность распределения вероятностей помехи р * . В том случае, если Д « 0, то р(Д) > р * в соответствии с аксиомами Маликова. Замена р(Д) меньшей величиной р * приведет к завышению в, что позволяет сохранить несмещенность и состоятельность оценки, но может вызвать незначительное увеличение дисперсии ошибки оценивания.
Размах помехи есть величина, обратная р*, и ее можно оценить через экстремальные порядковые статистики — х(1) и х(М) выборки М измерений, проведенных на холостом участке с помощью не требующих упорядочивания двух рекурсивных определений максимального и минимального значений
р*
М-1
М + 1)х(1) - х(М))'
(22)
В случае одномодальной симметричной или полимодальной (равномерной) помехи параметр в, рассчитанный на основе (22), будет завышенным, но с точки зрения выполнения условий состоятельности и несмещенности оценки сп этот эффект не является негативным. Возможность двухмодальных (типа арксинуса) распределений помех аналоговых и цифровых элементов измерительных схем подавляющим большинством исследователей, например [29, 32], опровергается.
Подбор начального приближения — оценки Сх
В том случае, если на первом шаге измерения реализуется помеха типа выброса, начальная ошибка оценивания будет весьма велика. Тем самым время сходимости оценки сп к истинному значению с* будет сильно увеличено.
Использовав информацию, полученную на холостом участке измерения, в частности экстремальные порядковые статистики, можно выбрать
начальное приближение как центр размаха (хтт + х тах)/2. В случае выброса на первом измерении начальная ошибка оценивания уменьшится по крайней мере вдвое. Кроме того, в качестве начального приближения с1 можно взять 1/М-усеченное среднее, исключив из всей суммы измерений потенциальные выбросы — хтт и хтах. Эта процедура также не требует упорядочивания, т. к. требуется простое суммирование измерений и вычитание экстремальных порядковых статистик. Очевидно, что в качестве первого приближения может быть взята и любая комбинация центра размаха и 1/М-усеченного среднего.
Методика определения разладки
в последовательности измерений
Критерием остановки оценивания величины с* может явиться либо достижение требуемой точности, либо выявление разладки в последовательности измерений. Под разладкой следует понимать произошедшее на т шаге измерения либо скачкообразное изменение величины сигнала с* вида
с. * = с * +Дс • 1 [ - т + 1/2], либо изменение его формы, например, из-за наложения линейного тренда с. * = с * +а • (/ - т) • 1 [ - т]. Здесь 1^) — функция Хевисайда, т — номер шага, на котором наступает разладка.
Первый вид разладки можно условно назвать "разладка типа скачка", второй — "разладка типа линейного слагаемого", или "дрейфа". При использовании алгоритма Я.З. Цыпкина второй вид разладки может интерпретироваться как систематическая линейная погрешность. В обоих случаях необходимо выявлять факт разладки и останавливать процесс оценивания с*.
Выявление разладки непосредственно связано с оценкой однородности выборки. Решения, основывающиеся на критерии Аббе, неприемлемы, поскольку, во-первых, закон распределения помехи может быть отличен от нормального; во-вторых, требуется сократить объем вычислений, отказавшись от расчета частичных дисперсий и т.д.; в-третьих, критерий должен обладать робастностью. Алгоритм кумулятивных сумм Пейджа и его робастные аналоги [44-46] неприемлемы по причине сложности алгоритма принятия решения, т. к. требуется устанавливать два пороговых значения, априорно задавая время задержки и вероятность ложной тревоги. Третий подход предполагает разрешение альтернатив типа "однородность—неоднородность" выборки, либо "принадлежность измерения семейству 1—семейству 2" на основе кумулятивных сумм или информации Фишера [47-49]. Очевидно, что и этот подход требует существенных вычислительных затрат. Кроме того, логично критерий разладки формировать на основе величин, уже используемых при оценивании сигнала с*.
В работе [50] предложен критерий разладки, базирующийся на анализе величины V п — после -довательности знаков поправок (сп+1 - сп ) . Заметим, что эти величины вычисляются на каждом шаге оценивания с*. В принципе данный критерий можно считать одним из робастных аналогов алгоритма кумулятивных сумм Пейджа. Поскольку помехи на каждом шаге формируются независимо, то вероятность наличия к положительных поправок в последовательности Ь измерений определяется в соответствии со схемой Бернулли (к успехов в Ь испытаниях при вероятности успеха 1/2 для симметричных помех). Вероятность успеха есть 1/2, поскольку для достаточно больших номеров п оценка сп близка к с*, и знак поправки практически определяется знаком помехи. С другой стороны, на начальных шагах оценивания, особенно в случае завышенных оценок параметра в, осуществляется практически детерминированное чередование знаков поправок. Последнее приведет к "утяжелению" центра гистограммы по сравнению с гистограммой, соответствующей схеме Бернулли.
Выдвинутые предположения о симметричности гистограммы в случае отсутствия разладки хорошо соотносятся с предложенными в работах [51, 52] правилами идентификации сигнала постоянного уровня, установленными эмпирически, без теоретического обоснования. В работе [51] предложено правило "1 из 5", примененное в условиях значительных уклонов. Так, признаком отсутствия линейного тренда первого порядка является наличие не менее одной перемены знака в последовательности из 5 измерений . То есть исключаются события "0 успехов" и "5 успехов", суммарная вероятность которых составляет 1/16, или около 6 %. В работе [52] сформулировано правило "3 из 10", применимое, по мнению автора, при условиях малости ошибки оценивания по сравнению с масштабом помехи. То есть приведенные выше рассуждения относительно применимости схемы Бернулли обоснованы. Признаком указанной ситуации является наличие не менее 3 перемен знака в 10 измерениях. В принятой нами терминологии 0, 1, 2 или 8, 9, 10 успехов
в 10 испытаниях считаются практически невозможными для симметричной схемы с р = q = = 1/2. На самом деле, исключаются события, суммарная вероятность которых есть 112/1024, или около 10 % . Естественно, что наличие гистограммы распределения V,, значительно более информативно по сравнению с одиночными измерениями
и, следовательно, позволяет делать статистически более достоверные выводы.
Таким образом, для решения указанной задачи используются значения знака поправки рп = sign(cn - хп+1), их сумма в скользящем окне
ширины Ь — Юп = и гистограмма значе-
к=п- Ь
ний юп — Ст[п]. Заметим, что пересчет Ст[п+1] из О0[п] представляет собой экономичную и простую рекурсивную процедуру, выполняемую в реальном времени и заключающуюся в добавлении единицы к содержимому одной из (Ь+1) ячеек памяти. Сами знаки поправок неизбежно вычисляются при реализации алгоритма (7), поскольку используются для пересчета оценки сп+1 из сп.
Достаточно быстро наступает момент, когда выполняется условие /3/п << s, где 5 — характерный масштаб помехи (стандартное отклонение, интерквартильная ширина, медиана абсолютных отклонений, в ряде случаев — размах). В этом случае величина фп, во-первых, определяется исключительно знаком помехи; во-вторых, практически не зависит ни от ошибки оценивания, ни от фп_1. То есть формирование юп происходит в соответствии со схемой Бернулли, а сама величина юп есть (2т - Ь), где т — число успехов в Ь испытаниях. Под успехом следует понимать, например, положительную помеху. Таким образом, при отсутствии разладки юп формируется в соответствии со схемой Бернулли с равными вероятностями успеха и неудачи, а, следовательно, гистограмма От должна быть симметрична. На начальных шагах
оценивания схема Бернулли нарушается, поскольку в/п становится соизмеримым и, может быть, превосходящим 5. В этом случае будет осуществляться практически детерминированное чередование знаков фп. По указанной причине не очень эффективным представляется использование
X -критерия согласия в качестве меры совпадения гистограммы с эталонной, соответствующей симметричной схеме Бернулли для Ь испытаний. Высокая чувствительность указанного критерия не может дать статистически значимых результатов в этих случаях.
При наличии разладки любого из двух типов симметрия схемы нарушается тем сильнее, чем более скачок Ас, величина а или "накопленные" ("аккумулированные") шаги измерения (п - т), поскольку вероятность успеха становится ^(±Ас) или Е(±аТ0(п - т)) вместо 1/2. Здесь Т0 — интервал дискретизации по времени, ^ — функция распределения помехи.
Таким образом, признаками разладки являются асимметрия гистограммы и смещение ее центра относительно нуля. Мерой может служить нормированный начальный статистический момент 3-го порядка. Примеры гистограмм при отсутствии разладки, а также при разладках двух типов представлены на рис. 6-8 соответственно.
еоо
-В -5 -4 -3 -2 -1 О 1 2 3 4 5 6
Рис. 6. Пример гистограммы знаков при отсутствии разладки
-6 -5 -4 -3 -2 -1 0 1 2 3 4 6 6
Рис. 7. Пример гистограммы знаков при разладке типа "скачок"
800
600
400
Рис. 8. Пример гистограммы знаков при разладке типа "дрейф"
СРАВНЕНИЕ РЕКУРРЕНТНОГО И ИНТЕРВАЛЬНОГО АЛГОРИТМОВ ПОЛУЧЕНИЯ ОЦЕНОК НА ИНФОРМАЦИОННОЙ ОСНОВЕ
Предположим, что оба алгоритма оценивания
(рекуррентный и интервальный) построены
на основе какого-либо критерия качества и соотносят оценку параметра положения сп с самим параметром с* так, что на определенном этапе обработки сигнала дисперсии ошибок оценивания
оказываются равными. Возникает вопрос, какой оценщик применять в приборах: рекуррентный или интервальный?
Критерием для выбора предлагается следующий [25, 53-55]:
Т(Д) = 1/ (А)/ (А), (23)
где Iу (А) — фишеровское информационное количество, получаемое на каждом шаге обработки при независимых выборках сигнала хг, г = 1,...,п ; I^(А) — шенноновское информационное количество, затрачиваемое на каждый шаг обработки.
Таким образом, согласно (23), преимущество имеет оценщик, использующий для получения оценки минимальное количество информации по Шеннону, или, иными словами, обладающий большей производительностью создания фишеровской информации.
Оценка параметра положения сигнала типа линейного тренда нулевого (постоянный сигнал) или первого порядка при интервальном и рекуррентном алгоритмах
Рассматриваемая модель предполагает, что имеются измерения сигнала хг = с * +£г неизвестного параметра с*, содержащие аддитивные случайные помехи . Помехи независимы и имеют известную плотность распределения вероятностей р(х). Нужно оценить с* по результатам наблюдений хг.
Одним из основных методов, применяемых в построении алгоритмов оценивания, является метод максимального правдоподобия (ММП). Суть ММП заключается в том, что в качестве оценки параметра сигнала с берется значение, максимизирующее функцию правдоподобия
п
Ь( х^..., Хп, с) = П р(хг - с) [56]. Переходя от
г =1
функции правдоподобия к ее логарифму, имеем
п
сп = а^шш 1п (с), 1п (с) = 1п р( х - с). (24)
г=1
Последнее выражение задает оценку максимального правдоподобия. Если оценочная функция ^ (х) = - 1п р( х) дифференцируема и выпукла, то задача минимизации в (24) эквивалентна реше-
п
нию уравнения ^ ^(хг - c)/dx = 0. Во многих
г =1
случаях это уравнение решено явно [56]. Так, при нормальном законе распределения помехи оптимальной в указанном выше смысле оценкой пара-
метра положения будет выборочное среднее, для распределения Лапласа — выборочная медиана, для равномерного распределения — полусумма экстремальных значений (центр размаха) и т.д.
Получение указанных оценок возможно как интервальными, так и рекуррентными алгоритмами. В работе [23] показано, что в случае выполнения условий регулярности оценки, полученные как рекуррентным, так и интервальным способами, по ММП эквивалентны и асимптотически эффективны, т. е. достигают границы Рао—Крамера. Однако оптимальность оценок ММП имеет свою оборотную сторону — они оказываются чувствительными к отклонениям распределения от предполагаемого [23]. В связи с этим возникает задача "огрубления" оценок ММП, или, пользуясь современной терминологией, необходимо найти робастные алгоритмы оценивания.
Обычно истинное распределение помех, или функция р(х), неизвестно, однако известна некоторая о ней информация. Так или иначе можно считать, что известен некоторый класс распределений, которому принадлежит истинное распределение. Известно [23], что в данном случае целесообразно найти в этом классе наименее благоприятное распределение, а затем применить соответствующую ему оценку ММП. При этом под наименее благоприятным распределением естественно полагать то, что дает наихудшую границу неравенства Рао—Крамера, минимизирует фишеровскую информацию в известном классе распределений.
В работе [23] показано, что для класса невырожденных распределений, для которых
р(0) > 1/(2а) > 0, т. е. в условиях практически полного отсутствия априорных сведений о распределении, наименее благоприятным является распределение Лапласа
р(х) = — ехр(-|Х / а). (25)
2а
Указанные выше рассуждения обосновывают использование оценщика вида
Сп+1 = Сп +Уп Э18П( Х„+1 - Сп).
Эта формула описывает сигнатурный алгоритм, совпадающий с алгоритмом Я.З. Цыпкина (7). То есть оценка, формируемая алгоритмом (7), обладает свойством минимаксности — минимальная дисперсия ошибки оценивания в случае наихудшего (в смысле максимума дисперсии) распределения помех в соответствии с лаплассовским законом распределения.
При использовании любого рекуррентного алгоритма существенную роль играет выбор начального приближения с. Для выбора хорошего на-
чального приближения можно применить сначала нерекуррентный алгоритм, например использовать алгоритм поразрядного уравновешивания или выбрать медиану первых пяти измерений [23].
Сравнение рекуррентного и интервального алгоритмов
Можно показать, что как рекурсивный, так и интервальный алгоритмы позволяют получать оптимальную оценку параметра положения — выборочную медиану, которая является асимптотически эффективной, т. е. достигает нижней границы Рао—Крамера.
Предположим, что в процессе обработки выборок сигнала дисперсии выборочных медиан, полученных интервально и рекуррентно, стали равными между собой. Для того чтобы сравнить алгоритмы по скорости создания фишеровской информации, т. е. чтобы воспользоваться критерием (23), вычислим необходимые объемы шенноновской информации. Тот алгоритм, который потребует меньше шенноновской информации будет в информационном смысле лучше. Найдем шенноновское количество информации 11п(. для интервального алгоритма на один шаг обработки при получении оценки. Количество информации, согласно Шеннону, определяется как разность энтропии [57] 11п = Н(хг) - Н(хг |хп), где первое
слагаемое — энтропия выборочного значения до измерения хп, второе слагаемое — энтропия действительного значения после использования элемента выборки хп .
Общее количество шенноновской информации, необходимое для получения оценки сп будет в п раз больше, т. к. производится п шагов оценивания. Последнее утверждение применимо и к рекурсивному алгоритму.
Энтропия выборочного значения до измерения вычисляется по известной формуле
н(х1) = -|р(х)1о§2(р(x))dx . (26)
Подставляя в выражение (26) значение р(х) для лапласовского распределения (25), получим, что Н(хг) = (4а). Если оценщик обладает по-
грешностью с равномерным распределением в диапазоне шириной 2Д, то энтропия результата измерений после получения показания хп есть 1о§2(2Д). Основание логарифма, равное двум, позволяет использовать в качестве единиц информации (энтропии) биты.
Итак, количество шенноновской информации для интервального оценщика при п шагах алгоритма = п (2а / А).
При рекуррентном алгоритме на каждый шаг получения оценки тратится только один бит информации, следовательно IЦ = п .
Заметим, что дисперсия оценки параметра положения по медиане с учетом эффекта квантования по уровню (поправка Шеппарда) есть
Б = (а2 + А2 /12)/ п .
Если оценщик спроектирован так, что А << а, то дисперсии оценки при интервальном и рекурсивном методе оценивания, отличающиеся наличием поправки Шеппарда, становятся практически одинаковыми, и выигрыш в скорости создания фишеровской информации для рекуррентного алгоритма по сравнению с интервальным составляет (2а / А) раз.
Указанный выигрыш рекуррентного алгоритма действителен при хорошем начальном приближении с. Реально, т.е. с учетом па начальных шагов адаптации (выхода на начальное приближение с0), количество шенноновской информации будет
12 = п - па + па (2а / А) . Соответственно кор-
ректируется и выигрыш в скорости создания фишеровской информации.
Эффект от использования рекуррентного алгоритма приведен в таблице для различных соотношений а/А и па/п .
Выигрыш в скорости создания фишеровской информации от применения рекуррентного алгоритма
Очевидно, что при оценке параметра положения сигнала рекуррентный алгоритм эффективнее в информационном смысле интервального. Наибольший эффект от использования рекуррентной процедуры получается при большом динамическом диапазоне изменения оцениваемого параметра сигнала. При этом уменьшение коэффициента загрузки канала ввода позволяет более эффективно использовать вычислитель для реализации остальных функций приборов.
Оценка параметра линейного тренда в режиме кинетического анализа
В предыдущем разделе было показано, что рекуррентное устройство эффективнее интервально-
го по скорости получения оценки параметра положения сигнала по крайней мере в 1.5 раза. Наибольший эффект от использования рекуррентной процедуры получается при большом динамическом диапазоне сигнала. Следствием этого преимущества является освобождение канала ввода информации в вычислитель оценщика, другими словами рекуррентная процедура более эффективна по критерию получения единицы фишеровской информации на единицу затрат шенноновской.
Структура любого современного оценщика включает в себя аналого-цифровой преобразователь и вычислитель с каналом ввода информации. В связи с этим модель модифицируется к виду
Уп =оп +£п А = Ь - ^, (27)
где п = 0,1,. — номера отсчетов.
Модель (27) показывает, что математическое ожидание М(у) = вп есть некоторая величина, линейно изменяющая свое значение. Необходимо найти алгоритм, согласно которому оценщик должен по наблюдениям величины Уп производить оценку вп в классе рекуррентных процедур, т. е. строить следующую зависимость
Ап+1 = / (Ап , Уп+1 ) , где Ап ,Ап+1 — оценки математического ожидания величины сигнала в моменты п + 1 и п соответственно, уп+1 — выборочное значение входного сигнала.
В условиях модели (27) изменение математического ожидания (или тренд величины у) может быть оценено с помощью модификации метода стохастической аппроксимации Роббинса—Монро, предложенной Дупачем [24]. Необходимые исходные положения и ограничения, накладываемые на параметры алгоритма, подробно разобраны в работах [53-55]. Одной из возможных модификаций такого алгоритма (аналогом сигнатурной модификации Я.З. Цыпкина (7)) может быть алгоритм в форме
хп+1=хп *+-ва э^^п+1- хп *]
па (28)
хп * = (1 + п _1)хп .
Здесь хп — оценка величины вп на п шаге оценивания линейного тренда (27).
Исходя из условий сходимости числовых рядов, параметр а должен выбираться из условия
0.5 < а < 1. Оптимальная скорость сходимости оценки для линейного тренда без смещения (Ь = 0) обеспечивается при а = 1, в общем случае — при а = 2/3.
Ввиду использования на входе алгоритма (28) лишь знака разницы прогноза (хп*) и реализации
Уп+1, эту процедуру также можно отнести к сигнатурным рекуррентным алгоритмам.
АППАРАТНАЯ СТРУКТУРА РЕКУРРЕНТНОГО ОЦЕНЩИКА
Использование микропроцессорных устройств в приборостроении позволяет применять достаточно сложные алгоритмы статистической обработки сигналов, поступающих с датчиков. При этом на аналоговую часть целесообразно возложить минимальные функции по обработке сигнала с целью их удешевления и уменьшения габаритных характеристик, особенно в многоканальных приборах. Аналоговый тракт содержит предварительный и нормирующий усилители, последний подключен к преобразователю. Вся дальнейшая обработка осуществляется в микропроцессорной среде. Очевидно, что наиболее сложным, дорогостоящим и "медленным" элементом типового тракта аналогового ввода сигналов в вычислитель является АЦП.
Полученные алгоритмы оценивания в качестве входной информации используют лишь знак разности аналогового входного сигнала и выработанный прогноз оценки. Это дало возможность предложить интересную и экономичную реализацию входного тракта [54, 55].
В основе подхода лежала идея использования компаратора и цифро-аналогового преобразователя (ЦАП) в качестве АЦП, что обеспечивало существенные преимущества в быстродействии, сходимости и экономичности. Комбинация ЦАП—компаратор превосходит АЦП на порядок по пропускной способности (бит/с) и более чем на порядок по потребляемой мощности (Дж/бит) [58]. Учитывая необходимость применения в комбинации с АЦП устройств выборки—хранения, преимущество обсуждаемой структуры становится еще заметнее.
Таким образом, предлагаемый рекуррентный сигнатурный оценщик имеет определенные преимущества перед классическими в части меньшего объема оперативного запоминающего устройства (ОЗУ), а также сходимости и быстродействия аналогового ввода.
Структура ЦАП—компаратор в качестве АЦП известна, изучена и, как показано в [58], хорошо сопрягается с ПЭВМ и микропроцессорами различных систем и является достаточно дешевой.
ЗАКЛЮЧЕНИЕ
Рассмотренный алгоритм (7) оценивания постоянного сигнала с* в условиях аддитивной помехи с априорно неизвестным законом распреде-
ления может эффективно использоваться в качестве элемента вычислительного устройства приборов химического и биологического (иммунного) экспресс-анализов.
Трудность описания закона распределения помехи связана со значительным разнообразием источников помехи. Источники помехи могут быть: химического (биологического) происхождения; связаны с неоднородностью чувствительного элемента и существенно стохастическим процессом установления равновесия между различными фазами реакции в начальные моменты времени; порождены оптическими эффектами (из-за различных засветок), аналоговыми и цифровыми элементами электрических схем (тепловые дрейфы, дро-бовый шум, погрешности квантования по уровню и т. д.). При существующих и признанных моделях [32] отдельных источников помех описать общий закон распределения для столь сложной их комбинации все равно не представляется возможным.
Сведение задачи определения информативного параметра к оцениванию сигнала постоянного уровня позволяет эффективно использовать алгоритм (7), отличающийся простотой аппаратной и программной реализаций. При этом, например, переход химической (биологической) реакции в стадию насыщения будет истолкован и выявлен алгоритмом как разладка типа скачок (т. к. оцениваемый сигнал вместо некоторой величины примет нулевое значение).
Таким образом, в работе рассмотрены различные аспекты разработки, исследования свойств и аппаратно-программной реализации устройства оценивания, ориентированного на применение алгоритма стохастической аппроксимации в робастной модификации Я.З. Цыпкина (7). Указанный алгоритм правомочен и эффективен в использовании для решения разнообразных задач химического и биологического (иммунного) экспресс-анализов. В частности, описанные алгоритмы применялись в разработанных в лаборатории "Информационно-измерительных био- и хемосен-сорных микросистем" Института аналитического приборостроения РАН спектрофотометрических приборах серии Sen и ^Sen.
Благодарности
Авторы благодарны Дмитрию Алексеевичу Бу-рылову, т. к. значительная часть описанных в данной работе результатов была получена нами совместно, а также Юрию Сергеевичу Шинакову — руководителю секции "Теория сигналов" на Международной конференции по цифровой обработке информации и ее применениям и на Научной сессии, посвященной Дню радио, чьи критические замечания безусловно способствовали более эффективной работе в данной области научных ис-
следований [59]. Работа поддерживалась ФЦП "Интеграция" (проект № А0141 "Оптика и научное приборостроение").
СПИСОК ЛИТЕРАТУРЫ
1. Тарасов В.В., Ягодин Г.А. Кинетика экстракции // Итоги науки и техники: Серия "Неорганическая химия". ВИНИТИ, 1974. Т. 4. 151 с.
2. Тарасов В.В., Ягодин Г.А., Пичугин А.А. Кинетика экстракции неорганических веществ // Итоги науки и техники: Серия "Неорганическая химия". ВИНИТИ, 1984. Т. 11. 187 с.
3. Тарасов В. В. Межфазные явления и кинетика экстракции неорганических веществ. Автореф. дис. ... д-ра хим. наук. 1980. 44 с.
4. Тарасов В.В., Ягодин Г.А., Пичугин А.А. Определение вкладов поверхностной и объемной реакций в системе жидкость—жидкость // Доклады АН СССР: Серия физич. химия. 1983. Т. 269, № 6. С. 1398-1404.
5. Иванов А.Б., Тарасов В.В., Ягодин Г.А. Об отличительных признаках объемных и поверхностных реакций, обусловливающих извлечение в экстракционных системах // Доклады АН СССР: Серия физич. химия. 1979. Т. 244, № 4. С.928-931.
6. Ягодин Г.А., Тарасов В.В. Межфазные явления в системах электролит—неэлектролит и их влияние на кинетику экстракции // Химия экстракции. Новосибирск, 1984. С. 35-52.
7. Курочкин В.Е., Теровский В.Б. Биосенсоры и иммуносенсоры: Обзор // Научное приборостроение. 1995. Т. 5, № 1-2. С. 3-12.
8. Stenberg M., Nygren H. Kinetics of antigen-antybody reactions at solid-liquid interfaces // J. of Immunogical Methods. 1988. V. 113. P. 3-15.
9. Евстрапов А.А., Курочкин В.Е., Макарова Е.Д. Отражательная фотометрия пластифицированных мембран в задачах обнаружения и оценки концентрации веществ в водных пробах // Научное приборостроение. 1991. Т. 1, № 4. С. 22-
35.
10. Буляница А.Л. Уточнение параметров моделей процессов детерминированного вида // Научное приборостроение. 1995. Т. 5, № 1-2. С. 38-
46.
11. Буляница А.Л. Массоперенос в коаксиальных элементах проточных аналитических микросистем. Автореф. дис. ... канд. физ.-мат. наук. СПб.: ИАнП РАН, 1997. 183 с.
12. Буляница А.Л., Курочкин В.Е., Макарова Е.Д. Оценивание диффузионного потока при конвективном массопереносе в тонком коаксиальном капилляре конечной длины // Научное приборостроение. 1997. Т. 7, № 1-2. С. 28-39.
13. Разработка электрохимических систем регистрации иммуноспецифических взаимодействий. Отчет о НИР / Рук. темы Ивницкий Д.М., исп.: Ситдыков Р.А., Юлаев М.Ф. и др. Самарканд,
1989. 196 с.
14. Henekamp H.B., Van Nieuwerk // Anal. Chim. Acta. 1980. V. 121. P. 13-22.
15. Kemula W., Kuther W. Modern trends // Anal. Chem. Budapest, 1984. P. A3-A31.
16. Поляк Б.Т., Цыпкин Я.З. Адаптивные алгоритмы оценивания (сходимость, оптимальность, стабильность) // Автоматика и телемеханика. 1979. № 3. С. 71-84.
17. Robbins H., Monro S. A stochastic approximation method // Anal Math. Stat. 1951. V. 22, N 3. P.400-407.
18. Cochran W., Davis M. Sequential experiments for estimating the median lethal dose // Colloques internationaux du centre National de la Recherche Scientifique. 1963. N 110. P. 181-194.
19. Kesten H. Accelerated stochastic approximation // Ann. Math. Stat. 1959. V. 29, N 1. P. 41-59.
20. Kiefer J., Wolfowitz J. Stochastic estimation of the maximum of a regression function // Ann. Math. Stat. 1952. V. 23, N 3. P. 462-466.
21. Fabian V. Stochastic approximation methods // Czech. Math. Journ. 1960. N 10 (85). P. 123-159.
22. Dvoretzky A. On stochastic approximation // Proc. of the 3rd Berkeley symposium on Mathematical Statistics and Probability. Berkeley, California: Univ. of Calif. Press., 1956. N 1. P. 39-55.
23. Цыпкин Я.З., Поляк Б. Т. Огрубленный метод
максимального правдоподобия // Динамика систем. Изд-во Горьковского ун-та, 1977.
№ 12. С. 22-46.
24. Dupac V. Dynamic stochastic approximation method // Anal. Math. Stat. 1965. V. 36, N 6. P.1695-1702.
25. Курочкин В.Е., Фельдман Б.Х. Оценка параметров тренда фотометрических анализаторов биосубстанций при кинетических исследованиях // Научное приборостроение: Методы и приборы биотехнологии. 1988. С. 80-86.
26. Цыпкин Я.З. Оптимизация в условиях неопределенности // Доклады АН СССР. 1976. Т. 228, № 6. С.1306-1309.
27. Бедельбаева А.А. Релейные алгоритмы оценивания // Автоматика и телемеханика. 1978. № 1.С. 87-95.
28. Красулина Т.П. О сходимости снизу модифицированного процесса Роббинса—Монро // Автоматика и телемеханика. 1992. № 4. С. 7376.
29. Маликов М.Ф. Основы метрологии. М.: Стан-дартгиз, 1949. 480 с.
30. Хампель Ф., Рончетти Э., Рассеу П., Шта-эль В. Робастность в статистике: Подход на
основе функций влияния (пер. с англ.). М.: Мир, 1989. 512 с.
31. Буляница А.Л., Курочкин В.Е. Свойство независимости точности оценивания от величины зоны нечувствительности в релейном алгоритме // Автоматика и телемеханика. 1999. № 9. С.187-189.
32. Новицкий П.В., Зограф И.А. Оценка погрешностей результатов измерений. Л.: Энергоатом-издат, Ленингр. отд-ние, 1991. 304 с.
33. Мусил Я., Новакова О., Кунц К. Современная биохимия в схемах. М.: Мир, 1984. 216 с.
34. Буляница А.Л., Бурылов Д.А. Частотные критерии устойчивости оценки рекуррентного алгоритма для сигнала постоянного уровня // Научное приборостроение. 1998. Т. 8. № 1-2. С.37-41.
35. Бурылов Д.А. Рекурсивный робастный оценщик информативных параметров сигналов в приборах химического экспресс-анализа. Дис. ... канд. техн. наук. СПб.: ИАнП РАН, 1999. 131 с.
36. Буляница А.Л., Курочкин В.Е., Бурылов Д.А. Анализ асимптотических свойств оценки рекурсивного алгоритма Я. З. Цыпкина с позиций теории автоматического управления // Доклады 3-й Межд. конференции "Цифровая обработка сигналов и ее применения". Москва, Россия, 29 ноября-1 декабря 2000 г. Т. I. С. 17-21.
37. Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1983. 780 с.
38. Аоки М. Оптимизация стохастических систем. М.: Наука, 1971. 216 с.
39. Буляница А.Л., Бурылов Д.А. Исследование сходимости оценки рекуррентного алгоритма для сигнала постоянного уровня // Научное приборостроение. 1998. Т. 8, № 1-2. С. 32-36.
40. Jeffreys H. Theory of probability (2nd eds). Oxford, 1948. 158 p.
41. Буляница А.Л., Бурылов Д.А. Модификация подхода М. Аоки для анализа сходимости оценки алгоритма стохастической аппроксимации в форме Я.З. Цыпкина // Доклады 4-й Межд. конференции "Цифровая обработка сигналов и ее применения". Москва, Россия, 27 февраля-1 марта 2002 г. Т. I. С. 17-19.
42. Буляница А.Л., Курочкин В.Е., Бурылов Д.А. Анализ и практическая реализация алгоритма стохастической аппроксимации в модификации Я.З. Цыпкина // Труды LVI Научной сессии, посвященной Дню радио. М., 2001. Т. 2. С. 323-325.
43. Бурылов Д.А., Курочкин В.Е., Буляница А.Л. Рекурсивный алгоритм оценивания параметров линейных трендов нулевого и первого порядков с контуром определения разладки // Доклады 2-й Межд. конференции "Цифровая
обработка сигналов и ее применения". Москва, Россия, 21-24 сент. 1999 г. Т. 1. С. 204-210.
44. Page E.S. Continious inspection scheme // Biometrika. 1954. V. 41, N 1. P. 141-154.
45. Page E.S. A test for change in a parameter occuring at an unknown // Biometrika. 1955. V. 42. N 4. P.523-525.
46. Никифоров И.В. Последовательное обнаружение изменения свойств временных рядов. М.: Наука, 1983. 254 с.
47. Воробейчиков С.Э. Об обнаружении изменения среднего в последовательности случайных величин // Автоматика и телемеханика. 1998. № 3. С. 50-56.
48. Дарховский Б.С., Бродский Б.Е. Непараметрический метод скорейшего обнаружения изменения среднего случайной последовательности // Теория вероятностей и ее применения. 1987. Т. 32, № 4. С. 703-711.
49. Бывайков М.Е., Ромащев А.А. О робастности в задаче обнаружения изменения параметра сдвига случайной последовательности // Автоматика и телемеханика. 1989. № 7. С. 138-143.
50. Буляница А.Л., Бурылов Д.А. Знаковый критерий определения разладки в последовательности измерений // Научное приборостроение. 1999. Т. 9, № 1. С. 60-64.
51. Ивницкий Д.М., Ситдыков Р.А., Курочкин В.Е., Рейфман Л.С. Высокоскоростной электрохимический анализатор с компьютерной обработкой данных для иммуноферментного анализа // Журнал аналитической химии. 1991. Т. 46, № 6. С. 1239-1244.
52. Поляк Б. Т. Новый метод типа стохастической аппроксимации // Автоматика и телемеханика.
1990. № 7. С. 98-107.
53. Курочкин В.Е., Фельдман Б.Х. Модель нерегулярного кусочно-детерминированного сигнала // Научное приборостроение: Автоматизация научных исследований. Л.: Наука, 1988. С. 6368.
54. Курочкин В.Е., Фельдман Б.Х. Оценка параметра положения сигнала // Научное приборостроение: Автоматизация научных исследований. Л.: Наука, 1988. С. 68-73.
55. Курочкин В.Е. Приборы иммунного и химического экспресс-анализа на основе гибридных методов. Дис. ... докт. техн. наук. СПб.: ИАнП РАН, 1994. 336 с.
56. Кендалл М., Стъюарт А. Статистические выводы и связи. М.: Наука, 1973. 899 с.
57. Левин Б.Р. Теоретические основы радиотехники (в 3-х кн.), кн. 2-я. М.: Советское радио, 1976. 392 с.
58. Алексеенко А.Г., Коломбет Е.А., Старо-дуб Г.И. Применение прецизионных аналоговых микросхем. М.: Радио и связь, 1985. 304 с.
59. Буляница А.Л., Курочкин В.Е., Бурылов Д.А. Реализация процедуры оценивания постоянного сигнала на основе метода стохастической аппроксимации в модификации Я.З. Цыпкина // Радиотехника и электроника. 2002. Т. 47, № 3. С. 343-346.
Институт аналитического приборостроения РАН, Санкт-Петербург
Материал поступил в редакцию 10.04.2002.
THE PROPERTIES AND HARDWARE—SOFTWARE IMPLEMENTATION OF YA. Z. TSYPKIN’S MODIFICATION OF THE STOCHASTIC APPROXIMATION ALGORITHM
A. L. Bulyanitsa, V. E. Kurochkin
Institute for Analytical Instrumentation RAS, Saint-Petersburg
The paper presents the analysis of Ya.Z. Tsypkin's modification of the Robbins—Monro stochastic approximation algorithm to evaluate constant signals in the presence of additive noise with an a priori unknown distribution law, including statistical outliers. The unbiasness property and the efficiency of estimation are evaluated and economical hardware realization of the algorithm is considered. The software version of the algorithm offered allows, in addition to signal estimation, algorithm parameter adjustment and initial approximation of the signal estimate, as well as detection of disorders in the measurement sequence.