一篇文档翻译，内有m文件，现在贡献出来，希望对你们有帮助
经典最大熵（ME）的问题包括根据已知函数的有限期望集确定概率分布函数（pdf）。该解决方案取决于$ N + 1 $拉格朗日乘数，该乘数是通过求解由$ N $数据约束和归一化约束形成的非线性方程组确定的。在这段简短的论文中，我们给出了三个Matlab程序来计算这些拉格朗日乘数。第一个考虑了功能可以是任何功能的一般情况。第二个考虑幂函数$ x ^ n $的特殊情况。在这种情况下，数据是$ p（x）$的几何矩。第三部分考虑傅立叶级数函数$ \ exp（-jn \ omega x）$的特殊情况。在这种情况下，数据为$ p（x）$的三角矩。
Shannon (1948)指出最大熵(ME)分布是如何通过变分法的直接应用而得到的。他定义概率密度函数p(x)的熵为
在文献中，已知求H的最大值是一种方法，它可以得到最小信息先验分布的形式;Jaynes(1968)和Zellner(1977)。Jaynes(1982)广泛地分析了离散情况下的例子，而在Lisman and Van Znylen(1972)、Rao(1973)和Gokhale(1975)中，考虑了Kagan、Linjik连续情况。在最后一种情况下，问题的一般形式如下
μ0 = 1,φ0 (x) = 1,φn (x), n = 0,……N, N是已知函数,μn N = 0,……，N是给定的期望数据。给出了该问题的经典解
(N + 1) Lagrangien参数λ=[λ0,……,λn] 通过 解决 following 组 (N + 1 ) 非线性方程
3式定义的分布形式大量已知的分布,通过选择适当的N和φn areobtained (x), N = 0,…, N。一般φn (x) areeither x或者x的对数的权力。看到Mukhrejee和赫斯特(1984),Zellner (1988), Mohammad-Djafari(1990)对许多其他的例子和讨论。许多作者广泛地分析和使用了一些特殊的案例。当φn (x) = xn, n = 0,……Nμn, N = 0,……，N是给定的N个分布矩。参见Zellner(1988)在N = 4的情况下的数值实现。
在这个部分我们提出三个用MATLAB编写的程序来解决方程组(4)。首先是一个通用的程序,φn (x)可以是任何功能。第二个是一个特殊情况φn (x) = xn, n = 0,…, n。在这种情况下,μn几何p (x)的时刻。第三个是特殊情况φn (x) = exp (?jnωx), n = 0,…, n。在这种情况下,μn的三角时刻(Fouriercomponents) p (x)。我们也给出了一些例子来说明这些程序的用处。 函数原理我们已经看到标准的解决我的问题是由(3)的拉格朗日乘数法λ是通过求解非线性方程(4)。一般情况下，这些方程都是用标准牛顿法求解的，该方法是将泰勒级数中的Gn（λ）展开到lambda的试值附近，去掉二次项和高阶项，迭代求解得到的线性系统。我们在这里给出了我们实现的数值方法的细节。当在试验λ0附近展开一阶泰勒级数的方程（4）中的Gn（λ）时，得到的线性方程由
注意向量δ和v 矩阵G 则式(5)变为
这个系统是求解δ的，我们从δ中驱动λ=λ0+δ，这将成为我们的新初始向量λ0，迭代将继续，直到δ变得适当小。注意矩阵G是对称的，我们有
因此，在每次迭代中，我们都要计算式(8)中的N(N?1)/2的积分，一般最大熵问题的算法如下: - 计算式(4)中的(N + 1)积分和矩阵G的N(N?1)/2个不同元素gnk，通过计算式(8)
- 解方程(7)找到δ
 - 计算λ=λ0 +δ,回到第三步,直到δ变得很小
方程(4)和(8)的积分可以用单变量森普森方法进行计算。 我们使用了这个方法的一个非常简单的版本
几何矩的情况现在考虑时间问题的特殊情况,φn (x) = xn, n = 0,…, n。
在这种情况下，方程(3)(4)(8) beco
这意味着 [(N +1)×(N +1)] 在方程矩阵 G (7) 成为 symmetric Hankel 矩阵完全由 2 N Gn(λ),n = 0 ,...,2N. + 1 的值因此，本例中的算法与前面的算法相同，只是做了两个简化。
- 在步骤2中我们不需要编写一个独立函数计算functionsφn (x) = xn, n = 0,…
- 在步骤4积分估值的减少,是因为元素gnkof矩阵G相关积分Gn(λ)方程(10)。这个矩阵完全由2N + 1个分量定义.
三角函数的例子另一个有趣的特殊情况是数据是p(x)的傅里叶部分
μn可能是复数的和存在值μ（下标?n）=μn。这意味着他们有以下关系 因此矩阵G的所有元素都与p(x)的离散傅里叶变换有关。注意，G是厄米特的托普利兹矩阵
实例和数值实验为了说明所提出的程序的有用性，我们首先考虑伽玛分布的情况
当约束条件为下列式子时，该分布可视为ME分布
因此很容易得到，方程(12)可以写成 现在考虑以下问题 鉴于μ1和μ2确定λ0,λ1和λ
这可以通过标准ME方法来完成。要做到这一点,首先我们必须定义rangeof x (xmin xmax, dx),和写一个函数fi_x计算functionsφ0 (x) = 1,φ1 (x) = x和φ2 (x) = lnx(参见附件中的函数fin1_x)。然后我们必须定义一个初始估计λ0λ,最后,让程序工作。 分析(α,β)和m = E {x}和均值方差σ2 = E {(x?m) 2}之间的关系你会发现他们，伽马分布的情况很有意思 或者反过来 这样我们就可以利用这些关系来确定m和σ。也要注意类似的熵的最终结果是函数的副产品表(1)给出了ME_DENS1程序得到的一些统计结果(见附件)
下一个例子是四次分布 当约束条件为此时，该分布可以认为是一个ME分布
现在考虑以下问题:给定的n,n = 1,…4计算出n,n = 0,…4。这可以由ME_DENS2程序完成。表(2)通过这个程序给出了一些数值结果 这些例子展示了如何使用建议的程序。第三个例子是在附件中,它展示了如何使用ME_DENS3程序来考虑三角矩的情况
总结本文中,我们先解决我的类分布当可用数据afinite组期望μn = E{φn (x)}一些已知函数φn (x), n = 0,…, n。我们提出三个Newton-Raphsonmethod Matlab程序来解决这一问题的一般情况下,对于几何数据,时刻φn (x) = xn和三角的时刻φn (x) = exp (?jnω0x)。最后，给出了一些具体算例的数值结果
附页
1 function [lambda,p,entr]=me_dens1(mu,x,lambda0)
2 %ME_DENS1
3 % [LAMBDA,P,ENTR]=ME_DENS1(MU,X,LAMBDA0)
4 % This program calculates the Lagrange Multipliers of the ME
5 % probability density functions p(x) from the knowledge of the
6 % N contstraints in the form:
7 % E{fin(x)}=MU(n) n=0:N with fi0(x)=1, MU(0)=1.
8 %
9 % MU is a table containing the constraints MU(n),n=1:N.
10 % X is a table defining the range of the variation of x.
11 % LAMBDA0 is a table containing the first estimate of the LAMBDAs.
12 % (This argument is optional.)
13 % LAMBDA is a table containing the resulting Lagrange parameters.
14 % P is a table containing the resulting pdf p(x).
15 % ENTR is a table containing the entropy values at each
16 % iteration.
17 %
18 % Author: A. Mohammad-Djafari
19 % Date : 10-01-1991
20 %
21 mu=mu(:); mu=[1;mu]; % add mu(0)=1
22 x=x(:); lx=length(x); % x axis
23 xmin=x(1); xmax=x(lx); dx=x(2)-x(1);
24 %
25 if(nargin == 2) % initialize LAMBDA
26 lambda=zeros(size(mu)); % This produces a uniform
27 lambda(1)=log(xmax-xmin); % distribution.
28 else
29 lambda=lambda0(:);
30 end
31 N=length(lambda);
32 %
33 fin=fin1_x(x); % fin1_x(x) is an external
34 % % function which provides fin(x).
35 iter=0;
36 while 1 % start iterations
37 iter=iter+1;
38 disp(’---------------’); disp([’iter=’,num2str(iter)]);
39 %
40 p=exp(-(fin*lambda)); % Calculate p(x)
41 plot(x,p); % plot it
42 %
43 G=zeros(N,1); % Calculate Gn
44 for n=1:N
45 G(n)=dx*sum(fin(:,n).*p);
46 end
47 %
48 entr(iter)=lambda’*G(1:N); % Calculate the entropy value
49 disp([’Entropy=’,num2str(entr(iter))])
50 %
51 gnk=zeros(N,N); % Calculate gnk
52 gnk(1,:)=-G’; gnk(:,1)=-G; % first line and first column
53 for i=2:N % lower triangle part of the
54 for j=2:i % matrix G
55 gnk(i,j)=-dx*sum(fin(:,j).*fin(:,i).*p);
56 end
57 end
58 for i=2:N % uper triangle part of the
59 for j=i+1:N % matrix G
60 gnk(i,j)=gnk(j,i);
61 end
62 end
63 %
64 v=mu-G; % Calculate v
65 delta=gnk\v; % Calculate delta
66 lambda=lambda+delta; % Calculate lambda
67 eps=1e-6; % Stopping rules
68 if(abs(delta./lambda)<eps), break, end
69 if(iter>2)
70 if(abs((entr(iter)-entr(iter-1))/entr(iter))<eps),break, end
71 end
72 end
73 %
74 p=exp(-(fin*lambda)); % Calculate the final p(x)
75 plot(x,p); % plot it
76 entr=entr(:);
77 disp(’----- END -------’)
1 %----------------------------------
2 %ME1
3 % This script shows how to use the function ME_DENS1
4 % in the case of the Gamma distribution. (see Example 1.)
5 xmin=0.0001; xmax=1; dx=0.01; % define the x axis
6 x=[xmin:dx:xmax]’;
7 mu=[0.3,-1.5]’; % define the mu values
8 [lambda,p,entr]=me_dens1(mu,x);
9 alpha=-lambda(3); beta=lambda(2);
10 m=(1+alpha)/beta; sigma=m/beta;
11 disp([mu’ alpha beta m sigma entr(length(entr))])
12 %----------------------------------
1 function fin=fin1_x(x);
2 % This is the external function which calculates
3 % the fin(x) in the special case of the Gamma distribution.
4 % This is to be used with ME_dens1.
5 M=3;
6 fin=zeros(length(x),M);
7 fin(:,1)=ones(size(x));
8 fin(:,2)=x;
9 fin(:,3)=log(x);
10 return
1 function [lambda,p,entr]=me_dens2(mu,x,lambda0)
2 %ME_DENS2
3 % [LAMBDA,P,ENTR]=ME_DENS2(MU,X,LAMBDA0)
4 % This program calculates the Lagrange Multipliers of the ME
5 % probability density functions p(x) from the knowledge of the
6 % N moment contstraints in the form:
7 % E{x^n}=mu(n) n=0:N with mu(0)=1.
8 %
9 % MU is a table containing the constraints MU(n),n=1:N.
10 % X is a table defining the range of the variation of x.
11 % LAMBDA0 is a table containing the first estimate of the LAMBDAs.
12 % (This argument is optional.)
13 % LAMBDA is a table containing the resulting Lagrange parameters.
14 % P is a table containing the resulting pdf p(x).
15 % ENTR is a table containing the entropy values at each
16 % iteration.
17 %
18 % Author: A. Mohammad-Djafari
19 % Date : 10-01-1991
20 %
21 mu=mu(:); mu=[1;mu]; % add mu(0)=1
22 x=x(:); lx=length(x); % x axis
23 xmin=x(1); xmax=x(lx); dx=x(2)-x(1);
24 %
25 if(nargin == 2) % initialize LAMBDA
26 lambda=zeros(size(mu)); % This produces a uniform
27 lambda(1)=log(xmax-xmin); % distribution.
28 else
29 lambda=lambda0(:);
30 end
31 N=length(lambda);
32 %
33 M=2*N-1; % Calcul de fin(x)=x.^n
34 fin=zeros(length(x),M); %
35 fin(:,1)=ones(size(x)); % fi0(x)=1
36 for n=2:M
37 fin(:,n)=x.*fin(:,n-1);
38 end
39 %
40 iter=0;
41 while 1 % start iterations
42 iter=iter+1;
43 disp(’---------------’); disp([’iter=’,num2str(iter)]);
44 %
45 p=exp(-(fin(:,1:N)*lambda)); % Calculate p(x)
46 plot(x,p); % plot it
47 %
48 G=zeros(M,1); % Calculate Gn
49 for n=1:M
50 G(n)=dx*sum(fin(:,n).*p);
51 end
52 %
53 entr(iter)=lambda’*G(1:N); % Calculate the entropy value
54 disp([’Entropy=’,num2str(entr(iter))])
55 %
56 gnk=zeros(N,N); % Calculate gnk
57 for i=1:N % Matrix G is a Hankel matrix
58 gnk(:,i)=-G(i:N+i-1);
59 end
60 %
61 v=mu-G(1:N); % Calculate v
62 delta=gnk\v; % Calculate delta
63 lambda=lambda+delta; % Calculate lambda
64 eps=1e-6; % Stopping rules
65 if(abs(delta./lambda)<eps), break, end
66 if(iter>2)
67 if(abs((entr(iter)-entr(iter-1))/entr(iter))<eps),break, end
68 end
69 end
70 %
71 p=exp(-(fin(:,1:N)*lambda)); % Calculate the final p(x)
72 plot(x,p); % plot it
73 entr=entr(:);
74 disp(’----- END -------’)
75 end
1 %ME2
2 % This script shows how to use the function ME_DENS2
3 % in the case of the quartic distribution. (see Example 2.)
4 xmin=-1; xmax=1; dx=0.01; % define the x axis
5 x=[xmin:dx:xmax]’;
6 mu=[0.1,.3,0.1,.15]’; % define the mu values
7 [lambda,p,entr]=me_dens2(mu,x);
8 disp([mu;lambda;entr(length(entr))]’)
1 function [lambda,p,entr]=me_dens3(mu,x,lambda0)
2 %ME_DENS3
3 % [LAMBDA,P,ENTR]=ME_DENS3(MU,X,LAMBDA0)
4 % This program calculates the Lagrange Multipliers of the ME
5 % probability density functions p(x) from the knowledge of the
6 % Fourier moments values :
7 % E{exp[-j n w0 x]}=mu(n) n=0:N with mu(0)=1.
8 %
9 % MU is a table containing the constraints MU(n),n=1:N.
10 % X is a table defining the range of the variation of x.
11 % LAMBDA0 is a table containing the first estimate of the LAMBDAs.
12 % (This argument is optional.)
13 % LAMBDA is a table containing the resulting Lagrange parameters.
14 % P is a table containing the resulting pdf p(x).
15 % ENTR is a table containing the entropy values at each
16 % iteration.
17 %
18 % Author: A. Mohammad-Djafari
19 % Date : 10-01-1991
20 %
21 mu=mu(:);mu=[1;mu]; % add mu(0)=1
22 x=x(:); lx=length(x); % x axis
23 xmin=x(1); xmax=x(lx); dx=x(2)-x(1);
24 if(nargin == 2) % initialize LAMBDA
25 lambda=zeros(size(mu)); % This produces a uniform
26 lambda(1)=log(xmax-xmin); % distribution.
27 else
28 lambda=lambda0(:);
29 end
30 N=length(lambda);
31 %
32 M=2*N-1; % Calculate fin(x)=exp[-jnw0x]
33 fin=fin3_x(x,M); % fin3_x(x) is an external
34 % % function which provides fin(x).
35 iter=0;
36 while 1 % start iterations
37 iter=iter+1;
38 disp(’---------------’); disp([’iter=’,num2str(iter)]);
39 %
40 % Calculate p(x)
41 p=exp(-real(fin(:,1:N))*real(lambda)+imag(fin(:,1:N))*imag(lambda));
42 plot(x,p); % plot it
43 %
44 G=zeros(M,1); % Calculate Gn
45 for n=1:M
46 G(n)=dx*sum(fin(:,n).*p);
47 end
48 %plot([real(G(1:N)),real(mu),imag(G(1:N)),imag(mu)])
49 %
50 entr(iter)=lambda’*G(1:N); % Calculate the entropy
51 disp([’Entropy=’,num2str(entr(iter))])
52 %
53 gnk=zeros(N,N); % Calculate gnk
54 for n=1:N % Matrix gnk is a Hermitian
55 for k=1:n % Toeplitz matrix.
56 gnk(n,k)=-G(n-k+1); % Lower triangle part
57 end
58 end
59 for n=1:N
60 for k=n+1:N
61 gnk(n,k)=-conj(G(k-n+1)); % Upper triangle part
62 end
63 end
64 %
65 v=mu-G(1:N); % Calculate v
66 delta=gnk\v; % Calculate delta
67 lambda=lambda+delta; % Calculate lambda
68 eps=1e-3; % Stopping rules
69 if(abs(delta)./abs(lambda)<eps), break, end
70 if(iter>2)
71 if(abs((entr(iter)-entr(iter-1))/entr(iter))<eps),break, end
72 end
73 end
74 % Calculate p(x)
75 p=exp(-real(fin(:,1:N))*real(lambda)+imag(fin(:,1:N))*imag(lambda));
76 plot(x,p); % plot it
77 entr=entr(:);
78 disp(’----- END -------’)
79 end
1 %ME3
2 % This scripts shows how to use the function ME_DENS3
3 % in the case of the trigonometric moments.
4 clear;clf
5 xmin=-5; xmax=5; dx=0.5; % define the x axis
6 x=[xmin:dx:xmax]’;lx=length(x);
7 p=(1/sqrt(2*pi))*exp(-.5*(x.*x));% Gaussian distribution
8 plot(x,p);title(’p(x)’)
9 %
10 M=3;fin=fin3_x(x,M); % Calculate fin(x)
11 %
12 mu=zeros(M,1); % Calculate mun
13 for n=1:M
14 mu(n)=dx*sum(fin(:,n).*p);
15 end
16 %
17 w0=2*pi/(xmax-xmin);w=w0*[0:M-1]’; % Define the w axis
18 %
19 mu=mu(2:M); % Attention : mu(0) is added
20 % in ME_DENS3
21 [lambda,p,entr]=me_dens3(mu,x);
22 disp([mu;lambda;entr(length(entr))]’)
1 function fin=fin3_x(x,M);
2 % This is the external function which calculates
3 % the fin(x) in the special case of the Fourier moments.
4 % This is to be used with ME_DENS3.
5 %
6 x=x(:); lx=length(x); % x axis
7 xmin=x(1); xmax=x(lx); dx=x(2)-x(1);
8 %
9 fin=zeros(lx,M); %
10 fin(:,1)=ones(size(x)); % fi0(x)=1
11 w0=2*pi/(xmax-xmin);jw0x=(sqrt(-1)*w0)*x;
12 for n=2:M
13 fin(:,n)=exp(-(n-1)*jw0x);
全部资料51hei下载地址： |
