Загребаев Методы обработки статистической информации в задачах контроля 2008

ση2	k	j)2 ,
	= 2 ∑ f j (		(3.5.45)

j=0

где – интервал разбиения.

Рис. 3.9. Преобразование непрерывной плотности распределения к кусочно-постоянной плотности распределения

Ограничения, накладываемые на f0 , ..., fk , определяются клас-

сом Φ, причем одним из ограничений обязательно является условие:

k
2 ∑ f j	=1,	(3.5.46)

j=0

которое представляет собой дискретный аналог условия (3.5.2).

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

Характерной особенностью задачи минимизации функции

ˆ

(3.5.44) является необходимость определения c ( f0 , ..., fk ) на каждом шаге итерационного процесса минимизации.

ˆ

Для нахождения c ( f0 , ..., fk ) могут быть использованы абсо-

лютно оптимальные рекуррентные алгоритмы (2.1.72), рассмотренные в предыдущей главе; с учетом кусочно-постоянного характера

функции f д(η) эти алгоритмы можно записать в виде:

251

ˆ	ˆ		−
c	(i) = c (i −1)	+ Γ(i)

fl − fl−1
			∂ψ(i,c )		, (3.5.47а)
fl			∂c
					ˆ
				c	=c (i−1)

Γ(i) = Γ(i −1) −

∂

т

Γ(i −1)

∂ψ(i,c )

ψ(i,c )

Γ(i −1)

∂c

∂c

−

c

=c (i−1)

, (3.5.47б)

IF−1( f0 ,..., fk ) +

∂тψ(i,c )

Γ(i −1) ∂ψ(i,c )

ˆ

∂c

ˆ

∂c

c

=c

(i−1)

c

=c

(i−1)

k −1

( f

j

+1

− f

j

)2

IF ( f0 ,..., fk ) 2 ∑

,

(3.5.47в)

f j

j=0

ˆ

ˆ

,

Γ(0) = λI ,

λ >>1.

c (0)

= c0

Число l в формуле (3.5.47а) определяется номером интервала, в

который попадает невязка процесса идентификации:

l = int

	ˆ	на i-м шаге рекуррентного
ε(i, c(i −1))
	ˆ
	ε(i,c (i −1))		+1.	(3.5.47г)

Параметр i изменяется от 1 до N, где N – число измерений входов и выхода объекта, используемых при определении оптимальной плотности распределения.

Общая схема процесса нахождения оптимальной функции плотности распределения (функции потерь) может быть представлена последовательностью действий.

1. Задаем начальные приближения:

ˆ		ˆ		Γ(0) = λI ,			λ >>1;
c (0) = c0 ;
			f00		, 0 ≤ η <		;
				0	,	≤ η < 2 ;
f д	0	(η) =	f
			L
			...
				0	, (k −1)		≤ η < k ,
			fk

где f д0 (η) удовлетворяет условию (3.5.46), f д0 (η) Φ ; l = 0 , l – номер итерации в процессе минимизация функции (3.5.44).

252

2. Используя рекуррентные формулы (3.5.47), находим оценку

ˆ l

при допущении, что

f (η) = f

дl

(η) .

параметра c

3. Вычисляем значение минимизируемой функции по формуле

(3.5.44):

( f jl+1 − f jl )2

k−1

1

J l = ∑

;

f jl

−1

ˆ l

2

l

j=0

tr

A

(c

,(ση )

)

k

(ση2 ) l = 2 ∑ f jl ( j)2 .

j=0

4. Рассчитываем

новое

значение

плотности

распределения

( f д)l+1 :

f0l

+ 0 ; 0 ≤ η ≤ ;

( f д)l+1(η) = ...

−1) ≤ η ≤ k ,

fkl + k ; (k

Тогда 0 , ...,

k определяются выбранным методом поиска экстре-

мума функции многих переменных (3.5.44) при ограничениях:

k

( f д)l+1(η) Φ ,

2 ∑ f jl+1 =1 .

j=0

5. Рассчитываем новое значение критерия:

k −1

f jl++11 − f jl+1

1

J l+1 = ∑

,

f jl+1

−1

ˆ l

2

l+1

j=0

tr

,

A

c

(ση )

(ση2 )l+1 =

k

2 ∑ f jl+1(

j)2 .

j=0

6.Проверяем условие окончания итерационного процесса:

J l+1 − J l
J l+1	≤ δ , δ <<1.

Если условие выполняется, то процесс поиска оптимальной плотности распределения заканчивается; в противном случае переходим к п. 2 итерационного процесса, присвоив l = l +1 .

253

Необходимо отметить, что минимизируемая функция (3.5.44) является чаще всего многоэкстремальной, что приводит к опасности нахождения локального минимума. Это необходимо иметь в виду при решении задач подобного типа.

3.5.5. Определение функции потерь для регрессионных объектов. Алгоритм Хубера

Как уже отмечалось выше, задача определения функции потерь в общем случае не может быть решена в явном виде. Однако для некоторых классов при идентификации линейного регрессионного объекта вида

y(i) = b1u1(i) +... + bnun (i) + η(i)

(3.5.48)

эта задача может быть получена в явном виде. Это связано с тем, что нормированная информационная матрица для линейных регрессионных объектов не зависит ни от оцениваемых параметров, ни от дисперсии помехи.

Учитывая это, минимизируемый функционал (3.5.43), а именно:

∞	f ′	2	(η)				1
J ( f (η)) = V( f (η)) = ∫				dη					,
					tr (A
−∞	f (η)					−1		2
							(c ,ση ))

f (η) Φ

может быть заменен более простым функционалом

∞	f ′	2	(η)
J ( f (η)) = ∫				dη , f (η) Φ ,	(3.5.49)
−∞	f (η)

который представляет собой информацию Фишера.

В том случае, когда помеха η имеет распределение, принадлежащее классу α-загрязненных распределений (см. классы Φ3 , Φ4 ,

Φ5 в разд. 3.2.1), для нахождения оптимальной плотности распре-

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

деления для Р-объектов является плотность распределения, которой соответствует минимальная фишеровская информация.

Интересно отметить, что работы Хубера в истории развития робастных оценок были пионерскими. Ему удалось обосновать принцип минимакса, опираясь только на выводы теории вероятностей и

254

математической статистики (в основном на центральную предельную теорему и ее следствия).

Итак, рассмотрим задачу идентификации коэффициентов линейного регрессионного объекта (3.5.48) в предположении, что распределение f (η) помех принадлежит классу α-загрязненных

распределений:

Φ ={ f (η) : f (η) = (1− α)h(η) + αg(μ)} ,	(3.5.50)
где h(α) – известная плотность распределения; g(η)	– произволь-

ная неизвестная плотность распределения; α – вероятность появления «выброса» с распределением g(η) , α удовлетворяет условию

0 ≤ α ≤1 .

Результат, полученный Хубером, получил название «Теорема Хубера».

Теорема Хубера. Пусть h(η) – дважды непрерывно дифференцируемая плотность распределения, такая, что (−ln h(η)) – выпуклая вниз функция. Тогда АМКО линейных регрессионных объек-

тов – V(F, f )	имеет седловую точку, т.е. существует плотность
распределения	f *(η) = (1− α)h(η) + αg *(η) и функция	F *(η) =
= −ln f *(η) , такие, что
	V(F*, f ) ≤ V(F*, f *) ≤ V(F, f *) .	(3.5.51)

Далее, пусть η0 и η1 ( η0 < η1 ) – концы интервала (один или оба конца могут быть бесконечными), где

			′
		h (η)				≤ k ,	η [η0 , η] ,				(3.5.52)
			h(η)
	и k , α, η0 и η1 связаны соотношением:
			η					h(η0 ) + h(η1)
	(1− α)−1 = ∫1 h(η)dη+									.	(3.5.53)
			η0					k
	Тогда плотность f
			* (η) имеет вид:
	(1− α)h(η0 )exp(k(η− η0 )),								η ≤ η0 ;
			− α)h(η),				η0 ≤ η ≤ η1;				(3.5.54)
	f *(η) = (1
	(1		− α)h(η )exp(−k(η− η )),						η ≥ η .
			1				1		1

255

Вначале покажем, что f * (η) , определяемая формулой (3.5.54),

удовлетворяет всем требованиям, предъявляемым к плотности распределения, и принадлежит классу (3.5.50). Затем докажем, что эта плотность f * (η) действительноявляетсяоптимальнойнаклассе(3.5.50).

Итак, рассмотрим условия второй части теоремы. Прежде пока-

жем, что	f * (η) удовлетворяет условию полноты интеграла:
∞	η				η
∫	f * (η)dη = ∫0	(1−α)h(η0 )ek(η−η0 ) + ∫1				(1−α)h(η)dη+
−∞	−∞				η0
	∞
	+ ∫ (1− α)h(η1)e−k (η−η1)dη =
	η1
		1		η				1
	= (1− α)h(η0 )		+ (1− α) ∫1 h(η)dη+ (1− α)h(η1)						=
								k
		k		η0
				η1	h(η ) + h(η )
	= (1−α)			∫ h(η)dη+	0	1	.
					k
			η0

По условию теоремы (3.5.53) выражение, стоящее в квадратных

скобках, равно	(1− α)−1 , и, следовательно, получаем
∞
∫ f * (η)dη =1, т.е.	f * (η) удовлетворяет условию полноты инте-
−∞

грала.

Покажем теперь, что и g * (η) , соответствующая f * (η) , также удовлетворяет условию полноты интеграла, а именно:

∞
∫ g * (η)dη =1.	(3.5.55)

−∞

Используя формулы (3.5.50) и (3.5.54), найдем выражения для g * (η) :

		(1− α)		h(η )ek(η−η0 ) − h(η)					,	η ≤ η ;
		α
					0					0
										η0 ≤ η ≤ η1; (3.5.56)
	g * (η) = 0,
		1− α				−k(η−η1)
		α	h(η1)e				− h(η)	,		η ≥ η1.

256

Интегрируя (3.5.56) на бесконечных пределах, нетрудно убедиться в справедливости выражения (3.5.55). Покажем, в заключение, что g * (η) неотрицательная функция. Докажем это для η< η0

(при η > η1 доказательство аналогичное).

Рис. 3.10. График функции –ln(f(η))

По условию теоремы − ln(h(η)) – выпуклая вниз функция. Следовательно, график этой функции лежит выше касательной, проведенной в любой точке, в том числе и в точке η0 (рис. 3.10).

Как известно, наклон касательной определяется производной функции в данной точке. Таким образом, учитывая условие выпуклости вниз функции −ln h(η) , можно записать:

−ln h(η) ≥ −ln h(η0 ) + (−ln h(η))′ η=η0 (η− η0 ) =

= −ln h(η						) −	h′(η0 )	(η− η	) ,	(3.5.57)
0							h(η0 )	0
						η≤ η0 .
По условию теоремы (3.5.52)
	′
	h (η)				≤ k , η [η0 , η] .
	h(η)

Учитывая, что h(η) – дважды непрерывно дифференцируемая

плотность распределения, можно заключить: h′(η0 ) = k h(η0 )

257

Тогда, подставляя последнее равенство в формулу (3.5.57), получим:

−ln h(η) ≥ −ln h(η0 ) − k(η− η0 ) ,

или

		′	≥ 0 .
	ln h(η0 ) + k(η− η0 ) − ln h (η)
Последнее неравенство эквивалентно выражению
	h(η0 )ek(η−η0 ) − h(η) > 0 .		(3.5.58)
Согласно формуле (3.5.56)
	1−α	(h(η0 )eh(η−η0 ) − h(η))= g * (η) .
	α

Сравнивая это выражение с неравенством (3.5.58), и так как 0 ≤ α ≤1 , то нетрудно заключить, что g * (η) ≥ 0 при η ≤ η0 , что и

требовалось показать.

Таким образом, g * (η) определена «вполне корректно» соот-

ношениями (3.5.56) в смысле требований, которым должна удовлетворять плотность распределения.

Покажем теперь, что f * (η) и соответствующая ей функция потерь F * (η) = −ln f * (η) действительно определяют седловую точку, т.е. удовлетворяют неравенству

	V(F * (η), f (η)) ≤ V( f * (η)) ≤ V(F(η), f * (η)) (3.5.59)
для	f (η) ,	~	~	~
		f (η) Φ , F(η) = −ln f (η) .

Правое неравенство выражения (3.5.59), очевидно, вытекает из свойств АМКО. Таким образом, необходимо доказать, что

V(F * (η), f (η)) ≤ V( f * (η))

(3.5.60)

для f (η) Φ .

Учитывая, что нормированная информационная матрица не зависит от ση2 и c , последнее неравенство можно переписать следующим образом:

∞

′

2

∫

f * (η)

f

′

(η)dη

∞

′2

(η)

f * (η)

∫

f *

dη ≤

−∞

.

(3.5.61)

f * (η)

∞

′

2

−∞

∫

f * (η)

f (η)dη

−∞

f * (η)

258

Подробное доказательство этого неравенства приведено в рабо-

те [4].

Таким образом, показали, что оптимальная на классе (3.5.50) плотность распределения существует и определяется по правилу

(3.5.54).

Запишем соответствующую ей оптимальную функцию потерь:

−ln[(1− α)h(η0 )] − k(ε − η0 ),			ε ≤ η0 ;
	− α)h(η0 )],	η0 ≤ ε ≤ η1;		(3.5.62)
F * (ε) = −ln[(1
−ln[(1	− α)h(η )] + k(ε − η ),		ε ≥ η .
	0	1	1

3.5.6. .Идентификация параметров регрессионного объекта при α-загрязненном нормальном распределении помехи

Пусть известно, что помеха принадлежит классу Φ3 – приближенно нормальных распределений:

	f (η) = (1− α) fN (η) + αg(η) ,									(3.5.63)
	где f N – нормальное распределение;
		1				1	2
	fN =								,	(3.5.64)
		2	exp		−	2	η
		2πση				2ση

α– вероятность появления «выброса» с распределением g(η) . Примерный вид α-загрязненных шумов изображен на рис. 3.11.

Рис. 3.11. Примерный вид α-загрязненных шумов

259

Учитывая результат, полученный в предыдущем разделе, запишем оптимальную на классе плотность распределения:

				1										1
																		η2
(1−α)								exp −												exp[k(η−η )], η<η ;
				2 1/2										2σ2
																			0						0		0
			(2πση)												η
				1										1
																η2 ,
f *(η) = 1	−α)						exp			−												η0 ≤η≤η1;					(3.5.65)
				2	1/2
										2σ2
		(2πση)												η
				1										1
(1−α)								exp −										η2		exp[−k(η−η )], η>η .
				2 1/2										2σ2
																			1						1		1
			(2πση)												η
Соответствующая функция потерь будет иметь вид
														1								1
			−ln		(1	− α)													+				η2		− k(ε − η	),	ε ≤ η ;
															2							2
																								0	0		0
											2πση										2ση
			−ln											1							1
F *(ε) =					(1	− α)													+			σ2	ε2 ,		η ≤ ε ≤ η		;
																					2
															2							η			0	1
											2πση
														1							1
			−ln		(1	− α)													+			σ2	η2		+ k(ε − η ),		ε ≥ η .
															2
																					2	η		1	1		1
											2πση
Очевидно, записанная функция потерь будет эквивалентна бо-
лее простой функции потерь:
										1				η2		− k(ε − η
																								),	ε ≤ η ;
								2ση2						0								0			0
										1
			F *(ε) =											ε2 ,				η			≤ ε ≤ η ;						(3.5.66)
									2ση2										0					1
										1				2
																+ k(ε − η1), ε ≥ η1.
									2σ2					η1
											η

На рис. 3.12 приведен примерный вид функции потерь (3.5.66). Видно, что на отрезке [η0 ,η1] функция потерь квадратичная (метод

260