Delphi 中的对数似然实现

Question

我正在尝试计算文本中单词对的对数似然得分，并且在我的Delphi实现中得到了与我从在线Java和Python源代码派生的同样异常的结果。Ted Dunning于1993年发表了这个源代码，为一个特定的单词对给出了以下结果：

• K11（AB，即联合频率）= 110，
• K12（没有B附近的A单词）= 2442，
• K21（没有A附近的B单词）= 111
• K22（除A或B以外的单词数）= 29114

并且给出了所需的结果为270.72。

Dunning还在R中提供了一种实现方法，网址为 http://tdunning.blogspot.co.uk/2008/03/surprise-and-coincidence.html

计算对数似然比得分（也称为G2）非常简单，LLR = 2 sum(k) (H(k) - H(rowSums(k)) - H(colSums(k)))
其中H是Shannon熵，计算为(k_ij / sum(k)) log (k_ij / sum(k))的总和。 在R中，此函数定义为 H = function(k) {N = sum(k) ; return (sum(k/N * log(k/N + (k==0)))}

但我不知道R，也不确定如何将其翻译为Pascal。

我的翻译尝试包括以下函数：

function LnOK(x : integer): extended;
begin
  if x<=0 then Result :=0
  else Result := Ln(x);
end;

function Entropy2(a, b: Integer): extended;
begin
  Result := LnOK(a + b) - LnOK(a) - LnOK(b);
end;

function Entropy4(a, b, c, d: Integer): extended;
begin
  Result := LnOK(a + b + c + d) - LnOK(a) - LnOK(b) - LnOK(c) - LnOK(d);
end;

function Log_likelihood_from_Java(f1, f2, joint, total_tokens: Integer): 
  single;
var
  k11, k12, k21, k22: Integer;
  matrixEntropy, rowEntropy, colEntropy: extended;
begin
  k11 := joint;
  k12 := f2 - joint;
  k21 := f1 - joint;
  k22 := total_tokens - f1 - f2 + joint;
  rowEntropy := Entropy2(k11 + k12, k21 + k22);
  colEntropy := Entropy2(k11 + k21, k12 + k22);
  matrixEntropy := Entropy4(k11, k12, k21, k22);
  if (rowEntropy + colEntropy < matrixEntropy) then
    Result := 0.0 // round off error
  else
   Result := 2.0 * (rowEntropy + colEntropy - matrixEntropy);
end;

以上代码返回的是7.9419，而非调用时期望的270.72。

Log_likelihood_from_Java(2552, 221, 110, 31777);

感谢您的帮助！

- Mike Scott

---

真正奇怪的是，你的代码与你链接的 R 代码没有明显的关系。 - David Heffernan

@DavidHeffernan 这似乎是从Java移植过来的 https://github.com/apache/mahout/blob/master/math/src/main/java/org/apache/mahout/math/stats/LogLikelihood.java - fantaghirocco

@fantaghirocco 看起来提问者选择了一个不太清晰的实现，因为他懂Java而不是R。但是R代码中的算法要清晰得多。 - David Heffernan

Answer 1

我发现在 LnOk 函数的翻译中存在问题，应该如下所示：

function LnOK(x: Integer): Extended;
begin
  if x = 0 then
    Result := 0
  else
    Result := x * Ln(x);
end;

---

离题

顺便说一句，为了改善编码风格，您可能更喜欢重载 Entropy 函数，而不是使用不同的名称调用它们：

function Entropy(a, b: Integer): Extended; overload;
begin
  Result := LnOK(a + b) - LnOK(a) - LnOK(b);
end;

function Entropy(a, b, c, d: Integer): Extended; overload;
begin
  Result := LnOK(a + b + c + d) - LnOK(a) - LnOK(b) - LnOK(c) - LnOK(d);
end;

- fantaghirocco

---

@MikeScott LnOK意味着如果Ln处理特定方式的域外情况，则LnOK实现受保护。但是此函数并未如此操作。是的，它会挑选出x = 0，但它忽略了x < 0。否则，它将返回x * Ln（x）而不是名称所示的Ln（x）。您应该将R源代码视为编写代码的更好基础。您的代码有些混乱。 - David Heffernan

Answer 2

我无法理解你写的代码与你提供的R代码毫无明显关联，我并没有尝试去协调这些差异。下面是R代码的逐字翻译，通过这种方式编写算法要简单得多，我相信你会同意这一点。

{$APPTYPE CONSOLE}

uses
  SysUtils, Math;

type
  TVector2 = array [1..2] of Double;
  TMatrix2 = array [1..2] of TVector2;

function rowSums(const M: TMatrix2): TVector2;
begin
  Result[1] := M[1,1] + M[1,2];
  Result[2] := M[2,1] + M[2,2];
end;

function colSums(const M: TMatrix2): TVector2;
begin
  Result[1] := M[1,1] + M[2,1];
  Result[2] := M[1,2] + M[2,2];
end;

function H(const k: array of Double): Double;
var
  i: Integer;
  N, kOverN: Double;
begin
  N := Sum(k);
  Result := 0.0;
  for i := low(k) to high(k) do begin
    kOverN := k[i]/N;
    if kOverN>0.0 then begin
      Result := Result + kOverN*Ln(kOverN);
    end;
  end;
end;

function LLR(const M: TMatrix2): Double;
var
  k: array [1..4] of Double absolute M; // this is a little sneaky I admit
  rs, cs: TVector2;
begin
  rs := rowSums(M);
  cs := colSums(M);
  Result := 2.0*Sum(k)*(H(k) - H(rs) - H(cs));
end;

var
  M: TMatrix2;

begin
  M[1,1] := 110;
  M[1,2] := 2442;
  M[2,1] := 111;
  M[2,2] := 29114;
  Writeln(LLR(M));
end.

Output

 2.70721876936232E+0002

- David Heffernan