% Tabella di contingenza di partenza 
%
% X= ricordo della pubblicità
% Y= acquisto del prodotto 
% 
% X\Y sì no
% 
% sì 87 188 
% no 42 406
% 

% N = tabella di contingenza di partenza
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 
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
Chi2=sum(((N-Ntheo).^2)./Ntheo,'all');


%% 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
Chi2chk1=0;
for i=1:nr
    for j=1:nc
        Chi2chk1=Chi2chk1+(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
% 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.
diff=N(1,1)-Ntheo(1,1);
diffstring=num2str(diff);
if N(1,1)-Ntheo(1,1)>0
    disp('Associazione positiva')
    disp(['La differenza tra n_{11} n e n_{11}^* è =' diffstring])
else
    disp('Associazione negativa')
    disp(['La differenza tra n_{11} n e n_{11}^* =' diffstring])
end

%% Calcolo dell'indice phi utilizzando il valore precalcolato dell'indice Chi2
if det(N)>0
    phi=sqrt(Chi2/n);
else
    phi=-sqrt(Chi2/n);
end

%% Indice \phi tramite la formula
% $\phi=\frac{n_{11}  n_{22} -n_{12}n_{21}}{\sqrt{n_{1.}n_{2.}n_{.1}n_{.2}}}$ 
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)])

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

%% Atro modo
% 6) Verificare che l'indice φ può essere ottenuto come coefficiente di
% correlazione tra i due fenomeni dicotomici X e Y. Data la tabella
% originale di contingenza ricostruire la matrice dei dati di partenza che
% ha prodotto la tabella di contingenza. Verificare che il coefficiente di
% correlazione calcolato sulla matrice dei dati originaria è uguale al
% valore dell'indice φ

% Osservazione: utilizziamo la funzione crosstab2datamatrix dell'FSDA
% toolbox. Per scaricare e installare FSDA, dal menu  Home|Add-Ons|Get
% Add-Ons.
% Nella casella di ricerca della finestra Add-on Explorer digitare FSDA
Xori=crosstab2datamatrix(N);
% Osservazione 
% crosstab(Xori(:,1),Xori(:,2)) mi produce la matrice di contingenza N
disp(corr(Xori))


%% Indice \theta (rapporto dei prodotti incrociati)
th=(N(1,1)*N(2,2))/(N(1,2)*N(2,1));
disp(['Il valore dell''indice \theta è ' num2str(th)])


%% 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
th=0:0.1:200;
U=(sqrt(th)-1)./(sqrt(th)+1);
plot(th,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')