%% 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] prima senza poi con il calcolo simbolico

%% 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)])

%% Modo alternativo per verificare la correlazione
% Osservazione: utilizziamo la funzione crosstab2datamatrix dell'FSDA
% toolbox 
% Creiamo Xori la matrice dei dati che ha generato la tabella di
% contingenza N
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 (faccio il grafico utilizzando una sequenza di valori)
figure
thx=0:0.1:200;
Q=(thx-1)./(thx+1);
plot(thx,Q)
xlabel('$\theta$','Interpreter','latex')
ylabel('$Q= \frac{\theta-1}{\theta+1}$','Interpreter','latex')
title('Relazione tra l''indice $Q$ e l''indice $\theta$','Interpreter','latex')

%% Relazione tra theta e Q (utilizzando il calcolo simbolico)
% Osservazione: è necessario scaricare ed installare il symbolic toolbox

% Definisco che la mia incognita è una variabile che si chiama th
syms th
Q=(th-1)./(th+1);
fplot(Q)
% fplot(Q) fa il grafico della funzione Q nell'intervallo scelto in
% automatico
xlabel('$\theta$','Interpreter','latex')
ylabel('$Q= \frac{\theta-1}{\theta+1}$','Interpreter','latex')
title('Relazione tra $\theta$ e $Q$','Interpreter','latex')


%% Relazione tra th e Q con intervallo x personalizzato
fplot(Q,[0 20])
% fplot(Q) fa il grafico della funzione Q nell'intervallo scelto in
% automatico
xlabel('$\theta$','Interpreter','latex')
ylabel('$Q= \frac{\theta-1}{\theta+1}$','Interpreter','latex')
title('Relazione tra $\theta$ e $Q$','Interpreter','latex')


%% Relazione tra theta e U
figure
thx=0:0.1:200;
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')
title('Relazione tra l''indice $U$ e l''indice $\theta$','Interpreter','latex')

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

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

%% Relazione tra theta e U con il calcolo simbolico
% Definisco che la mia incognita è una variabile che si chiama th
syms th
U=(sqrt(th)-1)/(sqrt(th)+1);
fplot(U,[0 200])
xlabel('$\theta$','Interpreter','latex')
ylabel('$U= \frac{\sqrt \theta-1}{\sqrt \theta+1}$','Interpreter','latex')
title('Relazione tra $\theta$ e $Q$','Interpreter','latex')