%% Testo
% Partendo dalla tabella di contingenza che segue:
% X= ricordo della pubblicità
% Y= acquisto del prodotto
% X\Y sì no
% sì 87 188
% no 42 406
% 1) Calcolare la tabella delle frequenze teoriche in caso di indipendenza. Implementare la procedura tramite moltiplicazione matriciale oppure tramite ciclo for.
% 2) Calcolare la statistica di \chi^2 di Pearson prima senza utilizzare cicli for e poi utilizzando i cicli for. Calcolare e commentare il p-value della statistica di Pearson.
% 3) Scrivere un programma che data una qualsiasi tabella di contingenza 2x2 scriva automaticamente "Associazione positiva" in presenza di associazione positiva e "Associazione negativa" in presenza di associazione negativa e che riporti la differenza tra la frequenza effettiva e quella attesa per cella (1,1) della tabella.
% 4) Calcolare l'indice \phi utilizzando la formula
% \phi= (\pm) \sqrt{ \frac{\chi^2}{n}} oppure la formula \phi=\frac{n_{11}  n_{22} -n_{12}n_{21}}{\sqrt(n_{1.}n_{2.}n_{.1}n_{.2})}
% 5) Verificare che l'indice \phi può essere ottenuto come coefficiente di correlazione tra i due fenomeni X e Y.
% 6) Calcolare e commentare il rapporto dei prodotti incrociati \theta Calcolare l'indice \theta prima utilizzando i 4 elementi della tabella di contingenza iniziale e poi utilizzando le funzioni MATLAB prod, diag e flip
% 6) Calcolare e commentare gli indici Q= \frac{\theta-1}{\theta+1} e U=  \frac{\sqrt\theta-1}{\sqrt \theta+1}
% 7) Studiare le funzioni  Q= \frac{\theta-1}{\theta+1} e U= \frac{\sqrt\theta-1}{\sqrt  \theta+1} nell'intervallo [0 20]
%% Soluzione

N= [87 188;
    42 406];
% n=numero di unità statistiche del campione.
n=sum(sum(N));

% sumc (1x2)= frequenze marginali di colonna
% sumc(1) = somma delle frequenze della prima colonna $n_{.1}$
% sumc(2) = somma delle frequenze della seconda colonna $n_{.2}$
sumc=sum(N,1);

% sumr (2x1)= frequenze marginali di riga
% sumr(1) = somma delle frequenze della prima riga $n_{1.}$
% sumr(2) = somma delle frequenze della seconda riga $n_{2.}$
sumr=sum(N,2);

%
%% Calcolo della tabella delle frequenze teoriche tramite moltiplicazione matriciale

% Calcolo delle frequenze teoriche
Ntheo=(sumr*sumc/n);
disp('Matrice delle frequenze teoriche')
disp(Ntheo)

%
%% Calcolo delle frequenze teoriche tramite doppio ciclo for

[nr,nc]=size(N);
Ntheochk=zeros(nr,nc);
for i=1:nr
    for j=1:nc
        Ntheochk(i,j)=sumr(i)*sumc(j)/n;
    end
end

%
%% Calcolo indice Chi2 senza cicli for
% Implementazione dell'indice Chi2 in maniera matriciale

Chi2=sum(sum(((N-Ntheo).^2)./Ntheo));
%
% Implementazione dell'indice Chi2 
% in maniera vettoriale
Chi2chk=sum(((N(:)-Ntheo(:)).^2)./Ntheo(:));
isequal(Chi2,Chi2chk)

%
%% Implementazione dell'indice Chi2 tramite doppio ciclo for

Chi2chk=0;
for i=1:nr
    for j=1:nc
        Chi2chk=Chi2chk+(N(i,j)-Ntheo(i,j))^2/Ntheo(i,j);
    end
end

%
%% Calcolo p-value della statistica di Pearson
% La statistica di Pearson si distriuisce come una v.c. Chi2 con (r-1)(c-1) gradi di libertà. In questo esempio dato che r=c=2;

gdl=1;
pval=1-chi2cdf(Chi2,gdl);
disp(['Il pvalue del test è: ' num2str(pval)])

% Commento: rifiuto decisamente l'ipotesi nulla di indipendenza tra 
% X= ricordo della pubblicità
% e
% Y= acquisto del prodotto

%
%% Mostrare il tipo di associazione
% 3) Scrivere un programma che data una qualsiasi tabella di contingenza 2x2 scriva automaticamente "Associazione positiva" in presenza di associazione positiva e "Associazione negativa" in presenza di associazione negativa e che riporti la differenza tra la frequenza effettiva e quella attesa per cella (1,1) della tabella.

if N(1,1)-Ntheo(1,1)>0
    disp('Associazione positiva')
    disp(['La differenza tra n_{11} n e n_{11}^* è =' num2str(N(1,1)-Ntheo(1,1))])
else
    disp('Associazione negativa')
    disp(['La differenza tra n_{11} n e n_{11}^* =' num2str(N(1,1)-Ntheo(1,1))])
end

%
%% Calcolo dell'indice phi utilizzando il valore precalcolato dell'indice Chi2

if N(1,1)*N(2,2)-N(1,2)*N(2,1)>0
    phi=sqrt(Chi2/n);
else
    phi=-sqrt(Chi2/n);
end

%
%% Indici \phi tramite la formula
% \phi=\frac{n_{11}  n_{22} -n_{12}n_{21}}{\sqrt{n_{1.}n_{2.}n_{.1}n_{.2}}} Calcolo dell'indice \phi senza passare attraverso il valore precalcolato dell'indice Chi2

phichk=(N(1,1)*N(2,2)-N(1,2)*N(2,1))/sqrt(prod(sumr)*prod(sumc));
disp(['Il valore dell''indice \phi è ' num2str(phi)])

%
%% Indice
% 4) Verificare che l'indice \phi può essere ottenuto come coefficiente di correlazione tra i due fenomeni X e Y.

x=[1;0];
y=x';
% Calcolo Media di X 
% (fenomeno dicotomico che assume modalità 0 e 1 con frequenze sumr(1) e sumr(2))
Mx=sum(x.*sumr)/n;
% Calcolo Varianza di X  
Varx=sum((x-Mx).^2.*sumr)/n;

% Calcolo Media di Y 
My=sum(y.*sumc)/n;
% Calcolo Varianza di Y  
Vary=sum((y-My).^2.*sumc)/n;

% Calcolo covarianza tra X e Y
covxy=sum((x*y).*N,'all')/n-Mx*My;
rxy=covxy/(sqrt(Varx*Vary));

disp(['Indice phi ottenuto come coefficiente di correlazione lineare ='  num2str(rxy)])

% Osservazione: utilizziamo la funzione crosstab2datamatrix dell'FSDA
% toolbox scaricabile da http://rosa.unipr.it
% Xori=crosstab2datamatrix(N);
%  disp(corr(Xori))

%
%% Indice \theta (rapporto dei prodotti incrociati)
% 6) Calcolare e commentare il rapporto dei prodotti incrociati \theta Calcolare l'indice \theta prima utilizzando i 4 elementi della tabella di contingenza iniziale e poi utilizzando le funzioni MATLAB prod, diag e flip

th=(N(1,1)*N(2,2))/(N(1,2)*N(2,1));
disp(['Il valore dell''indice \theta è ' num2str(th)])

% Calcolo dell'indice  $\theta$  utilizzando le funzioni MATLAB prod, diag e flip
thchk=prod(diag(N))/prod(diag(flip(N)));

%
%% Relazione tra theta e Q

thx=0:0.1:20;
Q=(thx-1)./(thx+1);
plot(thx,Q)
xlabel('$\theta$','Interpreter','latex')
ylabel('$Q= \frac{\theta-1}{\theta+1}$','Interpreter','latex')

%
%% Relazione tra theta e U

thx=0:0.1:20;
U=(sqrt(thx)-1)./(sqrt(thx)+1);
plot(thx,U)
xlabel('$\theta$','Interpreter','latex')
ylabel('$U= \frac{\sqrt \theta-1}{\sqrt \theta+1}$','Interpreter','latex')

%
%% Indice Q

Q=(th-1)/(th+1);
disp(['Il valore dell''indice Q è ' num2str(Q)])

%
%% Indice U

sqrtth=sqrt(th);
U=(sqrtth-1)/(sqrtth+1);
disp(['Il valore dell''indice U è ' num2str(U)])