隐式方程的matlab实现:酶动力学(Integrated Michaelis-Menten Equation)的无偏参数求解

神奇而强大的计算分析平台 Google搜索了一下,发现隐式方程的求参问题有不少解决方案,matlab, R都可以实现。基本思路与我用Excel求参的方法一样,先编写一个函数来实现隐式方程的显示化,matlab中要用到fzero或fsolve(Excel中用的是单变量求解插件),然后再用非线性求解函数进行参数拟合。下面是积分后的米曼方程(IMM)的matlab求参代码: 先构造IMM函数: function y=IMM(p,x) %输入参数 Vm=p(1); Km=p(2); y=zeros(size(x)); NN=length(x); %不输出提示 opt = optimset('display','off'); for i=1:NN %就是用这个方程来求出某个x(i)下的解,y的初值设定为x(i)*0.9 y(i)=fzero(@(y) x(i) -y-Vm*15+Km*log(x(i)/y), x(i)*0.9, opt); end end 上面的函数写好后,保存,最后只要用下面一行命令就可以求解参数了: X=lsqcurvefit(@(params, xdata) IMM(params, xdata),...…

你就是自己的建模师

Snoopy是谁的模型?模型是对事物本质的逼近,有时我们宁愿选择不那么逼近真实,因为这样我们才能更不受世俗约束的表达自己的真实本质。 ​        在高楼大厦林立的繁华城市中,一只小狗躺在自己的屋顶上,仰望天空,陪伴它的是一只忘记了季节的小鸟“糊涂塔克”。它幻想开飞机,开火箭,当超人,它从不担心是否需要考飞机驾照,是否要学航天知识。。。。。。这就是Snoopy! ​        我不知道为什么会选择药学这个专业,也不知道为什么读了硕士、博士、博士后。生活所迫?也许吧。那时候读书是唯一能让我生活不太依赖家人的方法,不用交 学费,还有生活补助。但内心那种冲动这么多年来一直与现实抗争着:本科时,我着迷QBasic的图形编程,泡在东大的计算机中心,忘记了时间;硕士时,我 迷上了汇编语言,抱着看雪的加密与解密在不能站直身体的春运列车上看得忘记一切;博士时,我学习了matlab, C++, perl,并用perl开发了一个专业网站,稍有盈利;博后时,我迷上了建模,把WinNonLin用的纯熟(有点吹牛啊),但我很少用它解决药学内的问 题(不务正业),如果有人让我帮忙分析药动数据,我建议他们还是去用3P87或DAS之类的软件去算吧。但是如果你问的是WinNonLin从没有人用过 的应用,我会比较感兴趣, 比如有个博士生让我模拟几个中药成分的化学反应动力学,或者如何解决有损失采样造成的系统偏差,或者养殖鱼苗中给药后药物多长时间可以清除等等,我会比较感兴趣。这就是有点叛逆的我吧! ​        昨天我把放下好久的NONMEM,WFN,XPOSE又配置了起来,因为这就是我现在的工作!这是一个结束,也是一个新的开始。 ​        分享一下自己在08年研究NONMEN的时候写的笔记,正因为有了这个,我很快就恢复了建模环境: ================我是可耐的分割线================ 用G95编译NONMEM6 ​        NONMEM6的源码可以用Compaq Visual Fortran 6(CVF)直接编译,没有问题,若是用其他编译器,比如G95,就需要对源代码稍微改动一下,因为不同的编译器对一些命令的使用有区别。 ​        G95是免费的Fortran编译器,安装后只有十几兆,在网上很易找到并下载。而CVF的安装文件近600兆,很难在网上找到下载,并且安装后,占硬盘空间也很大,至少200多兆。因此使用精巧、免费的G95编译器进行编译具有很实际的意义。 ​        整个编译的过程可以参考NONMEM6的CDSETUP6.BAT文件,至少稍微修改就可以对NONMEM目录中的*.for进行编译,G95编译后的文件后缀为*.o,而CVF编译的后缀为*.obj。我修改后的setupG95.bat如下(注意,先删除所有目录中的*.obj,以及nonmem.lib,并删除run目录): 1. 运行setupG95.bat前,需要将根目录下MAKEG95BAT中的第45和48行中的“.obj”改为“.o”; 2. NM\FLU.for中的第9行“CALL COMMITQQ(I)”,改为“CALL    FLUSH(I)”,否则编译后运行会出现错误“c:\nmvi1\NM\nonmem.lib(FLU.o):FLU.for:(.text+0xa): undefined reference to `commitqq_'    ”; 3. NM\BLKDAT.for中的第28行中的'con'改为'CONOUT$ ';…