BRNG: chi-quadro, ent, i timer e gli interrupt

Nel precedente post [1] abbiamo affrontato i concetti di base del generatore, lasciando aperto l’interrogativo del pessimo chi-quadro. Prima di addentrarci nelle ragioni del cattivo risultato, è necessario comprendere che cosa è questo valore e come la casualità viene “estratta” dall’evento di una radiazione. 

Partiamo dal capire cos’è e come funziona il chi-quadro (indicato anche con \chi^2). Si tratta di un valore utilizzato in statistica per verificare l’aderenza di un set di valori a una distribuzione teoricamente prevista. Affrontiamo prima l’utilizzo tradizionale in statistica, poi volgeremo lo sguardo all’uso del \chi^2 per questa applicazione.

Dato un set di dati, chiamiamo frequenza il numero di volte che si presenta un determinato dato, e gradi di libertà il numero di possibili valori meno uno. Perché meno 1? Quest’operazione, che a prima vista può sembrare abbastanza strana, è in realtà molto logica. Supponiamo di volere analizzare il lancio di una moneta. Abbiamo due possibili outcome: testa e croce. Se ci pensiamo, però, la percentuale di volte che esce testa, è direttamente determinato dalla percentuale di volte che esce croce. Se consideriamo un evento con tre possibili outcome, la percentuale del terzo outcome è direttamente determinata dagli altri due. In questo senso, parliamo di gradi di libertà. 

La formula del chi-quadro è:

    \[ \chi^2 = \sum_{i = 1}^{k} \frac{{(o_i - a_i)}^2}{a_i} \]

dove o_i è la frequenza osservata e a_i è la frequenza attesa teoricamente.

Prendiamo come esempio il classico lancio di uno dado a sei facce. Il lancio del dado può risultare in sei possibili risultati, questo ci da 5 gradi di libertà. Immaginiamo di lanciare il dado 1000 volte, vogliamo verificare quella che in statistica prende il nome di ipotesi zero, ovvero verificare che, entro una certa probabilità, i nostri risultati siano veramente casuali. 
Questi sono i dati che ottieniamo dal nostro esperimento:

Sommando i singoli contributi, otteniamo il valore di chi-quadro del nostro esperimento: 3.068. E ora? Prendiamoci un secondo per riflettere sui dati e sulla nostra formula, tenendo a mente che stiamo cercando di misurare, in questo caso, quanto i nostri dati aderiscono ad una distribuzione uniforme. Se tutto fosse perfetto (in un ipotetica dimensione in cui questo può succedere), avremmo le frequenze osservate o_i uguali alle frequenze teoriche a_i, per cui il nostro denominatore {(o_i - a_i)}^2 \longrightarrow 0. Di conseguenza \chi^2 \longrightarrow 0. Sappiamo però che la perfezione non appartiene al nostro mondo, altrimenti non avremmo nemmeno bisogno della statistica. In altre parole, è estremamente improbabile che i nostri dati rispecchino esattamente la distribuzione teorica, un valore di chi quadro troppo vicino a zero ci deve far insospettire. 

Dall’altro lato, più ci allontaniamo dalla distribuzione teorica, più il numeratore cresce, a parità di denominatore. Il che porta il valore di chi-quadro a crescere. Questo per l’uso normale del chi-quadro è un’ottima notizia, significa poter rifiutare l’ipotesi nulla, e quindi saper di stare lavorando con dati che non sono solo frutto del caso, ma che hanno un qualche significato. Per la nostra applicazione questa è invece una pessima notizia: vuol dire che i nostri dati non sono distribuiti in maniera uniforme.

Questo vuol dire che, per la nostra applicazione, dovremo cercare una via di mezzo. Ok, ma numericamente? Il valore accettabile di chi-quadro varia in base ai gradi di libertà del sistema. Qui inizia la parte difficile, fatta di integrali poco amichevoli e di equazioni non risolvibili algebricamente [2]. Fortunatamente qualcuno ha già fatto il lavoro per noi, disponendo i risultati in tabelle comodamente reperibili su internet.

Tavola distribuzione chi quadro

Le righe della tabella rappresentano i gradi di libertà del sistema, nel nostro caso del dado sono 5. Le colonne rappresentano i livelli di probabilità che il valore calcolato sia maggiore del valore tabulato. Esistono anche tabelle che indicano la probabilità che il valore calcolato sia minore, queste vengono chiamate left tail tables, mentre quella mostrata qua sopra è una right tail table. Questo perché in un caso viene considerata la parte destra del grafico, e nell’altro la parte sinistra. Nel nostro caso abbiamo \chi^2 = 3.068, che si colloca tra il 90% e il 25% dei casi. Questo è sufficiente per dire che non ci sono variazioni eccessive da un comportamento che possiamo classificare come random. [3]

 

Right tail del chi-quadro
Right tail del chi-quadro per d = 5 e chi-quadro = 3,068, da statdistributions.com

Tornando alle banane, prendiamo come percentuali di riferimento il 90% e il 10%. Purtroppo non sono riuscito a trovare tabelle per 255 gradi di libertà, fortunatamente però l’autore di ent mette a disposizione uno strumento per calcolare i valori per ogni probabilità [2].
Per 255 gradi di libertà (2^8 = 256\ possibili\ byte\ - 1) abbiamo i valori: P_{\chi^2,\ 90\%} = \textbf{226.5}P_{\chi^2,\ 10\%} =\textbf{284.3}. Dai test di ent sui valori registrati dal generatore, avevamo ottenuto un valore di 498.15, decisamente fuori dal range dell’accettabilità. La percentuale di probabilità restituita da ent era < 0.01\%.

Ent, tramite l’opzione -c ci restituisce la frequenza di ogni byte, che è riportata qui di seguito, omessi i byte poco significativi:

valerio@valerio ~/tests $ ent logbinario.txt -c
Value Char Occurrences Fraction
  0              739   0.004077
  1              402   0.002218
  2             1033   0.005699
  3              674   0.003718
...              ...      ...       
254   �         722   0.003983
255   �         691   0.003812

Total:        181256   1.000000

Entropy = 7.997995 bits per byte.

Optimum compression would reduce the size
of this 181256 byte file by 0 percent.

Chi square distribution for 181256 samples is 498.15, and randomly
would exceed this value less than 0.01 percent of the times.

Arithmetic mean value of data bytes is 127.4942 (127.5 = random).
Monte Carlo value for Pi is 3.138799695 (error 0.09 percent).
Serial correlation coefficient is 0.005408 (totally uncorrelated = 0.0).

Si nota subito che il byte 1 ha significativamente meno count degli altri, e che il byte 2 ne ha molti di più. A uno sguardo attento, i count che sembrano “mancare” a 1, sono assegnati a 2. Questo è sicuramente il nostro problema.

Dopo un po’ di prove infruttuose, ho deciso di dividere i byte di posizione pari dai byte di posizione dispari. Questo perché per ogni numero a 16 bit (2 byte) generato, vengono creati due byte, uno di posizione pari e uno di posizione dispari. Ho praticamente diviso il file in Most Significant Byte e Least Significant Byte. Questi sono i risultati:

valerio@valerio ~/tests $ ent logbinariomsb.txt -c
Value Char Occurrences Fraction
  0              358   0.003950
  1              393   0.004336
  2              347   0.003829
  3              333   0.003674
...              ...     ...
254   �          374   0.004127
255   �          330   0.003641

Total:         90628   1.000000
...
Chi square distribution for 90628 samples is 240.88, and randomly would exceed this value 72.82 percent of the times.
...
valerio@valerio ~/tests $ ent logbinariolsb.txt -c
Value Char Occurrences Fraction
  0              381   0.004204
  1                9   0.000099
  2              686   0.007569
  3              341   0.003763
...              ...     ...
254   �         348   0.003840
255   �         361   0.003983

Total:         90628   1.000000
...
Chi square distribution for 90628 samples is 891.23, and randomly would exceed this value less than 0.01 percent of the times.
...

Il gruppo dei MSB non riporta nessun problema significativo, mentre il gruppo dei LSB è dove risiede il problema.

Per capire da dove si origina questo problema, dobbiamo prima di tutto capire come i numeri sono generati internamente. Il tubo geiger, tramite un circuito di interfaccia, manda un segnale sul pin 2 (PB2/INT0) di arduino quando viene colpito da una radiazione. Il pin 2 è configurato in modo da generare un interrupt quando riceve un rising edge: attachInterrupt(digitalPinToInterrupt(2), randomCore, RISING);. L’interrupt chiamerà la funzione randomCore(), che è così definita:

void randomCore() {
  Serial.println((unsigned int)(( micros() >> 2 ) & 0xFFFF));
  return;
}

La funzione, quando chiamata, chiama a sua volta la funzione micros() di arduino. Questa funzione restituisce un numero a 32 bit che rappresenta il numero di microsecondi trascorsi dall’accensione del sistema. Essendo un numero unsigned a 32 bit, andrà in overflow dopo 4294.96 secondi, ovvero ogni circa 70 minuti. Essendo la velocità del microcontrollore insufficiente per ottenere un aggiornamento più preciso, la micros() si aggiorna a salti di 4 microsecondi, mantenendo sempre i due least significant bit a zero.

Per questo, eseguiamo uno shift del valore restituito dalla micros()di due bit a destra. In questo modo otteniamo un valore di 30 bit. Se usassimo anche i bit più signficativi, otterremmo dei numeri progressivi fino al successivo overflow del timer. Ogni numero sarebbe sicuramente maggiore del precedente e sicuramente minore del successivo, nell’arco dei 70 minuti necessari per l’overflow. Questo non è sicuramente random. 

Conserviamo quindi solo i primi 16 byte del valore della micros(). Questo valore avrà overflow ogni 262144 microsecondi, rendendo estremamente improbabile il verificarsi della situazione di cui sopra. (Il contatore geiger registra di fondo circa 30 pulsazioni al minuto).

Il valore problematico è quindi il seguente, dove il MSB è irrilevante:

16       8         0
XXXX XXXX 0000 0001 

Ha catturato la mia attenzione il fatto che questo valore si presenta ogni 4*2^8 = 1024 microsecondi, cioè circa 1 millisecondo, ed è il valore successivo a quello che genera un overflow interrupt. La mia attenzione si è quindi spostata sul codice della funzione millis() del core di arduino [4] [5].

La funzione millis() del core funziona impostando il prescaler del TIMER0 a 64, che per un timer a 8 bit con un clock di 16MHz comporta un overflow del timer ogni 1.024 millisecondi. L’overflow del timer genera un interrupt, con vettore TIMER0_OVF. Se l’impulso del tubo geiger arriva contemporaneamente all’overflow del TIMER0, avremo due interrupt concorrenti: TIMER0_OVF e INT0Questa situazione è gestita dal microcontrollore con un sistema di priorità degli interrupt, il cui ordine di priorità è indicato nel datasheet:

Priorità degli interrupt, dal datasheet dell'ATMega328/P
Priorità degli interrupt, dal datasheet dell’ATMega328/P [6]
L’overflow del TIMER0 ha priorità molto inferiore a quella dell’interrupt esterno, la mia ipotesi è quindi la seguente:

  1. Caso 1: l’interrupt INT0 arriva contemporaneamente all’interrupt TIMER0_OVF. Avendo priorità maggiore viene eseguito prima l’interrupt esterno a scapito della millis(), intaccando la precisione della funzione ma senza produrre effetti visibili sui numeri generati;
  2. Caso 2: l’interrupt INT0 arriva al ciclo di clock successivo rispetto all’interrupt TIMER0_OVF. Essendo già trascorso un ciclo di clock, l’interrupt TIMER0_OVF sarà già in esecuzione. Quando l’esecuzione sarà terminata, micros() sarà già al valore 2, per cui il numero generato sarà registrato con valore 2. 

È possibile che anche l’uso della seriale abbia un’influenza su questo ritardo, ma non ho approfondito. La soluzione a questo problema sarà esposta nel prossimo post, ed implica l’uso del TIMER1 per la generazione dei numeri, invece della funzione micros() di Arduino.

In chiusura, voglio ringraziare Mauro Mombelli e la community di StackExchange cryptography per il prezioso aiuto [7] [8].

Note, fonti e link:
[1] BRNG: generare numeri casuali a partire dalle banane
[2] Fourmilab – chi-squared calculator
[3] How to evaluate chi-squared results
[4] Millis in wiring.c su github
[5] µC eXperiment blog – Examination of the Arduino millis() Function
[6] Datasheet dell’ATMega328/P
[7] StackExchange Cryptography – Random number generator issue
[8] StackExchange Cryptography – How to evaluate chi squared result?
[9] AVRBeginners – Timers

CC BY-NC-SA 4.0 This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.

Lascia un commento

Il tuo indirizzo email non sarà pubblicato. I campi obbligatori sono contrassegnati *