Newsgroup: it.comp.lang.c
Titolo: WRGN thanx

Salve sto cercando di generare una sequenza di numeri casuale con 
probabilità gaussiana, come posso fare ?
rand() genera una sequenza di pseudo noise ovvero attraverso qualche 
polinomio generatore = con distribuzione uniforme.
E' chiaro che una gaussiana non è plausibile perchè oltre un certo valore 
non si può andare visto che i bit di ogni tipo di dato sono finiti, ma se 
volessi tipo una gaussiana troncata?
____________________
Se la so, rispondo sempre. 

 volessi tipo una gaussiana troncata?
Il metodo classico e' questo:
generi due numeri casuali tra 0 e 1, per esempio con
double r1,r2;
r1 = rand()/(1.0+RAND_MAX);
r2 = rand()/(1.0+RAND_MAX);
anzi, siccome come si vede dopo i numeri non devono
essere uguali a 0 (r1, altrimenti si ha errore numerico),
puoi definirti una funzioncina cosi'
double r01()
{
   double r;
   do r=rand()/(1.0+RAND_MAX); while (r==0.0);
   return r;
}
alternativamente, penso che senza introdurre errori rilevanti,
tu possa scrivere una cosa di questo tipo:
r1 = (1.0+rand())/(2.0+RAND_MAX));
Comunque tu faccia, devi estrarre due numeri casuali r1 e r2
in modo uniforme poi applichi questa trasformazione che ti
permette di generare due numeri distribuiti in modo normale
con media 0 e varianza 1:
x1 = sqrt(-2.0 * log(r1)) * cos(2.0*pi*r2);
x2 = sqrt(-2.0 * log(r1)) * sin(2.0*pi*r2);
dove log sqrt sin cos sono definite in math.h
e pi e' pigreco = 3.141592653589793238462643383 ...
ovviamente se ti serve un solo valore puoi scartare x2.
Se ti servono con altra media e varianza basta fare
qualcosa tipo y1 = M+sqrt(V)*x1 etc.etc.
Esistono anche altri modi equivalenti che a prezzo di
generare piu' numeri casuali possono fare a meno del
calcolo di sin e cos. (Potrebbero aver senso s
su certe architetture). Trovi queste cose per esempio
http://www.dspguru.com/howto/tech/wgn2.htm
http://www.taygeta.com/random/gaussian.html
ciao,
g.

Grazie mille 

> x1 = sqrt(-2.0 * log(r1)) * cos(2.0*pi*r2);
Ho solo una piccola domanda in proposito: dovendo fare questa operazione 
36*50*160 volte mi chiedevo se non vi fosse un metodo con complessità 
computazionale inferiore;
Io ho pensato (solo questa mattina) che un'esponenziale potesse ben 
approssimare (come fa comunque anche la funzione in alto essendo a valori 
finiti in uscita x forza) il mio caso, ovvero: supponiamo che il mio x1 sia 
uno signed short ovvero con ampiezza massima pari a 2^15 allora prendo un 
esponenziale  e^kx e lo tronco a 2^15 e scalo in modo che e^(k*1) = 2^15;
In questo modo dovrei ottenere qualcosa che si avvicina molto ad una 
gaussiana giusto?..forse posto nel posto sbagliato...

 computazionale inferiore;
Prima di tutto (prima ancora di leggere quello che hai scritto dopo),
bisognerebbe valutare se hai bisogno di romperti ulteriormente la testa.
Su un PC tradizionale calcolare 36*50*160  valori e' una inezia.
Per esempio, il seguente programmino che calcola questi 36*50*160 valori
e ne fa la media, sul mio ormai obsoleto Athlon 1400 richiede 
complessivamente circa 4 centesimi di secondo....
Chiaramente la cosa cambia se stai programmando un sistema embedded o 
qualcosa di simile in cui magari non hai il coprocessore matematico e
neanche i numeri in floating point..
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define Pi2 6.283185307179586476925286766559
#define QUANTI (36*50*160)
inline void coppiagauss(double *a, double *b)
{
   double t,r1,r2;
   r1=(1.0+rand())/(2.0+RAND_MAX);
   r2=(1.0+rand())/(2.0+RAND_MAX);
   t = sqrt(-2.0*log(r1));
   *a = t*cos(Pi2*r2);
   *b = t*sin(Pi2*r2);
}
int main()
{
   double a,b,acc=0.0;
   int i;
   for (i=0; i<(QUANTI/2); i+=2)
     {
       coppiagauss(&a,&b);
       acc+=(a+b);
     }
   /* se QUANTI e' dispari me ne manca 1 */
   if (QUANTI & 1)
     {
       coppiagauss(&a,&b);
       acc+=a;
     }
   printf("Media = %.7f\n",acc/QUANTI);
   return 0;
}

Giovanni Resta wrote:
Ah, dimenticavo.
Volendone generare di diversi ad ogni run si puo' aggiungere
un bel
#include<time.h>
all'inizio
e un bel
srand(time(0));
prima del ciclo for principale all'interno del main,
ciao,
g.

> srand(time(0));
si, questo l'avevo già messo


Un numero estratto da una distribuzione normale sara' un numero
compreso tra -infinito e +infinito.
Chiaramente un numero di questo tipo non ti serve e non fitta
soprattutto come mi sembra di capire nei tuoi 15 bit + segno.
Prima di tutto devi stabilire un limite superiore oltre al quale
non generi i numeri.
Per esempio puoi stabilire di limitare l'intervallo a -4 ... +4
che nel caso di gaussiane con media 0 e sigma 1 lascia fuori
una zona di probabilita' 0.0000633 circa.
stabilito l'intervallo puoi immaginare di scalare i valori
compresi tra -4 e 4 moltiplicandoli per 2^15/4 = 2^13=8192.
quindi se x e' un numero double generato tra -oo e +oo
potrai fare una cosa del tipo
se fabs(x)>=4.0 lo butti, altrimenti
short int s = (short int)(8192.0*x)
Quanto al fatto di generarli direttamente in quel range
senza usare operazioni matematica, come dicevo nell'altro
post, bisognerebbe prima vedere se nei hai proprio bisogno.
La cosa non e' tanto semplice.
ciao,
g.

> Su un PC tradizionale calcolare 36*50*160  valori e' una inezia.
Dunque forse ho fatto male i conti...spet..
#define pl 36
int pk = 27;
int v = (pl + 2)*pk; //=1026
while(v--){
//codice veloce
for ( c=0; c<160;c++){ //1026*160 volte
    wgn = 4*var*gauss();
    *(vett+c)=*(vett+c) + wgn;
}
//codice veloce
}
Bhe nn so perchè ma se decommento il for il tempo per realizzare il processo 
si allunga a qualcosa di tediosissimo
NB: non si impalla ne niente.. solo che ci mette tanto..
ho notato che gauss() contiene un do-while... mi domando se sia per 
quello... ma nn credo.. perche dovrebbe essere una condizione che viene 
rispettata al max un paio di volte... nn so..


Non sapendo cosa hai messo dentro la tua funzione gauss
e cosa sono tutte le varie variabili che usi,
non posso risponderti.
g.

 short int s = (short int)(8192.0*x)
Questo è chiaro, ma per modellare in numeri generati tra 0 e 4 con 
distribuzione gaussiana pensavo di rifarmi ad un'esponenziale crescente:
in una N[0,1] ho che la zona più probabile è quella nell'origine; ora prendo 
la funzione exp(x)-1 nell'intervallo x tra 0 e 1 per esempio e vado a vedere 
generando numeri uniformemente cosa mi sputa fuori: ebbene vedo che per come 
è fatta la pendenza della curva esponenziale ho buone probabilità di 
trovarmi vicino all'origine e scarse di trovarmi distante; questo perchè 
essendo perchè nell'intervallo [0,9;1]  (con prob. del 10%) mi trovo molto 
distante [5,04; 6,38] ed ho grande variazione (6,38-5,04) = 1,34 mentre 
prendendo l'intervallo[0; 0,1] mi ritrovo tra [0; 0,22] cioè assumo una 
gamma di valori molto + ristretta e quindi con maggior probabilità.
Ora provo in matlab 

 distribuzione gaussiana pensavo di rifarmi ad un'esponenziale crescente:
mi sembra di capire che sul tuo sistema tu possa usare l'esponenziale.
Quindi hai accesso alle varie funzioni matematiche.
Allora, scusa la curiosita', perche' vuoi complicarti la vita
con approssimazioni quando puoi calcolare i valori esatti ?
Tanto per curiosita', sulla tua macchina il programma completo
(che contiene la funzione coppiagauss()) e che ho postato
qualche messaggio fa, quanto tempo ci mette a girare ?
Come ti ho scritto sulla mia macchina piuttosto vecchia ci mette
circa 4 centesimi di secondo.
A meno che tu non stia facendo cose molto particolari o abbia 
disposizione hardware 10 volte piu' lento, 4 centesimi di secondo
mi sembrano un tempo piu' che dignitoso.
Se nell'altro messaggio con gauss() intendevi la versione
tipo questa, presa e modificata da uno dei link che ti ho segnalato:
void duegauss(double *a, double *b)
{
double x1, x2, w;
          do {
	   x1=(1.0+rand())/(2.0+RAND_MAX);
	   x2=(1.0+rand())/(2.0+RAND_MAX);
	   x1 = 2.0 * x1 - 1.0;
	   x2 = 2.0 * x2 - 1.0;
	   w = x1 * x1 + x2 * x2;
          } while ( w >= 1.0 );
          w = sqrt( (-2.0 * log( w ) ) / w );
          *a = x1 * w;
          *b = x2 * w;
}
che evita di usare sin e cos a prezzo di estrarre qualche
numero casuale in piu',
posso dirti, avendo provato a sostituirla all'altra,
che il risparmio in termini di tempo e' irrisorio
(almeno sul mio sistema).
Eseguendo 100 volte il numero di estrazioni che dicevi tu,
in un caso ci ha messo 4.19 secondi e nell'altro 3.83.
Quindi non e' significativamente piu' veloce.
Le cose cambierebbero se sulla tua CPU le operazioni
trigonometriche fossero molto piu' lente.
ciao,
g.

 in un caso ci ha messo 4.19 secondi e nell'altro 3.83.
Beato te.. a me sta talmente tanto tempo che mi scoccio e faccio ripartire..

 e pi e' pigreco = 3.141592653589793238462643383 ...
M_PI, definito in <math.h> ;)
Luca
-- 
Home: 
http://kronoz.cjb.net
Soltanto un istante...
tuffarsi nell'infinito...
annegare senza pensieri