美赛随机图代码包：randomgraph.m 与复杂网络建模实战

简介这份资源面向参加美国数学建模竞赛MCM/ICM的学生与复杂网络初学者聚焦随机图算法的代码实现帮助读者在建模中快速搭建网络模型并验证拓扑特性。压缩包内共1个文件为MATLAB脚本.m体积约2KB轻量易读可直接运行或按需修改参数。内容围绕Erdős-Rényi模型、Barabási-Albert无标度网络、Watts-Strogatz小世界模型等经典随机图生成方法展开并可能涉及生成树、社区检测、网络动力学与稳健性分析等延伸方向。目前已有149人学习适合作为赛前速查与代码参考。读者可借助脚本理解节点连接概率、幂律度分布、随机重连等机制将其嵌入自己的建模流程用于模拟传播过程、评估网络抗毁性或结合NetworkX、Gephi做可视化呈现从而提升模型构建效率与论文说服力。1. 美赛随机图代码包从 randomgraph.m 到可复现的复杂网络建模美赛赛题里只要出现「网络」「传播」「节点关系」这类字眼随机图基本是绕不开的建模底座。这个压缩包给了一份randomgraph.m属于美赛常见参考代码里偏底层的那一类——它不直接给你一张漂亮的结果图而是把随机图的生成逻辑摊开让你能改参数、换模型、接自己的动力学仿真。适合两类人一类是第一次打美赛、需要快速搭出网络模型骨架的参赛者另一类是平时用 NetworkX 或 igraph 做研究、想看看 MATLAB 侧怎么把 ER、BA、WS 这些模型用矩阵方式写干净的人。它解决的不是「画个图交差」而是「网络邻接矩阵怎么生成、度分布怎么验证、后续传播仿真怎么接上去」这条链路。下面按我拆包和跑代码的顺序讲重点放在参数含义、复现步骤和几个容易翻车的地方。2. 拆开 randomgraph.m邻接矩阵生成与三种经典模型2.1 先看清这份代码的输入输出约定拿到.m文件别急着 run先看函数签名和返回结构。这类美赛参考代码通常写成函数形式输入是节点数N、连边概率p或每步新增边数m输出是一个N×N的邻接矩阵A。判断约定是否合理看三点对角线是否为零无自环、矩阵是否对称无向图、元素是否只有 0/1无权图。如果代码里出现rand(N) p这种写法说明它用的是 ER 的 G(N,p) 版本而不是固定边数的 G(N,M) 版本两者在相变点附近行为不同选错会让你的连通性结论偏掉。我一般会先跑一个最小验证确认矩阵结构没问题再往下接% 最小验证生成 100 节点、p0.05 的 ER 随机图 N 100; p 0.05; A randomgraph(N, p); % 假设函数签名为 (N, p) % 结构检查 assert(isequal(A, A), 邻接矩阵不对称检查生成逻辑); assert(all(diag(A) 0), 存在自环对角线未清零); assert(all(A(:) 0 | A(:) 1), 存在非 0/1 元素); % 统计实际边数与平均度 E sum(A(:)) / 2; k_avg 2 * E / N; fprintf(边数 E%d, 平均度 k%.2f\n, E, k_avg);逻辑说明assert三连是防止代码内部用了randi或对称化处理时留下脏数据。参数说明N决定矩阵规模p直接控制稀疏程度p太小图会碎成多个连通片p太大就退化成近似完全图失去随机图的意义。经验上 ER 的连通相变阈值在p ≈ ln(N)/N附近100 节点对应约 0.046所以上面取 0.05 刚好在相变点上方能观察到巨片出现。2.2 ER、BA、WS 三种模型的参数差异这份代码包如果覆盖多个模型通常会用不同函数或一个model参数区分。ER 的核心参数是(N, p)BA 的核心是(N, m)m表示每个新节点带入的边数m1生成树状结构m≥2才有无标度特征WS 的核心是(N, K, beta)K是环形规则网络的近邻数beta是重连概率。三者不能混用参数我见过有人把p0.1直接塞给 BA结果代码按m0.1取整成 0生成一堆孤立节点还以为是算法错了。模型典型调用关键参数度分布特征适用赛题场景ERrandomgraph(N,p)连边概率 p泊松分布均匀随机接触、基线对照BArandomgraph(N,m)新增边数 m幂律分布社交网络、无标度传播WSrandomgraph(N,K,beta)近邻数 K、重连率 beta近似泊松高聚类小世界、规则到随机过渡选型理由如果赛题强调「少数节点影响全局」用 BA如果强调「短路径但局部抱团」用 WS如果只是要一个无偏的随机对照用 ER。别一上来就 BA无标度网络的鲁棒性结论和 ER 完全相反选错模型整篇论文的定性分析都会歪。2.3 把邻接矩阵接到度分布与连通性验证生成矩阵只是第一步真正决定论文能不能写下去的是验证环节。度分布用histcounts或hist统计每行求和连通性用广度优先或 MATLAB 自带的conncomp需要 graph 对象。我习惯把矩阵转成 graph 再算省得自己写遍历% 度分布与连通片统计 deg sum(A, 2); % 每行求和即节点度 figure; histogram(deg, Normalization, pdf); xlabel(度 k); ylabel(P(k)); title(度分布); G graph(A); comp conncomp(G); numComp max(comp); giantSize max(histcounts(comp, 1:numComp1)); fprintf(连通片数%d, 最大片占比%.2f\n, numComp, giantSize/N);逻辑说明sum(A,2)按行求和得到度序列比循环快一个量级。参数说明conncomp返回每个节点所属片编号max(histcounts(...))取最大片大小。判断标准ER 在相变点上方最大片占比应接近 1BA 通常天然连通WS 在beta较小时可能仍保持高聚类但连通片不多。如果最大片占比只有 0.3 还硬说网络连通审稿人一眼就能看出来。3. 从矩阵到仿真传播动力学与稳健性分析的接法3.1 在随机图上跑 SIR 传播的最小闭环美赛网络题十有八九要落到传播或扩散上SIR 是最常被拿来当骨架的模型。核心是把邻接矩阵当成接触通道每个时间步让感染节点按概率lambda感染邻居同时以概率mu恢复。下面这段可以直接接在randomgraph输出后面% SIR 在随机图上的离散时间仿真 lambda 0.3; % 感染率 mu 0.1; % 恢复率 T 100; % 仿真步数 N size(A, 1); state zeros(N, 1); % 0易感, 1感染, 2恢复 state(1) 1; % 初始感染源 S zeros(T,1); I zeros(T,1); R zeros(T,1); for t 1:T newInf false(N,1); infIdx find(state 1); for i infIdx neighbors find(A(i, :) 1); sus neighbors(state(neighbors) 0); newInf(sus) newInf(sus) | (rand(length(sus),1) lambda); end recIdx infIdx(rand(length(infIdx),1) mu); state(recIdx) 2; state(newInf) 1; S(t) sum(state0); I(t) sum(state1); R(t) sum(state2); end plot(1:T, S, 1:T, I, 1:T, R); legend(S,I,R); xlabel(时间步); ylabel(人数);逻辑说明先收集感染节点的邻居再对易感邻居做伯努利抽样避免同一节点在同一步被重复感染。参数说明lambda和mu的比值决定基本再生数R0 ≈ lambda/mu * kk是平均度。如果R0 1疫情会自行消退R0 1才出现峰值。这个关系是论文里解释「为什么某参数下传播不起来」的关键论据别只贴图不给数。3.2 稳健性分析随机删点和目标删点的差别网络稳健性常考「删掉多少节点网络会碎」。随机删点用randperm打乱后逐步移除目标删点按度从大到小移除。两种策略在 BA 网络上差异极大——BA 对随机删点极其鲁棒但对高度节点定向攻击非常脆弱这正是无标度网络「既稳健又脆弱」的来源。实现时每删一批就重算最大连通片占比% 随机删点 vs 目标删点 fractions 0:0.05:0.5; giantRand zeros(size(fractions)); giantTarget zeros(size(fractions)); for idx 1:length(fractions) f fractions(idx); nRemove round(f * N); % 随机删点 rmRand randperm(N, nRemove); A1 A; A1(rmRand,:) 0; A1(:,rmRand) 0; c1 conncomp(graph(A1)); giantRand(idx) max(histcounts(c1, 1:max(c1)1)) / N; % 目标删点按度降序 [~, order] sort(sum(A,2), descend); rmTgt order(1:nRemove); A2 A; A2(rmTgt,:) 0; A2(:,rmTgt) 0; c2 conncomp(graph(A2)); giantTarget(idx) max(histcounts(c2, 1:max(c2)1)) / N; end plot(fractions, giantRand, fractions, giantTarget); legend(随机删点,目标删点); xlabel(移除比例); ylabel(最大片占比);逻辑说明删点要把对应行和列同时清零只清行会留下不对称矩阵导致graph报错。参数说明fractions是移除比例步长 0.05 足够看出相变拐点。判断标准BA 的随机删点曲线在 0.3 之前几乎不降目标删点在 0.1 到 0.2 之间就急剧下降ER 两条曲线差别不大。这个对比图是稳健性章节的核心证据。3.3 社区检测与模块度Louvain 思路的 MATLAB 近似社区检测在美赛里常用来给网络做功能分区。MATLAB 没有内置 Louvain但可以用graph的conncomp做粗分或者自己写一个简化的标签传播。核心指标是模块度QQ 0.3通常认为有显著社区结构。如果代码包里带了社区检测函数重点看它返回的是标签向量还是社区元胞标签向量更方便算Q。没有的话用邻接矩阵和标签手算% 模块度 Q 计算给定社区标签 function Q modularity(A, labels) m sum(A(:)) / 2; k sum(A, 2); Q 0; communities unique(labels); for c communities idx find(labels c); Q Q (sum(sum(A(idx, idx))) / (2*m)) - (sum(k(idx)) / (2*m))^2; end end逻辑说明第一项是社区内部边占比第二项是随机情况下的期望占比差值即模块度。参数说明labels是N×1向量每个元素是社区编号。注意A(idx,idx)取子矩阵后求和要除以2m因为无向图每条边算两次。如果Q算出来是负数说明社区划分比随机还差得回去检查标签是不是随机生成的。4. 避坑与排查randomgraph.m 跑不通时的五个体检项4.1 现象矩阵不对称graph 报错「邻接矩阵必须对称」原因代码里用了randi([0,1], N)生成上三角后忘了对称化或者只对A(i,j)赋值没管A(j,i)。解决生成后强制A triu(A,1); A A A;再检查对角线清零。这个坑在自写 ER 时最常见尤其是想省内存只存上三角的情况。4.2 现象BA 模型生成一堆孤立节点平均度接近 0原因把概率p误传给了m参数m被取整成 0新节点没带任何边。解决确认 BA 调用时第二个参数是整数且m ≥ 1并在函数入口加assert(m 1 m round(m))。如果代码包里 ER 和 BA 共用一个入口参数名可能都是p这时候要看内部有没有model分支判断。4.3 现象WS 模型重连后出现重复边或自环原因重连时随机选目标节点没排除自身和已有邻居导致A(i,i)1或同一条边被连两次。解决重连前先取候选集setdiff(1:N, [i, find(A(i,:))])从候选集里抽抽完立即更新对称位置。WS 的重连逻辑是三个模型里最容易写错的血泪经验是每重连一条边就断言一次无自环。4.4 现象SIR 仿真曲线一直不降感染率明显偏高原因lambda和mu的量纲没对齐或者每个时间步对同一邻居重复抽样导致实际感染概率被放大。解决先算R0 lambda/mu * mean(sum(A,2))如果R0远大于 1 却期望疫情消退那是参数设错了再检查感染循环里是否对已感染节点也做了抽样。常见做法是每步只处理state1的节点且对每个易感邻居只抽一次。4.5 现象删点后最大连通片占比不降反升原因删点后conncomp的编号变了histcounts的边界用了旧的max(c)导致统计错位。解决每次删点后重新算max(c)再传给histcounts或者直接用max(accumarray(c, 1))统计各片大小避免边界依赖。这个坑很隐蔽图看起来正常但数值全错属于典型的黑匣子式翻车。5. 进阶技巧把 randomgraph.m 改成可复现的批量实验脚本单次生成随机图有个玄学问题换个随机种子度分布和连通性可能差很多论文里的结论就不稳。我一般会把randomgraph包一层批量循环固定种子集跑多次取均值和置信区间。MATLAB 里用rng(seed)控制可复现性把每次的giantSize/N、Q、R0存成表最后画误差棒。这样审稿人问「你是不是挑了一次好看的结果」时你有分布说话。% 批量实验20 个种子下的巨片占比统计 seeds 1:20; giantFrac zeros(size(seeds)); for s 1:length(seeds) rng(seeds(s)); A randomgraph(200, 0.05); c conncomp(graph(A)); giantFrac(s) max(accumarray(c, 1)) / 200; end fprintf(巨片占比均值%.3f, 标准差%.3f\n, mean(giantFrac), std(giantFrac)); errorbar(1, mean(giantFrac), std(giantFrac), o);逻辑说明rng(seeds(s))保证每个种子下的随机流独立且可复现accumarray统计各连通片大小比histcounts更稳。参数说明种子数 20 是经验值太少置信区间太宽太多跑得慢节点数 200、p0.05对应ln(200)/200 ≈ 0.026在相变点上方巨片应该稳定出现。如果标准差超过 0.1说明p太接近阈值得加大p或增加节点数。还有一个实用技巧把randomgraph.m里的随机数生成换成rand的显式调用并记录种子而不是依赖全局状态。这样即使别人拿到你的脚本只要种子一样结果就一样。我现在的习惯是每个实验脚本开头必写rng(42)中间不再动全局随机流需要独立流就用RandStream。从那以后我每次交论文前都强制走一遍「换三个种子重跑核心图」的流程确认结论不随种子翻转才敢定稿。希望帮到你。本文还有配套的精品资源点击获取