-
Notifications
You must be signed in to change notification settings - Fork 0
/
zipfGrignettiEntropy.m
100 lines (90 loc) · 3.56 KB
/
zipfGrignettiEntropy.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
## Copyright (C) 2017 Leonardo Araujo
##
## This program is free software; you can redistribute it
## and/or modify it under the terms of the GNU General Public License
## as published by the Free Software Foundation; either version 3 of
## the License, or (at your option) any later version.
##
## This program is distributed in the hope that it will be
## useful, but WITHOUT ANY WARRANTY; without even the implied warranty
## of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with the fuzzy-logic-toolkit; see the file COPYING. If not,
## see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{Hmean}, @var{Hdelta}] = } dpcmenco (@var{s}, @var{N}, @var{q}, @var{b})
## @deftypefnx {Function File} {[@var{Hmean}, @var{Hdelta}] = } dpcmenco (@var{s}, @var{N}), @var{q})
## @deftypefnx {Function File} {[@var{Hmean}, @var{Hdelta}] = } dpcmenco (@var{s}, @var{N}), @var{[]}, @var{b})
## @deftypefnx {Function File} {[@var{Hmean}, @var{Hdelta}] = } dpcmenco (@var{s}, @var{N})
##
## Estimate zipfian(-mandelbrot) entropy using Grignetti's approach.
## @table @code
## @item [Hmean, Hdelta] = zipfGrignettiEntropy(s,N,q,b)
## Estimate the entropy of a zipfian-mandelbrot distributted source with Zipf's exponnet s,
## Mandelbrot's flattening constant q, vocabulary size N and base b.
##
## @item [Hmean, Hdelta] = zipfGrignettiEntropy(s, N, q)
## Use base 2. The results are given in bits.
##
## @item [Hmean, Hdelta] = zipfGrignettiEntropy(s, N, [], b)
## Compute the Zipf's entropy in base b.
##
## @item [Hmean, Hdelta] = zipfGrignettiEntropy(s, N)
## Zipf's entropy is given in bits.
##
## @end table
## @end deftypefn
## @seealso{entropy}
function [Hmean, Hdelta] = zipfGrignettiEntropy(s,N,q,b)
narginchk(2,4);
if !isnumeric(s) || !isnumeric(N) || (exist("q") && !isnumeric(q)) || (exist("b") && !isnumeric(b)), error('Only numeric entries are allowed!'); endif;
if s < 0 || N < 0 || (exist("q") && q < 0) || (exist("b") && b < 0), error('You must provide positive values only!'); endif;
if N != floor(N), error('The vocabulary length must be integer!'); endif;
if nargin < 3, % Zipf's law
C = 1./sum( [1:N].^(-s) );
K = max( ceil(e^(1/s)) ,1);
SK = 0; for n=1:K-1, SK+=log(n)/(n^s); endfor;
lm = integrallogxx(N,s) - integrallogxx(K,s) + SK + log(N)/(N^s);
n=K; SK+=log(n)/(n^s);
LM = integrallogxx(N-1,s) - integrallogxx(K,s) + SK + log(N)/(N^s);
Hl = (s*C/log(2)) .* lm - log(C)/log(2);
Hh = (s*C/log(2)) .* LM - log(C)/log(2);
if nargin > 3, % change base from 2 (bits) to b
Hl = Hl / log2(b);
Hh = Hh / log2(b);
endif
Hmean = (Hh + Hl)/2;
Hdelta = Hh - Hl;
else, % Zipf-Mandelbrot
C = 1./sum( ([1:N]+q).^(-s) );
K = max( ceil(e^(1/s)-q) ,1);
SK = 0; for n=1:K-1, SK+=log(n+q)/((n+q)^s); endfor;
lm = integrallogxx(N,s,q) - integrallogxx(K,s,q) + SK + log(N+q)/((N+q)^s);
LM = integrallogxx(N-1,s) - integrallogxx(K,s) + SK + log(K+q)/((K+q)^s) + log(N+q)/((N+q)^s);
Hl = (s*C/log(2)) .* lm - log(C)/log(2);
Hh = (s*C/log(2)) .* LM - log(C)/log(2);
if nargin > 3, % change base from 2 (bits) to b
Hl = Hl / log2(b);
Hh = Hh / log2(b);
endif
Hmean = (Hh + Hl)/2;
Hdelta = Hh - Hl;
endif
endfunction
function r = integrallogxx(x,s,q)
if nargin < 3,
if s != 1,
r = ( (x^(1-s))/(1-s) ) * ( log(x) - 1/(1-s) );
else
r = 0.5*log(x).^2;
endif
else,
if s != 1,
r = ( ((x+q)^(1-s))/(1-s) ) * ( log(x+q) - 1/(1-s) );
else
r = 0.5*log(x+q).^2;
endif
endif
endfunction