function [W,p,chiqua,df,m] = kendall_W(X)

% [W,p,chiqua,df,m] = kendall_W(X) computes Kendall's Coefficient of
% (rater) Concordance W and its empirical significance p (for data 
% including tied ranks).

% Kendalls W mit Rangbindungskorrektur
%
% Input:
% X(N,m) : N - (zu beurteilende) Objekte, m - Urteiler
%
% Output:
% W: Kendalls Konkordanz-Koeffizient
% m: Anzahl der Urteiler
% chiqua: Chi-Quadrat-Wert zur Signifikanzprüfung
% df: Freiheitsgrade zur Signifikanzprüfung
% p: empirische Signifikanz

% Markus Brückl, 2010, makUs@brYkl.de

[Rdata, Vs] = bau_rang_matrix(X);
% Größe der (obligatorisch lückenlosen, sonst alle Rangstatistiken falsch) Matrix erfassen
N = size(Rdata,1); %Objektzahl
m = size(Rdata,2); %Urteilerzahl

T = sum(Rdata,2)'; % Rangsummen über die Urteiler, also pro Objekt
T_quer = sum(T)/N; % = (m/2) * (N+1); durchschnittliche Rangsumme

QSR = sum((T-T_quer).^2);  %Summe der quadrierten Differenzen
SUMV = sum(Vs); % Summe der Bindungskorrektur-Koeffizienten 
QSR_max = (m^2*(N^3-N)-m*SUMV) / 12; % mit Rangbindungskorrektur
% QSR_max = m^2*(N^3-N) / 12; % ohne Rangbindungskorrektur
W = QSR/QSR_max;

% Signifikanz berechnen
chiqua = m*(N-1)*W; % = (12*QSR) / (m*N*(N+1)-(1/(N-1))*SUMV);
df = N-1;
p = 1 - chi2cdf(chiqua,df);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%bau Rang-Matrix
function [R,Vs] = bau_rang_matrix(A)

% bau_rang_matrix transformiert eine (lï¿½ckenlose) Matrix pro Spalte (m) in eine
% Rangreihe der Lï¿½nge N und verwendet dazu die vektorbasierte Funktion bau_rang
% A(N,m) : N - Objekte, m - Urteiler 

 % N = size(A,1);
 m = size(A,2);
 
 for h=1:m
    [R(:,h),Vs(h)] = bau_rang(A(:,h));
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%bau Rang-Vektor
function [rang,V] = bau_rang(X)

% - transformiert einen Eingabevektor X in eine Rangreihe mit 
%   arithmetischen Rangmittelwerten bei Rangbindungen
% - berechnet die Bindungskorrektur pro Urteiler fï¿½r Kendalls W 

[Xsort,index] = sort(X); % Objekt-Werte pro Urteiler sortieren
hilf = 0;
V = 0;
for i=1:size(X,1)
  if hilf == 0
  ind = find(X==Xsort(i));
  vi = size(ind,1);
      if vi == 1
          rang(index(i)) = i;
          hilf = 1;
      elseif vi > 1
          for j=1:vi
              rang(index(i+j-1)) = mean(i:i+vi-1);
          end
          V = V + (vi^3-vi);    % Bindungskorrektur pro Urteiler berechnen
          hilf = vi;
      else
          rang(index(i)) = NaN;
          hilf = 1;
      end
  end
  hilf = hilf-1;
end

rang = rang';
