Имитируя сближение фамилия в популяции


Мне было любопытно несколько недель назад о том, как фамилии приходят и уходят (в основном на ходу). Я хотел бы знать, как быстро фамилий вымрем, потому что (в моем упрощенном мире моделирования) после того, как последний человек, с фамилией умирает, он ушел навсегда. Я начал писать небольшой симуляции, чтобы это произошло в "реальном времени". Я начал с простого "человека" представления, которые вели себя так же, как и все остальные "лица" в моделировании. Однако, как я попала в проект, я начала придумывать определенные черты характера, так как не каждый живет в том же возрасте, не у всех есть такое же количество детей, не у всех пар, даже у детей и т. д.

Несколько предположений:

  • Я не основываю это моделирование на каких-либо фактических статистических данных о населении, хотя я думаю, что это делает достойную работу производить некоторые интересные результаты.
  • Фамилии не восстанавливаются. Стартовая населения каждый получает уникальное "имя" (целое число, на единицу больше, чем последняя), что и формирует максимальное количество имен в симуляции.
  • Это отечески предвзято. Чтобы держать вещи простыми, самок населения берут на себя мужские фамилии в сочетании. Фамилии не бывают комбинированные (см. пункт 2), или изменены в любом случае. Это моделирование будет работать так же (в основном?) если по материнской линии предвзятым, хотя я предполагаю, что в этот момент мужчины нужны только для роста населения (тема не актуальна для этой модели, просто некоторые побочные мысли).
  • Фамилии, унаследованной от родителя. Я изначально не считаю это предположение, но, по словам Тоби спит, это не всегда верно в каждой культуре, Исландия, например.
  • Люди умирают, когда они достигают своей максимальной продолжительности жизни, но и может умереть рано из-за "несчастных случаев", будучи более склонен старые они.
  • Когда человек с партнером умер, оставшийся в живых партнер не получает нового партнера.
  • Экспоненциальное распределение используется здесь, чтобы помочь отойти от равномерного распределения, поскольку я обнаружил, что мундир населения умирают очень легко. Примером может быть, что Person более правоподобны для того чтобы иметь жизненный период 80, чем 20, но, что падает существенно за 80.

Поэтому я придумал следующее Person объект выступают в качестве отдельных единиц моего моделирования:


# Person.py

import random
import numpy as np

male = 1
female = 2

# these reduce performance contributions from lookups by half, and because of how often they are called, any bit helps
# # $ python -m timeit -s "import numpy as np" "np.random.randint(0, 1000)"
# # 10000000 loops, best of 3: 0.168 usec per loop
#
# # $ python -m timeit -s "import numpy as np; npRandInt = np.random.randint" "npRandInt(0, 1000)"
# # 10000000 loops, best of 3: 0.0814 usec per loop
npRandInt = np.random.randint
npRandUni = np.random.uniform
npRandEx = np.random.exponential


class Person(object):
    currentName = 0
    deaths = 0
    maxAgeDistribution = []
    maxChildDistribution = [0, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5]

    def __init__(self, parent=None):
        if parent is None:
            Person.currentName += 1
            self.surname = Person.currentName
        else:
            self.surname = parent.surname

        self.gender = random.choice([male, female])

        self.currentChildCount = 0
        self.currentAge = 0

        self.maxChildren = random.choice(Person.maxChildDistribution)  # cap how many children one person can have
        self.maxAge = random.choice(Person.maxAgeDistribution)

        self.childDifficulty = npRandInt(0, 100 + 1) / 100.0  # how difficult is it to have a child
        self.childYearsStart = npRandInt(10, 15 + 1)
        self.childYearsEnd = npRandInt(30, 50 + 1)

        self.partner = None  # some other person object
        self.withChild = False
        self.hasHadPartner = False  # give a way to check if they have had a partner before, so if one dies, they don't get another

        # determine age person is eligible to gain a partner. Having this too high (on average) leads to excessive population death
        if self.childYearsStart < self.maxAge:
            self.partnerEligibleAge = self.childYearsStart + int(self.childYearsStart * npRandEx(0.2))

            if self.partnerEligibleAge > self.maxAge:
                self.partnerEligibleAge = self.maxAge
        else:
            self.partnerEligibleAge = self.maxAge

    def __repr__(self):
        return 'Person(Gender: %6s, Surname: %5d, Age: %2d, MaxAge: %2d, Children: %1d, MaxChildren: %1d, Partner: %5s, PartnerEligibleAge: %2d)' % \
               ("Male" if self.isMale() else "Female", self.surname, self.currentAge, self.maxAge, self.currentChildCount, self.maxChildren,
                "True" if self.partner is not None else "False", self.partnerEligibleAge)

    def gainPartner(self, partner):  # gain a partner if both persons are eligible, female taking the surname of the other partner
        if self.isPartnerEligible():
            if self.gender == female:
                self.surname = partner.surname

            self.partner = partner
            self.hasHadPartner = True

    def age(self):
        self.currentAge += 1

        if self.currentAge >= self.maxAge:  # if they outlive their lifespan, they die
            return self.die()

        percentAge = 1.0 - ((self.maxAge - self.currentAge) / float(self.maxAge))  # get their current age as a percent of their max age

        # get a random [0, 1] with an exponential distribution.
        # Magic numbers were determined experimentally, seeing what allowed population growth, and then by checking the distribution using pyplot
        randDeathPercent = 1.0 - npRandEx(0.05) * 2.0

        # random accidents happen, so even if they haven't reached their max age,
        # there is some risk of death, especially as they get closer to their max age
        return self.die() if percentAge > randDeathPercent else False  # perform self.die, which returns True, if they have died

    def gainChild(self):  # gain a child if eligible and passes a "difficulty" check
        if self.isChildEligible():
            chance = 1.0 - (self.childDifficulty * self.partner.childDifficulty)
            if chance > npRandUni(0, 1):
                self.withChild = True

    def haveChild(self, modifier=2.0):  # produce child if they currently have one and pass a modified "difficulty" check
        if self.withChild is True:
            self.withChild = False  # pass or fail, they will no longer have a child

            chance = 1.0 - self.childDifficulty
            if chance * modifier > npRandUni(0, 1):  # if this modifier is set too low, not enough children are created and the population dies
                self.currentChildCount += 1
                return Person(parent=self)  # return a new person, giving it the last name of the parent

        return None  # if we get this far, a child fails to be created

    def die(self):
        # if they have a partner, remove themselves from that person
        if self.partner is not None:
            self.partner.partner = None
        Person.deaths += 1  # add to the death stats
        return True

    def hasPartner(self):
        return self.partner is not None

    def isChildEligible(self):  # check if person can have a child
        # Criteria:
            # only females can have children
            # they must have a partner
            # they haven't had more than their max children
            # they are within child bearing years, and so is their partner
            # they don't currently have a child
        return self.gender == female and \
               self.partner is not None and \
               self.currentChildCount < self.maxChildren and \
               self.childYearsStart <= self.currentAge <= self.childYearsEnd and \
               self.partner.childYearsStart <= self.partner.currentAge <= self.partner.childYearsEnd and \
               self.withChild is False

    def isPartnerEligible(self):  # check if person can gain a partner
        # Criteria:
            # they don't currently have a partner
            # they are old enough to have a partner
            # they haven't had a partner in the past
        return self.partner is None and \
               self.currentAge >= self.partnerEligibleAge and \
               self.hasHadPartner is False

    # helper functions to check for gender.
    # Performance profiling indicates these calls are expensive over time, so it is better to check this directly
    def isMale(self):
        return self.gender == male

    def isFemale(self):
        return self.gender == female

Теперь мне нужно на самом деле запустили симуляцию. После большой производительности проверяя, как профилирование и видя результаты моделирования, я смог придумать следующие основные приложения, с main обработка основная часть административной работы, и sim делать реальное моделирование:


# main.py

import random
import numpy as np
import cProfile
import operator
import datetime
import math
import matplotlib.pyplot as plt
import cPickle as pickle
import person
from person import Person
from stats import Stats
import sys
import uuid
import time

npRandInt = np.random.randint
npRandUni = np.random.uniform
npRandEx = np.random.exponential

fileId = str(uuid.uuid1())
SERIALIZE = False


def getTimeDelta(start, stop):
    delta = stop - start
    return (delta.microseconds + (delta.seconds + delta.days * 24 * 3600) * 10 ** 6) / (10 ** 6 * 1.0) * 1000.0


def setupAgeDistribution(distribution):
    for i in range(10000):
        distribution.append(int(getRandomAgeTo80()))

    for i in range(2000):
        distribution.append(int(getRandomAgeFrom80()))


# get a random age in range [0, 80]
# Magic numbers were determined experimentally,
# seeing what allowed population growth,
# and then by checking the distribution using pyplot
def getRandomAgeTo80():
    ran = 80.0 - npRandEx(50.0) * .4
    while ran < 0:
        ran = 80.0 - npRandEx(50.0) * .4
    return ran


# get a random age in range [80, 100]
# Magic numbers were determined experimentally,
# seeing what allowed population growth,
# and then by checking the distribution using pyplot
def getRandomAgeFrom80():
    ran = 80 + npRandEx(50.0) * .1
    while ran > 100.0:
        ran = 80 + npRandEx(50.0) * .1
    return ran


def countItems(aList):  # no access to collections.Counter, so use this instead
    array = np.bincount(np.array(aList))
    nonZero = np.nonzero(array)[0]
    return dict(zip(nonZero, array[nonZero]))


def gainChildren(people):  # have each person try to gain a child, letting the eligibility be determined by the object
    for p in people:
        p.gainChild()


def gainPartners(people):
    # get a list of eligible partners, split by gender
    mPartners = [x for x in people if x.gender == person.male and x.isPartnerEligible()]
    fPartners = [x for x in people if x.gender == person.female and x.isPartnerEligible()]

    numPartners = min([len(mPartners), len(fPartners)])  # get the smaller count of partners

    if numPartners > 0:
        numChoices = npRandInt(0, numPartners + 1)  # select a number of them to gain a partner

        if numChoices > 0:
            # shuffle both lists, splice to make lists the same length, zip them together, and pair up the zipped pairs.
            # Much faster than using random.choice.
            random.shuffle(mPartners)
            random.shuffle(fPartners)

            pairs = zip(mPartners[:numChoices], fPartners[:numChoices])

            for m, f in pairs:
                m.gainPartner(f)
                f.gainPartner(m)


def haveChildren(people):
    newCount = 0
    for p in people.copy():
        if p.withChild is True:
            newChild = p.haveChild(modifier=1.0)
            if newChild is not None:
                people.add(newChild)
                newCount += 1
    return newCount


def ageSim(people):
    for p in people.copy():
        if p.age():  # if they ded
            people.remove(p)


def sim(people):
    ageSim(people)  # age the whole sim, removing those that die
    gainPartners(people)  # pair up eligible people objects
    newCount = haveChildren(people)  # have children before gaining them to "develop" them for a year
    gainChildren(people)  # gain new children to attempt being created next year

    return newCount  # return how many total new people have been added


def main(people=None, stats=None):  # arguments allow resuming from serialized (saved) data
    setupAgeDistribution(Person.maxAgeDistribution)
    random.shuffle(Person.maxAgeDistribution)
    random.shuffle(Person.maxChildDistribution)

    if people is None and stats is None:
        # set up the stats object, setting the initial population count, and how often to print the stats
        stats = Stats(oldDate=datetime.datetime.now(), startingPopulation=10000, statYear=100.0)
        people = set([Person() for i in range(stats.startingPopulation)])  # create a new set of random people

    stats.oldDate = datetime.datetime.now()

    try:
        for year in range(stats.year, stats.year + 500 + 1):
            stats.year = year

            if len(people) == 0:
                print "Population has died! Final stats below:\n"
                printStats(people, stats)
                break

            if year % stats.statYear == 0:
                printStats(people, stats)
            stats.newPersonCount += sim(people)

    except KeyboardInterrupt:  # capture any ctrl-c, and print stats on exit
        printStats(people, stats)

        if len(people) > 0:
            sortedList = sorted(countItems([p.surname for p in people]).items(), key=operator.itemgetter(1))  # sort list by surname popularity
            if len(sortedList) < 1000:  # only print it if the sorted list isn't too long. Console has an issue if the line is too long.
                sortedList.reverse()  # display the list in descending order
                print sortedList

        sys.exit(0)

    # plot(people, stats)


def printStats(people, stats):
    statYear = stats.statYear
    if stats.year % statYear != 0:
        statYear = stats.year % stats.statYear

    delta = getTimeDelta(stats.oldDate, datetime.datetime.now())
    timeTaken = math.floor(delta / 1000.0)

    if timeTaken < 1:
        time.sleep(1)  # don't spam the console output
    stats.timeTaken.append(timeTaken)

    population = len(people)
    stats.population.append(population)

    meanAge = math.floor(np.mean([person.currentAge for person in people])) if len(people) > 0 else 0
    stats.meanAge.append(meanAge)

    stats.newPeople.append(stats.newPersonCount)
    peopleDelta = stats.newPersonCount - stats.oldPeople
    stats.oldPeople = stats.newPersonCount
    stats.deltaNewPeople.append(peopleDelta)

    surnameCount = len(countItems([p.surname for p in people])) if len(people) > 0 else 0
    stats.surnames.append(surnameCount)

    stats.deaths.append(Person.deaths)
    deathsDelta = Person.deaths - stats.oldDeaths
    stats.oldDeaths = Person.deaths
    stats.deltaDeaths.append(deathsDelta)

    netPop = (peopleDelta / statYear) - (deathsDelta / statYear)
    stats.netPopDelta.append(netPop)

    pairCount = len([True for person in people if person.partner is not None]) if len(people) > 0 else 0

    print datetime.datetime.now()
    print "Time taken: %d seconds (%d milliseconds / year)" % (timeTaken, math.floor(delta / statYear))
    print "Current population size at year %d: %d" % (stats.year, population)
    print "Current mean age: %d" % meanAge
    print "Current new people: %d" % stats.newPersonCount
    print "Current deaths: %d" % Person.deaths
    print "Current delta new people: %d (%d people / year)" % (peopleDelta, (peopleDelta / statYear))
    print "Current delta deaths: %d (%d deaths / year)" % (deathsDelta, (deathsDelta / statYear))
    print "Current surname count: %d" % surnameCount
    print "Current partner percentage: %d" % (pairCount / float(population) * 100) if len(people) > 0 else 0
    print "Net population per year: %d" % netPop
    print ""
    stats.oldDate = datetime.datetime.now()

    if SERIALIZE:  # controls whether serialization is enabled, useful to turn off for when debugging or testing
        serialize(people, stats, fileId)


def serialize(people, stats, id=""):
    with open("people%s.dat" % id, "wb") as f:
        pickle.dump(people, f, pickle.HIGHEST_PROTOCOL)

    with open("stats%s.dat" % id, "wb") as f:
        pickle.dump(stats, f, pickle.HIGHEST_PROTOCOL)


def deserialize(people=True, stats=True, id=""):
    peopleData = None
    statsData = None

    if people:
        with open("people%s.dat" % id, "rb") as f:
            peopleData = pickle.load(f)

    if stats:
        with open("stats%s.dat" % id, "rb") as f:
            statsData = pickle.load(f)

    return peopleData, statsData


def plot(people, stats):  # experimental plotting function, currently need a rework using subplots
    numplots = len(stats.getStats())
    colormap = plt.cm.gist_ncar
    plt.gca().set_color_cycle([colormap(i) for i in np.linspace(0, 0.9, numplots)])

    labels = []
    for stat, label in stats.getStats():
        plt.plot(stat)
        labels.append(label)

    plt.grid(True)

    plt.legend(labels, ncol=2, loc='upper center', bbox_to_anchor=[0.5, 1.1], columnspacing=1.0,
               labelspacing=0.0, handletextpad=0.0, handlelength=1.5, fancybox=True, shadow=True)

    # plt.yscale('log')
    plt.show()

    # Helpful diagnostic histogram for various stats

    # bins = np.arange(0, 100+1, 1)
    # plt.hist([person.maxAge for person in people], bins=bins)
    # plt.show()


if __name__ == '__main__':
    # deserializing usually happens here, but I have excluded that portion here
    main()

Я заметил, что мне нужен способ, чтобы отслеживать статистику с течением времени, так что я мог видеть это в виде графиков, а не просто напечатанные статистика. Я написал вспомогательный класс Stats для этого:


# stats.py

class Stats(object):
    def __init__(self, oldDate, startingPopulation, statYear):
        self.timeTaken = []
        self.population = []
        self.meanAge = []
        self.newPeople = []
        self.deaths = []
        self.deltaNewPeople = []
        self.surnames = []
        self.deltaDeaths = []
        self.netPopDelta = []
        self.startingPopulation = startingPopulation
        self.year = 0
        self.statYear = statYear
        self.oldDate = oldDate
        self.newPersonCount = 0
        self.oldDeaths = 0
        self.oldPeople = 0

    def getStats(self):  # convenience function to get all plot-able stats as an iterable with their labels attached
        return [
            (self.timeTaken, "Time Taken"),
            (self.population, "Pop Count"),
            (self.meanAge, "Mean Age"),
            (self.newPeople, "New People"),
            (self.deaths, "Deaths"),
            (self.deltaNewPeople, "Delta New People"),
            (self.deltaDeaths, "Delta Deaths"),
            (self.surnames, "Surnames"),
            (self.netPopDelta, "Net Pop Change")
        ]

После много проб и ошибок и диагностический участки, я чувствую, что я, наконец, получил это работает прилично. Ни в коем случае не хочу утверждать, что результаты являются статистически точный для любого реального населения планеты, но за то, что я хотел исследовать, результаты выглядят достаточно хорошо, особенно с Атлантического отмечает, что для Китая, с 1,3 млрд. человек, 87% населения разделяет только 100 уникальных фамилий.

Для страны с населением 1,3 миллиарда человек, есть удивительно небольшое число распространенные фамилии в Китае. Около 87 процентов акций населения один из 100 фамилий, и более чем один из пяти китайских граждан прозванного Ли, Ван или Чжан: более 275 миллионов человек.

Образец выполнения было сделано, имитируя 500 лет, печатая статистику каждые 100 лет:


2018-03-08 15:59:02.250965
Time taken: 0 seconds (0 milliseconds / year)
Current population size at year 0: 10000
Current mean age: 0
Current new people: 0
Current deaths: 0
Current delta new people: 0 (0 people / year)
Current delta deaths: 0 (0 deaths / year)
Current surname count: 10000
Current partner percentage: 0
Net population per year: 0

2018-03-08 15:59:11.135088
Time taken: 8 seconds (88 milliseconds / year)
Current population size at year 100: 26898
Current mean age: 25
Current new people: 48750
Current deaths: 31852
Current delta new people: 48750 (487 people / year)
Current delta deaths: 31852 (318 deaths / year)
Current surname count: 2044
Current partner percentage: 51
Net population per year: 168

2018-03-08 15:59:24.385102
Time taken: 13 seconds (132 milliseconds / year)
Current population size at year 200: 32645
Current mean age: 25
Current new people: 113218
Current deaths: 90573
Current delta new people: 64468 (644 people / year)
Current delta deaths: 58721 (587 deaths / year)
Current surname count: 1085
Current partner percentage: 52
Net population per year: 57

2018-03-08 15:59:41.181756
Time taken: 16 seconds (167 milliseconds / year)
Current population size at year 300: 39817
Current mean age: 25
Current new people: 191768
Current deaths: 161951
Current delta new people: 78550 (785 people / year)
Current delta deaths: 71378 (713 deaths / year)
Current surname count: 806
Current partner percentage: 53
Net population per year: 71

2018-03-08 16:00:02.352846
Time taken: 21 seconds (211 milliseconds / year)
Current population size at year 400: 47657
Current mean age: 25
Current new people: 286314
Current deaths: 248657
Current delta new people: 94546 (945 people / year)
Current delta deaths: 86706 (867 deaths / year)
Current surname count: 659
Current partner percentage: 50
Net population per year: 78

2018-03-08 16:00:28.200274
Time taken: 25 seconds (257 milliseconds / year)
Current population size at year 500: 58446
Current mean age: 25
Current new people: 400610
Current deaths: 352164
Current delta new people: 114296 (1142 people / year)
Current delta deaths: 103507 (1035 deaths / year)
Current surname count: 578
Current partner percentage: 53
Net population per year: 107

Теперь почему я здесь на код-ревью в первую очередь с этим чудовищем. Я чувствую, что сократили некоторые из основных проблем с производительностью, но свежим взглядом вижу то, что я скучаю. Я тоже ищу, чтобы сделать вещи более подходящие для Python и читаемым. Одна вещь, я забыл упомянуть, хотя оно появляется в коде, заключается в том, что я хочу иметь более надежные построение схемы, вероятно, с использованием сюжетных линий. Я еще не дошла до того, чтобы это произошло, так что если вы видите лучший способ, пожалуйста, поделитесь!

Я не претендую на гуру питона, но сайте StackOverflow проходит долгий путь, чтобы помочь выучить новый язык. Не стесняйтесь обсуждать то, что вы находите странным, и если что-то непонятно, что он делает, не стесняйтесь спрашивать.

С точки зрения читабельности, я знаю, что одна из вещей PEP8 указывает snake_case но я изначально разработчик и старые привычки умирают с трудом (и я ненавижу набрав подчеркивания). Я также добавил кучу комментариев, объясняющих мое обоснование, почему вещи появляются, как они делают, что я чувствую десны вещи немного, но я понимаю, что все это не в моей голове.


Некоторые побочные заметки про тайминги ниже:

$ python -m timeit -s "import numpy as np" "np.random.randint(0, 1000)"
10000000 loops, best of 3: 0.164 usec per loop
$ python -m timeit -s "import random" "random.randint(0, 1000)"
1000000 loops, best of 3: 1.52 usec per loop
$ python -m timeit -s "import random" "random.uniform(0, 1)"
1000000 loops, best of 3: 0.489 usec per loop
$ python -m timeit -s "import numpy as np" "np.random.uniform(0, 1)"
1000000 loops, best of 3: 0.219 usec per loop

timeit.timeit("rand(0, 1000)", setup="from numpy.random import uniform as rand")
0.1250929832458496
timeit.timeit("rand(0, 1000)", setup="from random import uniform as rand")
0.47645998001098633


Комментарии
2 ответа


С точки зрения читабельности, я знаю, что одна из вещей PEP8
указывает snake_case

Если вы знаете, PEP8, игнорируя именования, есть еще несколько моментов, чтобы рассмотреть, из руководства, такие как:


  1. import заказ модулей

  2. Писать комментарии

  3. Имея комментарии (я бы сильно рекомендуем это)


Во многих местах, я заметил, что у вас сниппеты как

if self.gender == female:
[x for x in people if x.gender == person.male and x.isPartnerEligible()]

и т. д., несмотря на свой self и x имея в наличии следующее определение:

def isMale(self):
return self.gender == male


Вместо определения male и female в качестве глобальных чисел, использовать Enum которая была портирована для старых версий Python. Или, вы можете просто True/False значение (так как это логическое в своем виртуальном мире). Это помогает с помощниками, а также:

def isMale(self):
return self.gender is MALE # MALE might be True or False

def isFemale(self):
return not self.isMale()


Имена переменных, как haveChild, gainChild немного смущают (может быть только мне), использовать более привычную терминологию; например. isPregnant, giveBirth и т. д.

4
ответ дан 9 марта 2018 в 11:03 Источник Поделиться

Одна вещь, которая намного ускорив дело в том, чтобы не использовать numpy без массивов. В частности, при использовании nprandintвы используете его для генерации одного случайного числа. Вместо этого, вы должны использовать random.randint, что примерно на 10% быстрее. Аналогичным образом, для равномерного random в 3 раза быстрее uniformи в 5 раз быстрее exponential.

Я получил эти тайминги следующим образом:

timeit.timeit('rand(0,1000)', setup='from numpy.random import uniform as rand')
0.7485988769913092
timeit.timeit('rand(0,1000)', setup='from random import uniform as rand')
0.2676771300029941

0
ответ дан 9 марта 2018 в 01:03 Источник Поделиться