在MIT的研究中，多次需要使用合成人口开展城市模拟，主要工具为python中的第三方模块synthpop。合成人口在国内的应用相对较少，synthpop更是非常缺少帮助文档，因此在这里记录一下它的具体算法和使用方式。
为什么要合成人口？
我们在做多代理人模拟时，经常需要个体级别的被模拟对象，例如一个一个的人。每个人都自带一系列的个体属性，如性别、年龄、收入等，这些属性毫无疑问会对模拟行为产生影响，因此不能随便乱设。然而，实际工作中我们很少有机会获得全部的个体数据，一般只有经过汇总统计后的普查数据，例如：已知某地总人口1000人，其中男女比例为6：4，低、中、高三组收入的比例为3：5：2。那么，合成人口的目的是要生成1000条带有性别、年龄等属性的个体记录，使得他们的属性汇总后与普查结果保持一致。
当然，上述过程并不是说要凭空制造出符合条件的1000条记录，而是基于另一套个体级别的微观数据抽样得到。这种微观数据就是以每个人为单位，详细列出了他/她的各项个体属性和其他指标（如意愿、行为等等），它们都是真实可靠的。在国内，这种微观数据的可获得性很差，绝大多数情况下需要自己去做抽样调查；在美国情况要好得多，有官方的数据源背书，叫PUMS（Public Use Microdata Samples），覆盖全美各地和各类人群，个体信息丰富，且随时间滚动更新（神奇吧……，羡慕……）。不过，不管是自己调查还是用官方数据，都必然只是抽样，与全人口相比抽样率小得可怜，而且空间范围不一致，个体属性分布往往有偏。所以，除非假装普查数据不存在，否则这样的微观数据是不宜直接用于模拟的。人口合成的过程，就是对这些微观数据样本进行重抽样，在保持数据真实性的同时，使最终样本的属性分布与普查结果保持一致。
听起来简单吗？如果人口属性很少，那还OK。比如说，普查数据告诉我们要有100个男女比为6:4的人群，而我们的微观数据有10条记录，2男8女，那么把所有男性样本复制30倍，所有女性样本复制5倍，不就有60个男的，40个女的了吗。可是，性别比虽然满足了，收入比、教育比等其他属性可能就乱套了，它们之间是相互制约的，得让所有属性的综合一致性达到最佳。
更麻烦的是，美国的微观数据PUMS中同时包含了家庭数据集和个人数据集，各自有家庭属性和个人属性，两个数据集之间也自然是关联的，一个家庭通常对应多个个人。于是，人口合成的任务就变成：如何在微观家庭样本中进行重抽样，使得合成后的家庭数据在各项家庭属性上的分布与普查数据一致，同时，使得与家庭对应的合成个人数据在各项个体属性上的分布也与普查数据一致。是不是麻烦多了？
于是我们需要用synthpop了。这个算法和工具包是美国亚利桑那大学的研究者提出和开发的，相关论文在这里。下面用一个真实的例子来说明它是怎么工作的。
数据
如上所述，我们需要2种数据：普查数据和微观数据。这里的普查数据为美国ACS（American Community Survey）数据，微观数据为美国PUMS（Public Use Microdata Samples）数据。每种数据又各分为家庭表和个人表，所以下面会有4张表，每张表都是一个pandas.DataFrame对象。需要指出：虽然美国数据的开放性好得多，但把最原始的可下载数据处理成下面4张表也得一番折腾，具体过程参见这篇文章。
普查家庭表
该表的变量名为h_acs_cat，其中，h代表household，acs代表数据源，cat代表category，即已将属性分类化。打印其前3行数据，如下所示。
| cat_name | cars | children | income | tenure | workers | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| cat_value | none | one | two or more | no | yes | gt100 | gt35-lt100 | lt35 | owned | rented | none | one | two or more |
| NAME | |||||||||||||
| Block Group 2, Census Tract 5172, Wayne County, Michigan | 51 | 94 | 20 | 158 | 10 | 52 | 100 | 16 | 168 | 0 | 28 | 106 | 33 |
| Block Group 1, Census Tract 5172, Wayne County, Michigan | 202 | 370 | 81 | 582 | 73 | 129 | 218 | 308 | 648 | 7 | 109 | 415 | 129 |
| Block Group 1, Census Tract 5207, Wayne County, Michigan | 319 | 278 | 125 | 657 | 67 | 55 | 179 | 490 | 724 | 0 | 275 | 336 | 111 |
这个数据来自美国密歇根州（Michigan）的Wayne县，县下面再分Census Tract，然后再下分Block Group。可以看到，DataFrame的行索引为NAME，具体为Block Group级别的一个一个地名。DataFrame的列为双重索引，外层索引为家庭属性的类名（cat_name），这里考虑5种属性：是否保有小汽车（cars），是否有小孩（children），家庭年总收入（income），住房类型（tenure），家庭中的劳动力数量（workers）。相应地，内层索引为属性的类别取值（cat_value），例如，住房类型（tenure）属性分括购买拥有房（owned）和租住房（rented）两种。于是，表中的每一行反映的是一个Block Group中各类家庭的数量。例如，第三行数据是说：“Block Group1, Census Tract 5207, Wayne County, Michigan”这个地方有319个家庭没车，278个家庭有1辆车，125个家庭有2辆及以上车；657个家庭没有孩子，67个家庭有孩子；55个家庭年收入在100万美元以上，179个家庭在35-100万美元之间，490个家庭年收入不足35万美元；724个家庭住在自己买的房子里，而租房住的家庭数量为0；275个家庭没有劳动力，336个家庭只有1个劳动力，111个家庭有2个及以上劳动力。
（PS：按不同属性求和得到的家庭总数是不完全一致的，例如，按小汽车保有情况计算为：319+278+125=722，而按有无孩子计算为：657+67=724。这种差异是由于缺失值和数据清洗程序造成的，没有关系。）
下面的示例中将只关注这第三行的数据，即要合成在“Block Group1, Census Tract 5207, Wayne County, Michigan”这个地方的人口，使家庭属性满足cars、children、income、tenure、workers的已知统计分布。
普查个人表
该表的变量名为p_acs_cat，其中，p代表person，acs代表数据源，cat代表category，即已将属性分类化。打印其前3行数据，如下所示。
| cat_name | age | race | sex | |||||||
|---|---|---|---|---|---|---|---|---|---|---|
| cat_value | 19 and under | 20 to 35 | 35 to 60 | above 60 | asian | black | other | white | female | male |
| NAME | ||||||||||
| Block Group 2, Census Tract 5172, Wayne County, Michigan | 170 | 884 | 488 | 38 | 33 | 963 | 98 | 486 | 232 | 1348 |
| Block Group 1, Census Tract 5172, Wayne County, Michigan | 7 | 417 | 251 | 142 | 111 | 305 | 30 | 371 | 410 | 407 |
| Block Group 1, Census Tract 5207, Wayne County, Michigan | 41 | 540 | 418 | 44 | 46 | 442 | 24 | 531 | 328 | 715 |
这个DataFrame与上面的h_acs_cat完全相同，只是从家庭属性换成了3个个人属性：年龄（age）、种族（race）、性别（sex）。对于拟合成人口的地方——“Block Group1, Census Tract 5207, Wayne County, Michigan”，这里按年龄分，有41位19岁以下的，540位20-35岁的，418位35-60岁的，以及44位60岁以上的；按种族分，有46位亚裔、442位黑人、531位白人、24位其他种族；按性别分，有715个男性，328个女性。我们合成的人口既要在家庭层面满足h_acs_cat中的家庭属性分布，也要在个人层面满足这些个人属性分布。
微观家庭表
该表的名称为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，而且行索引号不连续，因为经过了空间筛选。
| serialno | puma10 | cars | children | income | tenure | workers | |
|---|---|---|---|---|---|---|---|
| 0 | 2012000000359 | 3201 | none | no | lt35 | rented | none |
| 97 | 2012000223076 | 3201 | none | no | lt35 | rented | none |
| 117 | 2012000277533 | 3201 | none | no | lt35 | rented | none |
这样的微观数据其实包含了相当丰富的属性信息（有很多列），但为了方便，这里就只看与人口合成有关的几个属性——cars, children, income, tenure, workers，属性名称和分类都与2.1节中的普查数据属性相对应。每一行是一个家庭，serialno是唯一的家庭编号，用来与个人表做链接。可以看到，这里列出的3个家庭都没有车，没有孩子，家庭年收入都小于35万美元，都住在租来的房子，而且都没有劳动力。
h_pums表共有1155行，我们将从这1155个家庭中做抽样，以满足上述普查家庭表和普查个人表的属性分布
微观个人表
该表的名称为p_pums，其中，p代表person，pums代表数据源。这个表总共有2963行，每一行代表一个人。打印前3行，如下所示。
| serialno | puma10 | age | race | sex | |
|---|---|---|---|---|---|
| 0 | 2012000000359 | 3201 | above 60 | white | male |
| 3 | 2012000003126 | 3201 | above 60 | white | male |
| 6 | 2012000005738 | 3201 | above 60 | white | male |
个体属性很多，这里只显示了合成人口所用的3个属性——age, race, sex，它们的分类取值与2.2节的普查数据相对应。上图显示的3个人都是60岁以上的白人男性。此外，serialno用于链接到h_pums表，找到每个人所在的家庭。抽样并不发生在这个个人表上，而是发生在上面的家庭表上，家庭样本确定了，个人样本自然也就出来了。
于是，我们要做的人口合成任务可描述如下：
对于“Block Group1, Census Tract 5207, Wayne County, Michigan”这个空间单元，我们已经从普查数据中知道它在多个家庭和个体属性上的构成特征（h_acs_cat与p_acs_cat的第1行），同时，我们手头上还有与它空间上接近的“Michigan州3201号PUMA”的微观数据，包括1155个家庭和相应的2963个人（h_pums，p_pums），现在，我们要从这些家庭和个人中进行抽样，使抽出来的家庭和个人在总量及属性结构上尽可能接爱普查数据。由此，我们下一步喂给模拟系统的是包含大量真实属性信息的微观个体数据，同时它们又在汇总层面具有充分的代表性，经得起普查数据校核。
进行人口合成
这一节将说明如何通过synthpop的2~3个命令执行人口合成。如果不关心算法细节的话，直接按上面的数据格式准备好相应的4个DataFrame，然后按下面的命令操作就可以了。
首先当然是要导入模块。其中，categorizer是用来处理类别的辅助工具，synthesizer才是真正的合成工具。
获得联合分布
这一步的任务是获得统计微观数据在各属性之间的联合分布，代码如下，分别处理家庭和个人数据。我们要获得的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，即联合分布，因为它的每一类联合考虑了所有属性的组合。下面，就把这样的交叉细分类统称为联合类，样本在这些联合类上的分布为联合分布。
| cat_id | frequency | |||
|---|---|---|---|---|
| age | race | sex | ||
| 19 and under | asian | female | 0 | 47 |
| male | 1 | 46 | ||
| black | female | 2 | 34 | |
| male | 3 | 37 | ||
| other | female | 4 | 26 | |
| male | 5 | 29 | ||
| white | female | 6 | 240 | |
| male | 7 | 288 | ||
| 20 to 35 | asian | female | 8 | 36 |
| male | 9 | 24 | ||
| black | female | 10 | 23 | |
| male | 11 | 18 | ||
| other | female | 12 | 8 | |
| male | 13 | 7 | ||
| white | female | 14 | 206 | |
| male | 15 | 192 | ||
| 35 to 60 | asian | female | 16 | 56 |
| male | 17 | 57 | ||
| black | female | 18 | 58 | |
| male | 19 | 51 | ||
| other | female | 20 | 10 | |
| male | 21 | 15 | ||
| white | female | 22 | 449 | |
| male | 23 | 420 | ||
| above 60 | asian | female | 24 | 18 |
| male | 25 | 13 | ||
| black | female | 26 | 25 | |
| male | 27 | 17 | ||
| other | female | 28 | 5 | |
| male | 29 | 4 | ||
| white | female | 30 | 277 | |
| male | 31 | 227 |
然后，第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()函数之后，它才出现的。
| serialno | puma10 | age | race | sex | cat_id | |
|---|---|---|---|---|---|---|
| 0 | 2012000000359 | 3201 | above 60 | white | male | 31 |
| 3 | 2012000003126 | 3201 | above 60 | white | male | 31 |
| 6 | 2012000005738 | 3201 | above 60 | white | male | 31 |
家庭数据的处理类似，只是由于家庭属性较多，总共产生了108个联合类，h_jd的情况如下所示。cad_id编号从0到107，frequency统计了每个联合类所拥有的家庭样本数量。同时，h_pums中增加了一列cat_id，记录了每个家庭样本的联合类编码（0-107）。
| cat_id | frequency | |||||
|---|---|---|---|---|---|---|
| cars | children | income | tenure | workers | ||
| none | no | gt100 | owned | none | 0 | 1 |
| one | 1 | 0 | ||||
| two or more | 2 | 0 | ||||
| rented | none | 3 | 0 | |||
| one | 4 | 0 | ||||
| ... | ... | ... | ... | ... | ... | ... |
| two or more | yes | lt35 | owned | one | 103 | 7 |
| two or more | 104 | 5 | ||||
| rented | none | 105 | 0 | |||
| one | 106 | 3 | ||||
| two or more | 107 | 1 |
获得边缘分布
这一步的任务是获得普查数据的边缘分布，代码如下。
这一步其实非常简单，因为数据已经有了，只需要从普查家庭表和普查个人表中把对应的那一行切出来。我们要做人口中合成的空间单元是“Block Group1, Census Tract 5207, Wayne County, Michigan”，在普查表中的第3行，因此索引行号是2。把这一行切出来后做个转置，命名为h_marg和p_marg，其中，h_和p_代表household和person，后缀_marg代表marginal，即为边缘分布。它们的数据类型都是pandas.Series，打印出来长下面这个样子。
cat_name cat_value
cars none 319
one 278
two or more 125
children no 657
yes 67
income gt100 55
gt35-lt100 179
lt35 490
tenure owned 724
rented 0
workers none 275
one 336
two or more 111
Name: Block Group 1, Census Tract 5207,
Wayne County, Michigan, dtype: int32cat_name cat_value
age 19 and under 41
20 to 35 540
35 to 60 418
above 60 44
race asian 46
black 442
other 24
white 531
sex female 328
male 715
Name: Block Group 1, Census Tract 5207,
Wayne County, Michigan, dtype: int32学过统计学的肯定知道为什么叫边缘分布，这里再多解释一下。下表显示了性别和年龄的交叉表，2行代表性别的2分类，4列代表年龄的4分类。白色单元格区域即为上面提到的联合类，它们的分布即为联合分布，只是这里用了更直观的交叉表方式去表示。不过，这种表示方式只能用于二个维度的属性，更高维度就办不到了，所以还是得用上一节的表示方法。联合分布在普查数据中通常是未知的，而在微观数据中是可以统计获得的（上一节）。红色单元格区域是对行、列的求和，放在表的右边缘或下边缘，所以被形象地称为边缘分布。下表中，按行求和即得到性别分布（不管什么年龄组），按列求和即得到年龄分布（不管什么性别组）。普查数据里通常提供的正是这样的边缘分布。
| 19 and under | 20 to 35 | 35 to 60 | above 60 | Row SUM | |
|---|---|---|---|---|---|
| female | unknown | unknown | unknown | unknown | 328 |
| male | unknown | unknown | unknown | unknow | 715 |
| Column SUM | 41 | 540 | 418 | 44 |
执行synthesize()
在执行人口合成之前，再看看现在我们有哪些变量：
- 普查数据的边缘分布：h_marg、p_marg
- 微观数据的联合分布：h_jd、p_jd
- 微观数据的个体样本池：h_pums、p_pums
于是，人口合人的任务可再一次概括如下：
对于已知不同属性联合分布（_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。
| serialno | puma10 | cars | children | income | tenure | workers | cat_id | |
|---|---|---|---|---|---|---|---|---|
| 0 | 2013001372862 | 3201 | none | no | gt100 | owned | none | 0 |
| 1 | 2012000256072 | 3201 | none | no | gt35-lt100 | owned | none | 6 |
| 2 | 2012000256072 | 3201 | none | no | gt35-lt100 | owned | none | 6 |
| 3 | 2012000256072 | 3201 | none | no | gt35-lt100 | owned | none | 6 |
| 4 | 2012001223810 | 3201 | none | no | gt35-lt100 | owned | two or more | 8 |
| 5 | 2012001223810 | 3201 | none | no | gt35-lt100 | owned | two or more | 8 |
| 6 | 2012001260324 | 3201 | none | no | gt35-lt100 | owned | two or more | 8 |
| 7 | 2012001260324 | 3201 | none | no | gt35-lt100 | owned | two or more | 8 |
| 8 | 2012001260324 | 3201 | none | no | gt35-lt100 | owned | two or more | 8 |
| 9 | 2012001260324 | 3201 | none | no | gt35-lt100 | owned | two or more | 8 |
有了家庭样本，人口样本自然就导出来了。这里得到的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也能进行印证。
| serialno | puma10 | age | race | sex | cat_id | hh_id | |
|---|---|---|---|---|---|---|---|
| 0 | 2012000256072 | 3201 | above 60 | white | male | 139 | 1 |
| 1 | 2012000256072 | 3201 | above 60 | white | male | 139 | 2 |
| 2 | 2012000256072 | 3201 | above 60 | white | male | 139 | 3 |
| 3 | 2012000256072 | 3201 | above 60 | white | male | 139 | 713 |
| 4 | 2012001223810 | 3201 | above 60 | white | male | 139 | 4 |
| 5 | 2012001223810 | 3201 | above 60 | white | male | 139 | 5 |
| 6 | 2012001223810 | 3201 | above 60 | white | male | 139 | 13 |
| 7 | 2012001223810 | 3201 | above 60 | white | male | 139 | 15 |
| 8 | 2012001223810 | 3201 | above 60 | white | male | 139 | 20 |
| 9 | 2012001223810 | 3201 | above 60 | white | male | 139 | 21 |
synthesize()函数返回的后2个变量——people_chisq和people_p，反映的是合成人口的质量，下面会提到它们。
校验
best_households与best_people中的每个家庭、每个人都来自于h_pums和p_pums中的真实记录，那么问题来了：best_households与best_people在各家庭、个人属性上的边缘分布，是否与我们的目标——普查数据的边缘分布相一致呢？下面就做一个校验。通过如下代码可以生成best_households和best_people的边缘分布。
先比较一下家庭属性：下图中，左表为h_marg反映的基于普查数据的实际边缘分布，右表为best_households_marg反映的合成人口的边缘分布。可以看到，两者之间高度接近，差异微乎其微。
cat_name cat_value
cars none 319
one 278
two or more 125
children no 657
yes 67
income gt100 55
gt35-lt100 179
lt35 490
tenure owned 724
rented 0
workers none 275
one 336
two or more 111
cat_name cat_value
cars none 319
one 278
two or more 126
children no 655
yes 68
income gt100 56
gt35-lt100 180
lt35 487
tenure owned 723
rented 0
workers none 275
one 337
two or more 111
下面再比较一下个人属性：左表为p_marg反映的基于普查数据的实际边缘分布，右表为best_people_marg反映的合成人口的边缘分布。这个差异就大得多了。一方面，在总量上，合成人口总量明显高于实际人口；另一方面，在结构上，合成人口拥有更高比重的19岁以下人群和60岁以上人群，白人比重也偏高。
cat_name cat_value
age 19 and under 41
20 to 35 540
35 to 60 418
above 60 44
race asian 46
black 442
other 24
white 531
sex female 328
male 715
cat_name cat_value
age 19 and under 141
20 to 35 395
35 to 60 509
above 60 350
race asian 62
black 308
other 43
white 982
sex female 587
male 808
为什么会出现这种情况？简单来说就是抽样是直接作用在家庭层面上，所以会优先保证家庭层面的分布拟合。这种家庭优先的方式貌似也是人口合成领域的传统了。个人样本是随着家庭样本的抽出而自然确定的，是第二性的，分布拟合自然差了些，虽然如此，也已经比其他抽样方式要强得多。说到底，根本原因还是因为普查数据和微观数据是不同的数据源，在空间范围和样本选择上都有差异。理论上说，如果存在一种甚至多种抽样方案，能使家庭和个人的边缘分布都完美拟合，那么synthsize()经过足够的迭代一定能找到这些方案。遗憾的是，这种情况在实际中极少出现，那么我们就让算法把抽样做的尽可能好，然后接受不可避免的差异了。举2个例子：（1）如果普查数据上说买房和租房的家庭一半一半，而微观数据里只有买房的家庭，那么即使在家庭层面也不可能取得一致；（2）如果普查数据上说这里有100个家庭，300人，而微观数据里每个家庭都只有1个人，那么不管怎么抽这100个家庭，都只能出来100人。
如上所述，由于家庭属性的边缘分布拟合通常问题不大，因此个人属性的边缘拟合效果更值得关注。统计中常用卡方检验来衡量两个分布之间是否存在显著差异，那么正好用来比较普查数据和合成人口数据在个人属性上的边缘分布差异。synthesize()函数返回的后2个结果——people_chisq和people_p，正是个人属性卡方检验的卡方值和p值。卡方值越小，p值越不显著（>0.05），情况越理想。不过，由于数据源差异的客观存在，也不需要太在意这个这个问题。在普查数据不容质疑的前提下，人口合成的质量差，那只能说明我们的微观数据样本质量不高，对普查数据的代表性差，那恐怕只能重新去做微观抽样调查了。
更多空间单元
上面以一个Block Group为例，演示了用就近PUMA中的微观样本进行人口合成。实际中，我们的研究范围可能会覆盖多个PUMA，而每个PUMA又包含了多个Block Group，所以我们需要下面的循环结构。
技术细节
这一节细说一下synthsize()函数是怎样实现人口合成的。
概述：2个关系
synthsize()函数里的工作流其实并不复杂，主要是建立2个关系。
1. 通过constraints建立从“微观数据联合分布”到“普查数据边缘分布”的关系。
回到上面出现的那张性别*年龄的联合分布表，如下所示：红色单元格的数字是普查数据的边缘分布，白色单元格的数字是微观数据的联合分布。这两个分布肯定是对不上的，例如，对于19岁以下人口，微观数据共计有347+400=747个样本，而普查数据提示只需要41个样本。所以，在联合分布的半单元格中加了“→ ?”，表示需要调整。Constraints就是这一个一个的“？”，即联合分布中的每个细分类需要多少样本，才能符合边缘分布的要求。显而易见，Constraints的解不是唯一的。下表中，就有4*2=8个Constraints，而它们只需要满足4+2=6个边缘分布条件即可。但是，我们的Constraints是由微观数据的实际频率迭代导出的，这个算法被称为IPF（iterative proportional fitting)，因此可以使Constraints与实际频率相对接近。
| 19 and under | 20 to 35 | 35 to 60 | above 60 | Row SUM | |
|---|---|---|---|---|---|
| female | 347 → ? | 273 → ? | 573 → ? | 325 → ? | 328 |
| male | 400 → ？ | 241 → ？ | 543 → ？ | 261 → ？ | 715 |
| Column SUM | 41 | 540 | 418 | 44 |
2. 通过household weights建立从“微观数据个体样本”到“微观数据”的关系。
有了Constraints，是不是就可以直接从微观数据样本池中抽样了呢？是也不是。
如果只抽家庭样本，不关心个人；或者类似地，只抽个人，不关心家庭，那么确实是可以抽样了。举个例子，考虑“有车/无车”与“买房/租房”这2个家庭属性的联合分布，如果Contraints告诉我们：需要80个“无车+租房”的样本（实际样本池中有M个），同时需要有20个“有车+买房”的样本（实际样本池中有N个），那么我们只需要从M里随机抽80个，再从N里随机抽20个即可。或者，更一般的方式是：把那M个“无车+租房”的家庭样本无差别地赋予weight=80的权重，把那N个“有车+买房”的样本无差别地赋予weight=20的权重，然后从这（M+N）个样本中，依权重的大小，抽取（80+20=100）个样本即可。
但是，如果要同时考虑上层的家庭和下层的个人，而且两者之间不是独立的，那么像上面那样直接把Constraints作为抽样权重就不行了。因为很显然，尽管那M个“无车+租房”的家庭样本内部在家庭层面的属性上是无差别的（根据定义，只要有差别就会被分到不同的家庭属性联合类），但是在个人层面上的属性上是可以有差异的。比如，M0与M1都是那M个“无车+租房”家庭中的成员，但是M0家只有一个白人小伙，而M1家则是一对黑人老夫妻，那么如果抽样时M0被抽中，就会提高合成人口中年轻人和白人的比例，而如果M1被抽中，则会提高老年人和黑人的比例。考虑到个人属性也有需要瞄准的普查数据边缘分布，所以我们不应该让这M个“无车+租房”的家庭样本具有相同的权重，而是要设置合理的差别化权重，使之能够同时兼顾各家庭对于家庭层面联合分布和个人层面联合分布的贡献。
这就是我们所需的“household weights”。最终，算法将依照这个weights，从微观家庭样本池中按权重（或者说按概率）抽取家庭和对应的个人。
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=?”中。
| Detailed Classes with Joint Attributes | Updating Constraints for 1 Iteration | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| age | race | sex | starting=frequency | adjust for "19 and under", target=41 | adjust for "20 to 35", target=540 | adjust for "35 to 60", target=418 | adjust for "above 60", target=44 | adjust for "asian", target=46 | adjust for "black", target=442 | adjust for "other", target=24 | adjust for "white", target=531 | adjust for "female", target=328 | adjust for "male", target=715 | |
| 0 | 19 and under | asian | female | 47 | 2.580 | 2.580 | 2.580 | 2.580 | 1.052 | 1.052 | 1.052 | 1.052 | 0.624 | 0.624 |
| 1 | 19 and under | asian | male | 46 | 2.525 | 2.525 | 2.525 | 2.525 | 1.030 | 1.030 | 1.030 | 1.030 | 1.030 | 1.501 |
| 2 | 19 and under | black | female | 34 | 1.866 | 1.866 | 1.866 | 1.866 | 1.866 | 9.069 | 9.069 | 9.069 | 5.383 | 5.383 |
| 3 | 19 and under | black | male | 37 | 2.031 | 2.031 | 2.031 | 2.031 | 2.031 | 9.869 | 9.869 | 9.869 | 9.869 | 14.390 |
| 4 | 19 and under | other | female | 26 | 1.427 | 1.427 | 1.427 | 1.427 | 1.427 | 1.427 | 1.188 | 1.188 | 0.705 | 0.705 |
| 5 | 19 and under | other | male | 29 | 1.592 | 1.592 | 1.592 | 1.592 | 1.592 | 1.592 | 1.326 | 1.326 | 1.326 | 1.933 |
| 6 | 19 and under | white | female | 240 | 13.173 | 13.173 | 13.173 | 13.173 | 13.173 | 13.173 | 13.173 | 8.631 | 5.122 | 5.122 |
| 7 | 19 and under | white | male | 288 | 15.807 | 15.807 | 15.807 | 15.807 | 15.807 | 15.807 | 15.807 | 10.357 | 10.357 | 15.102 |
| 8 | 20 to 35 | asian | female | 36 | 36.000 | 37.821 | 37.821 | 37.821 | 15.425 | 15.425 | 15.425 | 15.425 | 9.155 | 9.155 |
| 9 | 20 to 35 | asian | male | 24 | 24.000 | 25.214 | 25.214 | 25.214 | 10.283 | 10.283 | 10.283 | 10.283 | 10.283 | 14.994 |
| 10 | 20 to 35 | black | female | 23 | 23.000 | 24.163 | 24.163 | 24.163 | 24.163 | 117.429 | 117.429 | 117.429 | 69.696 | 69.696 |
| 11 | 20 to 35 | black | male | 18 | 18.000 | 18.911 | 18.911 | 18.911 | 18.911 | 91.901 | 91.901 | 91.901 | 91.901 | 134.002 |
| 12 | 20 to 35 | other | female | 8 | 8.000 | 8.405 | 8.405 | 8.405 | 8.405 | 8.405 | 7.000 | 7.000 | 4.154 | 4.154 |
| 13 | 20 to 35 | other | male | 7 | 7.000 | 7.354 | 7.354 | 7.354 | 7.354 | 7.354 | 6.125 | 6.125 | 6.125 | 8.931 |
| 14 | 20 to 35 | white | female | 206 | 206.000 | 216.420 | 216.420 | 216.420 | 216.420 | 216.420 | 216.420 | 141.798 | 84.159 | 84.159 |
| 15 | 20 to 35 | white | male | 192 | 192.000 | 201.712 | 201.712 | 201.712 | 201.712 | 201.712 | 201.712 | 132.162 | 132.162 | 192.707 |
| 16 | 35 to 60 | asian | female | 56 | 56.000 | 56.000 | 20.975 | 20.975 | 8.554 | 8.554 | 8.554 | 8.554 | 5.077 | 5.077 |
| 17 | 35 to 60 | asian | male | 57 | 57.000 | 57.000 | 21.349 | 21.349 | 8.707 | 8.707 | 8.707 | 8.707 | 8.707 | 12.696 |
| 18 | 35 to 60 | black | female | 58 | 58.000 | 58.000 | 21.724 | 21.724 | 21.724 | 105.574 | 105.574 | 105.574 | 62.660 | 62.660 |
| 19 | 35 to 60 | black | male | 51 | 51.000 | 51.000 | 19.102 | 19.102 | 19.102 | 92.832 | 92.832 | 92.832 | 92.832 | 135.360 |
| 20 | 35 to 60 | other | female | 10 | 10.000 | 10.000 | 3.746 | 3.746 | 3.746 | 3.746 | 3.119 | 3.119 | 1.851 | 1.851 |
| 21 | 35 to 60 | other | male | 15 | 15.000 | 15.000 | 5.618 | 5.618 | 5.618 | 5.618 | 4.679 | 4.679 | 4.679 | 6.823 |
| 22 | 35 to 60 | white | female | 449 | 449.000 | 449.000 | 168.174 | 168.174 | 168.174 | 168.174 | 168.174 | 110.187 | 65.398 | 65.398 |
| 23 | 35 to 60 | white | male | 420 | 420.000 | 420.000 | 157.312 | 157.312 | 157.312 | 157.312 | 157.312 | 103.071 | 103.071 | 150.289 |
| 24 | above 60 | asian | female | 18 | 18.000 | 18.000 | 18.000 | 1.352 | 0.551 | 0.551 | 0.551 | 0.551 | 0.327 | 0.327 |
| 25 | above 60 | asian | male | 13 | 13.000 | 13.000 | 13.000 | 0.976 | 0.398 | 0.398 | 0.398 | 0.398 | 0.398 | 0.580 |
| 26 | above 60 | black | female | 25 | 25.000 | 25.000 | 25.000 | 1.877 | 1.877 | 9.122 | 9.122 | 9.122 | 5.414 | 5.414 |
| 27 | above 60 | black | male | 17 | 17.000 | 17.000 | 17.000 | 1.276 | 1.276 | 6.203 | 6.203 | 6.203 | 6.203 | 9.045 |
| 28 | above 60 | other | female | 5 | 5.000 | 5.000 | 5.000 | 0.375 | 0.375 | 0.375 | 0.313 | 0.313 | 0.186 | 0.186 |
| 29 | above 60 | other | male | 4 | 4.000 | 4.000 | 4.000 | 0.300 | 0.300 | 0.300 | 0.250 | 0.250 | 0.250 | 0.365 |
| 30 | above 60 | white | female | 277 | 277.000 | 277.000 | 277.000 | 20.799 | 20.799 | 20.799 | 20.799 | 13.627 | 8.088 | 8.088 |
| 31 | above 60 | white | male | 227 | 227.000 | 227.000 | 227.000 | 17.044 | 17.044 | 17.044 | 17.044 | 11.167 | 11.167 | 16.283 |
迭代过程很简单：每次处理一个边缘分布，在联合分布中找到与该边缘分布相关的行，将其当前的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如下所示。
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 | 31 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| age | 19 and under | 19 and under | 19 and under | 19 and under | 19 and under | 19 and under | 19 and under | 19 and under | 20 to 35 | 20 to 35 | 20 to 35 | 20 to 35 | 20 to 35 | 20 to 35 | 20 to 35 | 20 to 35 | 35 to 60 | 35 to 60 | 35 to 60 | 35 to 60 | 35 to 60 | 35 to 60 | 35 to 60 | 35 to 60 | above 60 | above 60 | above 60 | above 60 | above 60 | above 60 | above 60 | above 60 |
| race | asian | asian | black | black | other | other | white | white | asian | asian | black | black | other | other | white | white | asian | asian | black | black | other | other | white | white | asian | asian | black | black | other | other | white | white |
| sex | female | male | female | male | female | male | female | male | female | male | female | male | female | male | female | male | female | male | female | male | female | male | female | male | female | male | female | male | female | male | female | male |
| final constraints | 0.582794 | 1.40718 | 5.00093 | 13.4261 | 0.620585 | 1.70765 | 4.60929 | 13.6455 | 9.7389 | 16.0174 | 73.8058 | 142.498 | 4.16589 | 8.99271 | 86.3138 | 198.467 | 4.9096 | 12.3284 | 60.3172 | 130.845 | 1.6876 | 6.24502 | 60.9692 | 140.698 | 0.365124 | 0.650558 | 6.01539 | 10.0913 | 0.195231 | 0.385313 | 8.7027 | 17.5944 |
下面用groupby().sum()验证可知：对于age、race、sex的所有的10个边缘分布，与之相关的联合分布的Constraints的和与边缘分布的目标值均完全一致。如果手算的话，就对照上面那张表中各边缘分布列的白底黑字单元格，求和验证即可。
# final.group('age')
final constraints
age
19 and under 41.0
20 to 35 540.0
35 to 60 418.0
above 60 44.0
# final.group('race')
final constraints
race
asian 46.0
black 442.0
other 24.0
white 531.0
# final.group('sex')
final constraints
sex
female 328.0
male 715.0
家庭属性联合分布的Constraints求法完全相同，就不再演示了。至此，我们就获得了家庭和个人层次上的所有Constraints。家庭属性联合分布共有108个细分小类，所以有108个Constraints值，个人属性联合分布共有32个细分小类，所以有32个Constraints值。这些值告诉我们：对于每个家庭小类和个人小类，我们需要多少样本。当然，很多Constraints都是小数，甚至小于1，这些都没关系，最后抽样时会再处理。
Household Weights
下面再来看看同时考虑家庭和个人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。
| cat_id | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | ... | 99 | 100 | 101 | 102 | 103 | 104 | 105 | 106 | 107 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| hh_id | |||||||||||||||||||
| 0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 1 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 2 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 3 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 4 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 1150 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 1151 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 1152 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 1153 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 1154 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
下面的代码证明：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）的个人，应该是一个青年白人夫妇的两口之家。
| cat_id | 108 | 109 | 110 | 111 | 112 | 113 | 114 | 115 | 116 | 117 | 118 | 119 | 120 | 121 | 122 | 123 | 124 | 125 | 126 | 127 | 128 | 129 | 130 | 131 | 132 | 133 | 134 | 135 | 136 | 137 | 138 | 139 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| hh_id | ||||||||||||||||||||||||||||||||
| 0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 |
| 1 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 2 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 1.0 |
| 3 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 1.0 |
| 4 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 2.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 1150 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 1.0 |
| 1151 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 |
| 1152 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 1153 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 1154 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
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时类似，当一个家庭对多种家庭类/个人类有贡献时，它的权重在一轮更新中会被多次调整。
本例中，处理完全部140列即为完成了一轮更新，算法一般需要多轮更新才能达到收敛，收敛标准仍然是前一轮的weights和后一轮的weights的差异小于一定的阈值（tolerance）。由于也采用了迭代+等比例缩放的方式，该算法被称为IPU(iterative proportional updating)。每更新完一轮，可以计算一个优度指标——先计算每一列贡献值的加权和与constraint的偏离值，再对全部140列求平均。值得注意的是，weights在一轮一轮的迭代中并不能保证该优度指标单调下降，所以需要持续记录一个历史最佳的best_weights，只有当前轮次的weights表现更好，才把它作为新的best_weights。
综上可见，求取household weights的IPU算法与上一节中求取Constraints的IPF算法十分相似，下面对它们进行对照。
IPF：求联合类的constraints
- 需要多轮迭代，直接前后两轮的constraints无显著变化。
- 每一轮（iteration）迭代中，依次（step）循环每一个边缘分布中的大类。
- 每次迭代：
- 以边缘分布的取值为目标。
- 找到属于该边缘分布大类的联合分布细分小类。
- 等比例放缩这些联合类的constraints。
- 放纵后，这些联合类的constraints的和等于边缘分布的取值。
IPU：求家庭样本的weights
- 需要多轮迭代，直接前后两轮的weights无显著变化。
- 每一轮（iteration）迭代中，依次（step）循环每一个联合分布中的细分小类。
- 每次（step）迭代：
- 以联合分布的constraint为目标。
- 找到对该联合分布的细分小类有贡献的家庭样本。
- 等比例放缩这些家庭样本的weights。
- 放缩后，这些家庭样本的贡献值的加权和等于constraint。
利用IPU算法，household weights经过8轮迭代达到收敛。这里的收敛不是指所有的140个家庭类/个人类都实际了贡献值的加权和等于对应的constraint值，而是说weights已经停止更新。事实上，上述目标在实际中基本不可能实现，算法只是找到了它能找到的（可能也是理论上最优的）weights值。另外需要注意的是，weights之所以叫household weights，就是因为它只针对家庭样本，没有person weights这一说。我们在更新household weights时已经考察了个人类的影响，下面的抽样也将面向家庭开展。
抽样
最后一部就是抽样了，基本思路就是以上面所求的household weights为权重，对微观数据家庭样本池（h_pums）中的家庭进行随机抽样。权重越大，被抽到的概率越高。抽取的家庭总数由边缘分布中规定的本地家庭总数确定。此外有2个注意点。
第一，抽样不是面向全部家庭样本，直接按总数抽的，而是分不同的家庭属性联合类，一类一类抽的，比如：先抽“没车+没孩子+月收入小于35万美元+租房+没劳动力”的家庭，再抽“有车+有孩子+月收入小于35万美元+买房+有1个劳动力”的家庭，这样依次抽108类家庭。在抽每一类时，抽取的样本量是该联合类的constraint，备选的家庭样本是属于该类的子样本池，通过前面说的nz（none zero）位置索引变量中去找。
第二，前面也看到了，作为每一类抽样数量的constraints值都是小数，甚至是小于1的小数，怎么抽？答案是将其向下取整，按这个整数抽。例如，constraint=3.5就先抽3个，constraint=0.9就先一个都不抽。这样的话，最终各类的家庭样本总数会小于目标总数，回过头来再补抽一补。前面在向下取整的时候，还会记录取整后的差值，比如3.5取整为3会差0.5，0.9取整为0会差0.9，然后将这些差异按从高往低排列，总数还差几个，就从排列好的前N类中每类随机抽取1个。
第三，由于抽样的随机性，为了得到更好的结果，算法设定要抽20次，每次得到一组家庭样本和对应的个人样本，计算个人样本的边缘分布与实际边缘分布的差异，取卡方值最小（即差异最小）的那次抽样结果作为最终结果输出。
至此，就说完了synthpop这个包在人口合成上的全部技术细节，以后有需要的时候就用起来吧！
Pingback: 美国的ACS和PUMS数据 – 王灿的个人网站