sublimina.it

  • Full Screen
  • Wide Screen
  • Narrow Screen
  • Increase font size
  • Default font size
  • Decrease font size






  research gate


omino_logo_sublomina


La Conoscenza è come una linea, dove sono i confini non ci è dato di saperlo.

Sublimina.it è un viaggio personale nel mondo del pensiero umano. Per raggiungere ogni meta c'è una via ed ogni via ha un ingresso. Questa è la mia porta personale, l'ho appena aperta. Ognuno ha la sua porta, qualche volta si ha bisogno, però, di intravedere cosa c'è al di là della porta altrui per mirare l'altrove che sta dietro la propria.  Ispirato da: Franz Kafka, Il processo (1925)


Script Matlab

Il test di un programma può essere usato per mostrarela presenza di bug,ma mai per mostrare la loro assenza.

Edsger Dijkstra, Structured Programming, 1972

Ho selezionato una serie di script Matlab, alcuni elaborati durante i miei studi universitari, gli script sono editabili e liberamente modificabili

 

K-Means, B-SAS, problemi di Clustering con codice MATLAB

E-mail Stampa PDF

In questo articolo intendo presentare alcuni script in linguaggio MATLAB® per comprendere operativamente come funzionano gli algoritmi di "clustering", due in particolare: il "k-means" e il B-SAS. Prima di presentare gli script, come al solito un po di teoria.

1.0 Introduzione

1.1 Problema di k-clustering in logica crisp

1.2 K-means

1.3 Algoritmo B-SAS (Basic Sequential Algorithmic Scheme)

1.4 DOWNLOAD codice


1.0 Introduzione


Al presente, la ricerca scientifica, in particolare nel campo della Computer Science e dell'Atrificial Intelligence ha prodotto una serie di metodi che consentono alle macchine di "apprendere".  Tale termine andrebbe discusso e contestualizzato ampiamente, ciò non sarà fatto in questo articolo, ma il lettore attento potrà avere una idea, dopo aver compreso e provato gli algoritmi, su cosa si intende per "apprendimento automatico". In letteratura la disciplina che si occupa di studiare e sistematizzare le tecniche di apprendimento automatico è nota come "Machine Learning". In tale ambito è andata consolidandosi una serie di tecniche il cui scopo è modellare un qualsiasi processo fisico o astratto tramite i dati a disposizione. In particolare la Pattern Recognition, appartenente all'ambito generale del Datamining si prefigge lo scopo di classificare una serie di oggetti (pattern) in categorie o classi. La classificazione di tali discipline risulta ardua e non sarà discussa, si sottolinea che in letteratura anglosassone esse sono racchiuse nel generico termine di "Soft Computing".

Cioò che quì interessa è comprendere come il modello di un qualsivoglia processo possa essere ricavato dai dati a disposizione per tale processo (modellamento datadriven). Per inciso, la fisica presenta una serie di modelli matematici che spiegano la stragrande maggioranza dei fenomeni: le equazioni del moto sono un esempio. Tale tecica di modellamento può denominarsi analitica, poiché il modello del processo in questione è fornito analiticamente attraverso equazioni. Si possono pensare tuttavia tutta una serie di processi, la cui struttura analitica non è nota per varie ragioni: il processo è altamente complesso, il processo è altamente non lineare. Il modellamento analitico permette di indagare la natura delle leggi, ma ad oggi sappiamo che ciò non è sempre possibile. Al modellamneto analitico si contrappone il modellamento guidato dai dati (datadriven, appunto) dove il processo è considerato come una sorta diu scatola nera dove abbiamo una serie di ingressi ed una serie di uscite (questo accade se il processo è orientato, altrimenti si hanno solo le uscite ed il processo si definisce non orientato). Il processo è modellato tramite un sistema di calcolo capace di "apprendere un compito specifico tramite esempi".

Ricapitolando, il modellamento datadriven permette di sistetizzare il modello di un processo orientato o non tramite un campionamento sui dati, tale modello sostituisce qualsiasi sintesi analitica del processo considerato come una scatola nera (black box).

Nell'ambito del Machine Learning e del modellamento datadriven si distinguono due metodologie, la prima è nota come apprendimento o modellamento supervisionato, la seconda come come apprendimento o modellamento non supervisionato. La differenza è data dalla conoscenza a priori delle uscite del processo dati gli ingressi (quindi tale metodo si applica ai processi orientati, secondo la definizione accennata in precedenza). In altrer parole per alcuni particolari dati di ingresso (campionamento in ingresso) sono noti i dati di uscita (campionamento in uscita) così da poter sintetizzare una generica funzione test che permette di addestrare il modello per nuovi dati di ingresso (pattern) mai visti in precedenza, quindi differenti da quelli del campionamento iniziale. Per essere precisi questi due tipi di insiemi di campionamento hanno un nome specifico: training set e test set, con la proprietà che l'intersezione di tali insiemi sia vuota. Il modellamento non supervisionato non si avvale di tali conoscenze, ma nel modello stesso sono presenti gli ingredienti necessari per poter sintetizzare a dovere il modello di un opportuno processo.

Sintetizzando quanto detto la Pattern recognition attraverso il modellamento datadriven si prefigge lo scopo di operare una classificazione sui dati (classificazione quì è utilizzato come termine generico) ovvero di assegnare una etichetta ad ogni pattern generato dal processo e non presente quando il modello del processo è stato generato. A questo punto, brevemente si può introdurre un altro concetto legato a quanto detto: un sistema di modellamento siffatto possiede una capacità di generalizzazione ed è noto in letteratura come sistema di modellamento induttivo. Anche qui andrebbe aperta una discussione su cosa si intende per generalizzazione o modellamento induttivo in quanto si può finire nel paradosso secondo il quale un sistema computazionale ampiamente riconosciuto come deduttivo possa "creare" conoscenza induttiva grazie alla capacità di generalizzazione. Per approfondire qui ho tentato un approccio a tale discussione.

Dopo questa brevissima introduzione, possiamo definire cosa sia il clustering e quali istanze di problemi esso approccia.

La Cluster Analysis ha lo scopo di partizionare un insieme di dati in un certo numero di sottoinsiemi a seconda delle esigenze, le principali sono riassumibili in:

  • compressione dei dati;
  • generazione di ipotesi, ossia ricerca di possibili correlazioni tra le caratteristiche (feautures) che descrivono ciascun pattern;
  • verifica di ipotesi, ossia la validazione di ipotetiche correlazioni tra feautures;
  • sintesi di classificatori tramite tecniche di modellamento non supervisionato.

Intuitivamente, il numero di sottoinsiemi in cui un certo numero di pattern può essere clusterizzato può essere un dato di ingresso al problema ovvero già stabilito (problema di k-clustering) o da stabilire opportunamente (problema di free-clustering). Questa divisione rispecchia la tipolo9gia di modellamento che nel primo caso è supervisionata, nel secondo caso no, in quanto un opportuno criterio permette di stabilire il numero di cluster ottimale in cui l'insieme dei dati deve essere partizionato. Per completezza diamo la definizione di problema di k-clustering in logica crisp o aristotelica:

1.1 Problema di k-clustering in logica crisp

Dato un insieme di pattern «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»A«/mi»«mo»§#8834;«/mo»«mi»X«/mi»«/math» di cardinalità m ed un intero k>0, determinare una partizione «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»P«/mi»«mo»=«/mo»«mo»{«/mo»«msub»«mi»A«/mi»«mn»1«/mn»«/msub»«mo»,«/mo»«msub»«mi»A«/mi»«mn»2«/mn»«/msub»«mo»,«/mo»«mo».«/mo»«mo».«/mo»«mo».«/mo»«mo»,«/mo»«msub»«mi»A«/mi»«mi»k«/mi»«/msub»«mo»}«/mo»«/math» di «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»A«/mi»«/math» in k sottoinsiemi non vuoti e disgiunti: «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msub»«mi»A«/mi»«mi»i«/mi»«/msub»«mo»§#8800;«/mo»«mo»§#8709;«/mo»«mo»§nbsp;«/mo»«mo»§#8704;«/mo»«mi»i«/mi»«mo»;«/mo»«mo»§nbsp;«/mo»«msub»«mi»A«/mi»«mi»i«/mi»«/msub»«mo»§#8745;«/mo»«msub»«mi»A«/mi»«mi»j«/mi»«/msub»«mo»=«/mo»«mo»§#8709;«/mo»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mo»§#8704;«/mo»«mi»i«/mi»«mo»§#8800;«/mo»«mi»j«/mi»«mo»§nbsp;«/mo»«mi»i«/mi»«mo»,«/mo»«mi»j«/mi»«mo»=«/mo»«mn»1«/mn»«mo»,«/mo»«mn»2«/mn»«mo»,«/mo»«mo».«/mo»«mo».«/mo»«mo».«/mo»«mo»,«/mo»«mi»k«/mi»«/math» in modo che elementi appartenenti allo stesso sottoinsieme siano "simili" ed elementi appartenenti a sottoinsiemi diversi siano "differenti", secondo un predefinito criterio di similitudine. L'intero k, numero di sottoinsiemi in cui «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»A«/mi»«/math» è suddiviso è noto come ordine della partizione P.

Dalla definizione si evince che un problema di clustering ha bisogno di una qualsivoglia misura di similarità. In altre parole i pattern appartenenti all'insieme dei dati vivono in uno spazio topologico normato su cui è definita tale misura. Un esempio classico è la distanza euclidea tra gli elementi appartenenti ad uno spazio multidimensionale. Detto questo possiamo fornire una definizione di cluster:

"Un cluster può essere descritto come una regione connessa di uno spazio multidimensionale caratterizzata da una densità relativamente alta di punti, separata da altre regioni simili da una zona caratterizzata da una bassa densità di punti".

La definizione precedente è in logica crisp, ciò significa che ogni elemento è assegnato ad un solo cluster, tale ipotesi può essere rilassata assegnando ad ogni elemento una misura di appartenenza fuzzy.

Importante è definire come rappresentare un cluster: ciò si può ottenere attraverso il centroide (una sorta di eleemnto medio), o attraqverso un inviluppo. Dopo tale definizione bisogna altresì definire la distanza "punto-insieme", che può essere definita come distanza tra l'elemento e il centroide o tra l'elemento e l'inviluppo. Ciò dipende da scelte progettuali dell'algoritmo di clustering.

1.2 K-means

L'algoritmo k-means permette di trovarare la partizione P dato l'insime dei dati A e l'ordine della partizione k.

In pseudo codice si hanno i seguenti passi:

Initialize mi,  i = 1,…,k, for example, to k random xt

Repeat

For all xt in X

bit <--  1 if || xt - mi || = minj || xt - mj ||

bit <-- 0 otherwise

For all mi,  i = 1,…,k

mi <-- sum over t (bit xt) / sum over t (bit )

Until mi converge

Il vettore m contiene il riferimento al centroide di ogni cluster, x è il generico elemento (pattern di esempio) e b contiene l'etichetta stimata per ogni cluster. La misura di similarità è la distanza euclidea.

Si nota che l'inizializzazione dei centroidi può essere casuale.

La spiegazione dell'algoritmo è la seguente:

  1. scegli secondo qualche criterio (random ad esempio) i parametri di inizializzazione dei centroidi mi dei k cluster.
  2. per ogni pattern di esempio dell'insieme di partenza genera una assegnazione al centroide più vicino (mi).
  3. per ogni cluster ricalcola il centroide mi basandosi sull'assegnazione corrente;
  4. ripeti i punti 2. e 3. fino a convergenza o soddisfacimento di un opportuno criterio di arresto

Un possibile criterio di arresto è stabilire a priori il numero di "epoche" o iterazioni dell'algoritmo. Se tale valore è elevato si può assistere ad una non stabile convergenza, allora si può stabilire un nuovo criterio che si basa sulla distanza δ tra il centroide all'iterazione precedente ed il centroide all'iterazione attuale. In definitiva è usuale mettere le due condizioni in OR logico, al fine di assicurarsi l'arresto.

Il k-means è utilizzato in numerosi problemi come ad esempio la character recognition, speech recognition, e problemi di codifica decodifica dei dati.

Ecco lo script MATLAB® che consente facilmente di testare l'lgoritmo di clustering k-means. Lo script richiama una funzione, riportata in seguito, per la generazione dei cluster.



<p>%prova num.7 clastering data una matrice di pattern in R^2
%--simulazione al variare della dispersione dei punti
%--condizioni di partenza casuali
%--inserimento matrice U(k,n) appartenenza puttern i-esimo al j-esimo
%cluster
%--con condizioni di arresto
clear all;
clc;
%___________________________________________________________________
%inizializzazione problema di k-clustering
N=6
spar=[0,10,30,50,80,100]; %sparpagliamento
 
 
for i=1:N
figure(i)
n=300; %numero di pattern
dim=2;%dimensione pattern (numero di caratteristiche)
MAX=10;%numero massimo iterazioni
delta=0.1;
 
%matrice dei pattern casuale uniforme (n x dim)
% a=0.7.*rand(n,dim);
%  c=4+randn(n,2).*(rand(n,2));
% b=randn(n,2).*(rand(n,2));
%  a=c+b;
 
%generazione pattern tramite funzione
e1=8; %polarizzazione cluster 1
e2=0; %polarizzazione cluster 2
a=cluster_gen(n,e1,e2,spar(i));
 
 
%grafico relarivo alla matrice dei pattern
for s=1:n
plot(a(s,1),a(s,2),'.');
hold on
end
 
 
%inizializzazione randomica
%scelta randomica di k centroidi
k=2; %numero di centroidi (k-clustering)
mu=e1*rand(k,dim); %centroidi casuale (matrice)
 
%grafico centroidi casuali
for s=1:k
plot(mu(s,1),mu(s,2),'*');hold on
end
 
%matrice che colleziona l'appartenenza degli n pattern ai k cluster
U=zeros(k,n);
 
centroide=zeros(k,dim);%inizializzazione matrice centroidi
nk=zeros(1,k); %numero di pattern appartenente al centroide s-esimo
 
%inizializzazione ciclo principale
%_______________________________________________________________________
 
iterazioni=0;   
arresto=0; %variabile "binaria"
 
while iterazioni~=MAX &amp;&amp; arresto==0    
iterazioni=iterazioni+1;
 
for i=1:n
 
for s=1:k     
 
dist(s,i)= sqrt(sum((mu(s,:)-a(i,:)).^2));
provadist=dist';%trasposizione
end % fine ciclo sui k cluster
[valore(i),indice(i)]=min((provadist(i,:)));
 
%assegno ogni pattern ad un centroide in accordo al minimo della
%funzione costo            
U(indice(i),i)=1;
 
end %fine ciclo sugli i pattern
 
%calcolo nuovo centroide
for s=1:k
centroide(s,:)=(U(s,:)*a)/sum(U(s,:));
end
%-----------------------------------------------------------------
%calcolo la differenza tra centroidi a passo "iter.." e "iter..-1"
app=(mu-centroide).^2;%matrice differenza con componenti al quadrato
for s=1:k
intercentroide(s,iterazioni)=sqrt(sum(app(s,:)));
end
%calcolo condizione per criterio di arresto
soglia=(sum(intercentroide))./k;
 
 
%-----------------------------------------------
%criterio di arresto
if (soglia(iterazioni))<delta arresto=1; end;
%-----------------------------------------------
%aggiornamento nuovo centroide
mu=centroide;
 
 
%_________
%plotting centoide      
for s=1:k
plot(centroide(s,1),centroide(s,2),'+r');hold on
end      
%________
 
end   %end ciclo while
hold off </p>

Di seguito si riporta la funzione che genera, letti i parametri di ingresso che sono il numero di pattern n, il centro (e1,e2) la deviazione standard (vettore "spar"), due cluster in «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msup»«mi mathvariant=¨normal¨»§#8477;«/mi»«mn»2«/mn»«/msup»«/math» separati con distribuzione gaussiana.

<p>%funzione che genera  cluster spaziati con distribuzione gaussiana in R^2
 
function c=cluster_gen(s,e1,e2,spar)
 
n=s/2;
a=e1+randn(n,2);
b=e2+randn(n,2);
c=zeros(s,2);
cont=0;
for l=1:s
 
 
if l<s/2  c(l,:)=a(l,:); end;
if l>s/2 cont=cont+1; c(l,:)=b(cont,:); end;
end
 
c=c+spar*rand(s,2);
% for i=1:s
% plot(c(i,1),c(i,2),'.');
% hold on</p>


1.3 Algoritmo B-SAS (Basic Sequential Algorithmic Scheme)

Per il k-means l'insieme dei parametri cosiddetti di training è il numero di iterazioni (o da altro criterio di arresto) e dal fondale parametro k, che definisce il numero si sottoinsiemi in cui partizionare l'insieme dei pattern di partenza.

In molti casi k non è disponibile a priori come parametri di ingresso, ma deve essere calcolato. Algoritmi con tali caratteristiche sono noto come algoritmi di free clustering. Esso è definito come segue:

Dato un insieme di pattern «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»A«/mi»«mo»§#8834;«/mo»«mi»X«/mi»«/math» di cardinalità m determinare un intero k>0 e una partizione «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»P«/mi»«mo»=«/mo»«mo»{«/mo»«msub»«mi»A«/mi»«mn»1«/mn»«/msub»«mo»,«/mo»«msub»«mi»A«/mi»«mn»2«/mn»«/msub»«mo»,«/mo»«mo».«/mo»«mo».«/mo»«mo».«/mo»«mo»,«/mo»«msub»«mi»A«/mi»«mi»k«/mi»«/msub»«mo»}«/mo»«/math» di «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»A«/mi»«/math» in k sottoinsiemi non vuoti e disgiunti: «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msub»«mi»A«/mi»«mi»i«/mi»«/msub»«mo»§#8800;«/mo»«mo»§#8709;«/mo»«mo»§nbsp;«/mo»«mo»§#8704;«/mo»«mi»i«/mi»«mo»;«/mo»«mo»§nbsp;«/mo»«msub»«mi»A«/mi»«mi»i«/mi»«/msub»«mo»§#8745;«/mo»«msub»«mi»A«/mi»«mi»j«/mi»«/msub»«mo»=«/mo»«mo»§#8709;«/mo»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mo»§#8704;«/mo»«mi»i«/mi»«mo»§#8800;«/mo»«mi»j«/mi»«mo»§nbsp;«/mo»«mi»i«/mi»«mo»,«/mo»«mi»j«/mi»«mo»=«/mo»«mn»1«/mn»«mo»,«/mo»«mn»2«/mn»«mo»,«/mo»«mo».«/mo»«mo».«/mo»«mo».«/mo»«mo»,«/mo»«mi»k«/mi»«/math» in modo che elementi appartenenti allo stesso sottoinsieme siano "simili" ed elementi appartenenti a sottoinsiemi diversi siano "differenti", secondo un predefinito criterio di similitudine. L'intero k, numero di sottoinsiemi in cui «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»A«/mi»«/math» è suddiviso è noto come ordine della partizione P.

La definizione è simile al problema di k-clustering solo che qui k è una incognita del problema e non un dato fissato a priori.

Il parametro k spesso è calcolato in base ad opportuni criteri (specificatamente da un parametro strutturale «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»§#952;«/mi»«mo»§#8712;«/mo»«mi»§#920;«/mi»«/math»), ad esempio la massima dimensione consentita per ciascun cluster oppure una specifica legge di espansione del cluster stesso. In generale si dice che sia k che PA sono funzione dei parametri strutturali «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»§#920;«/mi»«/math».

Un seplice algoritmo di free clustering è il B-SAS (Basic Sequential Algorithmic Scheme) noto anche come "regola delta", che si basa su un parametro strutturale definito come una soglia di inclusione di un pattern in un determinato cluster.

Lo pseudo codice per il B-SAS è il seguente:

initialize: cluster list L represented by centroid

FOR i=1 TO m (for any vector xi in A)

{IF (L=Ø) THEN {generate first cluster and take his centroid as xi}

ELSE

{calculate centroid list Liδ which distance from xi is minor or equal to δ}

IF (Liδ=Ø)

THEN {create in L a new cluster end take relative centroid as xi}

ELSE

{put xi in Cj cluster (represented by cj) in Liδ which has minimum distance from xi and refresh the centroid «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msup»«msub»«mi mathvariant=¨bold-italic¨»c«/mi»«mi»j«/mi»«/msub»«mrow»«mo»(«/mo»«mi»n«/mi»«mi»e«/mi»«mi»w«/mi»«mo»)«/mo»«/mrow»«/msup»«mo»=«/mo»«mfrac»«mrow»«msup»«msub»«mi mathvariant=¨bold-italic¨»c«/mi»«mi»j«/mi»«/msub»«mrow»«mo»(«/mo»«mi»o«/mi»«mi»l«/mi»«mi»d«/mi»«mo»)«/mo»«/mrow»«/msup»«mo»§#183;«/mo»«msub»«mi»n«/mi»«mi»j«/mi»«/msub»«mo»*«/mo»«msub»«mi mathvariant=¨bold-italic¨»x«/mi»«mi»j«/mi»«/msub»«/mrow»«mrow»«msub»«mi»n«/mi»«mi»j«/mi»«/msub»«mo»+«/mo»«mn»1«/mn»«/mrow»«/mfrac»«/math» where nj is the cardinality of Cj before the insertion of xi}

}

}

L'algoritmo analizza i pattern in A una sola volta e a seconda del valore di δ, può clusterizzare in numerosi insiemi  (δ piccolo) o in pochi insiemi (δ grande).

In generale sia per il k-clustering che per il free clustering può essere effettuato un cosiddetto campionamento nello spazio dei poarametri «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»§#920;«/mi»«/math», ogni uno dei quali genera una specifica partizione P di A e scegliere quella che rispetta determinate condizioni di ottimalità. L'insermento di tali algoritmi di clustering in "contenitori" di ottimizzazione genera i cosiddetti algoritmi di clustering ottimizzati. Si accenna all'esistenza di numerosi parametri e misure di distanza tra partizioni per caratterizzare il probnlema di clustering.


[al momento è fornito lo script principale nel seguito l'articolo sarà aggiornato con gli indici di validazione e l'indice di compattezza e separabilità di Davis - Bouldin, nonchè il meccanismo di induzione nel clustering Winner Take All (WTA)]


%ALGORITMO BSAS con variazione parametro delta, calcolo numero di cluster
%ottimo
%(codice ripulito)
clear;clc;
 
%prova struct
clear;clc;
%__________________________________________________________________________
 
%Parametro free clustering
%valore limite pari alla diagonale ipercubo unitario se i dati sono 
%normalizzati
delta=[0.125:0.00675/2:sqrt(2)];
 
 
%--------------------------------------------------------------------------
%inizializzazione problema di free-clustering
n=200; %numero di pattern
dim=2;%dimensione pattern (numero di caratteristiche)
 
 
%__________________________________________________________________________
%si vuole dataset normalizzato?
norm=1; %normalizzazione dataset
 
%Parametri generazione cluster
%Numero di cluster
Nc=4;
 
%POLARIZZAZIONI
e11=10; %CL 1 
e21=10; %CL 1
 
e12=0;  %CL 2
e22=0;  %CL 2
 
e13=0;  %CL 3
e23=0;  %CL 3
 
e14=0; %CL 4
e24=0; %CL 4
 
var1=10;
var2=0;
var3=0;
var4=0;
spar=0; %sparpagliamento
 
POL={[e11;e21],[e12;e22],[e13;e23],[e14;e24]};
%__________________________________________________________________________
 
%%%%%%%GENERAZIONE CLUSTER
[a,nCeffett]=cluster_gen(n,Nc,POL,var1,var2,var3,var4,spar,norm);
 
 
n=nCeffett;
%__________________________________________________________________________
 
%grafico relarivo alla matrice dei pattern
figure(1)
for s=1:n
plot(a(s,1),a(s,2),'.');
hold on
end
%hold off
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%CICLO SUL PARAMETRO DELTA%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Lstore={};
Lstore2={};
for m=1:length(delta)
 
m%stampo il numero di ciclo
%inizializzazione LISTA centroidi come aray di celle
L={};
%inizializzo lista centroidi al passo iesimo (listadelta i-esima) come arry
%di celle
listadelta={};
 
%parametri cdi controllo non necessari
indicevett={};njvett=[];
cj=[];
%--------------------------------------------------------------------------
 
%%%%%%%%%%%%%%%%%%%%%%%%Algoritmo di clustering BSAS%%%%%%%%%%%%%%%%%%%%%%%
for i=1:n
%-----
%conmtrollo se L è vuota
if (isempty(L)) 
 
L(i,:)={[a(i,:)],[a(i,:)]}; 
%sprintf('ciclo %d,lista L vuota \ncreo primo centroide \n',i)
%------
else     
%inizializzazione parametri controllo distanza tra pattern e
%centroidi in L
distanza=0;
indice=0;
distanzaList=[];
listadelta={};
 
%parametri di controllo non necessari
contDistminDelta=0;
nj=0;
 
%ciclo sul numero di cluster in L al passo i-esimo
for k=1:size(L,1)  
 
%valutazione distanza tra pattern e centroide in L
distanza=sqrt(sum((a(i,:)-L{k}).^2));
 
%   sprintf('ciclo %d,distanza %d,indice %d \n',i,distanza,indice)          
 
 
if (distanza<delta(m)) 
%      sprintf('ciclo %d,distanza<delta \n',i)
 
contDistminDelta=contDistminDelta+1;
 
%listadelta(contDistminDelta,:)=L(k,:)
listadelta(k,:)=L(k,:);
end   
end
 
for u=1:size(listadelta,1)
if (~cellfun(@isempty,listadelta(u)))
distanzaList(u)=sqrt(sum((a(i,:)-listadelta{u}).^2));
 
else distanzaList(u)=1*10^10;
end
end
[distminima indice]=min(distanzaList);    
indicevett(i)={indice};
 
 
%controllo se la lista dei centroidi al passo i-esimo a distanza 
%minore di delta è vuota 
%------ 
if (isempty(listadelta)) %L=[L;a(i,:)]; 
%  sprintf('ciclo %d,listadelta vuota \naggiungo centroide a L \n',i)
 
%se vuota genero nuovo cluster in L con pattern attuale
%(i-esimo)
L=[L;{a(i,:),a(i,:)}];
 
%------       
else
%1)se non vuota lista delta passo i-esimo, calcolo il centroide più
%vicino al pattern attuale (i-esimo)
%2)inserisco il pattern attuale nella lista L corrispondente al
%3)centroide calcolato
%4)aggiorno il rappresentante attraverso la media dei pattern in lista 
 
%parametri non necessari
nj=size(L{indice,2},1);
njvett(i)=size(L{indice,2},1);
 
 
%aggiunta pattern a(i,:) alla matrice centroidi nella lista
%L(indice)
L{indice,2}=[L{indice,2};a(i,:)];
%calcolo il nuovo centroide come la media degli N=size(L{1,2},1)
%pattern i L{indice }
% cj=sum(L{1,2})/size(L{1,2},1);sbagliato
cj=sum(L{indice,2})/size(L{indice,2},1);
%aggiornamento centroide;
L{indice,1}=cj;     
 
end
end
 
end  %fine algoritmo BSAS (i)     
 
%calcolo il numero di centroidi al passo m (M=length(delta))        
NcentroidiDelta(m)=size(L,1);
% if NcentroidiDelta(m)==2 
%     figure(1)
%    plot(L{1,1}(1,1),L{1,1}(1,2),'r*')
%    plot(L{2,1}(1,1),L{2,1}(1,2),'r*')
% end
 
Lstore(m)={L(:,1)};
Lstore2(m)=[{L}];
 
end %fine ciclo principale su delta (m)
 
%lista per numero di centroidi (dimensione lista L) al variare del
%parametro delta
NumcentrDelta={};
for r=1:length(NcentroidiDelta)
NumcentrDelta(r,1)={NcentroidiDelta(r)};
NumcentrDelta(r,2)={delta(r)};
end
 
 
 
%Plotting Numero di centroidi al variare di delta
figure(2)
stem([NumcentrDelta{:,2}],[NumcentrDelta{:,1}])
title('Numero centroidi vs DELTA')
xlabel('DELTA')
ylabel('NUMCentroidi')
 
%calcolo il numero di occorrenze di uno stesso delta rispetto al numero di cluster
counter=[];
for j=1:max(NcentroidiDelta)
counter(j)=length(find(NcentroidiDelta==j));
 
%u=max(NcentroidiDelta)-j;
 
if counter(j)~=0;
 
%min(find(NcentroidiDelta==j))
prv(j)=min(find(NcentroidiDelta==j)); 
assex(j)=NcentroidiDelta(min(find(NcentroidiDelta==j)));
 
end
 
end
 
%plotting numero di occorrenze di uno stesso delta rispetto al numero di
%cluster
figure(3)
plot(assex,counter,'*')
title('occorrenze stesso delta vs Numero cluster')
xlabel('Numero cluster')
ylabel('occorrenze delta')
 
 
 
vettMass=counter;
vettMass(1)=0;
 
if var1==0 v1=0;
else v1=1;end
if var2==0 v2=0;
else v2=1;end
if var3==0 v3=0;
else v3=1;end
if var4==0 v4=0;
else v4=1;end
 
vsomma=v1+v2+v3+v4;
 
[max ind]=max(vettMass)
 
indelta=find([NumcentrDelta{:,1}]==ind);
 
fprintf('il numero di pattern è: %d \n',n)
fprintf('il numero di cluster effettivo è: %d \n',vsomma)
fprintf('il numero di cluester stimato è: %d \n',ind)
 
 
 
fprintf('per un intervallo del parametro delta pari a: %d -- %d \n',delta(min(indelta)),delta(-min(-indelta)))
fprintf('Cho azzeccato??!??!?')
 
for r=1:ind
figure(1)
plot(Lstore{min(indelta)}{r,1}(1),Lstore{min(indelta)}{r,1}(2),'r*')
 
end
hold off
 
%sistemazione Lista centroidi(eliminazione duplicati)
Lstore3=Lstore2(1);
cont=1;
 
 
for i=2:size(Lstore2,2)
 
if size(Lstore2{i},1)~=size(Lstore2{i-1},1)
cont=cont+1;
Lstore3(cont)=Lstore2(i);
end
 
end
save clusterdeltac Lstore NumcentrDelta Lstore3 a 




DOWNLOAD Script in formato compresso!



[*] Il presente materiale è una rielaborazione personale delle dispense del cosrso di "Soft Computing" tenuto dal prof. Antonello Rizzi presso la facoltà di Ingegneria dell'università La Sapeinza di Roma


[**] Lettura consigliata: Konstantinos Koutroumbas, Sergios Theodoridis, Pattern Recognition, 2006, Elsevier USA.

Ultimo aggiornamento Martedì 20 Settembre 2016 20:38

Discesa del gradiente con script MATLAB

E-mail Stampa PDF

steepest descendTra le tecniche di ottimizzazione di tipo locale si utilizza spesso il metodo della "discesa del gradiente" (steepest descend). Essa consiste in un metodo che, data una funzione multidimensionale permette di trovare un minimo locale.

Questo articolo, dopo una breve descrizione matematica del metodo, intende fornire uno script MATLAB® che tramite l'algoritmio della discesa del gradiente permette di calcolare il minimo locale di una funzione inserita come parametro (lo script ne presenta una versione semplificata).

Il metodo consiste, data una funzione a più variabili nello scegliere dei valori a caso per queste e calcolarne il valore della funzione. Perturbando tali valori iniziali di una quantità minima epsilon si ricalcola il valore della funzione. La differenza tra valore iniziale della funzione e tale nuovo valore fornisce la derivata discreta della funzione in ciascuna dimensione.

Data una funzione «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»f«/mi»«mo»(«/mo»«mover»«mi»x«/mi»«mo»§#175;«/mo»«/mover»«mo»)«/mo»«/math» dove «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«menclose notation=¨top¨»«mi»x«/mi»«/menclose»«/math» è un vettore appartenente a «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msup»«mi mathvariant=¨normal¨»§#8476;«/mi»«mi»n«/mi»«/msup»«/math», il problema è formulato come segue: «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mtable columnalign=¨left¨ rowspacing=¨0¨»«mtr»«mtd»«msub»«mover»«mi»x«/mi»«mo»§#175;«/mo»«/mover»«mrow»«mi»m«/mi»«mi»i«/mi»«mi»n«/mi»«/mrow»«/msub»«mo»=«/mo»«mi»a«/mi»«mi»r«/mi»«mi»g«/mi»«mo»§nbsp;«/mo»«mi»m«/mi»«mi»i«/mi»«mi»n«/mi»«mo»§nbsp;«/mo»«mi»f«/mi»«mo»(«/mo»«mover»«mi»x«/mi»«mo»§#175;«/mo»«/mover»«mo»)«/mo»«/mtd»«/mtr»«mtr»«mtd»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mi»x«/mi»«mo»§#8712;«/mo»«msup»«mi mathvariant=¨normal¨»§#8476;«/mi»«mi»n«/mi»«/msup»«/mtd»«/mtr»«/mtable»«/math» (1)

Scegliendo «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mo»-«/mo»«mi»f«/mi»«/math» in luogo di «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»f«/mi»«/math» si ottiene il duale problema di massimizzazione.

Si sceglie la direzione dove la  «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»f«/mi»«/math» decresce più velocemente che è la direzione opposta al gradiente «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mo»§#8711;«/mo»«mi»f«/mi»«mo»(«/mo»«mover»«mi»x«/mi»«mo»§#175;«/mo»«/mover»«mo»)«/mo»«/math» . Si sceglie un punto iniziale arbitrario «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msub»«mover»«mi»x«/mi»«mo»§#175;«/mo»«/mover»«mn»0«/mn»«/msub»«/math» e si discende lungo la funzione alla ricerca del minimo attraverso una procedura iterativa riportata in seguito:

«math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msub»«mover»«mi»x«/mi»«mo»§#175;«/mo»«/mover»«mrow»«mi»k«/mi»«mo»+«/mo»«mn»1«/mn»«/mrow»«/msub»«mo»=«/mo»«msub»«mover»«mi»x«/mi»«mo»§#175;«/mo»«/mover»«mi»k«/mi»«/msub»«mo»-«/mo»«msub»«mi»§#955;«/mi»«mi»k«/mi»«/msub»«mo»§#8711;«/mo»«mi»f«/mi»«mo»(«/mo»«mover»«msub»«mi»x«/mi»«mi»k«/mi»«/msub»«mo»§#175;«/mo»«/mover»«mo»)«/mo»«/math» (2).

«math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msub»«mi»§#955;«/mi»«mi»k«/mi»«/msub»«/math» è noto come parametro di learning (learning rate) in accordo con la letteratura che utilizza tale metodo nel machine learning.

Lo script presentato accetta una funzione, quella da minimizzare, un valore iniziale x=0.9, il passo di discesa che nello script è chiamato delta (il nostro «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»§#955;«/mi»«/math») ed una tolleranza tolerance=0.000001 (tali parametri sono settabili).

Inizialmente è computata la derivata prima della funzione (tramite il metodo simbolico "inline" di matlab) e data in pasto alla funzione "gradiente()" che ha come parametri di ingresso: il punto iniziale, la derivata,il delta e la tolleranza.

La funzione "gradiente()" restituisce il valore di x dove il gradiente è minimo.

Nello script principale tale valore è utilizzato per il calcolo del minimo e per il grafico dei risultati.

La funzione "gradiente()" utilizza un ciclo whille con condizione di uscita: finché il valore della derivata (in valore assoluto) nel punto x passato come parametro è maggiore della tolleranza passata come paramentro ricalcola il punto x mediante la (2).

Il metodo presentato non ha una convergenza molto veloce ma riesce comunque a trovare il minimo della funzione. Come si nota il learning rate è fisso, mentre altri metodi permettono di calcolarne il valore ottimale, di fatto nella (2) si nota che in generale esso può variare in finzione del passo di iterazione.

La tolleranza è un parametro che serve a definire il criterio di arresto. Il tuning dei parametri può essere complicato se la funzione è irregolare.

Quindi lo script si compone di un main principale e di una funzione.

L'esempio è pensato per una funzione monodiemnsionale tipo:«math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»f«/mi»«mo»(«/mo»«mi»x«/mi»«mo»)«/mo»«mo»=«/mo»«msup»«mi»e«/mi»«mrow»«mo»-«/mo»«mi»x«/mi»«/mrow»«/msup»«mo»+«/mo»«msup»«mi»x«/mi»«mn»4«/mn»«/msup»«/math» di cui si intende calcolare il minimo come mostra il grafico di output seguente dove l'asterisco '*' indica il minimo di tale funzione effettivamente calcolato.


diescesa gradiente

Questo è lo script principale:


%%%%inizia da qui lo script
 
%Esempio di applicazione dell'algoritmo di discesa del gradiente per il
%calcolo del minimo di una funzione:
%y(x)=exp(-x)+x^4
%la function richiamata è:
% [x,h]=gradiente(x0,fprime,delta,tolerance);
 
%Testata con Matlab 2015a
 
 
clear;clc;
%--------------------------------------------------------------------------
X=-1:0.01:2; %inizializzo vettore
f=exp(-X)+X.^4; %vettorizzo la funzione
figure(1)
plot(X,f,'r','LineWidth',1.8)%plotting funzuone di partenza (f)
legend('f')
hold on
 
 
syms x   %ambiente simbolico, x variabile simbolica
 
funct=exp(-x)+x^4;
fprime=diff(funct); %Calcolo la derivata simbolica
 
 
 
%%%plot funzione derivata
plot(X,eval(subs(fprime,X)),'--')
legend('f','fprime')
%__________________________________________________________________________
%definizione parametri algoritmo del gradiente
delta=0.1; %passo di discesa
tolerance=0.000001; %tolleranza(per la condizione di stop)
x0=0.9;%punto di partenza
%applico funzione
[x,h]=gradiente(x0,fprime,delta,tolerance);
%verifica risultati
 
 
%funct(x)
plot(x,eval(subs(funct,x)),'*','MarkerSize',10)
 
%%%%%finisce qui lo script

e questa è la funzione "gradiente()":

%%%inizia da qui la funzione
 
%funzione che applica il metodo del gradiente per il calcolo del minimo di
%una funzione passata come parametro
 
%x0=punto iniziale
%fprime=derivata funzione
%delta=passo di discesa del gradiente
%tollerance=tolleranza per il criterio di arresto
 
%h=numero di iterazioni
function [x,h]=gradiente(x0,fprime,delta,tolerance)
 
x=x0;
h=0
while (abs( eval(subs(fprime,x)) >tolerance))%condizione di arresto
x=x-delta*eval(subs(fprime,x))
h=h+1
end
xopt=x;iter=h;
 
%%%finisce qui la funzione
Ultimo aggiornamento Martedì 20 Settembre 2016 20:44

Campionamento Statistico: uno studio con MATLAB

E-mail Stampa PDF

In questo articolo sono considerate acquisite alcune conoscenze basilari di statistica matematica. Si intende altresì fornire un esempio al calcolatore, tramite uno script nel linguaggio di programmazione MATLAB®, di un campionamento statistico effettuato tramite il metodo Monte Carlo e studiarne le caratteristiche asintotiche.

Prima del codice un po' di teoria.

Sia dato un Fenomeno aleatorio F ed una variabile aleatoria (v.a.) x definita su quest'ultimo. Si può considerare un differente Fenomeno Aleatorio composito FN ottenuto ripetendo N volte F.

Le ripetizioni di F sono statisticamente indipendenti e alla generica ripetizione è possibile associare la v.a. xn con n = 1,...,N.

Ciascuna v.a. xn segue la distribuzione di probabilità di x, quindi la funzione di ripartizione di xn è pari a quella di x«math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msub»«mi»D«/mi»«msub»«mi»x«/mi»«mi»n«/mi»«/msub»«/msub»«mo»(«/mo»«mi»§#955;«/mi»«mo»)«/mo»«mo»=«/mo»«msub»«mi»D«/mi»«mi»x«/mi»«/msub»«mo»(«/mo»«mi»§#955;«/mi»«mo»)«/mo»«/math»

L'esperimento che genera il  fenomeno aleatorio composito FN prende il nome di campionamento statistico (semplice essendo gli esperimenti indipendenti e ripetuti in condizioni immutate) mentre la generica determinazione della v.a. N-dimensionale x1,...,xn si da il nome di campione statistico.

E' possibile mostrare come fissato un valore λ si definisce una partizione xnλ e i valori xn > λ. Se si indica con «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»§#947;«/mi»«mo»(«/mo»«mi»§#955;«/mi»«mo»)«/mo»«/math» il numeo di v.a. xn per le quali, nell'esperimento, si è ottenuto xnλ si ha che la quantità «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msub»«mover»«mi»D«/mi»«mo»^«/mo»«/mover»«mi»x«/mi»«/msub»«mo»(«/mo»«mi»§#955;«/mi»«mo»)«/mo»«mo»=«/mo»«mfrac»«mrow»«mi»§#947;«/mi»«mo»(«/mo»«mi»§#955;«/mi»«mo»)«/mo»«/mrow»«mi»N«/mi»«/mfrac»«/math» prende il nome di poligono somma ed è una stima della distribuzione «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msub»«mi»D«/mi»«mi»x«/mi»«/msub»«mo»(«/mo»«mi»§#955;«/mi»«mo»)«/mo»«/math» della v.a. x.

Si può dimostrare che il poligono somma converge in probabilità (cioè secondo l'espressione forte della Legge dei Grandi Numeri) alla funzione di distribuzione:

«math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»p«/mi»«mo».«/mo»«mo»§nbsp;«/mo»«munder»«mo largeop=¨true¨»lim«/mo»«mrow»«mi»n«/mi»«mo»§#8594;«/mo»«mo»§#8734;«/mo»«/mrow»«/munder»«mo»§nbsp;«/mo»«msub»«mover»«mi»D«/mi»«mo»^«/mo»«/mover»«mi»x«/mi»«/msub»«mo»(«/mo»«mi»§#955;«/mi»«mo»)«/mo»«mo»=«/mo»«msub»«mi»D«/mi»«mi»x«/mi»«/msub»«mo»(«/mo»«mi»§#955;«/mi»«mo»)«/mo»«/math».

Si può dimostrare anche che le quantità:

«math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msubsup»«mover»«mi»m«/mi»«mo»^«/mo»«/mover»«mi»x«/mi»«mrow»«mo»(«/mo»«mi»k«/mi»«mo»)«/mo»«/mrow»«/msubsup»«mo»=«/mo»«mfrac»«mn»1«/mn»«mi»N«/mi»«/mfrac»«msubsup»«mo»§#8721;«/mo»«mrow»«mi»n«/mi»«mo»=«/mo»«mn»1«/mn»«/mrow»«mi»n«/mi»«/msubsup»«msubsup»«mi»x«/mi»«mi»n«/mi»«mi»k«/mi»«/msubsup»«/math» (1)

(il secondo termine di questa espressione è noto come media campionaria o in lingua inglese statistical average)

sono una stima dei momenti «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msubsup»«mi»m«/mi»«mi»x«/mi»«mrow»«mo»(«/mo»«mi»k«/mi»«mo»)«/mo»«/mrow»«/msubsup»«/math» della v.a. x e sono note come caratteristiche del campione, in altre parole il valore atteso «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»E«/mi»«mo»{«/mo»«msubsup»«mover»«mi»m«/mi»«mo»^«/mo»«/mover»«mi»x«/mi»«mrow»«mo»(«/mo»«mi»k«/mi»«mo»)«/mo»«/mrow»«/msubsup»«mo»}«/mo»«mo»=«/mo»«msubsup»«mi»m«/mi»«mi»x«/mi»«mrow»«mo»(«/mo»«mi»k«/mi»«mo»)«/mo»«/mrow»«/msubsup»«/math»converge al valore vero «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msubsup»«mi»m«/mi»«mi»x«/mi»«mrow»«mo»(«/mo»«mi»k«/mi»«mo»)«/mo»«/mrow»«/msubsup»«/math» con probabilià 1.

A titolo di esempio, è noto  che la stima del  momento del primo ordine (k=1) corrisponde alla stima del valore medio della v.a. x: «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msubsup»«mover»«mi»m«/mi»«mo»^«/mo»«/mover»«mi»x«/mi»«mrow»«mo»(«/mo»«mn»1«/mn»«mo»)«/mo»«/mrow»«/msubsup»«mo»=«/mo»«msub»«mover»«mi»§#951;«/mi»«mo»^«/mo»«/mover»«mi»x«/mi»«/msub»«/math» (valor medio campionario). Per la varianza del valore medio stimato si ha:

«math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»v«/mi»«mi»a«/mi»«mi»r«/mi»«mo»{«/mo»«msub»«mover»«mi»§#951;«/mi»«mo»^«/mo»«/mover»«mi»x«/mi»«/msub»«mo»}«/mo»«mo»=«/mo»«mfrac»«msubsup»«mi»§#963;«/mi»«mi»x«/mi»«mn»2«/mn»«/msubsup»«mi»N«/mi»«/mfrac»«/math» (2)

che al limite per N tendente ad infinito, tende a 0.

E' possibile ottenere simili considerazini per la varianza campionaria definita come: «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msubsup»«mover»«mi»§#963;«/mi»«mo»^«/mo»«/mover»«mi»x«/mi»«mn»2«/mn»«/msubsup»«mo»=«/mo»«msubsup»«mover»«mi»m«/mi»«mo»^«/mo»«/mover»«mi»x«/mi»«mrow»«mo»(«/mo»«mn»2«/mn»«mo»)«/mo»«/mrow»«/msubsup»«mo»-«/mo»«mover»«msup»«msub»«mi»§#951;«/mi»«mi»x«/mi»«/msub»«mn»2«/mn»«/msup»«mo»^«/mo»«/mover»«/math», facendo attenzione alla polarizzazione.

Riassumendo, possiamo immaginare un certo numero di ripetizioni di uno stesso esperimento. Tale esperiemnto è caratterizzato da una v.a. con una propria distribuzione e propri parametri statistici. E' possibile mostrare che (1) è uno stimatore non distorto (la varianza dello stimatore tende a 0 nella (2)) del parametro (valore medio nel caso precedente) noto anche come strimatore efficiente.

Lo script presentato in questo articolo consiste in una generazione (sorteggio) secondo il metodo montecarlo di N valorii con distribuzione pseudo-random con una data media ed una data deviazione standard.

Il metodo montecarlo, per inciso, è un metodo statistico (non parametrico) per computare stime tramite simulazioni. Esso è utilizzato in ambito scientifico e ingegneristico data la capacità di risolvere problemi che altrimenti avrebbero un costo computazionale elevato. Esso consiste, in linea generale, in un algoritmo che genera tutte le realizzazioni possibili di un fenomeno aleatorio datane la densità di probabilità e di un metodo di misura sul campione casuale generato dei parametri statistici tanto più valido quanto più le medie statistiche si avvicinano ai valori veri.

Lo script è inizializzato con un numero pari a N=100 ed un numero pari a Nrip=10. N è la dimensione del vettore random, in altre parole l'esperimento aleatorio consiste nella generazione di 100 v.a. indipendenti a distribuzione pseudo-normale, e Nrip è il numero di ripetizioni dell'esperimento aleatorio che fornisce una idea dell'accuratezza del metodo. Nell'inizializzazione sono anche settabili deviazione standard e valore medio dellla densità di probabilità utilizzata per generare gli N valori casuali.

Lo script è strutturato in modo da calcolare la media e la varianza in maniera ricorsiva. Un  ciclo "for" sul numero di ripetizioni Nrip genera Nrip vettori pseudo-random e per ogni uno di questi è calcolata progressivamente la media conservando il vaore attuale in una variabile ('mean'), valore che sarà utilizzato per il prossimo ciclo, fino al termine. Stessa cosa con una formula leggermente più complicata per la varianza.

L'output finale è costituito da due grafici.

Il primo mostra la media campionaria dei 10 esperimenti aleatori costituiti da 100 sorteggi da una distribuzione pseudo-normale, il secondo mostra il valore quadratico medio campionario.


valor quadratico medio campionario metodo montecarlomedia campionaria metodo montecarlo

Accuratezza del metodo Monte Carlo. Andamenti del valor medio campionario e del valor quadratico medio campionario misurato fino a 100 sorteggi da una v.a. a distribuzione pseudo-normale

Si nota che gli andamenti tendono ad avvicinarsi man mano che il numero di esperimenti aumantano, ciò mostra graficamente la regolarità statistica dell'esperimento.

Si può notare, altresì, che gli andamenti dati dalle medie statistiche si stabilizzano intorno ai momenti che misurano al crescere del numero N di ripetizioni. La larghezza delle macchie (distanza tra le curve) rappresenta invece la dispersione degli errori di misura che tendono a diminuire al crescere di N.

%Esempio Metodo Montecarlo
%Generazione di "Nrip" vettori pseudorandom normali a media "niN" e deviazione
%standard "sigmaN"
%Calcolo della media e della varianza con algoritmo ricorsivo
%plotting del risultato della stima della media e della varianza per ogni
%ripetizione %Nrip" dell'esperimento.
%Libro Scarano, Teoria fenomeni aleatori
clc;clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%numero Nrip di ripetizioni esperimento: sorteggio di N valori da una
%distribuzione Pseudo-normale
Nrip=10;
%dati per distribuzione pseudo-Normale
N=100;
niN=0;
sigmaN=1;
 
for j=1:Nrip
X=niN+sigmaN*randn(1,N);
%-------------------------------------------------------------------------
 
%calcolo media ricorsiva
mean=0;
varianza=zeros(1,length(X));
for i=1:length(X)
%indice che parte da 0
l=i-1;
%accumulatore di di ciclo temporaneo
mean=(l*mean+X(i))/(l+1);
%accumulatore vettoriale (array di medie al passo i-esimo)
media(i)=mean;   
 
if i>1
varianza(i)=((1-(1/l))*varianza(l))+((1+l)*(media(i)-media(l))^2);
end
 
end
 
%plotting
figure(1)
plot(media)
title('media campionaria su Nrip')
xlabel('Nrip')
ylabel('media')
hold on
figure(2)
plot(varianza)
title('valor quadratico medio campionario su Nrip')
xlabel('Nrip')
ylabel('Varianza')
hold on
 
end
hold off
 
%codice di verifica che memorizza la media in un vettore
for i=1:length(X)
%indice che parte da 0
l=i-1;
 
 
%accumulatore di di ciclo temporaneo
mean=(l*mean+X(i))/(l+1);
%accumulatore vettoriale (array di medie al passo i-esimo)
media(i)=mean;



Il presente materiale è una personale rielaborazione di

Gaetano Scarano, "Segnali, Processi Aleatori, Stima", Università La Sapienza, Roma, 2009, edizione 1;

pertanto qualsiasi errore è soltanto mio.

Ultimo aggiornamento Martedì 20 Settembre 2016 20:55

Un armonica a 440 Hz con codice Matlab

E-mail Stampa PDF

Volete suonare attraverso la matematica?

In questo breve articolo mostrerò come sia semplice attraverso poche righe di codice MATLAB® riprodurre il “tu..tu” del segnale di “libero” al telefono (in Italia), cioè un LA a 440 Hz.

In verità se mandate in run il codice ci si accorge che il suono è “sporco”, di fatto all’armonica fondamentale sono aggiunte alcune armoniche con fase ed ampiezza randomiche.

Andiamo a vedere come funziona il codice.

La frequenza di campionamento è stata scelta pari a FS=8000 campioni al secondo.

Generiamo due vettori (a, phi) randomici con distribuzione di probabilità gaussiana e uniforme, rispettivamente, con Narm=20 componenti: il primo per le ampiezze ed il secondo per le fasi delle armoniche secondarie.

Il vettore a è costruito in modo che la prima componente a(1) contenga l’ampiezza (pari a 1) della frequenza fondamentale, cioè il LA a 440 Hz, le altre contengono le ampiezze casuali.

Il vettore della fase contiene 20 valori random a distribuzione uniforme,come appena detto.

Prima di effettuare il ciclo che computa i valori del segnale è generato un vettore x di 8000 componenti che sarà la base temporale: il valore iniziale è 1/8000 mentre quello finale è 1 con passo pari proprio ad 1/8000.

Il ciclo “for”  avanza sul numero di armoniche Narm e per ogni una computa un vettore che contiene il segnale sinusoidale a quell’armonica precisa. In altre parole per ogni indice di armonica i è computato il vettore contenente il segnale campionato dell’armonica.


Si nota che la funzione matematica utilizzata per generare il tono è la funzione trigonometrica “seno”:

m(i,:)=a(i)*sin(2*pi*f(i)*x+phi(i));%segnale sinusoidale


Il ciclo “for” restituisce quindi una matrice contenente per righe tali segnali.

Una semplice somma finale sulle righe restituisce l’andamento nel tempo x del segnale totale  “sporcato” con le armoniche con fase e ampiezza randomiche.

Lo script Matlab restituisce tramite la funzione nativa: wavplay(y,Fs) il segnale alla scheda audio del nostro PC, che provvederà alla conversione digitale/analogico per mandare in vibrazione gli altoparlanti.


Lo script inoltre restituisce due grafici:


L’analisi DFT del segnale e L’andamento temporale visualizzato come sequenza.


E’ possibile notare come l’armonica fondamentale sia a 440 HZ, le altre sono di ordine superiore (multipli) e sono state generate tramite  f(i)=i*440 nel ciclo “for”.


%generazione, suono e dft (FFT) di una armonica a 440Hz campionata a
%Fs=8000 sample
clear;
clc;
%____________________________________________________________________
Fs=8000;
Narm=20;
 
%casualizzazione fase e ampiezza
%ampiezza
a=0.02*randn(2,Narm);
a(1)=1;%prima armonica deterministica
%fase
phi=2*pi*rand(1,Narm);
 
x=(1/8000):1/8000:(1); %vettore
 
for i=1:Narm
f(i)=i*440;
 
m(i,:)=a(i)*sin(2*pi*f(i)*x+phi(i));%segnale sinusoidale
 
end
 
 
 
y=sum(m,1);
 
figure(1)
stem(y)
title('andamento temporale segnale x[n]')
xlabel('campioni n')
ylabel('x[n]')
 
 
Y=fft(y);
figure(2)
plot(real(Y))
title('DFT del segnale Y')
xlabel('campioni n')
ylabel('DFT  Y')
%Il grafico della DFT mostra, per un segnale sinusoidale campionato a
%Fs=8000 sample/sec, due impulsi centrati, uno a 440, l'altro a -440
%Mentre se il grafico permetteva le frequenze negative si sarebbero avuti 2
%impulsi, uno a 440, l'altro a -440
wavplay(y,Fs)
Ultimo aggiornamento Martedì 20 Settembre 2016 21:04

Studio di un Processo Aleatorio Bianco con la correlazione (con codice MATLAB)

E-mail Stampa PDF

Grafico correlazione in matlab di un segnale gaussianoIn quest’articolo ci occuperemo di un fatto squisitamente matematico, utile nelle applicazioni relative all’Ingegneria dell’Informazione. Si intende studiare brevemente la correlazione di due segnali, da un punto di vista qualitativo, senza entrare nei dettagli  matematici. Dopo aver fornito una breve descrizione di questo operatore, nel caso continuo e nel caso discreto, sarà fornito un esempio di codice nel linguaggio di programmazione matematico-scientifico MATLAB®, modificabile.

La correlazione, in linea di principio, è un operatore matematico che misura la “somiglianza” di due segnali. Per segnale si intende un qualcosa che vari nel tempo, in generale rispetto ad una variabile indipendente. Il segnale vocale scambiato tra due interlocutori coincide con la propagazione di un onda acustica opportunamente modulata dall’apparato fonatorio, onda che consiste nelle variazioni di pressione dell’aria interposta tra gli interlocutori. E’ possibile quindi rappresentare analiticamente tale segnale come una pressione dipendente dal tempo, i.e. s(t) = p(t).

Un tono puro, ad esempio, ha una espressione sinusoidale del tipo:

«math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»p«/mi»«mo»(«/mo»«mi»t«/mi»«mo»)«/mo»«mo»=«/mo»«mi»s«/mi»«mo»(«/mo»«mi»t«/mi»«mo»)«/mo»«mo»=«/mo»«mi mathvariant=¨normal¨»cos«/mi»«mo»(«/mo»«mn»2«/mn»«mi»§#960;«/mi»«msub»«mi»f«/mi»«mn»0«/mn»«/msub»«mi»t«/mi»«mo»+«/mo»«mi»§#966;«/mi»«mo»)«/mo»«/math»

dove «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msub»«mi»f«/mi»«mn»0«/mn»«/msub»«/math» è la frequenza, t è la variabile temporale e «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»§#966;«/mi»«/math» è una fase che serve ad identificare l’origine del sistema di riferimento.

Dati due segnali generici, se essi sono identici (cioè sono lo stesso segnale) allora si parla di autocorrelazione, altrimenti di crosscorrelazione o semplicemente di correlazione.

Nel caso cosiddetto tempo-continuo, dove la variabile indipendente t appartiene al campo dei numeri reali, quindi libera di variare con continuità su un dominio definito in questo campo, la correlazione possiede una espressione integrale del tipo:


«math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msub»«mi»R«/mi»«mrow»«mi»x«/mi»«mi»y«/mi»«/mrow»«/msub»«mo»(«/mo»«mi»§#964;«/mi»«mo»)«/mo»«mo»=«/mo»«msubsup»«mo»§#8747;«/mo»«mrow»«mo»+«/mo»«mo»§#8734;«/mo»«/mrow»«mrow»«mo»-«/mo»«mo»§#8734;«/mo»«/mrow»«/msubsup»«mi»x«/mi»«mo»(«/mo»«mi»t«/mi»«mo»)«/mo»«mi»y«/mi»«mo»(«/mo»«mi»t«/mi»«mo»-«/mo»«mi»§#964;«/mi»«mo»)«/mo»«mi»d«/mi»«mi»t«/mi»«/math»

dove i pedici x e y indicano che è una crosscorrelazione, e la variabile «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»§#964;«/mi»«/math» indica il ritardo temporale tra i due segnali. Si nota, per inciso, che questa è una versione semplificata per segnali appartnenti al campo dei Reali, in quanto se i essi sono complessi bisogna introdurre, per il segnale ritardato, l'operatore di coniugazione.

Per segnali discreti, meglio noti in letteratura tecnica come sequenze, l'espressione della correlazione è differente, pur conservando lo stesso  schema e significato (in letteratura si parla di dualità continuo-discreta):

.

«math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«msub»«mi»R«/mi»«mrow»«mi»x«/mi»«mi»y«/mi»«/mrow»«/msub»«mo»(«/mo»«mi»m«/mi»«mo»)«/mo»«mo»=«/mo»«munderover»«mrow»«mo»§#8721;«/mo»«mi»x«/mi»«mo»[«/mo»«mi»n«/mi»«mo»]«/mo»«mi»y«/mi»«mo»[«/mo»«mi»n«/mi»«mo»-«/mo»«mi»m«/mi»«mo»]«/mo»«/mrow»«mrow»«mi»n«/mi»«mo»=«/mo»«mo»-«/mo»«mo»§#8734;«/mo»«/mrow»«mrow»«mo»+«/mo»«mo»§#8734;«/mo»«/mrow»«/munderover»«/math»

Anche per questa relazione valgono le medesime considerazioni precedenti. Qui il ritardo è un valore discreto, m, di cui la correlazione ne è funzione.

In realtà queste formule, non sono applicabili a "tutti i segnali", ma solo ad una particolare classe, i cosiddetti segnali di energia,  in quanto per l'altra classe, i segnali di potenza, bisogna introdurre alcune modifiche (passaggio a limite). Per gli scopi di questo articolo, non si farà tale distinsione.

Queste formule provengono dall'interpretazione dei segnali come particolari stutture matematiche, i vettori, e le operazioni tra segnali, come operazione tra vettori, nel caso particolare: il prodotto scalare tra segnali.

La correlazione ha importanti proprietà di simmetria, consultabili su molti testi, ed è estremamente utile nel settore IT, in quanto è alla base dei metodi di elaborazione del segnale, qualsiasi esso sia. E' possibile interpretare graficamente la correlazione come l'operazione di far scivolare un segnale sull'altro, tenendo fisso il primo: x(t), e muovendo il secondo: y(t) di un fattore «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»§#964;«/mi»«/math» rispetto all'origine dei tempi di x(t). Per ogni «math  xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»§#964;«/mi»«/math», si computa l'area del prodotto dei due segnali (calcolo dell'integrale), il risultato è prorio la correlazione che è funzione del ritardo «math  xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»§#964;«/mi»«/math». Se i segnali sono identici, l'autocorrelazione computata per «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»§#964;«/mi»«mo»=«/mo»«mn»0«/mn»«/math», ci fornisce quella che è l'energia del segnale, o la potenza, se i segnali appartengono alla classe dei segnali di potenza [1].

Nella seconda parte di questo articolo, si introduce il concetto di processo aleatorio bianco, per poi terminare con un semplice script, nel linguaggio Matlab,  che consente di fare qualche esperimento.

E' noto che la Teoria matematica della Probabilità studia i Fenomeni Aleatori che producono realizzazioni di grandezze fisiche associabili alle cosiddette variabili aleatorie, sia nel contesto continuo che discreto. Essa oltre ad essere una potente branca della matematica, è un utile strumento per la tecnica, con numerosissime applicazioni, dall'elaborazione dei segnali, passando per le lotterie, fino alla teoria delle code che si formano negli uffici postali. Nel campo dell'elaborazione dell'informazione la Teoria della Probabilità permette la modellizazione dei sistemi, e consente di introdurre l'utile distinsione tra segnali deterministici e segnali aleatori.

Senza perdità di generalità, sacrificando il rigore in cambio della semplicità e della chiarezza, possiamo immaginare di misurare la tensione ai capi di un resistore R (reale). La tensione riscontrata è causata dall'agitazione termica degli elettroni nel conduttore che costituisce il resistore. L'agitazione è dovuta allo stato termico del resistore alla data temperatura T. La tensione misurata, sarà molto piccola, ma non costantemente nulla, con andamento temporale tale da riflettere la natura imprevedibile del fenomeno fisico. Importante è comprendere che in ogni esperimento di misurazione, quindi scegliendo intervalli temporali differenti, l'andamento temporale sarà sempre diverso, ma si conserverà un qualcosa relativo ad una descrizione globale del fenomeno stesso, come ad esempio, il valore medio nell'intervallo scelto. In letteratura tecnica si potrebbe trovare l'espressione: "i segnali hanno le stesse caratteristiche energetiche". Questo è un esempio di Fenomeno Aleatorio, dipendente dalla variabile temporale, ma se ne possono immaginare molti altri dipendenti da altri tipi di variabili, i. e. variabili spaziali. L'insieme di tutte le realizzazioni producibili dal Fenomeno Aleatorio è riferito, comunemente, con il termine Processo Aleatorio. Un altro possibile processo aleatorio, seguendo gli esempi precedenti, potrebbe, appunto, essere l'arrivo dei clienti in un ufficio postale, o delle automobili ad un semaforo. Tali fenomeni sono modellizzati da particolari processi  (i processi di Poisson) che prendono il nome dal matematico che li ha studiati: Siméon-Denis Poisson (1781-1840).

Tornando all'esperimento della misura della tensione ai capi del Resistore R, ciò che si riscontra è noto con il nome di "rumore termico", e in base a ciò che è stato affermato in precedenza, la natura dei valori di tensione è estremamente casuale. Per il lettore, con una quelche conoscenza sull'argomento, si riporta l'espressione dello spettro bilatero di densità di potenza, che vale:

«math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mi»P«/mi»«mi»n«/mi»«mo»(«/mo»«mi»f«/mi»«mo»)«/mo»«mo»=«/mo»«mn»2«/mn»«mfrac»«mrow»«mi»h«/mi»«mi»f«/mi»«/mrow»«mrow»«msup»«mi»e«/mi»«mfrac»«mrow»«mi»h«/mi»«mi»f«/mi»«/mrow»«mrow»«mi»K«/mi»«mi»T«/mi»«/mrow»«/mfrac»«/msup»«mo»-«/mo»«mn»1«/mn»«/mrow»«/mfrac»«mi»R«/mi»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«munder»«mrow»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mo»§#8773;«/mo»«/mrow»«mrow»«mi»f«/mi»«mo»§lt;«/mo»«mo»§lt;«/mo»«mfrac»«mrow»«mi»K«/mi»«mi»T«/mi»«/mrow»«mi»h«/mi»«/mfrac»«/mrow»«/munder»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mn»2«/mn»«mi»K«/mi»«mi»T«/mi»«mi»R«/mi»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mo»§nbsp;«/mo»«mo»[«/mo»«mfrac»«msup»«mi»V«/mi»«mn»2«/mn»«/msup»«mrow»«mi»H«/mi»«mi»z«/mi»«/mrow»«/mfrac»«mo»]«/mo»«/math»

Con:

  • K=costantte di Boltzman «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mo»(«/mo»«mn»1«/mn»«mo».«/mo»«mn»38«/mn»«mo»§#183;«/mo»«msup»«mn»10«/mn»«mrow»«mo»-«/mo»«mn»23«/mn»«/mrow»«/msup»«mo»§nbsp;«/mo»«mo»[«/mo»«mo»(«/mo»«mi»W«/mi»«mi»a«/mi»«mi»t«/mi»«mi»t«/mi»«mo»/«/mo»«mi»H«/mi»«mi»z«/mi»«mo»)«/mo»«mo»/«/mo»«mo»§#176;«/mo»«mi»K«/mi»«mo»]«/mo»«mo»§nbsp;«/mo»«mo»)«/mo»«/math»;
  • h=costante di Plank «math xmlns=¨http://www.w3.org/1998/Math/MathML¨»«mo»(«/mo»«mn»6«/mn»«mo».«/mo»«mn»62«/mn»«mo»§#183;«/mo»«msup»«mn»10«/mn»«mrow»«mo»-«/mo»«mn»34«/mn»«/mrow»«/msup»«mo»§nbsp;«/mo»«mo»[«/mo»«mi»W«/mi»«mi»a«/mi»«mi»t«/mi»«mi»t«/mi»«mo»§#183;«/mo»«mi»S«/mi»«mi»e«/mi»«mi»c«/mi»«mo»]«/mo»«mo»)«/mo»«/math»;

Considerazioni matematiche, basate sul teorema Centrale del Limite e la modellizzazione aleatoria della tensione in uscita al resistore portano a pensare il segnale misurato al resistore come appartenente ad un Processo Aleatorio Gaussiano e Bianco. Senza addentrarsi nella teoria, l'aggettivo gaussiano deriva dal fatto che il fenomeno aleatorio è modellizzato come un insieme di tante variabili aleatorie, la cui distribuzione di probabilità (Probability Distribution Function) è gaussiana, ovvero data dalla funzione a campana, nota anche come curva di Gauss dal matematico che maggiormente l'ha studiata: Carl Friedrich Gauss (1777-1855).gaussiana

Per inciso si nota che in un sistema di comunicazione a distanza, in cui la trasmissione dell'informazione è affidata a fenomeni elettrici reali, il rumore termico è ineliminabile, a meno che non si porti il sistema nei pressi dello zero assoluto (0 K , –273,15 °C, –459,67 °F). E anche se la sua dimensione in termini di valore assoluto è molto bassa, per tratte molto lunghe quando il livello del segnale informativo diventa molto basso, i due diventano comparabili, introducendo distorsione nell'informazione ricevuta dalla sorgente.

Il termine utilizzato in precedenza "bianco", è preso a prestito dal fenomeno elettromagnetico, secondo il quale una radiazione nello spettro del visibile, con energia distribuita uniformemente su tutte le lunghezze d'onda è percepita come luce bianca.

In ultimo si ricorda al lettore con un minimo di conoscenza sull'argomento, che se si opera la Trasformata di Fourier, sul segnale di correlazione si ottiene lo spettro di densità di potenza, funzione della frequenza. Nel caso tempo-discreto, la Trasformata di Fourier è sostrituita da una sua opportuna trasformazioni algoritmica nel campo discreto, velocizzata, nota come FFT (Fast Fourier Transform). Matlab possiede nel suo core algoritmi ottimizzati per tale computo.

Il codice inserito nell'articolo permette lo studio di un processo gaussiano, con media e deviazione standard modificabili. Il codice produce due grafici: assegnato il valore ai parametri (numero di campioni, valore medio, deviazione standard) il primo grafico rappresenta l'autocorrelazione del segnale, mentre il secondo la densità spettrale di potenza, calcolata come FFT dell'autocorrelazione, (il grafico rappresenta il modulo in quanto generalmente la FFT è a valori complessi).

La funzione matlab

randn(1,N)

produce campioni estratti da una distribuzione gaussiana.


Codice modificabile da inserire nell’editor di Matlab:


%studio di una sequanza pseudo-random a distribuzione pseudo-normale con%media mi (=0) e deviazione standard sigma (=1) -- Versione 1.0 --
%calcolo autocorrelazione
%calcolo DFT tramite FFT
%plotting
clc;
clear;
%scelta numero campioni del processo
N=100000;
 
%parametri statistici del processo
mi=0;
sighma=10;%nota:questa è la deviazione standard
 
%generazione del processo bianco (gaussiano)
x=mi + sighma*randn(1,N);
 
%calcolo autocorrelazione
corrx=xcorr(x,'biased');
 
%normalizzazione numero di campioni
corrxN=corrx;
 
figure(1)
plot(corrxN)
title('AUTOCORRELAZIONE')
xlabel('indice dei campioni')
ylabel('r[n]')
%%%%%%%%%%il valore dell'impulso è pari a sighma^2
 
%calcolo FFT della correlazione=spettro di densità di potenza (PSD)
Rxx=fft(corrxN);
%%%%%%%%%%la FFT è costante e pari a sighma^2
 
figure(2)
plot(abs(Rxx))%modulo, essendo valori complessi
title('modulo PSD')
xlabel('indice campioni')
ylabel('Rxx')
%assunzione di gaussianità---> campioni statisticamente indipendenti
%elaborazione adattativa dei segnali: Cap1, pag 75,76
<p>
</p>
<p>%studio di una sequanza pseudo-random a distribuzione pseudo-normale con
%media mi (=0) e deviazione standard sigma (=1) -- Versione 1.0 --
%calcolo autocorrelazione
%calcolo DFT tramite FFT
%plotting
 
clc;
clear;
%scelta numero campioni del processo
N=100000;
 
%parametri statistici del processo
mi=0;
sighma=10;%nota:questa è la deviazione standard
 
%generazione del processo bianco (gaussiano)
x=mi + sighma*randn(1,N);
 
%calcolo autocorrelazione
corrx=xcorr(x,'biased');
 
%normalizzazione numero di campioni
corrxN=corrx;
 
figure(1)
plot(corrxN) 
title('AUTOCORRELAZIONE')
xlabel('indice dei campioni')
ylabel('r[n]')
%%%%%%%%%%il valore dell'impulso è pari a sighma^2
 
%calcolo FFT della correlazione=spettro di densità di potenza (PSD)
Rxx=fft(corrxN);
%%%%%%%%%%la FFT è costante e pari a sighma^2
 
figure(2)
plot(abs(Rxx))%modulo, essendo valori complessi
title('modulo PSD')
xlabel('indice campioni')
ylabel('Rxx')
%assunzione di gaussianità---> campioni statisticamente indipendenti
%elaborazione adattativa dei segnali: Cap1, pag 75,76
 
</p>




[1] Gaetano Scarano, "Segnali, Processi Aleatori, Stima", Università La Sapienza, Roma, 2009, edizione 1;

[2] Aurelio Uncini, "Elaborazione adattativa dei segnali", ARACNE editrice S.r.l., Roma, 2010


Ultimo aggiornamento Martedì 20 Settembre 2016 21:25

Hanno detto..

You are here: Script Matlab