<jupyter_start><jupyter_text>Riferimenti bibliografici:* <NAME>. &  <NAME>. (2012), [Independent Component Analysis: Algorithms and Applications](https://www.cs.helsinki.fi/u/ahyvarin/papers/NN00new.pdf).* <NAME>. & <NAME>. (2008) [Visualizing Data using t-SNE](http://www.jmlr.org/papers/volume9/vandermaaten08a/vandermaaten08a.pdf).  Analisi delle componenti indipendenti (ICA)  Indice 1. [Problema del cocktail party](cocktail)    1.1 [Simulare i dati](simulare)    1.2 [Utilizzare ICA per ottenere i segnali originali](ica)2. [Statlog (Landsat Satellite) Data Set](landsat)    2.1 [Descrizione e lettura del data set](descrizione)    2.2 [Analisi esplorativa](esplorativa)3. [PCA, ICA e t-SNE a confronto](confronto)4. [Ridurre la dimensionalità come parte di una pipeline di classificazione](pipeline)    4.1 [Creare una baseline](baseline)    4.2 [Scegliere in numero di componenti principali](scegliere_principali)    4.3 [Scegliere in numero di componenti indipendenti](scegliere_indipendenti)<jupyter_code>from sklearn.decomposition import FastICA
from sklearn.pipeline import Pipeline
import inspect
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

%load_ext autoreload
%autoreload 2<jupyter_output><empty_output><jupyter_text>1. Problema del cocktail party    1.1 Simulare i dati<jupyter_code>from msbd.datasets import CocktailParty

print(inspect.getsource(CocktailParty))
A = np.array(
    [[0.5, 0.5], 
     [0.2, 0.8]]
)

cp = CocktailParty(A)
t = np.linspace(0, 100, 1000)
S = np.c_[cp.s1(t), cp.s2(t)]
X = np.c_[cp.x1(t), cp.x2(t)]

print("Dimensioni di S: {} X {}".format(*S.shape))
print("Dimensioni di X: {} X {}".format(*X.shape))
from msbd.grafici import grafico_cocktail_party

print(inspect.getsource(grafico_cocktail_party))
plt.figure(figsize=(8, 6))

grafico_cocktail_party(t, S, X)

plt.show()<jupyter_output><empty_output><jupyter_text>1.2 Utilizzare ICA per ottenere i segnali originali<jupyter_code>ica = FastICA(whiten=True, random_state=42)

S_hat = ica.fit_transform(X)
plt.figure(figsize=(8, 8))

grafico_cocktail_party(t, S, X, S_hat)

plt.show()<jupyter_output><empty_output><jupyter_text>EsercizioElencare tutte le differenze tra segnali veri e stimati.<jupyter_code>from msbd.preprocessamento import Sbiancare


print(inspect.getsource(Sbiancare))<jupyter_output><empty_output><jupyter_text>Esercizio1. Completare il metodo `inverse_transform()` della classe `Sbiancare` definita in `msbd/preprocessamento/sbiancare.py`;2. Verificare il corretto funzionamento del metodo appena definito.> Suggerimento: usare il fatto che la matrice $\mathbf{E}$ è una matrice [unitaria](https://it.wikipedia.org/wiki/Matrice_unitaria).  Esercizio1. Definire una pipeline `ica_pipe` con al primo step `Sbiancare()` e al secondo step `FastICA(whiten=False)`;2. Ottenere `S_hat`come `S_hat = ica_pipe.fit_transform(X)`;3. Fare nuovamente il grafico *cocktail party* utilizzando la stima di $\mathbf{S}$ ottenuta al passo precedente;4. Commentare i risultati.  2. [Statlog (Landsat Satellite) Data Set](https://archive.ics.uci.edu/ml/datasets/Statlog+(Landsat+Satellite))    2.1 Descrizione e lettura del data set    Data Set Information:The database consists of the multi-spectral values of pixels in 3x3 neighbourhoods in a satellite image, and the classification associated with the central pixel in each neighbourhood. The aim is to predict this classification, given the multi-spectral values. In the sample database, the class of a pixel is coded as a number. The Landsat satellite data is one of the many sources of information available for a scene. The interpretation of a scene by integrating spatial data of diverse types and resolutions including multispectral and radar data, maps indicating topography, land use etc. is expected to assume significant importance with the onset of an era characterised by integrative approaches to remote sensing (for example, NASA's Earth Observing System commencing this decade). Existing statistical methods are ill-equipped for handling such diverse data types. Note that this is not true for Landsat MSS data considered in isolation (as in this sample database). This data satisfies the important requirements of being numerical and at a single resolution, and standard maximum-likelihood classification performs very well. Consequently, for this data, it should be interesting to compare the performance of other methods against the statistical approach. One frame of Landsat MSS imagery consists of four digital images of the same scene in different spectral bands. Two of these are in the visible region (corresponding approximately to green and red regions of the visible spectrum) and two are in the (near) infra-red. Each pixel is a 8-bit binary word, with 0 corresponding to black and 255 to white. The spatial resolution of a pixel is about 80m x 80m. Each image contains 2340 x 3380 such pixels. The database is a (tiny) sub-area of a scene, consisting of 82 x 100 pixels. Each line of data corresponds to a 3x3 square neighbourhood of pixels completely contained within the 82x100 sub-area. Each line contains the pixel values in the four spectral bands (converted to ASCII) of each of the 9 pixels in the 3x3 neighbourhood and a number indicating the classification label of the central pixel. The number is a code for the following classes:  Number: Class* 1:  red soil * 2:  cotton crop * 3:  grey soil * 4:  damp grey soil * 5:  soil with vegetation stubble * 6:  mixture class (all types present) * 7:  very damp grey soil NB. There are no examples with class 6 in this dataset. The data is given in random order and certain lines of data have been removed so you cannot reconstruct the original image from this dataset. In each line of data the four spectral values for the top-left pixel are given first followed by the four spectral values for the top-middle pixel and then those for the top-right pixel, and so on with the pixels read out in sequence left-to-right and top-to-bottom. Thus, the four spectral values for the central pixel are given by attributes 17,18,19 and 20. If you like you can use only these four attributes, while ignoring the others. This avoids the problem which arises when a 3x3 neighbourhood straddles a boundary.  Attribute Information:The attributes are numerical, in the range 0 to 255.<jupyter_code>CANALI = ["R", "G", "B", "A"]
CLASSI = [
    "red soil", 
    "cotton crop", 
    "grey soil", 
    "damp grey soil", 
    "soil with vegetation stubble",
    "very damp grey soil",
]
NUMERO_CLASSE = dict(zip([1, 2, 3, 4, 5, 7], CLASSI))
NOMI = ["Pixel{}_{}".format(i, c) for i in range(9) for c in CANALI] + ["Classe"]

URL = "https://archive.ics.uci.edu/ml/machine-learning-databases/statlog/satimage/sat.trn"

dati = pd.read_table(URL, sep=" ", names=NOMI)
print("Dimensioni: {} X {}".format(*dati.shape))
dati.head(2)
from sklearn.model_selection import train_test_split

X_train, X_val, y_train, y_val = train_test_split(dati.iloc[:, :-1], dati.iloc[:, -1], 
                                                  test_size=0.25, random_state=42)<jupyter_output><empty_output><jupyter_text>2.2 Analisi esplorativa<jupyter_code>dati.iloc[:, :-1].describe()
y_train.replace(NUMERO_CLASSE).value_counts()
from msbd.grafici import grafico_immagine_satellitare

print(inspect.getsource(grafico_immagine_satellitare))
np.random.seed(42)

plt.figure(figsize=(8, 6))

grafico_immagine_satellitare(X_train, y_train)

plt.show()
from msbd.grafici import diagrammi_scatola_baffi_pixel

print(inspect.getsource(diagrammi_scatola_baffi_pixel))
plt.figure(figsize=(10, 10))

diagrammi_scatola_baffi_pixel(X_train, y_train, 4)

plt.show()
from msbd.grafici import diagrammi_scatola_baffi_classi

print(inspect.getsource(diagrammi_scatola_baffi_classi))<jupyter_output><empty_output><jupyter_text>Esercizio1. Stimare le prime due componenti ICA di `X_train` e assegnarle alla variabile `X_ica`;2. Utilizzare la funzione `diagrammi_scatola_baffi_classi` per fare i diagrammi a scatola e baffi delle classi sia per la prima che per la seconda componente ICA.> Nota: `x=y_train` e `y=X_ica[:, 0]` per la prima componente (pensare a cosa si vuole vedere sull'asse x del grafico).  3. PCA, ICA e t-SNE a confronto<jupyter_code>from sklearn.decomposition import PCA
from sklearn.manifold import TSNE

np.random.seed(42)

ica = FastICA(n_components=2, whiten=True)
pca = PCA(n_components=2, whiten=True)
tsne = TSNE(n_components=2)

X_ica = ica.fit_transform(X_train)
X_pca = pca.fit_transform(X_train)
X_tsne  = tsne.fit_transform(X_train)
from msbd.grafici import grafico_componenti

print(inspect.getsource(grafico_componenti))
plt.figure(figsize=(15, 5))

plt.subplot(131)
plt.title("PCA")
grafico_componenti(X_pca, y_train)

plt.subplot(132)
plt.title("ICA")
grafico_componenti(X_ica, y_train)

plt.subplot(133)
plt.title("t-SNE")
grafico_componenti(X_tsne, y_train)

# plt.tight_layout()
plt.subplots_adjust(top=0.85)
plt.show()<jupyter_output><empty_output><jupyter_text>EsercizioCommentare i risultati.  4. Ridurre la dimensionalità come parte di una pipeline di classificazione    4.1 Creare una baseline<jupyter_code>from sklearn.dummy import DummyClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import  accuracy_score
dc = DummyClassifier(strategy="most_frequent")

dc.fit(X_train, y_train)
y_pred = dc.predict(X_val)

dc_acc = accuracy_score(y_val, y_pred)

print("Accuratezza DummyClassifier(strategy=\"most_frequent\"): {:.1f}%".format(100 * dc_acc))
multinom = LogisticRegression(penalty="l2", multi_class="multinomial", C=float("inf"), solver="newton-cg")

multinom.fit(X_train, y_train)
y_pred = multinom.predict(X_val)

multinom_acc = accuracy_score(y_val, y_pred)

print("Accuratezza DummyClassifier(strategy=\"most_frequent\"): {:.1f}%".format(100 * dc_acc))
print("Accuratezza LogisticRegression(multi_class=\"multinomial\"): {:.1f}%".format(100 * multinom_acc))<jupyter_output><empty_output><jupyter_text>4.2 Scegliere in numero di componenti principali    Scegliere in numero di componenti principali studiando la varianza spiegata<jupyter_code>pca = PCA(whiten=True)
pca.fit(X_train)
from msbd.grafici import grafico_varianza_spiegata

print(inspect.getsource(grafico_varianza_spiegata))
plt.figure(figsize=(15, 3))
plt.subplot(121)
grafico_varianza_spiegata(pca.explained_variance_ratio_)

plt.subplot(122)
grafico_varianza_spiegata(pca.explained_variance_ratio_, cumulata=True)

plt.tight_layout()
plt.show()
N_COMPONENTS = 3

multinom = Pipeline([
    ("pca", PCA(n_components=N_COMPONENTS, whiten=True)),
    ("multinom", LogisticRegression(penalty="l2", multi_class="multinomial", C=float("inf"), solver="newton-cg")),
])

multinom.fit(X_train, y_train)
y_pred = multinom.predict(X_val)

print("Accuratezza DummyClassifier(strategy=\"most_frequent\"): {:.1f}%".format(100 * dc_acc))
print("Accuratezza LogisticRegression(multi_class=\"multinomial\"): {:.1f}%".format(100 * multinom_acc))
print("Accuratezza PCA(n_components={}) + LogisticRegression(multi_class=\"multinomial\"): {:.1f}%".format(N_COMPONENTS, 100 * accuracy_score(y_val, y_pred)))<jupyter_output><empty_output><jupyter_text>Scegliere il numero di componenti principali in base all'errore sull'insieme di validation<jupyter_code>multinom = Pipeline([
    ("pca", PCA(whiten=True)),
    ("multinom", LogisticRegression(penalty="l2", multi_class="multinomial", C=float("inf"), solver="newton-cg")),
])

acc = []

for n_components in range(1, 26):
    multinom.set_params(pca__n_components=n_components)
    multinom.fit(X_train, y_train)
    y_pred = multinom.predict(X_val)
    acc.append(accuracy_score(y_val, y_pred))
plt.figure(figsize=(15, 3))

plt.title("Accuratezza al variare del numero di componenti principali")
plt.plot(range(1, len(acc) + 1), acc, marker='o', ls='--')
plt.plot(np.argmax(acc) + 1, acc[np.argmax(acc)], c="tab:red", marker="o", 
         label="Miglior numero di componenti ({})".format(np.argmax(acc) + 1))
plt.xticks(range(1, len(acc) + 1), range(1, len(acc) + 1))
plt.xlabel("Numero componenti principali")
plt.ylabel("Accuratezza")
plt.legend()
plt.grid()

plt.show()<jupyter_output><empty_output><jupyter_text>EsercizioCambiare il `random_state` in `train_test_split` (nella cella dove abbiamo diviso i dati in `train` e `val`) e rieseguire le celle fino al grafico sopra. Il numero di componenti migliore è lo stesso? Commentare. > Nota: in questo caso sarebbe stato meglio utilizzare la convalida incrociata, in modo da ottenere stime più precise e i relativi errori standard. Una volta imparata la teoria relativa, in futuro considerare l'utilizzo della funzione di scikit-learn [validation_curve](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.validation_curve.htmlsklearn.model_selection.validation_curve).<jupyter_code>N_COMPONENTS = np.argmax(acc) + 1

multinom = Pipeline([
    ("pca", PCA(n_components=N_COMPONENTS, whiten=True)),
    ("multinom", LogisticRegression(penalty="l2", multi_class="multinomial", C=float("inf"), solver="newton-cg")),
])

multinom.fit(X_train, y_train)
y_pred = multinom.predict(X_val)

print("Accuratezza DummyClassifier(strategy=\"most_frequent\"): {:.1f}%".format(100 * dc_acc))
print("Accuratezza LogisticRegression(multi_class=\"multinomial\"): {:.1f}%".format(100 * multinom_acc))
print("Accuratezza PCA(n_components={}) + LogisticRegression(multi_class=\"multinomial\"): {:.1f}%".format(N_COMPONENTS, 100 * accuracy_score(y_val, y_pred)))
from msbd.grafici import grafico_matrice_confusione

print(inspect.getsource(grafico_matrice_confusione))
plt.figure(figsize=(6, 6))

grafico_matrice_confusione(y_val, y_pred, CLASSI)

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1091505650747357, 'fraction_numerical': 0.01844695588771418, 'mean_word_length': 6.044684129429892, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 52, 'https://': 5, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '342090', 'license': 'MIT', 'n_tokens_mistral': 5022, 'n_tokens_neox': 4717, 'n_words': 1305, 'path': '10_analisi_componenti_indipendenti.ipynb', 'repo_name': 'mnslarcher/metodi-statistici-big-data', 'size': '13745'}