Parkes误差网格MATLAB仿真：血糖仪临床风险评估与Type1边界实现

Parkes误差网格MATLAB仿真：血糖仪临床风险评估与Type1边界实现

发布时间：2026/9/16 22:18:31
Parkes误差网格MATLAB仿真：血糖仪临床风险评估与Type1边界实现 简介面向本硕博教研学习场景这份MATLAB仿真资源专注血糖样本的Parkes误差网格分析算法完整覆盖网格边界定义、样本点识别与结果可视化环节便于读者快速理解算法思想并对照编程实现。压缩包共6个文件以3个MATLAB脚本为核心分别承担主运行流程、第一类误差网格边界计算和样本自动识别功能另含1个操作录像视频、1个说明文档和1张结果示例图整体大小仅220KB轻量易用。目前已有430人学习/下载适合作为研究生课题或本科毕业设计的参考素材。配套操作录像针对MATLAB 2021a及更高版本演示了正确的环境配置与运行步骤清晰提示应运行主脚本而非子函数、当前文件夹必须指向工程目录等要点能有效避免新手常见失误说明文档和示例图片则可辅助核对仿真结果是否准确。1. Parkes误差网格为什么比相关系数更适合血糖样本评估在血糖仪算法验证里常有人拿相关系数R²说事但R²只能说明测量值与参考值的线性趋势无法回答最关键的临床问题当血糖仪报出一个错误数值时会不会让病人做出危险的处理决定同一个偏差在低血糖区和高血糖区的后果完全不同。Parkes误差网格正是为此设计的。它把参考血糖和测量血糖同时映射到二维平面用一组由临床共识构建的折线边界把结果划分为A到E五个风险等级。这套MATLAB仿真资源提供了Type1完整的边界定义、分区识别和可视化流程运行Runme_EGA.m即可得到类似论文里的散点网格图。对于做血糖仪数据后处理、连续血糖监测信号分析、或医学信息课程设计的本硕博学生它比那些只画散点图加拟合线的demo实用得多。而且Type1边界参数和Type2不同不能混用。2. Parkes_EGA_boundaries_Type1.mType1边界坐标构建与网格绘制2.1 边界数据的组织方式打开Parkes_EGA_boundaries_Type1.m会发现这个函数没有任何循环或判断它的全部作用就是返回一个结构体里面存着各分区上下边界的折线顶点。常见做法是用结构体字段区分区域每个字段是一个Nx2矩阵第一列是参考血糖浓度第二列是对应的测量血糖浓度。以Type1网格为例A区被对角线分成上下两部分upperA字段存A区上方边界的顶点序列lowerA存A区下方边界。B、C、D区域同理E区是这些边界之外的所有区域不需要单独定义。function B Parkes_EGA_boundaries_Type1() % 返回Type1 Parkes误差网格边界结构体 % 字段含义 % upperA: A区上边界顶点 [x, y] % lowerA: A区下边界顶点 % upperB/lowerB, upperC/lowerC, upperD/lowerD 同理 % % 以下顶点坐标仅为字段结构展示实际数值以本文件内完整定义为准 B.upperA [0, 0; 50, 25; 100, 80; 300, 280; 600, 600]; B.lowerA [0, 0; 60, -20; 180, -40; 600, -200]; B.upperB [0, 50; 70, 100; 200, 180; 600, 600]; % ... 其余边界按相同格式补充 end这段代码里的坐标是简化示例。真实Parkes网格的边界在低血糖区有水平折线在高血糖区有放大的锥形区域直接使用示例坐标会得到错误的临床分区结果。正确做法是打开资源里的Parkes_EGA_boundaries_Type1.m把里面实际定义好的顶点矩阵传给绘图或识别函数。为什么要用结构体而不是多个返回值因为调用时只需要一个变量B就能同时访问所有边界而且结构体字段名有自描述性后续代码里写B.upperA远比写boundaryA1、boundaryA2这种零散变量清晰。另一个设计细节是顶点顺序必须从左到右即x坐标递增排列。这样后续做线性插值的时候interp1不需要显式设置ascending选项也不容易出现因为顶点乱序导致的单调性错误。2.2 绘制网格的完整调用方式拿到边界结构体后画网格就是一个纯粹的plot过程。这里给出一个可直接运行的网格绘制脚本B Parkes_EGA_boundaries_Type1(); figure(Name, Parkes Error Grid Type1, Color, w, Position, [100 100 680 600]); hold on; % 对角线参考值等于测量值 plot([0 600], [0 600], --, Color, [0.6 0.6 0.6], LineWidth, 0.8); % 绘制各区分割线 plot(B.upperA(:,1), B.upperA(:,2), k-, LineWidth, 1.4); plot(B.lowerA(:,1), B.lowerA(:,2), k-, LineWidth, 1.4); plot(B.upperB(:,1), B.upperB(:,2), b-, LineWidth, 1.0); plot(B.lowerB(:,1), B.lowerB(:,2), b-, LineWidth, 1.0); plot(B.upperC(:,1), B.upperC(:,2), m-, LineWidth, 1.0); plot(B.lowerC(:,1), B.lowerC(:,2), m-, LineWidth, 1.0); plot(B.upperD(:,1), B.upperD(:,2), r-, LineWidth, 1.0); plot(B.lowerD(:,1), B.lowerD(:,2), r-, LineWidth, 1.0); axis([0 600 0 600]); xlabel(Reference Glucose (mg/dL), FontSize, 11); ylabel(Measured Glucose (mg/dL), FontSize, 11); title(Parkes Error Grid - Type 1 Diabetes, FontSize, 12); box on; grid off;执行这段代码后画布上会出现完整的Parkes网格框架。注意plot(B.upperA(:,1), B.upperA(:,2))这种写法它把矩阵的第一列作为x坐标、第二列作为y坐标连起来就是一条折线。由于顶点序列是连续的plot会把它们按顺序连接形成分段线性边界。参数说明坐标轴范围固定为0到600 mg/dL这个上限足以覆盖临床极端高血糖。线型与颜色建议区分边界等级A区用黑色和较粗的线宽B区用蓝色C区用品红D区用红色E区不画线用底色表达。这样叠加数据点之后读者能一眼看出点落在哪个风险区域。2.3 边界数据和识别逻辑分离的原因在这个仿真资源中Parkes_EGA_boundaries_Type1.m只负责“定义临床共识边界”Parkes_EGA_identify.m负责“用几何算法判断样本归属”。两者分离带来一个明显好处当需要评估不同版本网格时只要修改边界函数识别代码一行都不用动。例如后续做Type2网格的对照组实验只需复制一份边界函数替换坐标常量主脚本里改一个函数名即可。另一个工程层面的原因是可测试性。边界数据是静态的识别算法是动态的分离后可以单独为边界函数写自检脚本检查每条折线是否满足单调性、覆盖范围是否完整。如果混在一起每次跑仿真都要重新编译整个流程排错成本高。这个资源的文件命名也体现了这个思想函数名级别就标注了“boundaries”和“identify”从右侧操作录像里也能看到这种模块化设计。3. Parkes_EGA_identify.m样本分区判定与临床风险分级3.1 判定原理从折线边界到区域归属识别函数解决的是一个点(参考血糖, 测量血糖)到底落在哪个区域的问题。Parkes网格的边界不是封闭多边形而是几组分段的折线。对任意一个参考血糖值x通过边界函数可以在各条折线上插值出一个y的阈值。比如对A区插值得到upperA_bound和lowerA_bound如果测量值在这两者之间就属于A区否则继续判断B区、C区、D区最后落不到任何区域就是E区。function zone Parkes_EGA_identify(ref, meas, B) % ref: 参考血糖浓度标量或列向量 (mg/dL) % meas: 测量血糖浓度标量或列向量 (mg/dL) % B: Parkes_EGA_boundaries_Type1返回的结构体 % zone: 字符数组A ~ E % 将输入裁剪到网格定义域内避免插值越界 ref min(max(ref, 0), 600); meas min(max(meas, 0), 600); % 初始化zone为E后续逐级覆盖 zone repmat(E, size(ref)); % 判断是否在D区 zone((meas interp1(B.lowerD(:,1), B.lowerD(:,2), ref, linear, extrap)) ... (meas interp1(B.upperD(:,1), B.upperD(:,2), ref, linear, extrap))) D; % 判断是否在C区 zone((meas interp1(B.lowerC(:,1), B.lowerC(:,2), ref, linear, extrap)) ... (meas interp1(B.upperC(:,1), B.upperC(:,2), ref, linear, extrap))) C; % 判断是否在B区 zone((meas interp1(B.lowerB(:,1), B.lowerB(:,2), ref, linear, extrap)) ... (meas interp1(B.upperB(:,1), B.upperB(:,2), ref, linear, extrap))) B; % 判断是否在A区 zone((meas interp1(B.lowerA(:,1), B.lowerA(:,2), ref, linear, extrap)) ... (meas interp1(B.upperA(:,1), B.upperA(:,2), ref, linear, extrap))) A; end这段代码的关键点在于覆盖顺序先赋D再赋C、B、A。因为每个点理论上可能同时满足D区的上下界和B区的上下界MATLAB的字符串数组赋值是逐元素覆盖最后一次赋值最终生效。A区范围最小放在最后判断能保证所有满足A区条件的点最终落在A区。如果把A放最前面那么同样在B区范围内的A区点就会被错误覆盖成B。interp1的extrap参数允许外推这是因为部分边界的顶点并没有覆盖0到600的完整横轴范围。比如某些D区上边界可能只从150到600有定义当ref120时extrap会用最近的两个顶点线性延伸出一个值。这个外推虽然物理意义不强但能防止interp1返回NaN避免后续比较失败。3.2 边界插值方法的选型与边界条件Parkes网格的边界本身就是分段线性折线因此interp1使用linear方法就足够不需要smooth或pchip。使用spline反而会在折线拐点处产生过冲导致原本不存在的边界波动。这是很多新手会踩的坑看到边界是曲线形状就默认用更高阶插值实际上原定义就是直线段连接。另外要注意所有边界折线的x坐标不能有重复值。如果某个边界有一段垂直上升的线也就是同一个x对应多个y那interp1会崩溃。如果遇到这种情况解决方式是拆分成两条边或者对折线进行逆变换把垂直段转成极短的斜线段。在原始Parkes坐标系里边界函数通常是单调的所以该问题不常见但修改边界坐标时一定要检查。3.3 识别结果与临床风险的对应分区识别之后字母本身并不直接告诉用户“有多危险”需要对照Parkes误差网格的临床定义区域临床含义典型场景A区测量值在参考值允许偏差内治疗决策不受影响正常血糖监测调整胰岛素剂量B区偏差增加但不会导致错误治疗或遗漏风险轻微的传感器漂移C区可能导致不必要的修正如过度补糖低血糖值被高估引发不必要进食D区可能掩盖真实的低血糖或高血糖造成延误严重低血糖被误判为安全值E区测量值与参考值矛盾可能导致反向治疗低血糖被报成高血糖注射胰岛素做算法评估时通常要求A区占比不低于95%AB不低于98.5%D区和E区应为0。但不同论文和不同国家的标准略有差异。这套MATLAB仿真只负责计算“落在哪个区”最终是否通过标准由使用者自行设定阈值。3.4 向量化批量判断的注意事项如果样本量达到数万条逐点调用Parkes_EGA_identify会非常慢。推荐的做法是像上面代码那样把整个ref向量一次性传给interp1得到一组上边界阈值和下边界阈值向量然后做向量化比较。这里有个容易被忽略的问题repmat(E, size(ref))生成的是字符数组后续zone(...)D这种赋值要求三个逻辑索引矩阵大小完全一致。如果ref是列向量meas也是列向量一切正常但如果某个输入是行向量另一个是列向量MATLAB会触发广播语义导致结果维度变成矩阵后续统计就会出错。稳妥方式是在函数入口统一强制列向量ref ref(:); meas meas(:);这一行强制转换能避免80%的维度灾难。另一个问题是interp1要求输入x是单调递增。如果边界折线中出现了x递减段需要先按x排序[xs, idx] sort(vertices(:,1)); ys vertices(idx,2);。资源中的边界数据已经保证了单调但自己扩展Type2时一定要做这个检查。4. Runme_EGA.m 串联仿真流程与结果可视化4.1 主脚本的执行顺序Runme_EGA.m是这套仿真的入口。运行它MATLAB会依次执行加载数据、读取边界、批量分区、统计比例、绘制图像、保存output.png。它的设计意图是让使用者不需要懂内部算法也能跑通完整流程。下面是主脚本的核心结构%% 清空环境 clear; close all; clc; %% 1. 加载血糖样本 % 数据格式两列第一列参考血糖第二列测量血糖无表头 data load(sample_glucose_data.txt); ref data(:, 1); meas data(:, 2); %% 2. 获取Type1边界 B Parkes_EGA_boundaries_Type1(); %% 3. 调用识别函数 zones Parkes_EGA_identify(ref, meas, B); %% 4. 统计各区域百分比 n length(zones); fprintf(样本总数: %d\n, n); fprintf(A区: %.2f%%\n, sum(zones A) / n * 100); fprintf(B区: %.2f%%\n, sum(zones B) / n * 100); fprintf(C区: %.2f%%\n, sum(zones C) / n * 100); fprintf(D区: %.2f%%\n, sum(zones D) / n * 100); fprintf(E区: %.2f%%\n, sum(zones E) / n * 100); %% 5. 绘制散点与网格 figure(Name, Parkes EGA Simulation, Color, w); hold on; % 先画网格省略具体plot代码同第2章 Parkes_EGA_draw_grid(B); % 若有封装函数则调用 % 按区域着色 cmap { A, [0.2 0.7 0.2]; B, [0.2 0.4 0.9]; C, [0.9 0.8 0.1]; D, [1.0 0.4 0.1]; E, [0.8 0.0 0.0]; }; hold on; for i 1:n color cmap{strcmp(cmap(:,1), zones(i)), 2}; scatter(ref(i), meas(i), 24, color, filled, MarkerEdgeColor, k, LineWidth, 0.3); end xlim([0 600]); ylim([0 600]); xlabel(Reference Glucose (mg/dL)); ylabel(Measured Glucose (mg/dL)); title(Parkes Error Grid Type 1 - Sample Results); %% 6. 保存结果 saveas(gcf, output.png);代码中第5步是逐点绘制散点样本量达到上千时绘制会很慢。更高效的做法是用scatter的一次性调用并传入RGB三元组矩阵但需要提前构造颜色矩阵。这段代码的目的更偏向教学演示逐点循环更容易理解。注意load命令对文本文件的要求是纯数字矩阵如果数据文件带表头或逗号分隔要改用readmatrix。资源里提供的是txt格式数据组织方式以实际为准。4.2 output.png的解读要点运行成功后生成output.png它应该是网格背景加彩色散点的图。解读时先看整体分布散点应沿对角线散布颜色以绿色A区和蓝色B区为主。如果出现大面积黄色或红色说明样本误差大到会影响临床决策。再单独看低血糖区域即横坐标小于70 mg/dL的带状区域。这里最容易出现C区和D区点因为低血糖时血糖仪的相对偏差通常更大而Parkes网格在低值段有特殊的平行边界设计。比较隐蔽的问题是边界外推。比如ref650mg/dL在interp1的extrap作用下点能被分类但这个分类结果并没有真实临床依据。所以output.png上如果出现超出600的点首先要怀疑数据预处理没做异常值过滤。4.3 运行时的路径与版本注意事项资源内自述强调使用MATLAB 2021a或更高版本直接运行Runme_EGA.m不要单独运行子函数文件。原因是子函数依赖边界结构体的输入裸运行会报“输入参数不足”。另外运行时必须保证MATLAB左侧“当前文件夹”窗口指向工程目录。我遇到过几次报错“未定义函数或变量”多半是当前文件夹停在MATLAB默认路径。处理方式是cd到工程目录再运行。如果操作录像显示某一步界面不同不用紧张只要函数名和主脚本名称一致算法的输出不会因为MATLAB版本产生变化。2021a以后readmatrix和interp1行为保持一致绘图细节可能略有差异但output.png的核心内容不受影响。5. 边界参数调整与验证技巧让Parkes网格适配不同数据集5.1 切换Type1到Type2边界的标准化操作很多人在做完Type1仿真后会想对比Type2网格。Parkes误差网格最初针对1型和2型糖尿病分别定义了边界差别主要体现在低血糖区的宽容度。资源目前只给了Type1但通过合理扩展可以把它改写为Type2版本。标准做法是复制Parkes_EGA_boundaries_Type1.m为Parkes_EGA_boundaries_Type2.m把函数名和结构体中的顶点坐标替换成Type2临床共识中的数据。识别函数和主脚本不需要改动只需在主脚本里把B Parkes_EGA_boundaries_Type1();改成B Parkes_EGA_boundaries_Type2();此时zones的计算、统计和绘图都会自动适应新边界。这个操作能复用的前提就是前面坚持“边界与识别分离”如果当初把坐标写死在识别函数里现在就得重写整个算法。5.2 边界坐标合法性自检拿到任意一组边界坐标在正式做样本识别之前先跑一个自检脚本确认边界没有交叠或反向。自检逻辑很简单对每条边界折线检查x坐标是否严格单调递增再检查同一区域内upper边界的y值是否始终大于lower边界的y值。参考代码如下function check_boundaries(B) fields_upper {upperA,upperB,upperC,upperD}; fields_lower {lowerA,lowerB,lowerC,lowerD}; for i 1:4 up B.(fields_upper{i}); low B.(fields_lower{i}); if any(diff(up(:,1)) 0) || any(diff(low(:,1)) 0) error(%s或%s的x坐标未递增, fields_upper{i}, fields_lower{i}); end if any(up(:,2) low(:,2)) error(%s的上边界低于下边界, fields_upper{i}); end end disp(边界自检通过); end这个自检在课程设计答辩前很有用评审问“怎么保证边界数据可靠”时直接展示这个函数会很有说服力。另外它也能捕获复制粘贴时最容易犯的错——顶点顺序错位。5.3 分层统计技巧避免稀薄区域稀释结论最后给出一个我习惯用的分层验证技巧。当样本量较少或者数据集中在70到200 mg/dL时总体A区占比可能很高但低血糖区可能存在隐患却被平均掩盖。按参考血糖浓度分层统计能暴露问题。在Runme_EGA.m的统计部分之后追加edges [0 70 120 180 250 600]; labels {70, 70-120, 120-180, 180-250, 250}; for k 1:length(edges)-1 idx ref edges(k) ref edges(k1); if sum(idx) 0 local_A sum(zones(idx) A) / sum(idx) * 100; fprintf(区间 %s: n%d, A区%.1f%%\n, labels{k}, sum(idx), local_A); else fprintf(区间 %s: 无样本\n, labels{k}); end end这段代码把血糖分为五个临床相关区间分别计算A区比例。如果发现70 mg/dL区间A区比例急剧下降哪怕总体A区超过95%也说明设备在低血糖段存在系统性偏差。这种分层统计法多见于CGM传感器性能研究用来评价传感器在临床关键低值区的表现比单纯看总百分比可靠得多。配合资源里的output.png和操作录像可以快速验证自己的数据在不同区间下的表现从而决定是否需要调整边界参数或引入动态校准算法。本文还有配套的精品资源点击获取