Annonce

Bienvenue sur le site support de mes ouvrages d'introduction à SAS

La 4ème édition de mon ouvrage est disponible depuis le 11 avril 2019 !

Où trouver cet ouvrage ?


#1 20-06-2019 06:19:22

SAS-SR
Administrateur
Lieu: Université d'Orléans
Date d'inscription: 01-09-2008
Site web

Liliane, fais les valises, on rentre à Paris

ou le modèle de ségrégation de Schelling

Le titre de ce sujet des beaux mercredis a peu de chance de faire rire les moins de 50 ans....
https://www.youtube.com/watch?v=x_NihRC7W0s

bref...

Mon collègue de bureau a récemment attiré mon attention sur cette page du site de Thomas J. Sargent et John Stachurski "Quantitative Economics with Python".

https://lectures.quantecon.org/py/schelling.html

et là vous vous dites "quoi ! on fait du python sur www.sas-sr.com ? tout fout le camp..."

mais non, je vous rassure tout de suite mais le sujet que Sargent et Stachurski traite est tout à fait intéressant !

Schelling (prix nobel 1985) a développé un modèle très simple qui explique la formation des ghettos. Imaginons une ville dans laquelle vous avez deux types de population : les rouges et les bleus. Une personne est heureuse si, parmi ses 10 voisins les plus proches, il y a au moins 5 personnes du même type qu'elle.

S'il n'y a que 4 personnes ou moins qui lui ressemble, la personne va déménager et au bout de quelques périodes, on constate que des ghettos rouge et bleu se forment.

L'objectif de Sargent et Stachurski est de proposer un programme en python qui montre la formation de ces ghettos.

Si on peut le faire en python, on peut le faire en SAS... nous allons cependant légèrement modifier leur approche.

Ils vous disent : imaginons une ville contenant 250 verts et 250 oranges. Notre ville contiendra elle aussi 500 habitants mais le type (bleu/rouge) sera déterminé au hasard.

La ville de Sargent et Stachurski, tout comme la notre, est en fait un carré de 1*1 et la position de chaque habitant (x ; y) avec 0<x<1 et 0<y<1, est tirée au hasard.

la distance qui sépare chaque habitant est une distance euclidienne et, comme nous l'avons dit plus haut, un habitant est heureux s'il a au moins, parmi ses 10 voisins les plus proches, au moins 5 voisins du même type que lui : il est identiquement heureux si 10, 9, 8, 7, 6 ou 5 voisins sont identiques à lui (et donc peut être parfaitement heureux s'il a jusqu'à 5 voisins différents de lui) : il n'y a pas d'aversion à vivre dans un environnement mixte.

Quand l'individu n'est pas heureux, il déménage. Sargent et Stachurski indiquent qu'il va déménager vers un endroit où il pourrait être heureux, nous allons ici procéder différemment en considérant que l'individu déménage "au hasard" : lorsque je déménage, je ne choisis pas de moi même d'aller dans un ghetto.

Sargent et Stachurski montrent qu'en quelques cycles, on aboutit à un équilibre stable (dans lequel plus personne ne souhaite déménager) : la "carte" de la ville vous montre que c'est au prix de la formation de ghettos, zones peu mixtes.

comment montrer cela avec SAS ?

Construisons une ville et plaçons y 500 habitants rouge/bleu au hasard :

le nombre d'habitants de la ville est saisi dans une macro variable. Si l'envie nous prend de construire une ville avec un nombre différent d'habitants, nous n'aurons qu'à changer la valeur de la macro variable NB.

Code:

%let nb=500 ;

data loc(keep=x: y: type:) ;
   array type(&nb);
   array x(&nb);
   array y(&nb);

   do i=1 to &nb;
      type(i)=round(rand("uniform"),1);
      x(i)=rand("uniform");
      y(i)=rand("uniform");
   end;
run;

voilà, c'est fait : on dispose d'une table d'une observation et 1500 variables : X1 et Y1 vous donne la position du premier individu, TYPE1 vous donne son type. X2 et Y2 donne la position du second individu, son type est renseigné par TYPE2, etc.

Vous pouvez dès maintenant faire une première "carte de la ville". On va partir de la table LOC mais allons devoir la "transposer par partie" avant de réunir ces parties afin de disposer d'un table contenant trois variables pouvant être mobilisées par PROC SGPLOT et SCATTER.

Code:

proc transpose data=loc(keep=x:) prefix=X out=lac1;
run;

proc transpose data=loc(keep=y:) prefix=y out=lac2;
run;

proc transpose data=loc(keep=type:) prefix=t out=lac3;
run;

data lac4;
   merge lac1 lac2 lac3;
run;

ods graphics / width=600 height=600;

title "la carte initiale";

proc sgplot data=lac4 noautolegend ;
   scatter x=x1 y=y1 / group=t1 markerattrs=(symbol=circlefilled size=8pt);
   xaxis display=(nolabel);
   yaxis display=(nolabel);
run;

et on obtient :

http://www.sas-sr.com/img/ghetto1.png

La semaine prochaine, à partir de cet état initial de la ville, nous verrons comment construire trois nouvelles séries de variables :
- des variables DIST qui vont mesurer les distances qui me séparent de mes voisins
- des variables SAME qui indiqueront combien de voisins identiques à moi j'ai parmi mes 10 voisins les plus proches
- des variables MOVE qui vont indiquer si je vais déménager ou pas...

à suivre

Hors ligne

 

#2 26-06-2019 07:05:03

SAS-SR
Administrateur
Lieu: Université d'Orléans
Date d'inscription: 01-09-2008
Site web

Re: Liliane, fais les valises, on rentre à Paris

Poursuivons

Je vous l’avais annoncé la semaine dernière, nous avons trois ensembles de variables à construire : les variables DIST, les variables SAME et les variables MOVE.

Créons ces ensembles à la suite du programme qui a créé nos individus et leur a attribué un type :

Code:

%let nb=500 ;

data loc(keep=x: y: type:) ;
   array type(&nb);
   array x(&nb);
   array y(&nb);

   do i=1 to &nb;
      type(i)=round(rand("uniform"),1);
      x(i)=rand("uniform");
      y(i)=rand("uniform");
   end;

   array same(&nb);
   array dist(&nb);
   array move(&nb);

et voilà !

bon... il faut maintenant leur donner une modalité...

En ce qui concerne les variables DIST, vous remarquez qu’elles sont construites au moyen d’une instruction ARRAY et qu’il y aura autant de variables DIST qu’il y a d’individus dans ma « ville ». Il ne s’agit pas en effet de calculer toutes les distances entre tous les individus mais de calculer, individu par individu, la distance qui le sépare des autres. La différence peut paraitre subtile… mais la première approche nous obligerait à gérer un tableau double (500 , 500) (250 000 variables !), nous nous contenterons d’un tableau simple de 500 variables (puisque 500 individus).

Calculons ces distances :

Code:

do i=1 to &nb;
   same(i)=0;
   do j=1 to &nb ;
      dist(j)=sqrt((x(i)-x(j))**2+(y(i)-y(j))**2);
      if i=j then dist(j)=.;
   end;

Notre première boucle sur i nous dit : ce qui suit va être fait pour chacun de nos 500 individus. Si i=1, nous traitons le premier individu. On fixe au départ la variable qui nous dira son nombre de voisins identique à lui, SAME(1) à zéro. Ensuite, la boucle sur j va nous permettre de calculer la distance euclidienne qui sépare notre premier individu de tous les autres individus. Si i=j, on calcule une distance entre l’individu et lui-même, elle est forcément égale à zéro mais nous la mettons en valeur manquante.

La question suivante est maintenant de savoir quels sont les 10 voisins les plus proches de notre premier individu et aussi de savoir combien, parmi ces 10 voisins les plus proches, sont du même type que notre premier individu.

Code:

  do pic=1 to 10 ;
      do j=1 to &nb ;
         if dist(j)=smallest(pic,of dist:) then do ;
            same(i)+(type(i)=type(j));
            leave;
         end;
      end;
   end;

C’est la boucle sur PIC qui va permettre de trouver les 10 voisins les plus proches. Pour cela, il nous faut explorer, au moyen de la nouvelle boucle sur J, les variables DIST qui viennent juste d’être calculées.

On va regarder chacune des distances afin de voir laquelle est la plus faible (lorsque pic=1), laquelle est la seconde la plus faible (lorsque pic sera égal à 2), etc. Vous devez comprendre ici pourquoi nous fixons la distance de l'individu avec lui même à valeur manquante. Si nous ne le faisons pas, le voisin le plus proche du premier individu, c'est lui même... (et vous comprenez aussi que nous utilisons la fonction SMALLEST plutôt que ORDINAL (voir section 3.2.3 ED4)).

Si DIST(J) est la distance la plus faible, alors la variable SAME(1) (relative à notre premier individu) augmentera de 1 si le type de notre premier individu est identique au type de j.

Lorsqu’on aura examiné les types des 10 voisins les plus proches, on peut alors voir si notre premier individu va déménager ou pas :

Code:

   move(i)=(same(i)<5);
end;
totmov=sum(of move:);

Il déménage s’il n’a pas au moins 5 voisins du même type que lui. Le END clôt la boucle sur i, on remonte en haut de la boucle pour attaquer le second individu.

Lorsque tous les individus auront été analysés, la variable TOTMOV vous donnera le nombre de déménagements qu’il va falloir organiser.

Une dernière remarque : la boucle sur I débute par SAME(i)=0.

Pourquoi ?

Par mesure de sécurité… les variables SAME sont construites par incrémentation (same(i)+(type(i)=type(j)). Dans ce cas-là, (parce qu’on se fait régulièrement avoir…), mieux vaut partir d’une valeur dont on est certain…

Notre programme, in fine, construira une table contenant plusieurs observations. En ne fixant pas les variables SAME à zéro à l’intérieur de la boucle sur I, lorsque nous construirons la seconde observation, nous allons repartir avec les valeurs des variables SAME définies lors de la construction de la première observation et vous pourrez alors observer des cas dans lesquels, parmi les 10 voisins les plus proches d’un certain indique, il y a plus de 10 voisins qui sont du même type que cet individu… ce qui n’aurait aucun sens…

La suite ?

Et bien il faut faire déménager ceux qui doivent déménager puis voir si, suite à ces déménagements, certains ne souhaitent pas déménager.

Qui peuvent être les candidats à un déménagement dans ces phases ultérieures ?
-    Les personnes qui viennent de déménager et qui se retrouvent avec moins de 5 voisins du même type qu’eux
-    Les personnes qui, suite aux déménagements de leurs voisins, se retrouvent avec moins de 5 voisins du même type qu’eux.

on va donc devoir réexaminer la situation de chacun de nos 500 individus...

à suivre…

Hors ligne

 

#3 03-07-2019 09:57:12

SAS-SR
Administrateur
Lieu: Université d'Orléans
Date d'inscription: 01-09-2008
Site web

Re: Liliane, fais les valises, on rentre à Paris

Voici notre programme à ce stade :

Code:

%let nb=500;

data loc(keep=x: y: type:) ;
array type(&nb);
array x(&nb);
array y(&nb);

do i=1 to &nb;
   type(i)=round(rand("uniform"),1);
   x(i)=rand("uniform");
   y(i)=rand("uniform");
end;

array same(&nb);
array dist(&nb);
array move(&nb);

do i=1 to &nb;
   same(i)=0;
   do j=1 to &nb ;
      dist(j)=sqrt((x(i)-x(j))**2+(y(i)-y(j))**2);
      if i=j then dist(j)=.;
   end;
   do pic=1 to 10 ;
      do j=1 to &nb ;
         if dist(j)=smallest(pic,of dist:) then do ;
            same(i)+(type(i)=type(j));
            leave;
         end;
      end;
   end;
   move(i)=(same(i)<5);
end;

totmov=sum(of move:);
Run;

au terme du programme, nous avons analyser la situation initiale de nos individus, savons qui va déménager et qui restera en place. Avant de faire évoluer quoi que ce soit, il nous faut garder une trace de cette première analyse... et il suffit d'ajouter une instruction OUTPUT juste après l'instruction qui calcule la modalité de TOTMOVE (le nombre de déménagements qu'il va falloir organiser.

ensuite ?

et bien déménageons ceux qui souhaitent bouger !

Code:

output ;

do i=1 to &nb;
   if move(i)=1 then do ;
     x(i)=rand("uniform");
     y(i)=rand("uniform");
   end;
end;

et voilà, ils ont déménagé !

ce qu'il nous faut faire ensuite, c'est vérifier quels sont les individus, parmi l'ensemble de nos individus qui vont maintenant vouloir déménager. Il nous faut donc "remonter" dans notre programme, pour vérifier combien, parmi les 10 voisins les plus proches, sont du même type de l'individu pour ensuite faire déménager ceux qui le souhaite.

Nous sommes appelé à "remonter dans le programme jusqu'à ce que plus personne ne souhaite déménager. Nous allons pour cela utiliser une boucle DO UNTIL et la condition d'arrêt de la boucle, c'est TOTMOVE=0. Le programme devient :

Code:

data loc(keep=x: y: type:) ;
array type(&nb);
array x(&nb);
array y(&nb);

do i=1 to &nb;
   type(i)=round(rand("uniform"),1);
   x(i)=rand("uniform");
   y(i)=rand("uniform");
end;

array same(&nb);
array dist(&nb);
array move(&nb);

do until (totmov=0);
   do i=1 to &nb;
      same(i)=0;
      do j=1 to &nb ;
         dist(j)=sqrt((x(i)-x(j))**2+(y(i)-y(j))**2);
         if i=j then dist(j)=.;
      end;
      do pic=1 to 10 ;
         do j=1 to &nb ;
            if dist(j)=smallest(pic,of dist:) then do ;
               same(i)+(type(i)=type(j));
               leave;
            end;
         end;
      end;
      move(i)=(same(i)<5);
   end;
   totmov=sum(of move:);
   output;
   do i=1 to &nb;
      if move(i)=1 then do ;
         x(i)=rand("uniform");
         y(i)=rand("uniform");
      end;
   end;
end;
run;

un petit coup d'oeil à votre journal vous montrera en combien d'étapes vous êtes arrivé à un "équilibre".

Code:

NOTE: The data set WORK.LOC has 16 observations and 1500 variables.

16 observations = 15 déménagements jusqu'à l'équilibre dans lequel plus personne ne souhaite bouger.

il ne nous reste plus qu'à produire un ensemble de graphiques qui vont nous permettre de voir les évolutions successives de notre ville en train de ce "ghettoïser"

à suivre

Hors ligne

 

#4 Aujourd'hui 08:21:53

SAS-SR
Administrateur
Lieu: Université d'Orléans
Date d'inscription: 01-09-2008
Site web

Re: Liliane, fais les valises, on rentre à Paris

Avant de partir en vacances, il nous faut conclure ce sujet des beaux mercredis....

Nous disposons maintenant d'une table dans laquelle est présentée l'organisation de notre "ville". Une observation dans cette table représente notre ville à un instant t : la première observation est relative à votre ville au départ, avant tout déménagement, la dernière observation de la table représente notre ville "à l'équilibre" : plus aucune personne ne souhaite déménager. La table construite dans le précédent post contenait 16 observations : on doit donc produire, à partir ce cette table, 16 graphiques.

Nous allons pour cela rédiger un macro-programme mais commençons par modifier le programme rédigé jusqu'à maintenant pour "préparer le terrain" :

notre programme devient :

Code:

%let nb=500 ;

data loc(keep=x: y: type:) ;
array type(&nb);
array x(&nb);
array y(&nb);

do i=1 to &nb;
   type(i)=round(rand("uniform"),1);
   x(i)=rand("uniform");
   y(i)=rand("uniform");
end;

array same(&nb);
array dist(&nb);
array move(&nb);

do until (totmov=0);
   iter+1;
   do i=1 to &nb;
      same(i)=0;
      do j=1 to &nb ;
         dist(j)=sqrt((x(i)-x(j))**2+(y(i)-y(j))**2);
         if i=j then dist(j)=.;
      end;
      do pic=1 to 10 ;
         do j=1 to &nb ;
            if dist(j)=smallest(pic,of dist:) then do ;
               same(i)+(type(i)=type(j));
               leave;
            end;
         end;
      end;
      move(i)=(same(i)<5);
   end;
   totmov=sum(of move:);
   call symputx(compress("nbmov"||iter),totmov);
   output;
   do i=1 to &nb;
      if move(i)=1 then do ;
         x(i)=rand("uniform");
         y(i)=rand("uniform");
      end;
   end;
end;
call symputx("nbiter",iter);
run;

les nouveautés de ce programme ?

1- juste après le DO UNTIL, nous demandons la construction d'une variable ITER qui va augmenter de 1 à chaque tour dans la boucle
2- juste après la construction de TOTMOV, nous enregistrons dans une série de macro-variables, le nombre de déménagements à effectuer : la macro variable NBMOV1 indique le nombre de déménagements à faire dans la situation initiale, NBMOV8 donnera le nombre de déménagements à la 8ème étape.
3- juste avant le RUN final, nous enregistrons dans une macro variable le nombre d'itérations effectuées.

La suite ? proche dans l'esprit du programme présenté dans le premier post pour construire le graphique présentant la situation initiale.

Code:

%macro laville ;
ods html close;
ods html;
proc transpose data=loc(keep=x:) prefix=x out=lac1;
run;

proc transpose data=loc(keep=y:) prefix=y out=lac2;
run;

proc transpose data=loc(keep=type: obs=1) prefix=t out=lac3;
run;

data lac4;
merge lac1 lac2 lac3;
run;

Nous refaisons notre transposition "par partie" : la premier PROC TRANSPOSE nous permet d'obtenir une table de 500 observations (puisqu'on a 500 variables X) et qui aura un nombre de variables égal au nombre d'itérations effectuées. X1 donnera les coordonnées X des 500 individus à l'état initial, X2, les coordonnées X suite à la première série de déménagements etc.

Nous faisons la même chose avec les variables Y au moyen du second PROC TRANSPOSE. Le troisième est légèrement différents : les 500 variables TYPE ne connaissent pas d'évolution dans leur valeur au cours des itérations : on ne transpose donc que la première observation de la table LOC. Les trois tables ainsi créés sont réunies au moyen de MERGE. Pas besoin de variable BY ici...

ensuite ?
une jolie boucle macro pour créer nos graphiques !

Code:

ods graphics / width=600 height=600;
%do i=1 %to &nbiter;
   %if &&nbmov&i>1 %then %do ;
      title "itération n°&i - &&nbmov&i mouvements vont intervenir";
   %end;
   %if &&nbmov&i=1 %then %do ;
      title "itération n°&i - &&nbmov&i mouvement va intervenir";
   %end;
   %if &&nbmov&i=0 %then %do ;
      title "itération n°&i - équilibre - fin des déménagements";
   %end;
   proc sgplot data=lac4 noautolegend ;
      scatter x=x&i y=y&i / group=t1 markerattrs=(symbol=circlefilled size=8pt);
      xaxis display=(nolabel);
      yaxis display=(nolabel);
   run;
%end;

%mend;

%laville

et voilà !

le titre à produire au dessus des graphiques dépend de la valeur de la macro NBMOV associée et vous obtenez ainsi ce résultat :

http://www.sas-sr.com/img/bm719/sashtml5.htm

(cette fois ci, il y a eu 19 états jusqu'à l'équilibre - le fait que le nombre d'états changent régulièrement dans la description des programmes ne doit pas vous étonner : puisque nous ne fixons pas la graine au moyen de CALL STREAMINIT (voir page 170-171 ED4), c'est l'horloge qui est utilisée pour l'initialisation du générateur de base).

bref...

ce post clôt la 8ème saison des beaux mercredis de www.sas-sr.com - rendez vous en septembre pour une 9ème saison !

Hors ligne

 

Pied de page des forums

Propulsé par FluxBB
Traduction par FluxBB.fr
Flux RSS