Содержание


Моделирование стохастических процессов на языке Perl

Часть 2. Моделирование Марковских процессов на Perl

Comments

Серия контента:

Этот контент является частью # из серии # статей: Моделирование стохастических процессов на языке Perl

Следите за выходом новых статей этой серии.

Этот контент является частью серии:Моделирование стохастических процессов на языке Perl

Следите за выходом новых статей этой серии.

1. Моделирование дискретной марковской цепи при формировании цифровых изображений

Многие природные явления и технические процессы характеризуются как марковские стохастические процессы. Одним из таких процессов является формирование изображений в системах наблюдения и передачи видеоданных. Спектр мощности большинства реальных сюжетов сосредоточен в низкочастотной области1, что свидетельствует о существенной корреляционной зависимости значений яркости соседних элементов изображения в декартовом пространстве (x,y). Это свойство, в частности, используется для сжатия цифрового потока видео данных, когда в канал связи передается разностный сигнал яркости соседних элементов изображения. Или, иначе говоря, в последовательности значений яркости элементов реальных сюжетов, значения яркости предшествующих элементов оказывают влияние на значения яркости последующих.

Цифровые изображения дискретны по своей сущности и характеризуются дискретными последовательностями строк и столбцов. Квантование сигналов изображения по амплитуде, в аналого-цифровых преобразователях (АЦП), довершает картину представления цифрового изображения в виде марковской цепи. Случайный характер реальных изображений обусловлен многими факторами, начиная от сюжета до шумов квантования АЦП. Теория цепей Маркова является естественным инструментом для анализа и описания таких систем.

На рис.1 показано тестовое изображение "Lenna", изменения значений яркости B(x,y) в котором представляют собой случайный процесс Xn.

Рис. 1. Тестовое изображение
Рис. 1. Тестовое изображение

Процесс изменения случайной величины Xn в пространстве координат охарактеризуем как переходы из состояний яркости i в состояния яркости j. Диапазон изменения состояний яркости определяется числом уровней квантования M АЦП, т.е. 0<=B(x,y)<M. Как отмечалось в первой статье цикла, матрица переходных вероятностей P(i,j) является полным описанием Марковской цепи случайного процесса. Просканируем изображение B(x,y) и определим частоту появления элементов изображения с яркостью j, при условии, что предшествующий элемент имел яркость i. Полная совокупность переходов образует матрицу переходов или частость N(i,j):

Рис. 2. Матрица переходов N(i,j)
Рис. 2. Матрица переходов N(i,j)
Рис. 2. Матрица переходов N(i,j)

Условную вероятность перехода из состояния i в состояние j обозначим через pij = Ф(Xn=j|Xn-1=i). Совокупность всех переходов образует матрицу переходных вероятностей или цепь Маркова для дискретного процесса:

                     |p11 p12 ... p1m|
                     |    .          |
          P(i,j) =   |      .        |
                     |        .      |
                     |           .   |
                     |           pmm | .

Размерность матрицы P равна m2, где m - число уровней квантования. Для стохастической матрицы должно быть выполнено условие, когда все элементы матрицы pi,j≥0, и для каждой i строки сумма вероятностей pi,j равна единице:

           M-1
           Σpi,j  =  1.                               (1)
           j=0

Вариант программной реализации, процедуры вычисления матрицы условных вероятностей pi,j , представлена ниже. Для расчета матрицы переходов N(i,j) необходима предварительная процедура ввода в программу исходных данных, а именно файлов изображений. Для работы с графикой и изображениями в Perl хорошо зарекомендовала себя графическая библиотека GD. Рассмотрение модуля GD выходит за рамки настоящей статьи и поэтому в примере кода, для чтения файлов изображений, воспользуемся встроенными функциями. Второе упрощение касается формата файлов изображения - будем использовать формат растровой графики BMP grayscale с 8-ми разрядным представлением яркости. Размер изображения по горизонтали $Xw и по вертикали $Yh имеет формат двойного слова.

 	         open (LIST,"images/lenna_256.bmp");
	         binmode(LIST);
	         @buffer = <LIST>;
	         close (LIST);
	         $buf = join("",@buffer);
                     @BMP=unpack("C*", $buf);
                     $Xw = 256*$BMP[19]+$BMP[18];
	         $Yh = 256*$BMP[23]+ $BMP[22];  
## смещение данных относительно заголовка формата BMP grayscale 8 bits:
	         $cm = 1078;  
                     $wt = int($Xw/4);  
                     $wt = $wt*4;  
                     $wt = $Xw - $wt;
                     if($wt!=0){$wt = 4 - $wt;} 
                     $minB = 256;   ##макс. число уровней квантования АЦП
                     $maxB = 0;
                     for($y=0;$y<$Yh;$y++){
		               $y1=$y*($Xw + $wt);  
                                       $left=$Yh-$y-1;
		               for($x=0;$x<$Xw;$x++){
		                                 $B[$x][$left]=$BMP[$x + $y1 + $cm];   
                   if($B[$x][$left] > $maxB){ $maxB = $B[$x][$left];} 
                   if($B[$x][$left] < $minB){ $minB = $B[$x][$left];} 
                                                                                 }
}

Подготовка данных завершена - исходное изображение занесено в массив $B[$x][$y]. Одновременно были определены максимальное $maxB и минимальное $minB значения яркости, которые необходимы для вычисления диапазона уровней квантования тестового изображения. Для решения поставленной задачи, а именно моделирования цепи Маркова процесса формирования случайного изображения, реализован следующий код:

                     use PDL;
                     use PDL::Matrix;
                     $MM = $maxB - $minB + 1;          ## порядок матрицы цепи Маркова
                     for($i=0;$i<$MM;$i++){
                     for($j=0;$j<$MM;$j++){
	                                     $w[$i][$j] = 0;
                                          }
                     }
                     for($y=0;$y<$Yh;$y++){ 
                     for($x=0;$x<$Xw-1;$x++){
                                       $I = $B[$x][$y];
                                       $j = $B[$x+1][$y];       
                                       $w[$i][$j]  +=1;
                                            }
                     }
                     $zr = zeroes($MM,$MM);
                     $zr->diagonal(0,1) ++;
                     $dg = pdl[@w];
                     $dB = $dg + $zr;
                     $row = $dB->sumover;
                     $col = vpdl $row;
                                                   print "col=$col\n";
                     $p  = $dB/$col;
                     $check_mr = $p->sumover;
                     $check_col = vpdl $check_mr;
                                                   print "\$check_col=$check_col\n";

Матрица переходов N(i,j) сохраняется в промежуточном массиве $w. Диапазон уровней квантования $MM определяет порядок матрицы переходных вероятностей. Результирующая матрица P(i,j) сохранена в массиве $p, фрагмент которого показан ниже:

                            |                 ...                       |
                            |... .00 .08 .04 .00 .00 .00 .00 .00 .00 ...|
                            |... .02 .04 .02 .00 .00 .00 .01 .01 .00 ...| 
                            |... .00 .05 .00 .00 .00 .00 .00 .00 .00 ...| 
                            |... .00 .09 .00 .00 .00 .00 .00 .00 .00 ...| 
                            |... .06 .09 .01 .00 .00 .00 .00 .01 .00 ...| 
                            |... .08 .14 .01 .01 .00 .01 .01 .00 .00 ...| 
                            |... .10 .14 .00 .05 .00 .00 .05 .05 .00 ...| 
                            |... .18 .03 .11 .00 .00 .00 .03 .00 .00 ...| 
                            |... .21 .30 .04 .00 .04 .02 .00 .00 .00 ...| 
             P(i,j)  =      |... .05 .37 .11 .00 .08 .05 .01 .01 .01 ...| 
                            |... .05 .32 .13 .00 .03 .05 .00 .00 .03 ...| 
                            |... .00 .00 .00 .25 .25 .00 .00 .00 .00 ...| 
                            |... .06 .26 .03 .00 .24 .26 .00 .00 .00 ...| 
                            |... .00 .19 .04 .00 .33 .22 .00 .04 .04 ...| 
                            |... .00 .17 .17 .00 .00 .00 .17 .17 .00 ...| 
                            |... .00 .00 .00 .12 .12 .00 .12 .12 .00 ...| 
                            |... .00 .00 .20 .00 .00 .20 .00 .00 .40 ...| 
                            |... .14 .00 .14 .00 .14 .00 .00 .00 .00 ...| 
                            |... .00 .25 .00 .00 .00 .50 .00 .00 .00 ...| 
                            |... .00 .00 .00 .00 .00 .00 .01 .00 .00 ...| 
                            |                 ...                       |.

Суммирование элементов строк матрицы $p осуществляется с помощью подпрограммы PDL sumover(), а транспонирование векторов с помощью подпрограммы vpdl(). Для осуществления контроля, значений сумм элементов матрицы $w по строкам, предусмотрен промежуточный вывод на печать вектора $col (матрица nx1). Из-за конечного размера тестового изображения, совокупность переходов i->j может не содержать все возможные переходы. Во избежание нулевых строк матрицы предусмотрено добавление диагональной единичной матрицы $zr. При наличии последовательности изображений, например телевизионных кадров, эта проблема нивелируется. Вывод на печать вектора $check_col (матрица nx1) обеспечивает контроль условия (1).

2. Анализ практических приложений использования Марковской модели

Благодаря наглядности и относительной простоте математического аппарата с одной стороны, и высокой достоверности получаемых решений с другой стороны, марковские процессы привлекли большое внимание исследователей в области экономики и бизнеса. В силу инерционности случайных процессов в экономике и их существенной зависимости от деятельности участников рынка и государственных структур, отличительной чертой таких процессов является возможность управлять (воздействовать) значениями условных вероятностных переходов. Известными примерами могут служить торговые операции, в которых сбыт продукции зависит от проведения (или не провидения) рекламных кампаний, изменения номенклатуры и рынка сбыта товаров, ценовой политики. Биржевые процессы чувствительны к высказываниям политических деятелей и позиции государственных органов и тп. В марковских моделях такие процессы относятся к управляемым марковским процессам. Рассмотрим несколько не сложных, но интересных примеров.

2.1. Анализ рынка авиалиний бизнес-класса

Ниже представлены данные, которые были получены Британскими Авиалиниями (БА) после анализа переходов клиентов бизнес-класса между авиалиниями БА и их конкурентами (Competition):

                               Следующий рейс:
                               БA      Competition     
Прошедший рейс:  БA          | 0.85          0.15 |
                 Competition | 0.10          0.90 |.

Из матрицы условных вероятностей следует, что если последний рейс клиентов бизнес-класса был связан с БА, вероятность, что их следующий рейс в БА будет 0.85. Клиенты бизнес-класса делают в среднем 2 рейса в год. В настоящее время БА имеет 30 % рынка авиалиний бизнес-класса.

Задача прогноза: какова будет доля БА на рынке авиалиний бизнес-класса после двух лет? Исходными данными являются матрица переходных вероятностей:

                       P(i,j) =   |  0.85    0.15  |
                                  |  0.10    0.90  |

и вектор (матрица 1xn) начального состояния системы s1 = [0.30, 0.70].

Решение. Сначала определим матрицу переходных вероятностей в конце года:

                       P2 = P x P = |  0.7375    0.2625  |
                                    |  0.1750    0.8250  |.

Возведение в квадрат матрицы P возникает, поскольку клиенты бизнес-класса делают 2 рейса в год в среднем. Следовательно, после того, как один год прошел, состояние системы будет иметь вид:

                       S2 = S1 x P = [0.34375, 0.65625].

После того, как прошли два года, состояние системы S3 = S2 x P = [0.368, 0.632]. Таким образом, после двух лет, доля БА на рынке авиалиний бизнес-класса составит 36.8 %. Отметим, что для векторов (матриц nx1) s2 и s3 справедливо условие (1).

Программа на Perl сводиться к перемножению матрицы P и векторов (матриц 1xn) состояний S1 и S2:

                     use PDL;
                     use PDL::Matrix; 
                     $p  =  pdl[[0.85,0.15],[0.1,0.9]];
                     $s1 =  pdl[[0.30, 0.70]]; 
                     $p = $p x $p; 
                     $s2 = $p x $s1;
                     $s3 = $p x $s2;

2.2. Прогнозирование долей рынка

Компания по продаже компакт дисков использовала теорию Маркова для анализа рекламной кампании при смене брендов компакт дисков между тремя различными брендами. Собранные данные чередования трех типов (I, II, III) брендов использовались для оценки переходных вероятностей за каждый месяц:

                         переход к: 
                         I      II     III
     переход от:   I  | 0.80   0.10   0.10 | 
                  II  | 0.03   0.95   0.02 |
                 III  | 0.20   0.05   0.75 |.

Текущие доли на рынке (в первый месяц), которые занимают эти бренды, составляет 45%, 25% и 30% для брендов I, II и III соответственно. Поставлены следующие вопросы. Каковы будут ожидаемые доли на рынке по истечении двух месяцев, то есть в 3-ем месяце? Каково отдаленное предсказание ожидаемой доли на рынке для каждой из трех марок?

Решение. Исходное состояние системы задано вектором S1 = [0.45, 0.25, 0.30] и матрицей переходных вероятностей:

                         |0.80   0.10   0.10|
                     P = |0.03   0.95    0.02|                       (2)
                         |0.20   0.05   0.75|

Следовательно, после истечения одного месяца, состояние системы S2 определяется произведением S2 = S1 x P = [0.4275, 0.2975, 0.2750], и после того, как пройдет уже два месяца, состояние системы будет равно S3 = S2 x P = [0.4059, 0.3391, 0.2550]. Отметим, что для векторов (матриц nx1) S2 и S3 справедливо условие (1). Таким образом, доли на рынке после двух месяцев составляют 40.59%, 33.91% и 25.50% для брендов I, II и III соответственно. Предполагая, что в отдаленном будущем система достигает равновесия [x1, x2, x3] = [x1, x2, x3] x P, где матрица [x1, x2, x3] = 1, а ее элементы определяются как:

                     x1 = 0.80 x1 + 0.03 x2 + 0.20 x3,
                     x2 = 0.10 x1 + 0.95 x2 + 0.05 x3,      
                     x3 = 0.10 x1 + 0.02 x2 + 0.75 x3,

где коэффициенты в системе уравнений заданы столбцами матрицы переходных вероятностей (2). Преобразование относительно x1, x2 и x3 дает:

                      0.20 x1 = 0.03 x2 + 0.20 x3                     (3) 
                      0.05 x2 = 0.10 x1 + 0.05 x3                     (4) 
                      0.25 x3 = 0.10 x1 + 0.02 x2.                    (5)

Вычитание уравнения (5) из уравнения (4), приводит к уравнению относительно x2:

                      0.05 x2 - 0.25 x3  =  0.05 x3 - 0.02 x2,
                      x2 = (0.30/0.07) x3.                            (6)

Подстановка x1 = 1 – x2 – x3 в уравнение (3) дает

                      0.20(1 - x2 – x3) = 0.03 x2 + 0.20 x3,
                      0.20 = 0.23 x2 + 0.40 x3.

Подстановка в последнее уравнение, уравнения (6) приводит к определению x3:

                       0.20 = 0.23(0.30/0.07) x3 + 0.40 x3,
                       x3 = 0.1443.

Подстановка x3 в уравнение (6) дает решение относительно x2:

                       x2 = (0.30/0.07) x3 = 0.6184,

а соответствующие подстановки x2 и x3 в x1 = 1 - x2 – x3 дает решение относительно x1 = 0.2373. Полученные значения x1, x2 и x3, в пределах ошибок округления, удовлетворяют уравнениям (3)-(5). Следовательно, отдаленные прогнозы долей на рынке составляют 23.73%, 61.84% и 14.43% для брендов I, II и III соответственно. Из сравнения оценок состояния системы S3 после двух месяцев и отдаленного распределения долей рынка видно их существенное отличие. Таким образом, можно сделать вывод о необходимости изменить обстоятельства, которые будут способствовать получению более работоспособной начальной матрицы переходных вероятностей (2).

Программная реализация решений линейных уравнений с помощью модуля PDL выглядит лаконичнее. Отметим сначала, что уравнения (3)-(5) являются однородными и, с учетом подстановки x1 = 1 – x2 – x3, сводятся к двум независимым уравнениям:

                      0.20  =  0.23 x2  +  0.40 x3,
                     -0.1    = -.15 x2    -  0.05 x3.

Запишем эту систему в виде известного матричного уравнения AX = B, решением которого будет X = A-1B, где A-1 - обратная матрица, A - матрица коэффициентов, B - свободные члены уравнения.

                     use PDL;
                     use PDL::Matrix;
                     use PDL::MatrixOps;
                     $a = pdl [[.23,.4],[-.15,-.05]];
                     $inv_a = $a->inv;
                     $b = pdl [[.2],[-.1]];
                     $result = $inv_a x $vec;
                     print $result;

На печать будут выведен столбец значений x2 и x3:

                     [
                       [ 0.6185567]
                       [ 0.1443299]
                     ],

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

3. Выводы

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

Литература

1 Цифровое кодирование телевизионных изображений - Под редакцией И.И. Цуккермана.


Ресурсы для скачивания


Комментарии

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

static.content.url=http://www.ibm.com/developerworks/js/artrating/
SITE_ID=40
Zone=Linux, Open source
ArticleID=521521
ArticleTitle=Моделирование стохастических процессов на языке Perl: Часть 2. Моделирование Марковских процессов на Perl
publish-date=09142010