PDA

Visualizza la versione completa : [C++] Generazione numeri casuali GSL


Follyer
18-12-2007, 15:44
Ciao a tutti!!!

Dato che oramai sono stato abbandonato dal mio prof che non mi risponde alle email e addirittura non fa neanche il ricevimento mi affido al forum.

Devo fare un esercizio dove devo creare una generazione di numeri casuali che seguano una distribuzione chi-quadrato.

Chiedo a chiunque sappiamo utilizzare la libreria GSL di provare questo programma.
Se al posto di chisq metto tdist mi genera i numeri casuali e fa il test chi-quadrato se metto chisq per creare la sequenza l'output errato...

Grazie a chiunque riesca ad aiutarmi!




#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <float.h>
#include <math.h>
#include <limits.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_rng.h> // Per l'utilizzo dei numeri casuali
#include <gsl/gsl_cdf.h> // Per uso delle distribuzioni cumulative
#include <gsl\gsl_math.h>
#include <gsl/gsl_randist.h> // per generatori di distribuzioni non uniformi
#define CAT 10
#include <gsl/gsl_qrng.h> //numero categorie

char *orario()
{
time_t t;
t=time(&t);
return asctime(localtime(&t));
}

void intestazioneFile(FILE *fp)

{
fprintf(fp,"\n * Decima esercitazione * \n\n");
/*nome dell'algoritmo*/
fprintf(fp, " Nome algoritmo: Esercizio10.c \n\n");
/*autore*/
fprintf(fp, " Variabile casuale scelta: Chi-Quadrato \n\n");
/*v.c*/
fprintf(fp," Casella Alberto 055941 \n\n");
/*data*/
fprintf(fp, " Data: %s\n\n\n\n",orario());




}
void calcolaFreqAss(double num,double limInt[],long dimVt,long VI[])
{
long j;
for(j=0;j<dimVt;j++)
if (num<limInt[j])
{VI[j]++;
return;
}
}
int main(void)
{
FILE *fu;
int i, j, N_Iter, k; int s;
double v; int n=20;
double *Vettore_Frequenze_Teoriche, limitiInt[CAT];
double Vettore_Frequenze_Assolute[10000];
double Xquadrato, xTeoriche;
long vFrA[CAT];
gsl_rng *rg;

if((fu=fopen("Esercizio10.txt","w")) == NULL) //creazione del file contenente l'output
return -1;
intestazioneFile(fu);

printf ("Inserisci il numero di iterazioni: ");
scanf ("%d", &N_Iter);



printf ("Inserisci gradi di liberta': ");
scanf ("%lf", &v);


fprintf(fu,"Generazione di %d numeri casuali estratti da una distribuzione Chi-Quadrato con i seguenti parametri:\n", N_Iter);


fprintf(fu,"Numero gradi di liberta': %lf\n\n", v);

rg = gsl_rng_alloc(gsl_rng_taus);
gsl_rng_set(rg,time(NULL));



Vettore_Frequenze_Teoriche = (double*) calloc (n+1, sizeof (double));
if(Vettore_Frequenze_Teoriche == NULL)
{ printf("Non possibile allocare memoria!\n");
exit(-1);
}



for(j=0;j<N_Iter;j++) { s=(int)gsl_ran_chisq (rg, v); //genera nc e aggiorna freq ass
Vettore_Frequenze_Assolute[s]++;}

for (i=0;i<=n;i++)
Vettore_Frequenze_Teoriche[i]=gsl_ran_chisq_pdf (i, v)*N_Iter;

for(i=0;i<=n;i++)
{
fprintf(fu,"x=%d\tp(%d)=%lf\tF(%d)=%lf\n",i,i,gsl_ran_chisq_pdf(i,v),
i,gsl_cdf_chisq_P(i,v));
}
fprintf(fu,"\n");

for(i=0;i<CAT;i++) //inizializza vettore limiti intervalli
limitiInt[i]= n*i/(double)CAT ;

for(j=0;j<CAT;j++) //azzera frequenze assolute
{vFrA[j]=0;
}

for(j=0;j<N_Iter;j++) //genera nc e aggiorna freq ass
{
calcolaFreqAss(gsl_ran_chisq (rg, v),limitiInt,CAT,vFrA);
}

fprintf(fu,"\nFrequenze osservate:\n\n");
for(k=0;k<CAT;k++)
fprintf(fu,"%d^ classe: %d\n",k+1,vFrA[k]);

//Verifica H0 su X^2

double risultato=0.00;

for (i=0;i<=n;i++)

risultato = risultato + (((Vettore_Frequenze_Assolute[i]-Vettore_Frequenze_Teoriche[i])*(Vettore_Frequenze_Assolute[i]-Vettore_Frequenze_Teoriche[i]))/Vettore_Frequenze_Teoriche[i]);

Xquadrato = risultato;
xTeoriche = gsl_cdf_chisq_Pinv (0.97, (double)n-1);

if (Xquadrato <= xTeoriche)
fprintf(fu, "\nPoiche' %4.4f <= %f: test superato\n", Xquadrato, xTeoriche);
else
fprintf(fu, "\nPoiche' %4.4f > %f: test non superato\n", Xquadrato, xTeoriche);

fclose(fu);
gsl_rng_free(rg);
system("pause");
return 0;
}

alka
22-12-2007, 15:50
Il linguaggio va indicato anche nel titolo, come da Regolamento (http://forum.html.it/forum/showthread.php?s=&threadid=973887).

Qui l'ho aggiunto io. Mi dirai eventualmente se va modificato.

Ciao! :ciauz:

Loading