// EVOGOD // SciLab program to simulate the evolution of religion // Version: 1.0 // Modification date: 1/28/2008 // Copyright 2008 James W. Dow // This program is released under the GNU General Public License Version 2, June 1991. See http://www.gnu.org/licenses/gpl.txt // SYSTEM PARAMETERS ns = 140; ip1 = 10; le = 65; mc = 5; dcr = 0.01; dcu = 0.01; dcur = 0.01; dfir = 0.010; dfiu = -0.005; dfisu = -0.02; dp = 0; lr = 0; ur = 1; lv = 0; uu = 1; dist = 0; popgrow = 0; greenbeard = 1; mult =1; cmr = 0.5; cmu = 0.5; icp = 0.2; // MODEL // Build the initial population matrix (a) of agents . We don't know the size yet. It is easier to build it and then compute the size. Each agent is a row. nva = 6; // Number of variables (columns) for each agent vector // The column variables for each agent are // 1. age (ag) // 2. initial probability of communicating a real experiences (ipr). A real experience is an experience that gives the holder knowledge that can enhance the fitness of others. This knowledge can be communicated. // 3. initial probability of communicating unreal experiences (ipu) A number representing the probability that the agent will personally have an unreal experience. An unreal experience is something that can be communicated but contains no knowledge that can enhance fitness. The probability variables pr and pu represent biological evolution in the brain. // 4. developing probability of communicating a real experience (cpr). This changes during the simulation // 5. developing probability of communicating an unreal experience (cpu). // 6.fitness (fi) The number of offspring produced over a lifetime. // INITIALIZATION of ages and number of agents // The number of agents in cohort k is ip1 - (k-1)*dp. // Function for the number of agents in cohort k function y = na(k), y = ip1 - (k-1) * dp, endfunction; // Set up ages // i indexes the cohort, j indexes the column of the a matrix. // Do the first cohort. a(1) contains the age. The maximum number of cohorts is le, the initial life expectancy. However there are ilkley to be fewer if the population in each cohort decrements (dp) by a large amount a = ones(na(1),nva); // Do remaining cohorts for i = 2:le; if na(i) < 0 then break; end; b = ones(na(i),nva) * i ; a = cat(1,a,b); end; // Determine the size of the population // Function for the number of rows in a matrix. function y = sz(x), y=size(x,'r'); endfunction; // Set initial population size ips = sz(a); // INITIAL COMMUNICATION PROBABILITIES // Set up initial probalities of communicating a real experience, a(:,2) and, an unreal experience, a(:,3). // Set up a random vector of probabilities. If dist = 0 use a uniform distribution. If dist = 1 use a binomial distribution. select dist ; case 0 then; // Uniform distribution // Function to set up a random vector of uniformly distributed probabilities between l and u. function y = rv(l,u), y = (u - l) * rand(ips,1) + l, endfunction; a(:,2) = rv(lr,ur); a(:,3) = rv(lv,uu); case 1 then; // Binomial distribution a(:,2) = grand(ips,1,'bin',ip1,ur) / ip1; a(:,3) = grand(ips,1,'bin',ip1,uu) / ip1; end; // Learned communication probabilities. Start with icp and add according to inherited fitness. a(:,4:5) = ones(ips,2)*icp; // INITIAL FITNESS a(:,6) = ones(ips, 1); // Set up recording for the steps. rps = ones(ns,1); // Population size mcp = zeros(ns,2); // Mean of the probability of commnicating real and unreal information // STEP THROUGH THE SIMULATION // Note that in this simulation it is possible for an agent to communicate with himself. This does not seem to be a conceptual problem and occurs in reality. // Functions to provide a set of x random indices to agents in a. greenbeard determines the distribution from which the random indices are selected. select greenbeard; case 0 then ; // The probabilities of selection are distributed uniformly over the population function y=ria(x), y1=grand(1,'prm',(1:sz(a))'), y=y1(1:x), endfunction; case 1 then ; // The probabilities of selection are distributed according to the current probabilities, cpu, of generating unreal communications. The higher the cpu of an agent, the greater the chance of it being selected. function y=ria(x); ps = size(a,'r') // Normalize the cpu so that they sum to 1 y1 = a(:,5) / sum(a(:,5)) // Form a cumulative distribution y2 = cumsum(y1) y = [] // Initialize y as a blank vector. for i = 1:x y3 = y2 < rand() // A T-F vector for lower values y =[y, sum(y3)+1] // Counts agents lower in the cumulative distribution and adds the count to the vector, thus creating the index of the area within which the random number lies. end endfunction; end; // START STEPPING for s = 1:ns; // Do the step // Calculate population size for this step ps = sz(a); nrc = 0; // monitor total real communications nuc = 0; // monitor unreal communications for j = 1:ps; // Do the agents // Communicate REAL experiences // The number of agents communicated with (nca) will depend on the communnicators probablity of communicating real information. It will be a proportion of mc nca = round(a(j,4) * mc) ; // Get the number of agents to receive real communications from this one. // Select a set of agents to receive real communications. Note that i is a vector of indices. i = ria(nca); // select a subset of agents in a at random. i is a column vector. It can be empty // Have them learn to communicate real experiences if ~ (i == []) then; // Change the probability of the receivers communicating real information. This will depend on the inherited probability of communicating unreal real information a(i,4) = min(ones(size(i,'r'),1),a(i,4) + a(i,2)*dcr); // Add to their fitness. if mult == 0 then; a(i,6) = a(i,6) + (cmr * dfir); // Multiply the amount by a constant, cmr else; // Multiply by the receivers learned probability of communicating real information a(i,2) a(i,6) = a(i,6) + (a(i,4) * dfir); end; end; // Communicate UNREAL experiences. The number of agents contacted will depend on the probability of communicating unreal information nca = a(j,5) * mc ; // Calculate number nuc = nuc + nca; // Add to unreal communications // Select a set of agents randomly i = ria(nca); // note that i is a vector and possibly null if there are less than .5 agents to contact // Are there agents to contact if ~ (i == []) then; // Change the probability of the receivers communicating unreal information a(i,5) = min(ones(size(i,'r'),1),a(i,5) + a(i,3)*dcu); // Change the probability of the receiver's communicating real information a(i,4) = min(ones(size(i,'r'),1),a(i,4) + a(i,2)*dcur); //Change the fitness of the receiver. dfiu will normally be negative if mult == 0 then; a(i,6) = a(i,6) + (cmu * dfiu); // Multiply by a constant else; a(i,6) = a(i,6) + (a(i,5) * dfiu); //Multiply by the recivers learned probability of making unreal communications end; end; // Make the sender pay a cost for sending unreal information. This is the cost of costly signaling, and can occur even if there are no receivers a(j,6) = a(j,6) + dfisu; // Note dfisu is normally negative end; // Agents are done. j is free // EVOLVE THE POPULATION OF AGENTS // Advance the ages a(:,1) = a(:,1) + 1; //Find agents who will die ida = []; //Start with a blank index vector for i = 1:sz(a); if a(i,1) > le then; ida = [ida ; i]; // Add the index to ida end; end; // Index of dead agents (ida) is ready // Add offspring of dead agents // Calculate the number of offspring of each dead agents no = round(a(ida,6)); // vector of the number of offspring of agents indexed by ida // Form new offspring with age 1 and fitness 1 (not inherited) // Maintain the probabilities of communication of the parents for i = 1:sz(ida); // Scan dead agents j = ida(i) ; // index of the dead agent in a for k=1:no(i); // Scan number of offspring new = [1, a(j,2), a(j,3), icp, icp, 1]; // Catenate new agent at the end of a a($+1,:) = new; end; end; a(ida,:) = []; // Remove the dead ones from a. // POPULATION GROWTH CONTROL // Function to randomly select x agents selected from a uniform distribution. function y=riau(x), y1=grand(1,'prm',(1:sz(a))'), y=y1(1:x), endfunction; select popgrow ; case 0 then; // Population can grow or decline. Do nothing case 1 then;// Control the increase or the decline dps = sz(a) - ips; // Calculate the change in population size if dps > 0 then ; // There is an increase in population // Remove dps agents at random i = riau(dps); // Get a uniform set at random a(i,:) = [] ; // Kill them off elseif dps < 0 then ; // There is a decrease in population i = riau(-dps) ; // Pick a set of agents at random // Have them reproduce // Maintain the probabilities of communication of the parents new = [ones(-dps,1), a(i,2:3), a(i,2:3), ones(-dps,1)]; // Catenate new agents at the end of a a = [a ; new]; end; // End if // There is neither an increase of a decrease end; // End of population control // save results for this step ps = sz(a); // Calculate population size rps(s,1) = ps; // Record the population size mcp(s,:) = [mean(a(:,2)), mean(a(:,3))]; // record means of pr and pu // immediate display printf("step %i population %i \tmean pr %f mean pu %f \n",s,ps,mcp(s,1),mcp(s,2)); end; // Move to next step // End of simulation // PLOT RESULTS xbasc(0); scf(0); // Open first plot window. plot2d(1:ns,rps); // Population size per step ax = get("current_axes"); ax.x_label.text = "Step"; ax.x_label.font_size = 3; ax.y_label.font_size = 3; ax.sub_ticks = [3,1]; xbasc(1); scf(1); // Open second plot window. plot2d(1:ns,mcp); // Probability of real and unreal communications per step ax = get("current_axes"); ax.x_label.text = "Step"; ax.x_label.font_size = 3; ax.y_label.font_size = 3; ax.sub_ticks = [3,1]; // Plot histograms of pr and pu xbasc(2); scf(2); // Open third plot window title("Final Distribution of ipr", 'fontsize', 4); x = nfreq(int(a(:,2)*10));bar (x(:,1),x(:,2)) xbasc(3); scf(3); // Open fourth plot window title("Final Distribution of ipu", 'fontsize', 4); x = nfreq(int(a(:,3)*10));bar (x(:,1),x(:,2)); //