Файл: Цель работы Изучить способы аппроксимации функции плотности вероятности по данным наблюдений. Задачи.docx

Аппроксимация плотности вероятности

Цель работы

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

Задачи

1. Сгенерировать выборку случайных чисел.2. Сформировать локальные оценки логарифма плотности вероятности по полученной выборке значений случайной величины.3. По новой выборке локальных оценок логарифма плотности выполнить аппроксимацию функции плотности;Исходные параметры для функции аппроксимации заданы в таблице 1.

n	A	a	lampda	g
200	[0.66; 0.33]	[0; 0]	[1; 0.1]	'N'

Ход работы

М-файл для работы представлен на рис. 2.

n = 200 A = [0.66 0.33]; a = [0 0]; lambda = [1 0.1]; g = 'N'; Approxpdf(n, A, a, lambda, g)

Рис. 2. М-файл для работыИспользуя функцию Approxpdf(n,A,a,lampda,g), были получены локальные оценки логарифма плотности (см. рис. 3).Рис. 3. Формирование локальных оценок логарифма плотностиНа рис. 4 показаны результаты аппроксимации по полученным локальным оценкам.Рис. 4. Результаты аппроксимации

Вывод

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

---

functionApproxpdf(n,A,a,lampda,g)X=GenRandom(n,A,a,lampda,'N');[xv,fv,fz,m]=approx(X,1);F=Truepdf(xv,A,a,lampda,'N');subplot(211);plot(xv,log(fv),'ro-',xv,log(F),'b');xlabel('X');ylabel('logf');legend('LF','LF*');

% n	объем-выборка
% R	Равномерное	N	Нормальное
% L	Логистическое	E	Экспоненциальное
% G	Гамма (хи-квадрат)

subplot(212); plot(xv,F,'k',xv,fz,'b-'); xlabel('X'); ylabel('f'); legend('F','F*'); function X=GenRandom(n,A,a,lampda,g) %гдеA=[A1A2...Am]%a=[a1a2...am]lampda=[lampda1lampda2...lampdam]%еслиf=0.5N(5,1)+0.5N(11,2)% A=[0.5 0.5] a=[5 11] lampda=[1 2]X=[];for i=1:length(A)ni=n*A(i);switchgcase'R'Xi=a(i)+ lampda(i)*rand(ni,1);case'N';Xi=a(i)+lampda(i)*(sqrt(-2*log(rand(ni,1)))).*sin(2*pi*rand(ni,1));case'E'Xi=a(i)- lampda(i) *log(rand(ni,1));case'L'Xi=a(i)+lampda(i)*log((1-rand(ni,1))./rand(ni,1));case'G'ng=input(['n',num2str(i),'\n']);forj=1:niXi(j)=a(i)+ lampda(i)*sum((randn(ng,1)).^2)/2;endendif size(Xi,2)>size(Xi,1)Xi=Xi';endX=[X;Xi];endX=sort(X);function F=Truepdf(X,A,a,lampda,g)F=0;for i=1:length(A)switchgcase'N';F=F+A(i)*exp(-0.5*((X-a(i))/lampda(i)).^2)/(sqrt(2*pi)*lampda(i)); case'E'F=F+A(i)*exp(-(X-a(i))/lampda(i))/lampda(i);case'L'z=(X-a(i))/lampda(i);F=F+A(i)*1/lampda(i)*exp(-z)./(1+exp(-z)).^2;case'G'n=input(['n',num2str(i),'\n']);F=F+A(i)*(X-a(i)).^(n/2-1)/(lampda(i)^(n/2)*gamma(n/2)).*exp(-(X-a(i))/lampda(i)); endendfunction[xv,fv,fz,m]=approx(x,criteria)%ДубовИ.Р.Аппроксимациявероятностныхзаконоввсистемахмониторинга:%Владимир,2001.–137с.%criteria=0или1x=sort(x); %упорядоченнаявыборкаn=length(x);fori=1:(n-1)xv(i)=(x(i)+x(i+1))/2;endxv=xv';k=1:n;d=sum(k.^-2); %ДисперсияошибокM=-1*sum(k.^-1); % пси-функция Эйлера (дигамма – функция)nr=n-1;for i=1:nrfv(i)=exp(M)/(x(i+1)-x(i));end

---

fv=fv'

lfv=log(fv); % lfv локальные оценок логарифма плотности вероятности

% Аппроксимация логарифма плотности

if criteria==0 [m,fest,lfz,a]=Malloca(xv,lfv,d);

else

for i=1:n-1 [yest,fest]=fit(xv,lfv,'x',i-1,1);

f(i)=abs(mean((lfv-yest).^2)-d);

j=i+1; [yest,fest]=fit(xv,lfv,'x',j-1,1);

c=abs(mean((lfv-yest).^2)-d); if c>f(i)

break; end

end m=i-1;

[lfz,fest]=fit(xv,lfv,'x',m,1); end

fz=exp(lfz);
function [yest,fest,a]=fit(x,y,f,m,s)

% y=X1(x)+X2(x)+X3(x) … +Xm(x)+error of std s

% where Xj(x) bases function

% A=[X1(x) X2(x) X3(x) … Xm(x)]

% % size of A is N*M

% for example if y(x)=a_0+a_1*x+a_2*x^2 … a_m*x^m

% A=[1 X X^2 … X^m]

% a=[a_0 a_1 … a_m]

% where X=[x_1 x_2 … x_n]'
if size(x,2)>size(x,1) x=x';

end

if size(y,2)>size(y,1) y=y';

end

if size(s,2)>size(s,1) s=s';

end n=length(x); s=s.*ones(n,1); si=1./s;

sg=sqrt(1/mean(si.^2)); sb=si*sg; Af=zeros(n,m+1);

for i=0:m Af(:,i+1)=(((eval(f)).^i));

end sa=sb(:,ones(1,m+1)); A=Af.*sa; alpha=A'*A;

b=y.*sb; Beta=A'*b; a=inv(alpha)*Beta; yest=Af*a;

fest=[num2str(a(1))]; for i=2:length(a)

fest=[fest,'+',num2str(a(i)),'*',f,'.^',num2str(i-1)]; end

%fest=['y=',fest]; temp=fest; n=length(temp)-1; j=1;

i=1;

while (j
if (temp(j)=='+')&(temp(j+1)=='-') fest(i)=[];

i=i-1; end i=i+1; j=j+1;

end
function [degree,fest,LF,a]=Malloca(xf,lfv,sd2)

%==================================================

% malloca criteria

% c=abs(mean((lfv-yest).^2)-d*(1-2*(j-1)/n))

%==================================================

n=length(xf); sd=sd2.^0.5; for i=2:n

[yest,fest]=fit(xf,lfv,'x',i-1,sd); f(i)=sum(((yest-lfv).^2)./sd2)-(n-2*(i-1)); j=i+1;

[yest,fest,a]=fit(xf,lfv,'x',j-1,sd);

c=sum(((yest-lfv).^2)./sd2)-(n-2*(j-1)); if c>f(i)

break; end

end degree=i-1;

[yest,fest,a]=fit(xf,lfv,'x',degree,sd); LF=yest;

---