| 查看: 1504 | 回复: 4 | |||
| 【悬赏金币】回答本帖问题，作者happy_rolly将赠送您 5 个金币 | |||
happy_rolly新虫 (初入文坛)
|
[求助]
用matlab求动力学参数代码 求指教已有1人参与
|
||
|
附件中的文章 跟我的研究方向比较相关 所以想重现一下作者的结果, 首先建立传染病动力学SEIQR模型（分别表示易感者、潜伏期、感染未住院、感染并住院、康复），然后写出了模型的常微分方程组，共有8个参数，能给出数值的有μ和δ，另外六个参数给出了初始值范围，并且需要估计 S E I Q R五个未知数中，S Q R的初始值可以给出，E I是给定范围 ，要估计 最优参数的估计采用的是fminsearch，目标函数是MSS=(log2 (1+data per month) 2-log2 (1+ infectious and hospitalized on day 1+...+on day 28(or 29,30,31)))2 目标函数表达式中data per month就是表1中 每月发病人数/10, 另外 infectious and hospitalized on day 1+...+on day 28(or 29,30,31)应该是来源于常微分方程组的解 目前为止我只写出了常微分方程组的代码，见li.m 至于如何构建MSS这样的目标函数，以及如何利用常微分方程组的解和每月发病数据来求MSS最小值 并且要得出另外六个参数值和E I的初始值 我一点头绪也没有 求大神指教 必有重谢 最好能给出matlab代码 医学生学matlab好困难，整个课题组就我一个人在做这个 *li.m function f=li(t,y) u=3.9139*10^-5; f=[u*(y(1)+y(2)+y(3)+y(4)+y(5))-b*(y(3)+y(4))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5))+r*y(5)-u*y(1); b*(y(3)+y(4))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5))-a*y(2)-u*y(2); a*(1-p)*y(2)-m*y(3)-u*y(3); a*p*y(2)-n*y(4)-u*y(4); m*y(3)+n*y(4)-r*y(5)-u*y(5)];     常微分方程组.png  表1 每月发病数据.png |
回复此楼» 本帖附件资源列表
-
欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：libolin3@tal.com - 附件 1 : Li-2014-Modeling_and_Prevent.pdf
2018-04-18 23:20:19, 315.08 K
» 收录本帖的淘贴专辑推荐
 常微分方程 |
» 猜你喜欢
求助将cif文件转化为jade
已经有1人回复
-
金属氧化物薄膜空气退火和真空退火对表面粗糙度的影响
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有58人回复
-
Apex2单晶测试
已经有0人回复
-
结构建模
已经有0人回复
-
声子谱计算问题
已经有2人回复
-
已结束
已经有23人回复
-
Error in routine lambda (100): wrong # or too many modes
已经有0人回复
-
祈福国自然面上
已经有55人回复
-
请问castep计算出的能带图横轴数据缩放因子怎么计算
已经有1人回复
happy_rolly
新虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 57.5
- 帖子: 15
- 在线: 6.1小时
- 虫号: 4336034
- 注册: 2016-01-05
- 性别: MM
- 专业: 流行病学方法与卫生统计
2楼2018-04-20 14:51:49
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15356.8
- 散金: 101
- 红花: 234
- 帖子: 3409
- 在线: 1222.7小时
- 虫号: 291104
- 注册: 2006-10-28
3楼2018-04-20 17:09:41
happy_rolly
新虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 57.5
- 帖子: 15
- 在线: 6.1小时
- 虫号: 4336034
- 注册: 2016-01-05
- 性别: MM
- 专业: 流行病学方法与卫生统计
|
程序我仿照写出来了 用的是2010年数据 用fmincon优化（文章中用fminsearch,但是我不知道用这个该怎么加lb和ub） 但是估计的ｋ跟文章中２０１０年相差很大 k0 = [0.5 0.5 0.5 0.5 0.5 0.5]; % 参数初值我是在lb和ub范围内随便选的,是不是会对结果有很大影响？ 文章中E (0) in [1.8 × 104,1.4 × 107] I(0) in [1.8 × 103,1.8 × 105], 但是我不知道怎么将ｙ的初值设成一个范围，于是用2.043*10^4和 1.8001*10^3代替了，有什么办法可以改一下吗 %6个参数 运行不报错 但是参数值和文章相差很大 结果拟合也不好 %function k1k2k3 format long clear all clc tspan = [1:12]; k0 = [0.5 0.5 0.5 0.5 0.5 0.5]; % 参数初值 lb = [0.0001 0.0001 0.0001 0.1 0.01 0.1]; % 参数下限 ub = [50 1 1 1 1 1]; % 参数上限 x0 = [1.4*10^8 2.043*10^4 1.8001*10^3 1212 0]; %y初值 ExpData=[ 1 3756.7 2 2386.2 3 7775.6 4 24860.9 5 35434.7 6 34310 7 26126.3 8 11909.6 9 10165.4 10 8761.2 11 7959.1 12 6087.9 ]; t=ExpData(:,1); yexp = ExpData(:,2); % yexp: CDC数据 % 使用函数fmincon()进行参数估计 [k,fval,flag] = fmincon(@MSS,k0,[],[],[],[],lb,ub,[],[],x0,yexp); fprintf('\n使用函数fmincon()估计得到的参数值为:\n') fprintf('\tk1 = %.4f\n',k(1)) fprintf('\tk2 = %.4f\n',k(2)) fprintf('\tk3 = %.4f\n',k(3)) fprintf('\tk4 = %.4f\n',k(4)) fprintf('\tk5 = %.4f\n',k(5)) fprintf('\tk6 = %.4f\n',k(6)) fprintf(' The sum of the squares is: %.1e\n\n',fval) %% =======绘图显示结果======= [tt xx] = ode45(@li,tspan,x0,[],k); plot(t,yexp(:,1),'r.-') hold on plot(tt,xx(:,4),'ro-') legend('exp y1', 'pre y1') xlabel('t/月') ylabel('人数') function f = MSS(k,x0,yexp) % fmincon（） tspan = [1:12]; [t x] = ode45(@li,tspan,x0,[],k); y(:,1) = x(:,4); f = sum((log2(1+ yexp(:,1))-log2(1+ y(:,1))).^2); %微分方程组 function dydt=li(t,y,k) u=3.9139*10^-5; dydt=[(u*(y(1)+y(2)+y(3)+y(4)+y(5))-k(1)*(y(3)+y(4))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5))+k(2)*y(5)-u*y(1)) (k(1)*(y(3)+y(4))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5))-k(3)*y(2)-u*y(2)) (k(3)*(1-k(5))*y(2)-k(4)*y(3)-u*y(3)) (k(3)*k(5)*y(2)-k(6)*y(4)-u*y(4)) (k(4)*y(3)+k(6)*y(4)-k(2)*y(5)-u*y(5))]; 使用函数fmincon()估计得到的参数值为: k1 = 0.0001 k2 = 0.3354 k3 = 0.4065 k4 = 0.8597 k5 = 1.0000 k6 = 0.1000 The sum of the squares is: 1.2e+01 |
4楼2018-04-25 05:39:35
happy_rolly
新虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 57.5
- 帖子: 15
- 在线: 6.1小时
- 虫号: 4336034
- 注册: 2016-01-05
- 性别: MM
- 专业: 流行病学方法与卫生统计
5楼2018-04-25 05:43:18
10