Доклады БГУИР
2012 № 5 (67)
УДК 519.246
ОБОБЩЕННОЕ РАСПРЕДЕЛЕНИЕ SECHk
А.В. ОВСЯННИКОВ
Белорусский государственный университет пр. Независимости, 4, Минск, 220030, Беларусь
Поступила в редакцию 14 марта 2012
Исследуется модель обобщенного закона распределения Sechk и ее применение в статистической радиотехнике. Приводятся основные свойства и статистики распределения. Получены алгоритмы генерирования случайных чисел и алгоритмы моделирования стохастических процессов с одномерной плотностью Sechk Приведены алгоритмы точечной оценки параметра сдвига при известных и неизвестных параметрах масштаба, а также алгоритмы последовательной интервальной оценки параметра сдвига. Результаты моделирования подтверждают теоретические расчеты.
Ключевые слова: распределение, нелинейное преобразование, алгоритм, робастность, генерация, точечная и интервальная оценка, параметр сдвига, параметры масштаба.
Введение
В литературе, затрагивающей вопросы статистической обработки данных (результатов наблюдений) для математического описания разнообразных помеховых процессов в каналах связи, используются стохастические дифференциальные уравнения (СДУ) в которых детерминированная функция /(Ъ) определяется одномерной плотностью вероятности (ПВ)
dЪ(t) / dt + / (Ъ) = уп^), / (Ъ) = (Ъ) / 2, (1)
где 2(Ъ) = ^ 1пР(Ъ)/ dЪ - нелинейное преобразование над процессом ) с плотностью [1]
Р(Ъ) = (С / )ехр(-2 [Ъ/(х^х / Ъ), у = у/2Ъ / N , N- односторонняя спектральная плотность
белого шума п(0 с нулевым математическим ожиданием и дельтообразной корреляционной функцией, коэффициент диффузии Ъ=сош^ С - константа.
Для описания широкого класса помех, распределенных по частоте, времени или пространству (флуктуационные), а также сосредоточенных во времени или пространстве (импульсные), используется обобщенная экспоненциальная модель [2]:
Р(Ъ) = 2е (у)г(1/ V) ехр(-1 Ъ / С (V)Г), С (V) = а [Г(! / V) / Г(3 / V)]12, (2)
где с2 - дисперсия (мощность) помехи, V > 1/2 - параметр, принимаемый в различных ситуациях за V = 1/2, V = 1 (лапласовская ПВ) или V = 2 (гауссовская ПВ).
Модель вида (2) при несомненных достоинствах экспоненциальной формы содержит изменяющийся параметр V, что не всегда удобно при построении алгоритмов и практических схем демодуляторов, поскольку требует их структурной перестройки. Например, нелинейное преобразование с V = 2 имеет вид 2 (Ъ) = Ъ / с2 (линейный вход демодулятора), а в случае V = 1 - жесткий ограничитель в составе демодулятора 2(Ъ) = V2sign(Ъ) / с .
В практике статистических расчетов в ряде случаев для вероятностного описания помехи используют модель распределения с в -загрязнением (приближенно нормальное распределение, допускающее удлиненные хвосты) Р(£) = (1 -в)+ вР1(£), где - гауссовская ПВ, Р1(^) - произвольная ПВ. Если Р1(^) лапласовская ПВ, то нелинейное преобразование имеет вид 2= {Дsign(^)/с2, при |^|>Д; / с2, при |^|<д}, где Д - ширина зоны линейной части нелинейного преобразования.
В этой связи в [3] получено распределение имеющее ПВ
P(^) = аС(k)Sechk(а£),k = ß/а ,а,ß> 0,
(3)
где С^) = Г((1 + k)/2)Г1(k /2)/\/л , Г(х) - гамма-функция.
Такая модель помехи (одномодальная, двухпараметрическая) позволяет учесть медленные нестационарные изменения интенсивности и амплитуды процесса ^). При параметре
а ^ 0 распределение стремится в пределе к гауссовскому, а при а ^ да к лапласовскому (рис. 1). Нелинейное преобразование для ПВ (3) имеет вид мягкого ограничения 2 = рШ(а^) (рис. 2), это позволяет получать алгоритмы обработки без структурной перестройки схемы, используя только адаптивную настройку параметров 2.
Щ)
........ <1-1(1 1
/
<(=11.5 7
_____________ „-•. 1
_1
Рис. 1. Распределение (3), 1 - а=10; 2 - а=0,5; 3 - а=0,1
Рис. 2. Нелинейное преобразование, 1 - а=10; 2 - а=0,5; 3 - а=0,1
Цель работы - обобщить сведения о свойствах, статистике и информационных характеристиках распределения с ПВ (3); получить алгоритмы генерации случайной величины и процесса ^), оценки параметра сдвига при известных и неизвестных параметрах масштаба, интервальной оценки параметра сдвига.
Основные свойства и статистики распределения Sechk
Распределение с ПВ (3) с точки зрения построения робастных алгоритмов и устройств обработки является распределением в классе [3] = |p : j^2qP(^)d C2q}, где
C2q = f (У2,У45 -,У2q)®2q, У2q - безразмерные кумулянтные коэффициенты (у2 = 1, у4 - коэффициент эксцесса), q = 1,2... - индекс. Класс распределений шире, чем классы: ={P:P(0) >5> 0} («наихудшая» ПВ-лапласовская) и = |p : P(^)d^<с2} («наихудшая» ПВ-гауссовская), которые являются предельными для .
Распределение с ПВ (3) включает в себя, как частный случай, распределение Чампер-науна [4] P(^) = aSech(a^) / % (здесь в (3) к = 1) и т.н. распределение «Sech2» с
ПВP(^) = aSech2(a^)/2 (здесь в (3) к = 2). Таким образом, распределение (3) можно определить как обобщенное «Sechk» распределение. Функция распределения имеет вид
>-\к
F © = — С (к )e к
ак^ F 2-М
к, к 2
1+к
2
,-e
2а^
(4)
где 2Fl - обобщенная гипергеометрическая функция. Для практических расчетов удобно воспользоваться аппроксимацией (4) в виде Fa = + th + Ь1ксс / 2 . При такой аппроксимации на интервале у = а^е[-102;102], к е[0,05;1] с коэффициентами а1 =-4,10 -10-3, Ь1=64,23-10-2, с1 = 92,76-10-2 среднеквадратическая ошибка D(AF = ¥а -F) не превосходит величины 5 -10-5, а в случае к е[1;20], а1 =-18,98 -10-2, Ь1=79,37 -10~2, с1 = 56,23 -10~2 эта ошибка D(AF) < 5 -КГ6.
Характеристическая функция ПВ (3)
Ф(^) = 2к -1 С(к )е
ж(г + jаk)
2а
В
-1Лк -а ),1 - к
2 а
жг
+ е а В
-^(к),1 - к 2а
(5)
где В(z, а,Ь) = [ ta 1(1 -1)Ь 1 dt - неполная бета-функция. » 0
В табл. 1 в качестве примера приведены ПВ (3), функции распределения (4) и характеристические функции (5) для параметра к = 3,4.
Таблица 1. ПВ (3), функции распределения и характеристические функции для &=3,4
Параметр к к = 3 = 4
ПВ Р(%) 2а8еЛ3(а^)/ ж 3аБеЛ4(а^)/ 4
Функция распределения F© - 0,5 1 (2аг^(Ш(а£, / 2)) + ЗесКа^)Ш(а^)) ж 1 2 4(2 + Sech2(а^))th(а^)
Характеристическая функция Ф© t2 жг (1 + 0- ) а 2а жГ ,, жГ —(4 + 2 )csch(— ) 8а а 2а
Математическое ожидание и асимметрия ПВ (3) равны нулю, а дисперсия и эксцесс (рис. 3) определяются функциями, зависящими от параметров р и к :
0(к, р) = ^кр 4 Fз
к к к к
к к к — +1— +1— +1 2 2 2
-1
(6)
у (к) =
, 6 F5
3к2-к 6 5
к к к к к к 2'2'2'2'2'
к л к л к л к л к л — +1— + 1,— +1, — +1, — + 1 22222
-1
С (к)
- 3
4 Fз
к к к к 2'2'2'
к к к — +1— + 1, — +1 222
,-1
(7)
где 4 F3, 6 F5 - обобщенные гипергеометрические функции. Функция эксцесса у (к) при значениях (к ^ 0) стремится к трем (лапласовская ПВ), а при а ^ 0 (к ^ да) к нулю (гаус-совская ПВ), что соответствует отмеченному выше свойству предельных распределений для обобщенного распределения. Также можно отметить, что не менее, чем у 50 % всех используемых в практических задачах распределений [4], коэффициент эксцесса находится в пределах [0,...,3]. Ввиду громоздкости точных формул (6), (7) их можно аппроксимировать упрощенными зависимостями Da (к, Р) = (а2 + Ь2кС2) / р2, уа (к) = 3 / (а3 + Ь3кс), где для интервала значений ке[0,1] параметры а2 = 2, Ь2 = 0,47, с2 = 1,65, а3 = 1, Ь3 = 0,5, с3 = 1,7 и среднеквадратическая ошибка не превосходит величины 10-5 (для Da) и 3• 10-5 (для уа). В интервале значений к е [1,10] параметры а2 = 1,29, Ь2 = 0,99, с2 = 0,99, а3 = 0,8, Ь3 = 0,7, с3 = 1,3 и среднеквадратическая ошибка не превосходит 3 -10-3 (для Da) и 10-4 (для уа). Заметим, функции (6), (7) имеют точные значения D(0, Р) = 2, у(0) = 3 , у(1) = 2, у(2) = 1,2 .
Количество информации Фишера рассматриваемой ПВ
2
1. (к, Р) = | 2 2(.)Р(.)ё . = р2/(1 + k).
Коэффициент подавления помехи имеющей ПВ (3), ц(к) = D(k,Р) 1. (к,Р) и находится в пределах ц(к) е[1,...,2] (рис. 4).
определяется
(8) формулой
>т «
1 1
; ; 1 1
>; | |
1 1 ; ;
1 1 ; ;
:...............;...............
т
............
............ ............ ............
Рис. 3. Коэффициент экцесса
Рис. 4. Коэффициент подавления
Генерация случайных величин и процессов
Аналоговый алгоритм генерирования процесса с одномерной ПВ (3) непосредственно следует из (1) ё .(О/ Ш + Ьр^(а.)/2 = уп(0 и легко реализуется, например, в системе Simulink МаАаЬ. При реализации аналогового алгоритма следует обеспечить выполнение соотношения
Т. «М(.)/ = 2/ (Ь1.) >> Тп, где Т. - постоянная времени моделируемой системы, Тп -постоянная времени процесса п^).
В случае дискретного варианта алгоритма численного моделирования СДУ (1) для симметризированной формы уравнения можно использовать метод Рунге-Кута второго или четвертого порядка, а для формы Ито - метод Эйлера. Так оптимальная разностная схема, при понимании СДУ (1) в форме Ито, имеет вид (А << 1):
.+1 -АЬр&(а.,)/2 + Ауи., (9)
где А - шаг квантования, п1 - стандартный дискретный белый шум (независимые гауссовские случайные величины с нулевым средним и М(ninj) = 8.., где 8... - символ Кронекера). Глобальная среднеквадратическая погрешность схемы (9) определяется формулой [5]:
с< А|
.10
(]>[(£К»2]*) = (А)3'2(ь/2)2т , (10)
здесь Т = Ат , т - объем выборки. Неравенство (10) позволяет определить шаг квантования А .
Для генерирования случайных чисел с ПВ (3) можно применить стандартный подход с использованием обратной функции распределения = F) (определяемой из ¥а (.)) и генератора равномерно распределенных чисел ц1.
Более точный результат можно получить с использованием современных вычислительных средств на основе модифицированного алгоритма ступенчатой аппроксимации функции распределения. Алгоритм состоит в следующем.
1. Задать параметры k, а.
2. Сгенерировать равномерную сетку по У на интервале [-с / k,с / k], где параметр с = 10...20 обеспечивает покрытие сеткой не менее 99,995 % всей области значения функции F(.) (4) при к < 10 . В этом случае ширина ступеньки аппроксимации функции распределения будет определяться шагом сетки АУ .
3. Вычислить F(У..) по формуле (4) - = 1,I, I = 1е^Ш(У).
4. Сгенерировать равномерно распределенное число ц1, где i = 1, m , m > l.
5. Вычислить индекс j из условия j = arg min I ^ - F(Y) I.
ми] L J
6. Вычислить для индекса i число . = Yj / а .
7. Вернуться к шагу 4 если i < m .
Абсолютная погрешность метода Д. = max IF(Yj+1) - F(Y,) I / 2а = Д YF' (0) / 2а , где ве-
je[1,l ] ^ j j ' "
личины F"(0) = C(k), ДY = 2c / kl = (20...40)/ kl. Заметим, при равномерной сетке по zj = F(Yj), zj е [0,1], j = 1, J, (неравномерная сетка по Yj) абсолютная погрешность определяется величиной Д. = Дг / 2а = (2аJ)—1. Для получения случайной величины ... = Yj / а необходимо решать одним из численных методов обратную задачу Yj = Fzj) для каждого конкретного значения k.
Оценка параметра сдвига при известных параметрах масштаба
Рассмотрим задачу оценки параметра сдвига 0 по имеющейся однородной независимой выборке Г = [r1 ... rm ] , принадлежащей распределению с плотностью P0(4)=P(r 10), где 0 е ® и ® - интервал на действительной оси. При фиксированном объеме выборки r выбо-
Rr>m
i r _ _ _ _ = R , на котором зада-
на плотность P0 (4), r е R . В этом случае в (3) следует положить = r — 0, i = 1, m .
Частная функция потерь имеет вид В(.) = — lnP(.i) = C — kln^ech^.)). Метод максимального правдоподобия и использование выборочного эмпирического функционала
W3 (.,) = (1/ m)^ В(. ) приводит к оценке 0* = arg max W3 (ri — 0) и, следовательно, получаем i r неявный интегральный робастный алгоритм
m ) Л = fi Z(.(i)) ^ Jt| 0=0* J^fi ^^ —9))Jt|0=0- = 0. (11)
Применение стохастической градиентной процедуры к последнему уравнению позволяет перейти к последовательным робастным оценкам
0*, =0* — 1+1 i
M
( 52ln P(.)1
Ш2
= 0* + KZ) Л =0* + ^th^r —0*)), (12)
d 0 d0 iI%
= —1/ il..
где К = [м(<521пРф/502)]_1 =- М((51пРф/50)2)
В табл. 2 приведены качественные характеристики алгоритма (12) для случая 0 = 0,1, 0О = 0, k = 1, ОСПвх, ОСПвых - отношение сигнал/помеха до и после применения алгоритма (11).
Таблица 2. Качественные характеристики алгоритма (8)
Длина выборки Параметр а 0,2 0,5 1 2 5
ОСПвх, дБ -37,8 -29,9 -24,0 -17,7 -10,1
n=50 ОСПвых, дБ -22,6 -18,3 -8,4 -4,6 8,1
n=500 ОСПвых, дБ -13,5 -9,7 -0,1 2,0 15,5
n=5000 ОСПвых, дБ -7,9 -1,2 8,6 10,3 22,9
Оценка параметра сдвига при неизвестных параметрах масштаба
Общая постановка задачи аналогична рассмотренной выше, однако параметры масштаба а и в (а и к) являются неизвестными. Интегральный алгоритм оценки получаем из уравнения )/ йЛ = 0, где Л* = [а*,к ]. Тогда система уравнений интегральных оценок примет вид
1
|а=а*,к=к*
= 1,
к
-I ¥(-) — 21 2 2
)1 — - | Ы^ес^а^))^ = 0,
Т о |а=а*,к=к*
(13)
где ¥(х) = й1пГ(х) / йх - логарифмическая производная гамма-функции. Полученная система
* *
не позволяет найти явные выражения для оценок а , к , однако процедурами стохастической аппроксимации можно получить систему алгоритмов последовательных оценок. В сочетании с уравнением (11) система (13) принимает замкнутый вид.
Для приближенного решения систему (13) можно упростить, используя аппроксимацию сложных нелинейных функций более простыми. Так, функция /п^еЛ(х)) аппроксимируется
зависимостью — | х|, при этом среднеквадратическая ошибка такой аппроксимации на интервале х е [—103,103] не превосходит единицы. Функцию ^((1 + к)/2) — ^(к /2) на интервале к е[102,102] можно аппроксимировать зависимостью 1,94/ к со среднеквадратической ошибкой аппроксимации не превышающей 10—2.
Более простой и надежный (при достаточном объеме выборки, времени наблюдения), легко реализуемый программными средствами способ идентификации состоит в следующем.
1. Оценить выборочную дисперсию D .
2. Оценить выборочный коэффициент эксцесса у*. Если у(к )ё[1,...,3], то необходимо увеличить объем выборки.
3. Вычислить оценку параметра к* из формулы у (к ) = у*.
4. Вычислить оценку параметра в* из формулы D(k*, Р) = D*.
5. Вычислить величину а* = Р* / к *.
В качестве функций D(k,Р), у (к) можно использовать их аппроксимации Da, у а. В табл. 3 приведены качественные характеристики алгоритма (11) с оценкой параметров масштаба по приведенному выше алгоритму. Исходные условия соответствуют табл. 2.
Таблица 3. Качественные характеристики алгоритма (8) с оценкой параметров масштаба
Длина выборки Параметр а 0,2 0,5 1 2 5
ОСПвх, дБ -37,8 -29,9 -24,0 -17,7 -10,1
п=50 ОСПвых, дБ -23,8 -18,7 -9,2 -4,9 7,9
п=500 ОСПвых, дБ -13,7 -9,8 -0,9 1,1 15,4
п=5000 ОСПвых, дБ -8,1 -1,3 8,2 9,4 22,3
Интервальная оценка параметра сдвига
Для получения интервальной оценки используем тот факт, что величина —58(4) / 59 асимптотически нормальна с параметрами (0, . I^) (Г. Крамер, 1946), т.е. если . то
—(. /^)—1/2 68(4)/ 59 => N {0,1} . Алгоритм последовательных интервальных оценок
е: -агс^^^ + G;-1)/а<9, <е* + агс^^^ -а, (14)
G; = + Ш(а(г,— е*)),
где . - квантиль распределения Стьюдента. Минимальный объем выборки, который обеспе-
2 — 1 —2
чивает заданную точность оценки Де , определяется неравенством т < 1р (к +1) Л (аДе).
2 — 1 —2
Если величина а Де << 1, неравенство можно упростить т < 1р (к +1) (аДе) .
Заключение
В заключении необходимо отметить теоретическую и практическую ценность рассматриваемого распределения, которая заключается в следующем. Модель ПВ (3) может использоваться для описания импульсных промеховых процессов с изменяющейся интенсивностью и амплитудой выбросов, при этом вместо двух различных распределений описывающих поведение процесса вблизи математического ожидания и на удлиненных хвостах (хьюберовский класс приближенно нормальных распределений, распределения с 8 -загрязнением) используется одно. В связи с этим модель (3) интересна и может найти применение не только в области статистической радиотехники, но и в разнообразных задачах статистической обработки результатов наблюдений, возникающих в финансово-экономической сфере, метрологии, социологии и т.д.
Для модели (3) в отличии от (2) и модели с 8 -загрязнением всегда существует дР(г|е)/Шг и непрерывна V/ = 1,2,... . ПВ (3) является «наихудшей» в обобщенном классе распределений и дает возможность получать робастные алгоритмы оценивания, фильтрации. Модель (3) позволяет, не меняя структуру демодулятора, использовать адаптивную процедуру настройки параметров а, в и, следовательно, применять робастно-адаптивные алгоритмы работы демодулятора. Практическая реализация нелинейного преобразования X(.) может быть осуществлена на базе дифференциального усилителя-ограничителя.
GENERALIZED DISTRIBUTION SECHk
A.V. AUSIANNIKAU
Abstract
Model of generalized law of distribution of Sechk and its researched application in statistical radio engineering. Base attributes and statisticians of distribution are given. Algorithms of point estimation of parameter of shift for known and unknown parameters of scale, and also algorithms of a sequential interval estimation of parameter of shift are resulted.
Список литературы
1. Тихонов В.И., Кульман Н.К. Нелинейная фильтрация и квазикогерентный прием сигналов. М., 1975.
2. Левин Б.Р. Теоретические основы статистической радиотехники. М., 1989.
3. ОвсянниковА.В. // Радиотехника. 2011. №3. С. 85-89.
4. Вадзинский Р.Н. Справочник по вероятностным распределениям. СПб., 2001.
5. Никитин Н.Н., Разевиг В.Д. // Журнал вычислительной математики и математической физики. 1978. Т 18, №1. С. 106-117.