自我质疑机制下公共物品博弈模型的相变特性∗

杨波 范敏 刘文奇 陈晓松

1)(昆明理工大学数据科学研究中心,昆明 650500)

2)(昆明理工大学理学院,昆明 650500)

3)(中国科学院理论物理研究所,理论物理前沿重点实验室,北京 100190)

4)(中国科学院大学物理科学学院,北京 100049)

自我质疑机制下公共物品博弈模型的相变特性∗

杨波1)2)†范敏1)2)刘文奇1)2)陈晓松3)4)

1)(昆明理工大学数据科学研究中心,昆明 650500)

2)(昆明理工大学理学院,昆明 650500)

3)(中国科学院理论物理研究所,理论物理前沿重点实验室,北京 100190)

4)(中国科学院大学物理科学学院,北京 100049)

公共物品博弈是研究群体相互作用的经典模型,广泛用于解释自私个体间合作的涌现和保持.本文从理论分析和蒙特卡罗模拟两个方面研究了二维正方格子上一个有偿惩罚机制下随自我质疑更新规则演化的公共物品博弈模型的相变特性.理论分析方面,将公共物品博弈模型转化为一个外场不为零的铁磁Ising模型.通过有效能量发现:不存在惩罚时,个体间的耦合强度为零,体系只有外场作用;存在惩罚时,个体间包含最近邻、次近邻和第三近邻相互作用且外场不为零.蒙特卡罗模拟方面,首先验证了理论分析的正确性,然后对公共物品博弈模型相关的一级相变和二级相变进行了有限尺度标度分析.研究发现:1)蒙特卡罗模拟所得结果与类Ising模型分析结果完全吻合;2)相比二维Ising模型,公共物品博弈的二级相变临界指数发生了变化;3)公共物品博弈的一级相变与二维Ising模型相同.

Ising模型,有限尺度标度理论,蒙特卡罗模拟,自我质疑更新规则

1 引 言

演化博弈理论(evolutionary game theory)突破了传统博弈论关于“理性人”和“完全信息”的限制,它强调有限理性的博弈个体在重复博弈过程中根据自身掌握的局部信息,通过自适应学习做出尽可能占优的策略［1,2］.在人类社会和自然界中,自私个体之间产生合作是一个“惊人”的现象,得到了许多学者的重视和研究［3,4］.空间演化博弈理论(将多个个体放置在社会网络上进行重复博弈)的提出从空间结构和社会网络的角度很好地解释了合作的涌现,同时也开启了演化博弈研究的新篇章［1,2,5,6］.

公共物品博弈(public goods game)是一种典型的多人空间演化博弈模型,在博弈过程中,N个参与者独立决定向公共池子投资(合作)或不投资(不合作),初始投资总额通过增益系数放大r(1＜r＜N)倍后,平均分给群组内所有个体,无论初始投资与否.显然,采用不合作策略可以在事先不进行任何投资的情况下获得收益,搭上合作者的便车.公共物品博弈生动再现了个人与集体间的矛盾,有着广泛的应用,如:在班级大扫除中,每个学生付出的劳动是不对等的,但干净明亮的学习环境却是大家一起共享的.尽管所有成员都合作时,群组可以实现利益最大化,然而对一个自私且理性的个体而言,最优策略是选择背叛.从而在个体最优策略与群组最优策略之间形成了社会困境,即所谓的公共品悲剧(the tragedy of the commons)［1,7−9］.然而,现实生活中这种悲剧并未出现.为了解释合作的广泛传播与个体的自私性之间的矛盾,许多机制被相继提出,如:期望诱导重连机制、删边机制、惩罚机制、奖励机制、社会多样性机制、马太效应、志愿者参与机制、迁徙机制、调整投资机制、自适应有限投资反馈机制和条件策略［10,11］等.本文基于前人对惩罚机制的研究［12−18］,采用一类简化的有偿惩罚机制作用下的公共物品博弈模型进行研究.

空间演化博弈的两个研究重点是策略更新规则和网络的空间结构［1］.网络演化博弈的策略更新规则可以划分为两类:学习机制（模仿）和自我质疑机制（自省）.学习机制主要是向网络中的最近邻学习,包括学习邻居中收益最大的［5］;以一定的概率学习那些收益比自己高的个体们［19,20］;任意选择最近邻中的一个,比较收益差,以较大的概率模仿收益比自己高的个体［21−23］.自我质疑机制的研究相对较少［24−27］,它是指任意博弈个体先计算当前策略所得总收益,然后采用自身反策略进行一次虚拟博弈,并计算虚拟总收益,通过比较当前博弈和虚拟博弈收益变化的多少,来决定自己下一局所采用的策略.自我质疑机制与统计物理中Ising模型的单自旋翻转蒙特卡罗模拟方法十分类似.基于此,本文研究自我质疑更新规则下的公共物品博弈模型与Ising模型的等价关系［28］.之前,演化博弈模型的相关研究主要依靠计算机模拟方法,本文提出的类Ising方法将为该领域的研究提供理论基础.值得一提的是,经济学中广泛研究的随机最优反应均衡(quantal response equilibrium,QRE)与自我质疑机制十分类似,它用随机反应代替传统的最优反应,将纳什均衡作为它的一个特例进行处理［29−31］.

本文提出一个有偿惩罚机制下按自我质疑机制演化的公共物品博弈模型,从以下三个方面展开研究.理论分析方面,找到与该演化博弈模型等价的Ising模型并运用相关理论进行分析预测;蒙特卡罗模拟方面,验证理论分析的结果;相变特性方面,任选一组惩罚参数,研究博弈模型中伴随的相变现象,分析相变的种类,测定相关临界指数.

2 自我质疑更新规则下的公共物品博弈模型

2.1 模型的引入

在由N个个体参与的公共物品博弈中,每个个体都有两种策略可供选择,投资或不投资,假设原始投资额为1.投资总额放大r倍后(可认为是投资后获利),平均分给群组内所有参与博弈的个体,而不管个体是否向公共池内投过资,其中r为强化因子(enhancement factor),满足1＜r＜N.假设群组内有N+个合作者,则投资者的收益为gC=rN+/N−1,而不投资者的收益为gD=rN+/N.显然,不投资者的收益永远大于投资者的收益.所以,不投资者可以在投资者投资的基础上“搭便车”.从而,公共物品博弈的纳什均衡是,所有人都不投资.为了解释自私个体间仍有合作行为产生,本文提出一种简化的有偿惩罚机制,即惩罚的代价为C(cost),被惩罚者的处罚额度为F( fi ne).此时,投资者的收益为gC=rN+/N−1−C(N−N+),不投资者获得收益gD=rN+/N−FN+.

采用一个固定的网络来刻画个体之间的博弈关系.在网络模型确定之后,种群中的个体占据网络中的节点.每个个体都是纯策略者并只能选择合作或者不合作策略.每个个体与自己的邻居组成一个群体进行公共品博弈.在每一轮博弈中,任选一个个体i,然后根据自我质疑机制决定是否改变当前策略.即选定个体i后,个体i采用当前策略与周围邻居博弈获得当前收益Gi,采用当前策略的反策略与周围邻居进行一次虚拟博弈,获得虚拟收益个体i通过比较当前收益和虚拟收益,以一定概率决定是否改变当前策略,概率的选择为

2.2 公共物品博弈模型与Ising模型间的转化关系

在空间网络中,任意个体i,不仅与周围最近邻个体博弈,还与次近邻进行博弈.即个体i不仅参与以自己为中心节点进行的群体博弈,而且参与以邻居为中心节点的群体博弈,即总共经历ki+1次博弈,其中ki为节点i的度.当个体i参与以自己为中心节点的群体博弈时,参与博弈的个体总数为ki+1,当参与邻居j为中心节点的群体博弈时,参与人数为kj+1.图1显示了二维格点上个体i参与的五次公共物品博弈,深色格子表示博弈过程中的中心节点,浅色格子表示中心节点的最近邻节点,每轮博弈在由深色和浅色节点组成的群体中进行,白色格子表示其余节点.

图1 (网刊彩色)二维正方格子上任意个体i每轮博弈中参与以不同节点为中心的五次公共物品博弈Fig.1.(color online)The public goods games which involve player i on two-dimension square lattice each round.

个体i每一轮博弈的收益来自于ki+1个群体博弈所得收益的总和,为便于推导,用Si表示个体i的策略,有两种取值,Si=+1表示合作,Si=−1表示不合作,在以i为中心结点的群体博弈中,当Si=+1时,

其中ni+表∑示i的邻居中采用合作策略的个体总和,ni+=j(1+Sj/2);ni−表示∑i的邻居中采用不合作∑策略的个体总和,ni−=j(1−Sj)/2=(ki−jSj)/2,j是i的最近邻.当Si=−1时,

个体i翻转所带来的收益差为

将ni+和ni−的表达式代入(4)式,并用新记号∆gii描述个体i在以i为中心节点进行博弈的过程中所带来的收益变化,

同理可得,∆gij,个体i在以j为中心节点进行博弈时,改变策略带来的收益差为

其中,j是i的最近邻,表示对i的最近邻求和;表示在以j为中心的博弈群体中,对除i以外的节点进行求和.从而,一轮博弈中i改变策略带来的总收益差∆Gi为

将(5)式和(6)式代入(7)式,并整理可得

从上面的表达式可以看出,总收益的改变由两部分组成:前面部分只与博弈个体i所处的状态和空间结构有关,后面部分则包含个体i与最近邻j以及j的邻居k间的相互作用.显然,当不存在惩罚时,总收益差只与个体i和网络结构有关.也就是说,惩罚的引入直接导致个体间产生相互作用.为了便于研究和分析,本文以二维正方格子上的公共物品博弈为例,展开进一步讨论.在二维方格子中,ki=kj=4且∑j=4,表达式(8)中第一项可化简为Si［5−r−10(F−C)］,第二项中求和项是对所有与个体i相关的ki+1个集团中所有的参与者进行求和.图2是对所有求和节点的求和次数的示意图,可以将这样的求和规则改写成:两次最近邻求和加上两次次近邻求和再加上一次第三近邻求和,即其中j表示最近邻,j′表示次近邻,j′′表示第三近邻.代入到(8)式的第二项可得:

因此,表达式(8)可写成

公共物品博弈中追求利益最大化,而Ising模型则希望系统能量达到最小.所以,定义有效能量差为∆Ei=−∆Gi,相应的,

考虑最近邻、次近邻和第三近邻的Ising模型的表达式为

其中J表示最近邻相互作用强度,J′表示次近邻相互作用强度,J′′表示第三近邻相互作用强度,〈ij〉表示对最近邻求和,〈ij′〉表示对次近邻求和,〈ij′′〉表示对第三近邻求和,H表示外场强度.

图2 (网刊彩色)任意个体i参与的公共物品博弈等价于与最近次、次近邻分别进行两次两两博弈和与第三近邻进行一次两两博弈Fig.2.(color online)The number of interaction between player i and it’s neighbours.

Ising模型的局部能量变化关系∆E=E′−E可表示为

对比表达式(10)和表达式(12)可得公共物品博弈对应类Ising模型的耦合强度和外场为:

显然,耦合强度是惩罚费用和惩罚代价的函数,因为F＞0且C＞0,所以耦合强度J＞0.由Ising模型相关理论可知,该模型为铁磁Ising模型.此外,耦合强度的大小与空间结构和强化因子无关,只与惩罚的取值相关.外场由两部分构成:前一项表示一轮博弈中,均分后的利益是否大于原始投资;后一项来源于惩罚的引入,只要惩罚的费用大于惩罚的代价,外场都为正,即促进合作的产生.特别指出,运用矩阵分解法也可得到相似的等价关系［2］.

综合写出有效能量为

2.3 蒙特卡罗模拟

本部分采用Metropolis算法对周期边界条件下二维空间正方格子上的公共物品博弈模型进行蒙特卡罗模拟.空间格子的大小为N=L×L,用伪随机数生成任意初始构型后按照自我质疑机制进行演化.首先,随机选择群体中的任意博弈个体i,然后用自我质疑机制决定个体i是否改变当前策略,重复以上操作,直至系统趋于稳定状态为止.为了保证计算的准确性,扔掉前10000个蒙特卡罗模步(N次随机试验称为一个蒙特卡罗步).计算接下来的2×105个蒙特卡罗步的平均值.为消除初始条件对系统演化的影响,再选择100个不同的初始构型进行系综平均.

在演化博弈模型的社会学研究中,通常选择合作者占比f作为观测量,其表达式为

其中N代表参与博弈的总人数,N+代表总人数中采用合作策略的人数.

由表达式(17)可知,当不存在惩罚时,耦合强度为零,体系的有效能量为

显然,当r＞5时为正向外场,整个体系趋向于合作,而r＜5时为负向外场,体系趋向于不合作.如图3所示,合作与不合作以r=5为分界线,温度较低时,在强化因子逐渐增大的过程中,体系在r=5处,从完全不合作跳变到完全合作,随着温度的升高,这种跳变逐渐变得平缓.

图3 (网刊彩色)不存在惩罚时合作频率随强化因子的变化Fig.3.(color online)Frequency of cooperators vs.r for various T.

当存在惩罚时,体系表现为铁磁相互作用,系统状态随外场而变.显然,临界外场Hc=0,即1−rc/5−2(F−C)=0.解得rc=5−10(F−C),则当r＞rc时,系统为完全合作状态,r＜rc时系统为完全不合作状态.图4显示不同的惩罚额度F和惩罚的代价C下,合作频率f随强化因子r的变化,固定温度T=0.1.显然,蒙特卡罗模拟的结果与Ising模型理论分析的结果很好地吻合.合作频率以rc为突变点发生跳变.从社会学角度看,惩罚的引入增强了个体间的联系,相应地也增强了体系对不确定性因素T的鲁棒性.

图4 (网刊彩色)不同惩罚条件下合作频率随强化因子的变化Fig.4.(color online)Frequency of cooperators vs.r at T=0.1 for various punishment.

3 公共物品博弈模型的相变特性

从与公共物品博弈模型对应的Ising模型中看到,耦合强度J恒大于零,体系表现为铁磁相互作用.铁磁Ising的相变和临界现象在统计物理中已经广泛研究,基于判断有效能量推导的正确性和解释公共物品博弈诱发的相变现象两个原因,我们采用蒙特卡罗模拟和有限尺度标度理论详细分析公共物品博弈模型的相变特性.在Ising模型中,零场下随着温度的变化,存在铁磁序到顺磁序的二级相变(即由完全合作或完全不合作态向混乱状态的转变),相变温度称为临界温度,在临界温度附近表现出标度特性;外场非零时,只要温度小于临界温度,随着外场的变化会引发一级相变(由完全合作态到完全不合作态的转变或其逆过程).为了便于与Ising模型进行对比,本部分仍采用Ising模型的相关观测量展开研究,即选择磁化强度(magnetization)作为序参量.其表达式为

磁化强度与合作者占比的关系为m=2f−1.此外,定义磁化率为

一级相变和二级相变具有完全不同的相变特性,如二级相变由于关联长度的发散引发一系列奇异特性而一级相变则表现出δ函数奇异性.为了讨论方便,先简单介绍一级相变和二级相变的标度函数,然后再进行蒙特卡罗模拟.

二级相变的相变点通常称为临界点,在临界点附近,系统的关联长度发散,从而引起系统相关热力学量的奇异性,如比热和磁化率在临界点处发散.通常无穷大系统才表现出临界现象,然而无论是实际系统还是计算机模拟都不可能无穷大.Fisher于20世纪70年代提出的有限尺度标度理论通过对小系统的计算机模拟实现了对临界现象的研究［32］.具体而言,在临界点附近磁化强度和磁化率存在以下标度关系:

其中,fm和fχ为普适函数,t=(T−Tc)/Tc和h=(H−Hc)/Hc分别称为约化温度和约化外场,L是有限系统的特征长度,ν和νh分别代表关联长度对温度和外场发散的临界指数,1/νh=(γ+β)/ν,β和γ是描述热力学函数m和χ临界行为的临界指数.为确定临界点的位置,定义四阶Binder累积量为

Binder累计量的标度形式满足

其中U为普适函数,在临界点处,它与系统尺寸无关,即不同尺度下U的曲线交于临界点.

在无穷大系统中,一级相变的关联长度不发散且相变表现为δ函数奇异性,即磁化率可用一个δ函数表示.但在有限系统模拟中,δ函数的奇异性将被平滑掉.文献［33—37］运用热力学涨落理论获得了一级相变磁化强度、磁化率和四阶累积量的有限尺度标度函数.磁化强度为

其中Msp是热力学极限下的自发磁化强度,χ0是系统处于单相时的磁化率,d为系统维度,kB表示玻尔兹曼常数.

磁化率为

上式表明,在L→∞极限下发生在H=0处的δ函数奇异性,当L有限时被平滑为一个峰,其高度正比于Ld,其宽度正比于L−d.

四阶累计量的表达式为

从磁化强度、磁化率和四阶累计量的表达式看出,空间维度d是一级相变的唯一标度指数.

3.1 外场为零时的相变特性

公共物品博弈对应的有效能量包括:最近邻、次近邻和第三近邻相互作用,精确解未知.外场为零时,体系随温度变化表现为二级相变,通过四阶Binder累积量可以确定临界温度的位置,通过有限尺度标度理论可以获得相应的临界指数.模拟过程中,固定惩罚额度F=0.2,惩罚费用C=0.1,此时外场h=(r−4)/2.显然,当r＞4时,h＞0为正向外场,r＜4时,为负向外场,r=4时,外场为零.为研究的方便,先研究外场为零时系统随温度变化表现出的二级相变.

图5和图6是零外场时磁化强度和磁化率随温度变化的曲线,显然,不同尺寸下存在显著的有限尺度效应.磁化强度随着系统尺寸的增大,临界点附近的曲线变得越来越陡,越来越接近系统无穷大时的变化规律.热力学极限下,磁化率在临界点处发散,但在有限系统下为有限值,图6正好说明了这一点,随着系统尺度的不断增大,曲线的峰值变得越来越大.

图5 (网刊彩色)外场为零时不同尺度下磁化强度随温度的变化Fig.5.(color online)Magnetization m plotted vs temperature T at h=0 for various L.

图6 (网刊彩色)外场为零时不同尺度下磁化率随温度的变化Fig.6.(color online)Susceptibility χLplotted vs temperature T at h=0 for various L.

不同尺寸下的Binder累积量曲线相交于临界温度. 如图7所示,可确定临界温度Tc=1.065±0.0067.通过有限尺度标度理论,可以测定临界指数β/ν=0.154±0.0042,γ/ν=1.71±0.023,1/ν=1.085±0.027.将模拟结果与二维Ising模型的精确值对比,发现二者的临界指数存在一定差异,有可能是普适类发生了变化.为了进一步验证我们的设想,图8和图9用二维Ising模型的临界指数去坍塌模拟所得数据.可以明显看到临界点右半边曲线的坍塌效果较差,导致这一现象的原因是博弈模型中不仅包含最近邻相互作用,还包括次近邻和第三近邻相互作用.

图7 (网刊彩色)外场为零时不同系统尺度下四阶累积量随温度的变化Fig.7.(color online)Reduced cumulant U plotted vs temperature T at h=0 for various L.

图8 (网刊彩色)外场为零时临界温度附近磁化强度的有限尺度标度关系Fig.8. (color online)Scaled magnetization mL1/8 plotted vs scaled temperature εL1/1at T=1.065 and various L.

图9 (网刊彩色)外场为零时临界温度附近磁化率的有限尺度标度关系Fig.9.(color online)Scaled susceptibility χL/L7/4 plotted vs scaled temperature εL1/1at T=1.065 and various L.

3.2 外场非零时的相变特性

当温度小于临界温度时,沿外场演化模型表现为一级相变.当温度等于临界温度时,表现为二级相变.为更好地对比和区分两种相变,分别讨论T=1(T＜Tc,一级相变)和T=1.065(T=Tc,二级相变)两种情况.

图10和图11分别为T=1和T=1.065时磁化强度随外场的变化曲线.当系统无穷大时,在相变点r=5处,磁化强度不连续,发生跳变.当系统尺寸有限时,磁化强度曲线的奇异性被平滑掉.随着系统尺寸的不断增大,曲线变得越来越陡峭,系统无穷大时变得不连续.对比两图发现:在rc附近,T=Tc处磁化强度的变化比T＜Tc处的变化要缓慢一些.

图10 (网刊彩色)T＜Tc时不同尺度下磁化强度随强化因子的变化Fig.10.(color online)Magnetization m plotted vs the factor r at T=1 for various L.

图11 (网刊彩色)T=Tc时不同尺度下磁化强度随强化因子的变化Fig.11.(color online)Magnetization m plotted vs the factor r at T=1.065 for various L.

图12和图13分别展示了T=1和T=1.065时磁化率随外场的变化曲线.从(27)式知道,磁化率在相变点处随系统尺寸幂率发散,在两图中可以明显看到,随着系统尺度的增大,磁化率曲线的尖峰不断增大,当系统尺度为无穷大时发散,T=Tc处的峰值没有T＜Tc处的峰值高.这是由于系统无穷大时,一级相变在临界点附近δ函数发散而二级相变则按临界指数γ幂率发散.

图12 (网刊彩色)T＜Tc时不同尺度下磁化率随强化因子的变化Fig.12.(color online)Susceptibility χLplotted vs the factor r at T=1 for various L.

图13 (网刊彩色)T=Tc时不同尺度下磁化率随强化因子的变化Fig.13.(color online)Susceptibility χLplotted vs the factor r at T=1.065 for various L.

图14 (网刊彩色)T＜Tc时不同尺度下四阶累积量随强化因子的变化Fig.14.(color online)Reduced cumulant U plotted vs the factor r at T=1 for various L.

图14和图15得到了公共物品博弈相变点的位置,可明显看到不同尺度下所有曲线都交于相变点r=4处.

图15 (网刊彩色)T=Tc时不同尺度下四阶累积量随强化因子的变化Fig.15.(color online)Reduced cumulant U plotted vs the factor r at T=1.065 for various L.

图16 (网刊彩色)T＜Tc时临界外场处磁化率随尺度变化的双对数图Fig.16.(color online)Log-log plot of susceptibility maximum vs linear dimension at T=1.

图16和图17为磁化率随系统变化的双对数曲线,拟合曲线的斜率可以获得磁化率随尺度变化的临界指数,当T＜Tc时临界指数d=1.935±0.0093;当T=Tc时,临界指数γ/ν=1.714±0.0026.二维Ising模型临界指数的精确结果为d=2和γ/ν=1.75.显然二者间仍然有微小的差别,为了进一步验证模拟结果的正确性,图18和图19用二维Ising模型的临界指数对模拟所得数据进行坍塌,其中ε=(r−rc)/rc为约化增益系数,γ/ν的理论值为7/4,(γ+β)/ν的理论值为15/8.图16基于表达式(27)进行数据坍塌,在坍塌过程中忽略修正项χ0.从图中可以看出,公共物品博弈模型随外场变化的相变特性和二维Ising模型完全相同.

图17 (网刊彩色)T=Tc时临界外场处磁化率随尺度变化的双对数图Fig.17.(color online)Log-log plot of susceptibility maximum vs linear dimension at T=1.065.

图18 (网刊彩色)T＜Tc时临界外场附近磁化率的有限尺度标度关系Fig.18. (color online)Scaled susceptibility χL/L2 plotted vs scaled fi eld εL2at T=1 and various L.

图19 (网刊彩色)T=Tc时临界外场附近磁化率的有限尺度标度关系Fig.19.(color online)Scaled susceptibility χL/L7/4plotted vs scaled fi eld εL15/8at T=1.065 and various L.

图20和图21用临界指数的精确解对模拟得到的磁化强度的数据进行了坍塌.T＜Tc时,图20显示,在临界点rc附近曲线的坍塌效果较好,随着r逐渐远离rc,坍塌效果变差.其原因来源于表达式(26),在rc附近m≈hMspL2kBT.图21显示了较好的坍塌效果.

图20 (网刊彩色)T＜Tc时临界外场附近磁化强度的有限尺度标度关系Fig.20.(color online)Scaled magnetization m plotted vs scaled fi eld εL2at T=1 and various L.

图21 (网刊彩色)T=Tc时临界外场附近磁化强度的有限尺度标度关系Fig.21.(color online)Scaled magnetization mL1/8 plotted vs scaled fi eld εL15/8at T=1.065 and various L.

4 结 论

本文研究了有偿惩罚机制下随自我质疑更新规则演化的公共物品博弈模型.首先采用博弈模型收益差类比Ising模型能量差的方法,获得了二维正方格子上博弈模型的有效能量.有效能量显示:没有惩罚时,个体间不存在相互作用,合作策略的选取只与外场有关,外场则取决于原始投资和最终获得收益分红间的关系;存在惩罚时,个体间包括最近邻、次近邻和第三近邻相互作用,同时也减小了外场的作用(惩罚需要付出代价).也就是说,惩罚的引入增强了个体间的关联,促使博弈群体具有更强的鲁棒性,微小的扰动很难引起个体策略的改变.其次,计算机模拟博弈过程获得了不同理性程度和不同惩罚额度下合作者占比的变化曲线.结果显示:模拟结果和类Ising分析结果完全吻合.最后,以Ising模型相关的热力学量为基础研究了公共物品博弈模型的相变和临界现象.分别对二级相变和一级相变进行了讨论,结果显示,公共物品博弈模型随理性参数演化的二级相变与铁磁Ising模型随温度变化的二级相变不同,而一级相变则具有相同的有限尺度效应.需要特别说明的是,公共物品博弈通过集体分红获得收益,而传统的囚徒困境和雪堆博弈通过两两相互博弈获得收益.虽然二者的博弈形式不同,但都可以转化为与之等价的Ising模型,囚徒困境和雪堆博弈对应的Ising模型已经在参考文献［28］中给出.

［1］Szabó G,Fath G 2007Phys.Rep.446 97

［2］Szabó G,Borsos I 2016Phys.Rep.624 1

［3］He D R,Liu Z H,Wang B H 2009Complex Systems and Complex Networks(Beijing:Higher Education Press)p73(in Chinese)［何大韧,刘宗华,汪秉宏 2009复杂系统与复杂网络(北京:高等教育出版社)第73页］

［4］Nowak M A 2006Science314 1560

［5］Nowak M A,May R M 1992Nature359 826

［6］Nowak M A(translated by Li Z Q,Wang S C)2010Evolutionary Dynamics:Exploring the Equations of Life(Beijing:Higher Education Press)pp108–125(in Chinese)［Martin A.Nowak著 (李振清,王世畅 译)2010进化动力学-探索生命的方程(北京：高等教育出版社)第108—125页］

［7］Hardin G 1968Science162 1243

［8］Feeny D,Berkes F,McCay B J,Acheson J M 1990Hum.Ecol.18 1

［9］Palfrey T R,Rosenthal H 1984J.Public Econ.24 171

［10］Wang X W 2015M.D.Dissertation(Anhui:University of Science and Technology of China)(in Chinese)［王旭文2015博士学位论文(安徽:中国科学技术大学)］

［11］Wang B K 2013M.S.Dissertation(Xi’an:Xidian University)(in Chinese)［王保魁 2013硕士学位论文 (西安:西安电子科技大学)］

［12］Gao J 2012Ph.D.Dissertation(Xi’an:Xidian University)(in Chinese)［高佳 2012博士学位论文 (西安:西安电子科技大学)］

［13］Hauert C,Traulsen A,Brandt H,Nowak M A,Sigmund K 2007Science316 1905

［14］Dreber A,Rand D G,Fudenberg D,Nowak M A 2008Nature452 348

［15］Helbing D,Szolnoki A,Perc M,Szabó G 2010New J.Phys.12 083005

［16］Helbing D,Szolnoki A,Perc M,Szabó G 2010Phys.Rev.E81 057104

［17］Ohtsuki H,Iwasa Y,Nowak M A 2009Nature457 79

［18］Szolnoki A,Szabó G,Perc M2011Phys.Rev.E83 036101

［19］Hauert C,Doebeli M 2004Nature428 643

［20］Santos F C,Santos M D,Pacheco J M 2008Nature454 213

［21］Szabó G,Tőke C 1998Phys.Rev.E58 69

［22］Szabó G,Hauert C 2002Phys.Rev.Lett.89 118101

［23］Szabó G,Vukov J,Szolnoki A 2005Phys.Rev.E72 047107

［24］Wang W X,Ren J,Chen G,Wang B H 2006Phys.Rev.E74 056113

［25］Gao K,Wang W X,Wang B H 2007Physica A380 528

［26］Liu Y K,Li Z,Chen X J,Wang L 2009Chin.Phys.Lett.26 377

［27］Liu Y K 2010Ph.D.Dissertation(Xi’an:Xidian University)(in Chinese)［刘永奎 2010博士学位论文 (西安:西安电子科技大学)］

［28］Yang B,Li X T,Chen W,Liu J,Chen X S 2016Commun.Theor.Phys.66 439

［29］McKelvey R D,Palfrey T R 1995Games Econ.Behav.10 6

［30］He J T 2006Behavioral Economics The Economic Behavior of China(Beijing:China Financial Economic Publishing House)p79(in Chinese)［贺京同 2006行为经济学与中国经济行为(北京:中国财政经济出版社)第79页］

［31］Zhuang Q,Di Z R,Wu J S 2014PloS One9 e105391

［32］Fisher M E 1967Rep.Prog.Phys.30 615

［33］Privman V 1990Finite Size Scaling and Numerical Simulation of Statistical Systems(Singapore:World Scienti fi c)p223

［34］Landau D P,Binder K 2014A Guide to Monte Carlo Simulations in Statistical Physics(United Kingdom:Cambridge University Press)p82

［35］Binder K,Heermann D 2010Monte Carlo Simulation in Statistical Physics An Introduction(Berlin Heidelberg:Springer-Verlag)p56

［36］Binder K,Landau D P 1984Phys.Rev.B30 1477

［37］Lee J,Kosterlitz J M 1991Phys.Rev.B43 3265

Phase transition properties for the spatial public goods game with self-questioning mechanism∗

Yang Bo1)2)†Fan Min1)2)Liu Wen-Qi1)2)Chen Xiao-Song3)4)
1)(Data Science Research Center,Kunming University of Science and Technology,Kunming 650500,China)
2)(Faculty of Science,Kunming University of Science and Technology,Kunming 650500,China)
3)(Institute of Theoretical Physics,Key Laboratory of Theoretical Physics,Chinese Academy of Sciences,Beijing 100190,China)
4)(School of Physical Sciences,University of Chinese Academy of Sciences,Beijing 100049,China)

The spatial public goods game is one of the most popular models for studying the emergence and maintenance of cooperation among sel fi sh individuals.A public goods game with costly punishment and self-questioning updating mechanism is studied in this paper.The theoretical analysis and Monte Carlo simulation are involved to analyze this model.This game model can be transformed into Ising model with an external fi eld by theoretical analysis.When the costly punishment exists,the e ff ective Hamiltonian includes the nearest-,the next-nearest-and the third-nearestneighbor interactions and non-zero external fi eld.The interactions are only determined by costly punishment.The sign of the interaction is always greater than zero,so it has the properties of ferromagnetic Ising.The external fi eld is determined by the factorrof the public goods game,the fi neFon each defector within the group,and the relevant punishment costC.The Monte Carlo simulation results are consistent with the theoretical analysis results.In addition,the phase transitions and critical behaviors of the public goods game are also studied using the fi nite size scaling theory.The results show that the discontinuous phase transition has the same fi nite size e ff ects as the two-dimensional Ising model,but the continuous phase transitions is inconsistent with Ising model.

ising model, fi nite size scaling theory,Monte Carlo simulations,self-questioning update rules

27 May 2017;revised manuscript

4 July 2017)

(2017年5月27日收到;2017年7月4日收到修改稿)

10.7498/aps.66.196401

∗昆明理工大学引进人才科研启动基金项目(批准号：KKSY201607047)和国家自然科学基金(批准号:61573173)资助的课题.

†通信作者.E-mail:yangbo@kmust.edu.cn

©2017中国物理学会Chinese Physical Society

PACS:64.60.De,87.55.K–,02.50.Le,87.23.Ge

10.7498/aps.66.196401

*Scienti fi c Research Foundation for Introduced Scholars,Kunming University of Science and Technology(Grant No.KKSY201607047)and the National Natural Science Foundation of China(Grant No.61573173).

†Corresponding author.E-mail:yangbo@kmust.edu.cn