Marche aléatoire (1j) – PROGRAMME C

Marche aléatoire (1j) – PROGRAMME C

Dans le dernier article, j'ai écrit sur la façon de simuler un lancer/retournement de pièce en utilisant des nombres aléatoires générés dans la plage :.

Nous pouvons utiliser ce code pour simuler un processus stochastique populaire, appelé la marche aléatoire .
REMARQUE : Cela servira également de test pour notre générateur de nombres aléatoires.

Considérons un exemple élémentaire de marche aléatoire unidimensionnelle sur une droite numérique. Le marcheur commence à 0 et peut faire un pas en avant (incrément positif) ou un pas en arrière (incrément -ve), les deux avec des probabilités égales.

Nous savons que pour une pièce non biaisée, la probabilité d'obtenir un pile ou face est égale. J'en ai déjà parlé dans le dernier article. Nous n'utiliserons donc ce code que pour le tirage au sort, qui décidera si notre marcheur aléatoire avance ou recule.

Alors, écrivons un programme qui simule une marche aléatoire et traçons la distance parcourue en fonction du nombre de pas effectués. Ce tracé nous permettra de vérifier/confirmer si notre programme représente vraiment une marche aléatoire ou non.

CODE :

/*********************************
RANDOM WALK 1-d
Plot the path of a 1-d random walker and print out the final displacement
*********************************/
#include<stdio.h>
#include<math.h>
/**Function that generates a random number.
Parameters: 
r0: initial (first) seed
a: scale factor , so that a*r0 give the first random number
m: gives the max. value of random numbers that can be generated (m-1)
c: additional displacement factor
**/
int rand(int r0, int a, int m, int c){
	int r1=(a*r0+c)%m;
	return r1;
}
/**Function that generates random numbers given a seed, and stores them in an array that is passed as an argument.
Parameters:
r0: initial (first) seed
a: scale factor , so that a*r0 give the first random number
m: gives the max. value of random numbers that can be generated (m-1)
c: additional displacement factor
n: no. of random numbers to be generated
x[n]: array that will store the random numbers
**/
void randomNos(int r0, int a, int m, int c, int n, int x[n]){
	int r1=rand(r0,a,m,c);;
	int i;
	for(i=0;i<n;i++){
		x[i]=r1;
		r1=rand(r1,a,m,c);
	}
}
/**Function that results the result of a coin toss:
Parameters:
r: a random number between 0 and 1
Returns 1 for Heads and 0 for tails
**/
int coinTossSingle(double r){
	if(r>0.5){
		return 1;
	} else if(r<0.5){
		return 0;	
	}
}
/**Function that generates n coin tosses results, given a seed and other starting conditions, and stores them in an array that is passed as an argument.
Parameters:
r0: initial (first) seed
a: scale factor , so that a*r0+c give the first random number
m: gives the max. value of random numbers that can be generated (m-1)
c: additional displacement factor
n: no. of coin tosses to be generated
x[n]: array that will store the random numbers
**/
void coinToss(int r0, int a, int m, int c, int n, int results[n]){
	int randNos[n];
	randomNos(r0, a, m, c, n, randNos);
	//Renormalize to 0 to 1
	int i;
	double randNosNew[n];
	for(i=0;i<n;i++){
		randNosNew[i]=(double)randNos[i]/(m-1);
	}
	for(i=0;i<n;i++){
		results[i]=coinTossSingle(randNosNew[i]);	
	}
	
}
main(){
	int a, m, c, r0, n;
	printf("Enter the value of a:\n");
	scanf("%d",&a);
	printf("Enter the value of m:\n");
	scanf("%d",&m);
	printf("Enter the value of c:\n");
	scanf("%d",&c);
	printf("Enter the value of r0(initial):\n");
	scanf("%d",&r0);
	printf("Enter the no. of steps require:\n");
	scanf("%d",&n);
	int tossResults[n];
	coinToss(r0, a, m, c, n, tossResults);
	int i;
	//Step-size
	double h=1;
	//Origin (Start of random walk)
	double x0=0,origin=x0;
	double x1;
	//Array to store the position of the random walker at the ith step
	double x[n];
	for(i=0;i<n;i++){
		if(tossResults[i]==1){
			//Heads=>Move right
			x1=x0+h;
		} else{
			//Tails=>Move left
			x1=x0-h;	
		}
		//Store the position at the ith step in array x[i]
		x[i]=x1;
		x0=x1;
	}
	//Plot the random Walk (Trajectory)
	FILE *fp=NULL;
	fp=fopen("randomWalk1.txt","w");
	for(i=0;i<n;i++){
		fprintf(fp,"%d\t%lf\n",i+1,x[i]);
	}
	double dist=x1-origin;
	printf("\nThe distance travelled is:\n%lf",dist);
}

SORTIE :

Les graphiques ci-dessus semblent un bon exemple pour un marcheur aléatoire d'un point de vue naïf, nous pouvons donc maintenant aller plus loin et travailler sur plus de problèmes sur la marche aléatoire 1-d.

Maintenant, utilisons le programme ci-dessus pour vérifier certaines propriétés communément connues d'un marcheur aléatoire, qui sont,

  1. La valeur d'espérance de la distance parcourue par un marcheur aléatoire 1-d est 0.
    Considérons un marcheur aléatoire qui part de l'origine, et nous le laissons prendre étapes et notez la distance parcourue ( ) de l'origine après étapes. Répétez ce processus fois, et prendre la moyenne des valeurs que vous obtenez. Pour l'infini vous obtiendrez .

  2. La valeur moyenne du carré de la distance parcourue par un marcheur aléatoire 1-d après étapes, est

    ou
    La quantité ci-dessus est appelée la distance quadratique moyenne, et c'est à peu près la distance que nous pouvons nous attendre à ce que notre marcheur aléatoire ait marché après N pas.

Donc, modifions le programme ci-dessus et ajoutons quelques lignes supplémentaires pour effectuer les calculs pour un et .

Ce que je vais faire, c'est exécuter la simulation de marche aléatoire ci-dessus, pour un nombre différent d'étapes de 0 à 1 000 par pas de 100. Pour chaque valeur de , la simulation de marche aléatoire est exécutée fois. Il y a donc la variable M dans le code initialisé comme 100000 , pour exécuter la simulation de marche aléatoire fois.

Pour faire chacun des M exécutions de la simulation différentes les unes des autres, nous aurons besoin d'une graine différente et complètement choisie au hasard pour chaque exécution. Ainsi, à la fin de chacune des M exécutions, je génère un nouveau nombre aléatoire à partir de la dernière graine de l'itération précédente. Ensuite j'ai les tableaux d et d2 qui stockera la valeur de et pour chacun des M court. Enfin je viens de calculer les moyennes des valeurs stockées dans d et d2 , et les stocke dans un .txt fichier avec la valeur de n. Pour que nous puissions tracer et vs. .

CODE :

/*********************************
RANDOM WALK 1-d
Plot <d(N)> and <d^2(N)> vs N
*********************************/
#include<stdio.h>
#include<math.h>
/**Function that generates a random number.
Parameters: 
r0: initial (first) seed
a: scale factor , so that a*r0 give the first random number
m: gives the max. value of random numbers that can be generated (m-1)
c: additional displacement factor
**/
int rand(int r0, int a, int m, int c){
	int r1=(a*r0+c)%m;
	return r1;
}
/**Function that generates random numbers given a seed, and stores them in an array that is passed as an argument.
Parameters:
r0: initial (first) seed
a: scale factor , so that a*r0 give the first random number
m: gives the max. value of random numbers that can be generated (m-1)
c: additional displacement factor
n: no. of random numbers to be generated
x[n]: array that will store the random numbers
**/
void randomNos(int r0, int a, int m, int c, int n, int x[n]){
	int r1=rand(r0,a,m,c);;
	int i;
	for(i=0;i<n;i++){
		x[i]=r1;
		r1=rand(r1,a,m,c);
	}
}
/**Function that results the result of a coin toss:
Parameters:
r: a random number between 0 and 1
Returns 1 for Heads and 0 for tails
**/
int coinTossSingle(double r){
	if(r>0.5){
		return 1;
	} else{
		return 0;	
	}
}
/**Function that generates n coin tosses results, given a seed and other starting conditions, and stores them in an array that is passed as an argument.
Parameters:
r0: initial (first) seed
a: scale factor , so that a*r0+c give the first random number
m: gives the max. value of random numbers that can be generated (m-1)
c: additional displacement factor
n: no. of coin tosses to be generated
x[n]: array that will store the random numbers
**/
void coinToss(int r0, int a, int m, int c, int n, int results[n]){
	int randNos[n];
	randomNos(r0, a, m, c, n, randNos);
	//Renormalize to 0 to 1
	int i;
	double randNosNew[n];
	for(i=0;i<n;i++){
		randNosNew[i]=(double)randNos[i]/(m-1);
	}
	for(i=0;i<n;i++){
		results[i]=coinTossSingle(randNosNew[i]);	
	}
	
}
main(){
	int a=1093, m=86436, c=18257, n, r0=43, M=100000, stepCount=0, N=1000;
	//int m=121500, a=1021,c=25673, n, r0=51,M=100000, stepCount=0, N=1000;
	//int m=259200, a=421, c=54773, n, r0=12, M=100000, stepCount=0, N=1000;
	//int m=121500, a=2041, c=25673, n, r0=25, M=100000, stepCount=0, N=1000;
	
	/*printf("Enter the value of a:\n");
	scanf("%d",&a);
	printf("Enter the value of m:\n");
	scanf("%d",&m);
	printf("Enter the value of c:\n");
	scanf("%d",&c);
	printf("Enter the value of r0(initial):\n");
	scanf("%d",&r);*/
	
	
	FILE *fp="NULL";
	fp=fopen("randomWalk4.txt","w");
	double d[M];
	double d2[M];
	//Run the random-walk simulation for n steps
	for(n=0;n<=N;n=n+100){
		printf("%d\n",stepCount); //To keep trak of where we are in the execution
		stepCount++; //To keep trak of where we are in the execution
		int j;
		//Run the same simulation M times
		for(j=0;j<M;j++){
			int tossResults[n];
			//use the coin toss/flip result to define forward or backward movement 
			coinToss(r0, a, m, c, n, tossResults);
			int i;
			double h=1;
			double x0=0,origin=0;
			double x1;
			int count[2];
			count[0]=0;
			count[1]=0;
			for(i=0;i<n;i++){
				if(tossResults[i]==1){
					//x1=x0+h;
					count[0]++;
				} else{
				//	x1=x0-h;	
				count[1]++;
				}
				//x0=x1;
			}
			//find the distance from origin
			//d[j]=x1-origin;
			d[j]=count[0]-count[1];
			//square of the distance
			d2[j]=pow(d[j],2);
			//generate a new seed at each of the M runs
			r0=rand(r0,a,m,c);	
			
		}
		//find out the averages of the d and d^2 after M runs
		double sum1=0,sum2=0;
		for(j=0;j<M;j++){
			sum1=sum1+d[j];
			sum2=sum2+d2[j];
		}
		double dav=sum1/M;	// <d>
		double dav2=sum2/M; // <d^2>
		//store the value of n, <d> and <d^2> in .txt file for each n
		fprintf(fp,"%d\t%lf\t%lf\n",n,dav,dav2);
	}
}

SORTIE :

—>Pour a=1093, m=86436, c=18257,

lorsque nous traçons les données du fichier texte généré après l'exécution, nous obtenons les tracés suivants.
Pour N=1000, M=100 000

Mais, pour N=5000, M=100 000

Ici, quelque chose d'intéressant se passe.
La valeur de pour N(le nombre de pas effectués)>1000 nous n'obtenons pas les résultats attendus. Cela implique que notre générateur de nombres aléatoires n'est pas idéal. Il existe des corrélations.

N'oubliez pas que nous utilisons l'algorithme suivant pour générer des nombres aléatoires,

appelé générateur congruentiel linéaire
Cet algorithme génère un maximum de nombres aléatoires avec la valeur maximale de (Essayez de voir pourquoi).
Ici, est la graine.
Les valeurs de , et sont des valeurs soigneusement choisies.

Ici, nous avons une option pour changer les valeurs de a, m et c.
Les mathématiciens ont testé de nombreuses valeurs pour celles-ci et ici, j'écris quelques-unes de celles tirées de Numerical Recipes in C.
1. m=86436, a=1093, c=18257
2. m=121500, a=1021,c=25673
3. m=259200, a=421, c=54773
4. m=121500, a=2041, c=25673

Exécutons à nouveau le code ci-dessus pour différentes valeurs de a, m et c et voyons les résultats.

—> Pour m=121500, a=1021, c=25673, on obtient

Pour N=2000, M=100 000

Pour N=10 000, M=10 000

—>Pour m=259200, a=421, c=54773, on obtient

N=5000, M=100 000

N=10 000, M=10 000

—>Pour m=121500, a=2041, c=25673, on obtient

Pour , M(essais)=100000

Pour , M(essais)=10000

Nous voyons que pour différentes valeurs, nous obtenons la valeur attendue de pour un plus grand N.

Mais nous n'obtenons toujours pas le comportement attendu pour N supérieur à 2000.

Alors, peut-être avons-nous besoin d'un meilleur générateur de nombres aléatoires. Ou peut-être pouvons-nous essayer d'utiliser une graine encore plus aléatoire en utilisant l'horloge système.

Essayons ça.

Références :

http://www.mit.edu/~kardar/teaching/projects/chemotaxis(AndreaSchmidt)/random.htm

http://mathworld.wolfram.com/RandomWalk1-Dimensional.html

https://en.wikipedia.org/wiki/Random_walk