该表的变量名为h_acs_cat，其中，h代表household，acs代表数据源，cat代表category，即已将属性分类化。打印其前3行数据，如下所示。

该表的变量名为p_acs_cat，其中，p代表person，acs代表数据源，cat代表category，即已将属性分类化。打印其前3行数据，如下所示。

该表的名称为h_pums， 其中，h代表household，pums代表数据源。PUMS数据有自己的一套空间分区体系，其基本单元为PUMA（Public Use Microdata Area）。PUMA比普查数据中的Block Group要大得多，这也正常：本来就是抽样率较低的数据，如果按Block Group来抽样，那每个Block Group能分配的样本量就非常少了。在为特定的Block Group做人口合成时，首先要明确该Block Group所对应的PUMA是哪个，然后从这个PUMA中的微观样本池中抽样。虽然由于PUMA较大，必然会抽到该Block Group以外的样本，但已经是空间层面最接近的了。就像在上海，如果我们的普查区是同济大学，而我们收集了五角场和崇明岛两个区域的微观个体样本，那么从前者抽样肯定比后者来得准确。

本例中要合成人口的Block Group是“Block Group1, Census Tract 5207, Wayne County, Michigan”，它所对应的PUMA是Michigan州的3201号PUMA，因此，这里的家庭表和下面的个人表，都是该PUMA小区内的微观样本。打印h_pums的前3行记录，如下所示，可以看到PUMA10=3201，而且行索引号不连续，因为经过了空间筛选。

该表的名称为p_pums，其中，p代表person，pums代表数据源。这个表总共有2963行，每一行代表一个人。打印前3行，如下所示。

于是，我们要做的人口合成任务可描述如下：

对于“Block Group1, Census Tract 5207, Wayne County, Michigan”这个空间单元，我们已经从普查数据中知道它在多个家庭和个体属性上的构成特征（h_acs_cat与p_acs_cat的第1行），同时，我们手头上还有与它空间上接近的“Michigan州3201号PUMA”的微观数据，包括1155个家庭和相应的2963个人（h_pums，p_pums），现在，我们要从这些家庭和个人中进行抽样，使抽出来的家庭和个人在总量及属性结构上尽可能接近普查数据。由此，我们下一步喂给模拟系统的是包含大量真实属性信息的微观个体数据，同时它们又在汇总层面具有充分的代表性，经得起普查数据校核。

这一步的任务是获得统计微观数据在各属性之间的联合分布，代码如下，分别处理家庭和个人数据。我们要获得的2个新变量——h_jd和p_jd，h_和p_分别代表household和person，_jd后缀代表joint distribution，即联合分布。

第2、6行中，cat.category_combinations()函数的功能是获得所有可能的类别组合，它接收的是一个二维索引，就像上面普查表里的列索引一样。例如，在普查个人表中，外层索引是3种个人属性——AGE, RACE, SEX，内层索引是它们的各自分类，其中，age分有4类，race有4类，sex有2类。cat.category_combinations()在获得这个二维索引（p_acs_cat.columns）以后，会把3种个人属性的类别交叉起来，得到4*4*2=32个交叉细分类，并加以编号，由此返回一个32行1列的DataFrame，如下表所示。注意，唯一的1列就是cat_id，目前还没有frequency列，它是后面加上的。cat_id从0编到31，是每一个交叉细分类的索引，例如：cat_id=0就代表了“19岁及以下、亚裔、女性”这一细分类别，cat_id=11就代表了“20-35岁、黑人、男性”这一类别。cat.category_combinations()返回的这个表被命名为h_jd / p_jd。这里的_jd后缀代表了joint distrubtion，即联合分布，因为它的每一类联合考虑了所有属性的组合。下面，就把这样的交叉细分类统称为联合类，样本在这些联合类上的分布为联合分布。

然后，第3、7行的cat.joint_distribution()函数的功能是的基于上面细分出来的联合类，对微观个体数据进行统计。它接收2个参数：第一个参数是微观数据样本表（h_pums和p_pums），第二个参数就是上面的联合类编码（h_jd和p_jd）；它的返回值同样是这2个变量，只是分别做了点修改，各自增加了一列信息。仍然以个人数据为例，该函数读取p_pums中的每一行对应的个人，根据其属性的组合情况，确定他/她在p_jd中的联合类编码，记录在新增的cat_id列中，如下表所示。与原始的p_pums相比，只是增加了cat_id这一列，该表的前3个样本都属于cat_id=31的联合类，即“60岁以上的白人男性”。同时，该函数还统计了每一个联合类所拥有的个人样本的数量，并记录在新增的frequency列中，这就是上表中frequency列的由来——本来不存在，在执行了cat.joint_distribution()函数之后，它才出现的。

家庭数据的处理类似，只是由于家庭属性较多，总共产生了108个联合类，h_jd的情况如下所示。cad_id编号从0到107，frequency统计了每个联合类所拥有的家庭样本数量。同时，h_pums中增加了一列cat_id，记录了每个家庭样本的联合类编码（0-107）。

这一步的任务是获得普查数据的边缘分布，代码如下。

这一步其实非常简单，因为数据已经有了，只需要从普查家庭表和普查个人表中把对应的那一行切出来。我们要做人口中合成的空间单元是“Block Group1, Census Tract 5207, Wayne County, Michigan”，在普查表中的第3行，因此索引行号是2。把这一行切出来后做个转置，命名为h_marg和p_marg，其中，h_和p_代表household和person，后缀_marg代表marginal，即为边缘分布。它们的数据类型都是pandas.Series，打印出来长下面这个样子。

学过统计学的肯定知道为什么叫边缘分布，这里再多解释一下。下表显示了性别和年龄的交叉表，2行代表性别的2分类，4列代表年龄的4分类。白色单元格区域即为上面提到的联合类，它们的分布即为联合分布，只是这里用了更直观的交叉表方式去表示。不过，这种表示方式只能用于二个维度的属性，更高维度就办不到了，所以还是得用上一节的表示方法。联合分布在普查数据中通常是未知的，而在微观数据中是可以统计获得的（上一节）。红色单元格区域是对行、列的求和，放在表的右边缘或下边缘，所以被形象地称为边缘分布。下表中，按行求和即得到性别分布（不管什么年龄组），按列求和即得到年龄分布（不管什么性别组）。普查数据里通常提供的正是这样的边缘分布。

于是，人口合成的任务可再一次概括如下：

对于已知不同属性联合分布（_jd）的家庭及其对应的个人微观样本池（_pums），从中进行适当的重抽样，使结果分别满足各属性预设的边缘分布（_marg）。

下面就是最终执行人口合成的代码，它将调用synthesizer模块中的synthesize()函数。该函数接收的前2个参数是普查数据在家庭和个人级别上的边缘分布，即h_marg和p_marg；后面2个参数是微观数据在家庭和个人级别上的联合分布，即h_jd和p_jd；再后面2个参数是微观数据在家庭和个人级别上的个体样本池，即h_pums和p_pums；再后面的参数都是一些算法和格式的设定，不太重要，用默认值即可。

该函数会输出4个变量，最重要的是前2个——best_households和best_people。它们都是pandas.DataFrame格式，与h_pums和p_pums的形式几乎完全一样。best_households就是抽样后得到的家庭样本集，best_people是抽样后得到的个人样本集。

我们得到的best_households的维度为723*8，表明共抽取得到了723个家庭。打印best_households的前10行，如下所示。通过serialno可以看到，原样本被不同程度地重抽样了。在h_pums中原编号为2013001372862的样本被抽样了1次（第0行），原编号为2012000256072的样本被抽样了3次（第1-3行），原编号为2012001223810的样本被抽样了2将 （第4-5行），原编号为2012001260324的样本被抽样了4次（第6-9行）。由于有重复，所以原来的serialno肯定不能再当作唯一标识的索引了，于是best_housholds中自动重置了索引，新索引采用自然数编号，从0到722。

有了家庭样本，人口样本自然就导出来了。这里得到的best_pepole的维度是1395*7，表明共抽出了1395个人。打印best_pepole的前10行，如下所示。best_people同样重置了索引，采用了0-1394的新索引，serialno反映的是这个人在原先的p_pums数据中的家庭编号，现在肯定也有了重复。最右侧新加的hh_id反映了这个人来自于best_households中的哪个家庭。例如，前4行的个人其实是一个样本重抽4次得到的，他们分别来自best_households中索引为1、2、3、713的家庭，检查serialno也能进行印证。

best_households与best_people中的每个家庭、每个人都来自于h_pums和p_pums中的真实记录，那么问题来了：best_households与best_people在各家庭、个人属性上的边缘分布，是否与我们的目标——普查数据的边缘分布相一致呢？下面就做一个校验。通过如下代码可以生成best_households和best_people的边缘分布。

先比较一下家庭属性：下图中，左表为h_marg反映的基于普查数据的实际边缘分布，右表为best_households_marg反映的合成人口的边缘分布。可以看到，两者之间高度接近，差异微乎其微。

下面再比较一下个人属性：左表为p_marg反映的基于普查数据的实际边缘分布，右表为best_people_marg反映的合成人口的边缘分布。这个差异就大得多了。一方面，在总量上，合成人口总量明显高于实际人口；另一方面，在结构上，合成人口拥有更高比重的19岁以下人群和60岁以上人群，白人比重也偏高。

上面以一个Block Group为例，演示了用就近PUMA中的微观样本进行人口合成。实际中，我们的研究范围可能会覆盖多个PUMA，而每个PUMA又包含了多个Block Group，所以我们需要下面的循环结构。

synthsize()函数里的工作流其实并不复杂，主要是建立2个关系。

1. 通过constraints建立从“微观数据联合分布”到“普查数据边缘分布”的关系。

回到上面出现的那张性别*年龄的联合分布表，如下所示：红色单元格的数字是普查数据的边缘分布，白色单元格的数字是微观数据的联合分布。这两个分布肯定是对不上的，例如，对于19岁以下人口，微观数据共计有347+400=747个样本，而普查数据提示只需要41个样本。所以，在联合分布的半单元格中加了“→ ?”，表示需要调整。Constraints就是这一个一个的“？”，即联合分布中的每个细分类需要多少样本，才能符合边缘分布的要求。显而易见，Constraints的解不是唯一的。下表中，就有4*2=8个Constraints，而它们只需要满足4+2=6个边缘分布条件即可。但是，我们的Constraints是由微观数据的实际频率迭代导出的，这个算法被称为IPF（iterative proportional fitting)，因此可以使Constraints与实际频率相对接近。

synthesize()通过调用ipf.ipf中的calculate_constraints()函数计算Contraints。该函数主要有2个输入参数：一个作为拟合目标的（普查数据的）边缘分布，一个作为出发点的（微观数据的）联合分布，它的输出是与联合分布相匹配的Constraints。

下面回到我们之前的例子：对“Block Group1, Census Tract 5207, Wayne County, Michigan”这个空间单元进行人口合成。我们需要分别计算家庭层面和个人层面各个联合类的Constraints，这里就以联合类较少的个人层面为例。下表中，左侧“Detailed Classes with Joint Attributes”中就是由age（4分类）、race（4分类）、sex（2分类）交叉得到的32（=4*4*2）种联合分布的细分类；右侧“Updating Constraints for 1 Iteration”展示了Constraints的一轮迭代过程，其中，最左边的“starting=frequency”是Constraints的初始化值，采用的就是联合分布的观测频率（参见获得联合分布）。同时，我们有10个（4+4+2）边缘分布的目标值（参见获得边缘分布），它们被显示在表头首行的“target=?”中。

迭代过程很简单：每次处理一个边缘分布，在联合分布中找到与该边缘分布相关的行，将其当前的Constraints值加总，与边缘分布的目标值进行比较，按照差异对Constraints进行等比例缩放，使缩放后的Constraints值之和等于边缘分布的目标值。例如：我们来看第一次调整“adjust for 19 and under, targt=41”，这是说，我们现在处理的边缘分布是age=”19 and under”，已知其目标值为41个人。在联合分布中，与之相关的联合类有8个，它们的类名和当前的Constraints分别是：“19 and under + asian + female (47)”、“19 and under + asian + male (46)”、”19 and under + black + female (34)”、“19 and under + black + male (37)”、“19 and under + other + female (26)”、“19 and under + other + male (29)”、“19 and under + white + female (240)”、”19 and under + white + male (288)”，由于是第一步，所以括号里的Constraints就是每个联合类的观测频率。把这些Constraints加总，其和为47+46+34+37+26+29+240+288=747，表明我们的联合分布中，符合“age = 19 and under”的总共有747人，而我们只需要41人，于是把这8个Constraints同时乘以（41/747=0.05489），所得到的新的Constraints就显示在“adjust for 19 and under, target=41”这一列的白底黑字单元格中，可以自行验证一下，它们的和正好等于41。至于这一列下的黑底白字单元格，它们不属于“age = 19 and under”下面的细分类，所以它们的Constraints在这一步不更新。就这样，算法对10个边缘分布依次处理，每次处理的结果显示在“adjust for ***, target=***”中，被更新Constraints值以白底黑字突出显示，其他Constraints则保持不变。这一过程中，每个联合类的Constraint值都将被多次更新，例如：对于第0行的“19 and under + asian + female”联合类，它的Constraint起始值为47，是该类的观测频率，在处理“age = 19 and under”这一边缘分布时被更新为2.580，在处理“race = asian”这一边缘分布时被更新为1.052，在处理“sex = female”时被更新为0.624。

所有的边缘分布都被处理后，即完成了一轮（iteration）更新。算法一般需要多轮更新才能使Constraints收敛，收敛标准为前后两轮的Constraints之前差异小于一定的阈值（tolerance）。由于采用了这种“不断迭代+每次等比例缩放”的启发式方法，所以被称为IPF（iterative proportional fitting）。本例的个人分布经过5次迭代收敛，最后的Constraits如下所示。

下面用groupby().sum()验证可知：对于age、race、sex的所有的10个边缘分布，与之相关的联合分布的Constraints的和与边缘分布的目标值均完全一致。如果手算的话，就对照上面那张表中各边缘分布列的白底黑字单元格，求和验证即可。

下面再来看看同时考虑家庭和个人Constraints的家庭权重是怎么生成的。synthesize()将调用ipu.ipu中的household_weights()函数计算家庭权重。它的基本思路是这样的：对任意一个家庭样本，考虑它对于联合分布中每个家庭类和每个个人类的贡献，如果一个家庭类或个人类的Constraint较高（本类别样本的需求量多），那么对它贡献大的家庭样本的权重将被调高；反之，如果一个家庭类或个人类的Constraint较低（本类别样本的需求量少），那么对它贡献大的家庭样本的权重将被调低。依次考虑每一个家庭类和个人类的Constraint，不断调整各个家庭的权重，直至算法收敛。

那么怎么衡量家庭样本对于家庭类/个人类的贡献呢？很简单，用频数。synthesize()将会首先利用cat.frequency_tables()函数生成2个频数表：household_freq为家庭类的频数表，person_freq为个人类的频数表。这两个表的行数是一样的，都是微观数据家庭样本的数量；而对于列数，household_freq的列数为家庭属性联合类的数量，person_freq为个人属性联合类的数量。两张频数表中单元格的意义是：行所在的家庭样本中，包含了几个属于列所在的家庭类/个人类的样本。

这样听起来可能有点怪，特别是对于household_freq，岂不是在说：行所在的家庭样本中，包含了属于列所在的家庭类的样本？那还能有几个，如果行所在的家庭样本属于列所在的家庭类，那就是1个；反之，如果行所在的家庭样本不属于列所在的家庭类，那么就是0个。没错，这张household_freq的取值就是0或1，而且由于每个家庭样本必然属于一个且仅能属于一个家庭类（类间互斥），所以每一行中有且仅有一个元素为1，其他都为0。

下表就是本例中的household_freq，其维度为1155行*108列。行索引“hh_id”即为家庭样本的ID，这里为了方便已进行了重编号，使其值落在0~1154之间。列索引”cat_id”的取值为0~107，代表了108个家庭属性联合类，其各自代表的类可参照获得联合分布。这是一个0/1二值表，1就代表了行所在的家庭样本属于列所在的家庭类，或者换个拗口但更通用的说法：行所在的家庭样本中包含了属于列所在的家庭类的1个样本，贡献值为1。

下面的代码证明：household_freq中的单元格取值要么是0，要么是1；同时，每一行的和均为1，即每个家庭样本必须且只能为一个家庭类贡献一个样本。

然后再看个人类的频数表person_freq，如下所示，它会显得更有意义一些。该表的维度为1155行*32列，其中，行索引“hh_id”的意义与上面的household_freq相同，列索引“cat_id”的取值为108~139，它是由原来取值为0~31的cat_id各自加上108构成的，各自代表了一种个人属性联合类，具体对应的类名可参照获得联合分布。这张表里的单元格取值的意义是：行所在的家庭样本中，包含了几个属于列所在的个人类的个人样本。虽然看上去仍然以0为主，但是它可以是0，可以是1，也可以是更大的数，同时，每一行的和等于这个家庭的个人成员总数，当然可以比1大了。例如：对于hh_id=1的家庭样本，它包含了1个cat_id=122（20-35, white, female）的个人，还包含了1个cat_id=123（20-35, white, male）的个人，应该是一个青年白人夫妇的两口之家。

household_weights()函数的主要输入参数有4个：前两个是家庭联合类和个人联合类的频数表（household_freq, person_freq），后两个是家庭联合类和个人联合类的Constraints（household_constraints, person_constraints）。现在，这4个变量我们都已经有了。household_weights()函数将基于这两个频数表对家庭权重做迭代更新，具体的更新机制和逻辑与上面的Constraints迭代更新非常相似，可以参照类比。

首先，household_weights()函数会初始化一个全为1的weights向量，其长度等于微观家庭样本的数量。本例中有1155个家庭样本，那么weights就是1155个1，代表在最初状态下，所有家庭的权重相同。

然后，函数会把person_freq放在household_freq的右侧，将两张频数表合并成一个大频数表。该表的维度为1155行*140列（140=180+32），其意义是：对于1155个家庭样本和140个家庭类/个人类，行所在的家庭样本中为列所在的家庭类/个人类做出了多少贡献。PS：由于要合并，为了避免家庭类的cat_id索引与个人类的cat_id索引出现重复，所以才将个人类的cat_id统一加上108，即家庭类的总数。

接下来，函数将会生成一个名为freq_wrap的迭代器实例，它可以调用内部方法iter_columns()去循环大频数表中140列的每一列，每次循环返回4个变量：cat_id, col, constraint, nz。

• 第一个元素cat_id是本列的索引（列编号），可以帮助我们知道本列是个什么样的家庭类/个人类，不过基本用不到。
• 第二个元素col是本列中所有不为0的元素的取值列表，即找到那些对本列家庭类/个人类的贡献大于0的家庭样本，取出它们的贡献值。可以参照上面Constraints计算表中每一列的白底黑字单元格，把它们的取值作为一个列表。
• 第三个元素constraint是一个标量，是本列的家庭类/个人类所对应的constraint值。
• 第四个元素nz代表none zeros，是本列中所有不为0的元素对应的位置列表，与col相对应，是col各元素的位置索引。可以参照上面Constraints计算表中每一列的白底黑字单元格，把它们的位置索引（行编号）作为一个列表。

下面的代码是household_weights()函数的核心，它利用freq_wrap迭代器依次循环大频数表的每一列（即每一个家庭类/个人类），不断对weights进行更新，更新函数为_update_weights()。具体而言，在每次更新时，利用非零元素位置索引（nz）找到weights向量中不为零的值，将它与相同长度的col做对应元素相乘，翻译一下，就是找到那些对当前家庭类/个人类有贡献的家庭样本，用它们各自的贡献值乘以此时此刻的权重值，然后求和，相当于对贡献值做加权和。用这个加权和与当前家庭类/个人类的constraint做比较，以constraint为目标，求一个比例adj，然后对上面这些有贡献的家庭样本的权重做等比例放缩（同时乘以adj），放缩后的加权重将正好等于constraint。这样说绕吗？那再简单一点：假设当前循环到的列是个人类“above 60, white, female”，该个人类的constraint是100，表示我们希望能抽出100个这样的个人；而我们的家庭样本中，只有3个家庭A、B、C有这种类别的个人，其中A家庭2个，B、C家庭各1个，此时，所有家庭的权重均为初始值，那么在_update_weights()函数中，贡献值colum=[2,1,1]，权重值weights=[1,1,1]，目标值constraint=100，可求得调整系数adj=100 / (2*1 + 1*1 + 1*1)=25，然后将权重调整为原来的25倍，即新的weights=[25,25,25]，表示我们希望这3个家庭能各被重复抽取25次。与计算Constraints时类似，当一个家庭对多种家庭类/个人类有贡献时，它的权重在一轮更新中会被多次调整。

    IPF：求联合类的constraints

• 需要多轮迭代，直接前后两轮的constraints无显著变化。
• 每一轮（iteration）迭代中，依次（step）循环每一个边缘分布中的大类。
• 每次迭代：
  • 以边缘分布的取值为目标。
  • 找到属于该边缘分布大类的联合分布细分小类。
  • 等比例放缩这些联合类的constraints。
  •  放纵后，这些联合类的constraints的和等于边缘分布的取值。

    IPU：求家庭样本的weights

• 需要多轮迭代，直接前后两轮的weights无显著变化。
• 每一轮（iteration）迭代中，依次（step）循环每一个联合分布中的细分小类。
• 每次（step）迭代：
  • 以联合分布的constraint为目标。
  • 找到对该联合分布的细分小类有贡献的家庭样本。
  • 等比例放缩这些家庭样本的weights。
  • 放缩后，这些家庭样本的贡献值的加权和等于constraint。