Полная версия

Главная arrow Информатика arrow ВВЕДЕНИЕ В АНАЛИЗ ДАННЫХ

  • Увеличить шрифт
  • Уменьшить шрифт


<<   СОДЕРЖАНИЕ ПОСМОТРЕТЬ ОРИГИНАЛ   >>

Две программы па МатЛабе

АЛЛ. Эволюционный процесс для отыскания центра Минковского

%cm.m, computing Minkowski p-distance central point c of a series x %along with the average distance d and its proportion pe in the sum

function [c,d,pe]=cm(x,p)

n=length(x);

lb=min(x);

rb=max(x); %.................lb, rb are boundaries of the area (i)—

de=0; for ik=l:n

end

de=de/n;%...........................average p-th power of the data

%.............population setting (ii) and setting the limit, iter, to iterations

de=de+(abs(x(ik)))Ap;

pp=15; %popuIation size

feas=(rb-lb)*rand(pp,l)+lb; %generated population of p c values within the range

flag=l;

count=0;

iter=5000;

%..........evaluation of the initially generated population (iii) —

funp=0; for ii=l:pp

vv(ii)=mink(p,x,feas(ii));

end

[funi, ini]=min(vv); soli=feas(ini) %initial best c value funi %initial error si=l;%0.5; %step of change

%.............evolution of the population (iv).................

while flag==l

count=count+l;

feas=fcas+si*randn(pp,l); % Gaussian mutation added with step si for ii=l:pp

feas(ii)=max(lb, feas(ii));

feas(ii,:)=min(rb,feas(ii));% keeping the population within the range vec(ii)=mink(p,x,feas(ii)); %evaluation end

%..............elite maintenance (v)................

[fun, in]=min(vec); %best distance value sol=feas(in,:);%corresponding c value [wf,wi]=max(vec); wun=fcas(wi); %worst c if wf>funi

feas(wi)=soli;

vec(wi)=funi; % changing the worst for the elite end

if fun < funi

soli=sol;

funi=lun;

end

if (count>=iter) flag=0; end

pe=funi/de;

%............screen the results of every 1000th iteration

if rem(count,1000)=0 %funp=funi; disp([soli pe]); end end c=soli; d=funi; pe=d/de;

return

%........computing the quality of cc, the average deviation in p-th power

function dis=mink(p,x,ce)

nn=length(x);

dis=0;

for ik=l:nn

dis=dis+(abs(x(ik)-ce))Ap;

end

dis=dis/nn;

return

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

А.4.2. Оценка степенного закона: эволюционный метод и линеаризация

% plan.m, анализ степенного закона в предположении, что предиктор х и цель у % уже находятся в рабочей памяти МатЛаба % степенной закон: у=ахАЬ (1)

% линеаризованная версия: log(y)=Iog(a)+b*log(x) (2)

function plan(x,y)

%.....linear analysis of log(x) and log(v)/ линейный анализ логарифмов.....

for ii=l:length(x);xc(ii)=max(.05,x(ii));yc(ii)=max(0.05,y(ii));end;

%0.05 instead of 0 to make logarithms possible

xll=log(xc);

vll=log(vc);

[all,bll,cll, rvll]=lr(xll,yll);

%all — the slope, bll — the intercept of the linear regression

%cll the correlation coefficient, rvll — the residual variance of the linear regression yle=all*xll+bll;% yll, оцененная по линейной регрессии cd=cllA2;%determination coefficient, it should be cd=l-rvll cd rvll

%figure(l);plot(xll,yll/k.,,xll,yle,,rp’);

%.....линеаризация: оценка уравнения (1) через оценку уравнения (2)

[al,bl, rl]=llr(x,y);

% al the estimate of а, bl the estimate of b and % rl the proportion of the residual variance in the variance of у

% ylr — the linearized rule estimate for the power law for ii=l:length(x);ylr(ii)=al*x(ii)Abl;end;

%.......as is: fitting equation (1) by straightforwardly minimizing the

%.......residual variance with an evolutionary algorithm/непосредственная оценка (1) (эволюционно)

[an,bn,f, rn]=nlr(x,y);

% an the estimate of a, bn the estimate of b and % rn the proportion of the residual variance in the variance of у for ii=l:length(x);yn(ii)=an*x(ii)Abn;end; %estimatcd power law

%...........output: two-plot figure, real on the left, log on the right/двухоконная фигура,

%........слева в реальных шкапах, справа — в логарифмах

%figure(2); subplot( 1,2,1);

plot(x,y/k.x,ylr,’b.’,x,yn,’r.’);%data scatter with two estimated power laws,

% blue-linearized, red- as is / две оценки: синий — линеаризованная, красный — как есть subplot( 1,2,2);plot(xll,yll,'k.,,xll,yle/rp,>;

%...........output: text file of the results/ выдача в виде текстового файла (см. ниже)

saveplan(‘rep ell, al, bl, rl, an, bn, rn,cd);

return

% llr.m, fitting a nonlinear regression function y=axAb % using linearization/ оценка через линеаризованную версию % х is predictor, у is target, a,b -regression parameters to be fitted

% regression is power law у=а*хлЬ as reflected in the procedure % residvar is the average square error’s proportion to the variance of y; % xt, yt are predictor and target

%.....an elementary check of length compatibility........

lMength(xt); if ll~=length(yt)

disp(‘Something wrong is with the data’); pause; end

%.........calculating a and b using the linearization

for ii=l:ll;xc(ii)=max(.05,xt(ii));yc(ii)=max(0.05,yt(ii));end; %putting 0.05 instead of zero to make possible logarithms of the data xl=log(xc); %taking log of x and у yl=log(yc);

[al,bl,dl]=lr(xl,yl);

b=al;

a=exp(bl); ab=[a b];

residvar=delta(ab,xt,yt)/var(yt,l);

return

%........computing the quality of the approximation y=a*(xAb)

%which is the residual variance

function esq=delta(tt,x,y)%tt=[a, b]; x predictor, у target

a=tt(l);

b=tt(2);

esq=0;

for ii=l:length(x)

vp(ii)=a*(x(ii)Ab); %this power law function can be changed/можно взять и другую функцию esq=esq+(y(ii)-yp(ii))A2; end

esq=esq/length(x);

return;

% nlr.m, evolutionary fitting of a nonlinear regression function y=f(x,a,b) / эволюционный метод % x is predictor, у is target, a,b -regression prameters to be fitted

function [a,b, funi,residvar]=nlr(xt,yt);

% in this version the regression equation is power law y=a*xAb which is % reflected only in the subroutine ‘delta’ in the bottom for computing the % value of the average error squared;/ вид минимизируемой функции «зашит» в ‘delta’ внизу % funi is the average square error’s best value;

% residvar is its proportion to the variance of y;

% xt, yt are predictor and target

%.....an elementary check of length compatibility........

ll=length(xt); if ll~=length(yt)

disp(‘Something is wrong with the data’); pause; end

%— determine rectangle at which (a.b)-populations fluctuate / оценка области изменения а и b

[ab,bb]=ddr(xt,vt);

lb=[ab(l) bb(l)j;

rb=|ab(2) bb(2)|;

lb rb

% интервалы изменения оцениваемых величин disp(‘Hit ENTER if you wish to proceed. ‘); pause;

%.............organisation of iterations, iter the limit to their number/ число итераций равно iter

p=15; %population size/ размер популяции for ii=l:p;feas(ii,:)=(rb-lb).*rand(l,2)+lb;end;

% generated population of p pairs coefficients within the range

flag=l;

count=0;

iter=10000;%5000;

%..........evaluation of the initially generated population/ оценка начальной популяции

funp=0; for ii=l:p

vv(ii)=delta(feas(ii,:),xt,yt);

end

[funi, ini]=min(vv); soli=feas(ini,:) %initial coeffts funi %initial error si=l;%0.5; %step of change

%.............evolution of the population/эволюция популяции

while flag=1

count=count+l;

fcas=fcas+si*randn(p,2); %mutation added with step si/ мутация for ii=l:p

feas(ii,:)=max([lb;feas(ii,:)]);

feas(ii,:)=min([rb;feas(ii,:)));% keeping the population within the range vec(ii)=delta(feas(ii,:),xt,yt); %evaluation/ оценка end

[fun, in]=min(vec); %best approximation value / наилучший член популяции sol=feas(in,:);%corresponding parameters [ wf, wi j=max(vec);

wun=fcas(wi,:); %worst case / наихудший член популяции if wf>funi

feas(wi,:)=soli;

vec(wi)=funi;

%changing the worst for the best of the previous generation end

if fun < funi

soli=sol;

funi=fun;

end

if (count>=iter) flag=0;

end

residvar=funi/var(vt, 1);

%.....screen the results of every 500th iteration / вывод результатов на экран каждые 500 итераций

if rem(count,500)=0 %funp=funi; disp([soli residvar]); end end

a=soli(l);

b=soli(2);

return

а*(хЛЬ)/ оценка качества оценки

%........computing the quality of the approximation y==

function esq=delta(tt,x,y)%tt=[a, b]; x predictor, у target

a=tt(l);

b=tt(2);

esq=0;

for ii=l:length(x)

yp(ii)=a*(x(ii)Ab); %this is a power law function esq=esq+(y (ii )-y p( ii)) A2; end

esq=esq/length(x); %the average difference squared return;

% ddr.m, determination of the domain for power law y=a*xAb with Ь/ определение допустимой области % restricted

function [ab,bb]=ddr(x,y) n=length(x);

bm=(log(y(l))-log(y(2)))/(log(x(l))-log(x(2))); am=y( 1 )/(x( 1 )Abm); ab=[am am]; bb=|bin bm];

%.............finding extreme values for a and b using pairwise equations

bs=0;as=0; bsq=0;asq=0;

count=0;

for ii=l;(n-l);

if min(x(ii),y(ii))>.25 forjj=(ii+l):n

if min(x(jj),y(jj))>.25

if (x(ii)/x(jj)<0.75)|(x(ii)/x(jj)>1.25)

count=count+l;

bt=(log(y(ii))-log(y(jj)))/(log(x(ii))-log(x(jj))>; aij=y(ii)/(x(ii)Abt); aij=min(aij,100);%restriction %if (aij> 100)

% disp([ii jj]); aij

%end;

bs=bs+bt;

bsq=bsq+bt*bt;

as=as+aij;

asq=asq+aij*aij;

if bt>bb(2)

bb(2)=bt;

end;

if bt

bb( 1 )=bt;

end;

if aij>ab(2)

ab(2)=aij;

end;

ifaij

ab(l)=aij;

end;

end;

end;

end;

end;

end;

as=as/count

asq=asq/count;

sas=sqrt(asq-asA2)

bs=bs/count

bsq=bsq/count;

sbs=sqrt(bsq-bsA2)

ab(l)=as-4*sas;ab(2)=as+4*sas;

bb( 1 )=bs-4*sbs;bb(2)=bs+4"'sbs;

count

return

% saveplan.m, saving results of the power-law analysis in plan.m/ сохранение отчета о результатах ct =num2str(cc);

first=[‘Results of the power-law analysis у=ахл1У ]; alla=[ ‘On the level of logarithms, the correlation is ‘ num2str(cc) |; alex=[‘Explained proportion of log(y)-variance is ‘ num2str(100*cd) nt=[ ];

ltl=[‘Lincarizcd estimate parameter values are a= ‘ num2str(al) b= ‘ num2str(bl)]; lt2=[‘Explained proportion of у-variance is r=‘ num2str(100*(l-rl))'%’ ]; ntl=[‘”As is” estimate parameter values are a=' num2str(an) b= ‘ num2str(bn)|; nt2=[‘Explained proportion ofy-variance is r= ‘ num2str(100*(l-rn)) alltext=strvcat(alla. It 1 ,lt2,nt 1 ,nt2);

oul=[‘ These are visualized on the Figure produced:’!

our=|' The power-law estimates on the left, the logarithms, on the right’];

alltt=strvcat(alltext, oul, our);

alltt

Filename=[ file ‘.out’]; fid= fopen( Filename, ‘at’); if fid-=-l

fprintfffid, ‘%s ’, first); fprintfffid, ‘%s ’, ‘ ‘); fprintfffid, ‘%s alia); fprintfffid, ‘%s ’, alex); fprintfffid, ‘%s ’, ‘ ‘); fprintfffid, ‘%s ’, Itl); fprintfffid, ‘%s ’, It2); fprintfffid, ‘%s ‘'); fprintfffid, ‘%s ’, ntl); fprintlflid, ‘%s ’, nt2); fprintfffld, ‘%s ’, ‘ ‘); fprintfffid, ‘%s ’, oul); fprintfffid, ‘%s ’, our); fprintf(fid, ‘%s ’, ‘'); fprintlflid, ‘%s ’, ‘ ‘); fclose(fid);

end;

return

Две случайные выборки для экспериментов

П5.1. Short.dat — случайные выборки из трех разных распределений согласно данным табл. П.2.

Таблица П.2

Три столбца из грех распределений

8

20

1512

12

21

50

11

23

48

10

21

206

9

9

12

7

20

199

10

22

51

12

18

50

9

20

198

13

21

843

9

5

12

13

13

8

10

10

7

11

14

9

9

18

39

9

13

12

7

21

51

11

20

46

11

21

50

9

18

54

8

20

1391

10

19

49

10

19

41

13

24

35

12

23

45

10

13

11

12

9

9

10

21

49

7

10

10

8

17

52

12

8

8

11

20

48

12

17

199

8

11

9

8

11

13

9

20

978

12

17

51

9

20

6233

13

19

23

10

21

47

11

11

8

11

20

973

11

7

43

13

20

201

9

18

200

10

19

49

9

10

7

14

20

36

9

10

8

11

21

203

П5.2. Выборка 280 значений из нормального распределения

Ниже представлена выборка из Гауссова распределения N(0, 10) с нулевым математическим ожиданием и стандартным отклонением, равным 10.

Таблица П.З

Выборка 280 значений из Лг(0,10), отсортированных в порядке возрастания

-30,29

-12,48

-7,01

-2,99

1,76

5,58

10,35

-25,57

-12,29

-6,94

-2,91

1,98

5,59

10,50

-25,34

-12,27

-6,83

-2,83

1.98

5,63

10,94

-23,79

-11,89

-6,79

-2,78

2,07

5,65

10,98

-23,34

-11,61

-6,65

-2,75

2,08

5,65

11,08

-22,38

-11,50

-6,64

-2,66

2,14

5,74

11,13

-22,37

-11,33

-6,11

-2,66

2,14

5,74

11,64

-21,78

-11,10

-6,02

-2,58

2,18

5,81

12,28

-21,05

-10,78

-5,98

-2,52

2,21

5,82

12,33

-20,89

-10,57

-5,87

-2,23

2,27

5,89

12,59

-20,65

-10,52

-5,53

-2,07

2,28

6,13

12,79

-19,10

-10,44

-5,35

-2,06

2,29

6,26

12,93

-18,16

-10,13

-5,33

-1,91

2,36

6,29

13,15

-17,95

-10,09

-5,22

-1,90

2,37

6,51

13,24

-17,79

-10,08

-5,17

-1,74

2,56

6,55

13,42

-17,58

-10,06

-4,91

-1,60

2,71

6,59

13,44

-16,47

-9,79

-4,82

-1,51

2,79

6,59

13,48

-16,43

-9,11

-4,62

-1,44

2,85

6,65

13,56

-16,31

-9,08

-4,58

-1,42

2,91

7,00

13,99

-16,19

-9,01

-4,53

-1,28

2,94

7,09

14,27

-16,15

-8,95

-4,43

-1,26

2,98

7,16

14,69

-16,14

-8,93

-4,26

-0,80

3,16

7,30

14,95

-15,90

-8,71

-4,18

-0,79

3,21

7,58

15,35

-15,89

-8,53

-4,17

-0,73

3,27

7,99

15,74

-15,67

-8,49

-4,08

-0,50

3,27

8,34

15,82

-15,56

-8,01

-4,01

-0,49

3,46

8,57

15,84

-15,50

-7,98

-3,98

-0,23

3,66

8,58

15,99

-15,04

-7,97

-3,95

-0,21

3,74

8,70

16,03

-15,00

-7,75

-3,84

-0,08

3,80

8,85

16,84

-14,91

-7,67

-3,78

-0,02

4,29

8,87

16,87

-14,16

-7,48

-3,74

0,03

4,39

8,97

17,29

-14,14

-7,46

-3,65

0,33

4,41

9,02

17,62

-14,04

-7,44

-3,61

0,65

4,42

9,08

18,43

-13,88

-7,37

-3,59

0,70

4,48

9,12

19,57

-13,84

-7,37

-3,47

0,78

4,60

9,39

19,58

-13,72

-7,35

-3,46

0,80

4,78

9,57

20,80

-13,58

-7,27

-3,39

1,10

4,94

9,83

22,38

-13,33

-7,24

-3,14

1,20

5,28

10,02

22,66

-12,98

-7,20

-3,02

1,38

5,41

10,08

29,50

-12,68

-7,03

-3,01

1,58

5,54

10,09

32,03

 
<<   СОДЕРЖАНИЕ ПОСМОТРЕТЬ ОРИГИНАЛ   >>