Version pour Octave / MATLAB

Voici une version pour Octave / MATLAB. À la fin, on affiche le pourcentage de la richesse totale détenu par les 1% les plus riches...

J'ai regardé les temps de calcul (pour Python avec le code publié plus haut) :

Python => 10.5 s
Octave 5.2.0 => 320 s
MATLAB 2020a => 1 s

Je suis quand même étonné de la mauvaise performance d'Octave, même si je sais qu'il s'en sort mal quand il y a des boucles...

-------- Fichier richesse.m ----------
% 27/07/2021
% Modèle de marché libre du type "vide grenier"
% Cf. article de Pour la Science par B. Boghosian
close all;
clear all;
tic;
 %% Paramètres
 % Nombre de transactions
 nb_transactions = 5000000;
 % Pourcentage de la richesse du plus pauvre gagné par le plus pauvre à chaque transaction gagnante
 pct_gain_pauvre = 20;
 % Pourcentage de la richesse du plus pauvre gagné par le plus riche à chaque transaction gagnante
 pct_gain_riche = 17;
 % Nombre d'agents
 nb_agents = 1000;
 % Richesse initiale (identique pour tous les agents)
 richesse_ini = 1000;
 % Seed des générateurs aléatoires
 rand('state', 35291);
 %% Simulation
 % Initialisation
 Richesse = richesse_ini * ones(1, nb_agents);
 % Boucle sur les transactions
 for nt = 1:nb_transactions
 if (mod(nt, 100000) == 0)
 fprintf(1, 'Transaction %d\n', nt);
 end
 %% Tirages aléatoires
 % Tire aléatoirement deux agents distincts
 [A, B] = randi2(nb_agents);
 % Tire à pile ou face si A perd ou gagne
 if (rand < 0.5)
 A_perd = true;
 else
 A_perd = false;
 end
 %% Transaction
 % A perd
 if A_perd
 % A et B ont la même richesse
 if (Richesse(A) == Richesse(B))
 % Même gain (personne n'est avantagé)
 gain = pct_gain_pauvre/100 * Richesse(A);
 Richesse(B) = Richesse(B) + gain;
 Richesse(A) = Richesse(A) - gain;
 % A pauvre, B riche (A est avantagé)
 elseif (Richesse(A) < Richesse(B))
 gain = pct_gain_riche/100 * Richesse(A);
 Richesse(B) = Richesse(B) + gain;
 Richesse(A) = Richesse(A) - gain;
 % A riche, B pauvre (B est avantagé)
 else
 gain = pct_gain_pauvre/100 * Richesse(B);
 Richesse(B) = Richesse(B) + gain;
 Richesse(A) = Richesse(A) - gain;
 end
 % A gagne
 else
 % A et B ont la même richesse
 if (Richesse(A) == Richesse(B))
 % Même gain (personne n'est avantagé)
 gain = pct_gain_pauvre/100 * Richesse(A);
 Richesse(A) = Richesse(A) + gain;
 Richesse(B) = Richesse(B) - gain;
 % A pauvre, B riche (A est avantagé)
 elseif (Richesse(A) < Richesse(B))
 gain = pct_gain_pauvre/100 * Richesse(A);
 Richesse(A) = Richesse(A) + gain;
 Richesse(B) = Richesse(B) - gain;
 % A riche, B pauvre (B est avantagé)
 else
 gain = pct_gain_riche/100 * Richesse(B);
 Richesse(A) = Richesse(A) + gain;
 Richesse(B) = Richesse(B) - gain;
 end
 end
 end
 %% Pourcentage de la richesse totale détenue par les 1 % les plus riches
 R = quantile(Richesse, 0.99);
 k = find(Richesse >= R);
 pct = sum(Richesse(k)) / sum(Richesse);
 fprintf(1, '\nLes 1 %% les plus riches détiennent %.1f %% de la richesse totale\n\n', pct * 100);
 toc;
 --------- Fichier randi2.m --------
 % Tire aléatoirement deux entiers distincts entre 1 et N
 function [m, n] = randi2(N)
 m = floor(1 + rand * N);
 n = m;
 while n == m
 n = floor(1 + rand * N); 
 end
 end