24小时热门版块排行榜    

CyRhmU.jpeg
查看: 646  |  回复: 11
当前主题已经存档。

hht666

金虫 (正式写手)

[交流] 【求助】用matlab中的高斯牛顿法或单纯形法进行模型嵌合

各位高手, 我最近碰到一问题,需要分析人体内的药物含量随时间变化,需要用上述方法进行模型嵌合, 只知道matlab 软件可以计算,想自己来时间不够,只能求助各位高手。
t        X
0.5        29701602
1        28952373
2        27578765
3        26202931
4        24616875
5        23031635
6        21218045
7        19869353
8        18343096
9        16833890
10        15504818
11        14335500
12        13110579
13        12048990
14        10987643
模型公式如下:X=A×e∧(-a×t)+B× e∧(-b×t)+C×e∧(-c×t)+D× e∧(-d×t)
  A,a, B,b,C,c, D,d 需满足如下条件。
A+B+C+D=29792400, A×a+B×b+C×c+D×d=0,
a,b,,c,d 均大于零,且 b>a,b>c,b>d;
请高手解决,30金币敬请笑纳。

[ Last edited by 点点风 on 2009-3-21 at 12:44 ]
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

fspdlh

金虫 (正式写手)

★ ★ ★ ★ ★
hht666(金币+5,VIP+0):多谢啊.有一个条件不满足,请帮忙再算算 3-21 12:02
以下是我用非线性拟合做的,结果与几个量的初值关系比较大,你能不能说一下各值大约在什么范围内啊,或者数据再多些,因为待定参数太多,相对来说数据有点少
abcd =

    0.0631    0.9828    0.0593    0.0650


ABCD =

  1.0e+009 *

    9.3546   -0.0019   -3.1108   -6.2121


r =

    0.9999



abcd =

    0.0949  184.3219   37.3425    0.2386


ABCD =

  1.0e+007 *

    4.3306    0.0086   -0.0452   -1.3148


r =

    0.9999

[ Last edited by fspdlh on 2009-3-20 at 19:53 ]
2楼2009-03-20 19:50:23
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

fspdlh

金虫 (正式写手)

★ ★
Doctorcbw(金币+2,VIP+0):谢谢参与 12-8 15:07
这是程序
datafit.m中:
function [abcd ABCD r ]=datafit(beta)
load data.txt
t=data(:,1);
X=data(:,2);
beta=nlinfit(t,X,'fun',beta);
a1=beta(1);
delta_b=beta(2);
c1=beta(3);
d1=beta(4);
A=beta(5);
B=beta(6);
a=a1^2;
c=c1^2;
d=d1^2;
b=max(max(a,c),d)+delta_b^2;
C=-(-B*d+A*a+B*b+29792400*d-A*d)/(c-d);
D=(-c*B+A*a+29792400*c-c*A+B*b)/(c-d);
X1=fun(beta,t);
r=corrcoef(X,X1);
r=r(1,2);
abcd=[a b c d];
ABCD=[A B C D];

fun.m中:
function X=fun(beta,t)
a1=beta(1);
delta_b=beta(2);
c1=beta(3);
d1=beta(4);
A=beta(5);
B=beta(6);
a=a1^2;
c=c1^2;
d=d1^2;
b=max(max(a,c),d)+delta_b^2;
C=-(-B*d+A*a+B*b+29792400*d-A*d)/(c-d);
D=(-c*B+A*a+29792400*c-c*A+B*b)/(c-d);
X=A.*exp(-a.*t)+B.*exp(-b.*t)+C.*exp(-c.*t)+D.*exp(-d.*t);

data.txt中:
0.5        29701602
1        28952373
2        27578765
3        26202931
4        24616875
5        23031635
6        21218045
7        19869353
8        18343096
9        16833890
10        15504818
11        14335500
12        13110579
13        12048990
14        10987643
3楼2009-03-20 19:58:10
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

fspdlh

金虫 (正式写手)


点点风(金币+1,VIP+0):感谢参与讨论,加分鼓励! 3-21 12:45
我编的函数的参数为
sqrt(a) sqrt(b-max(a,c,d)) sqrt(c) sqrt(d) A B 六个
由这六个参数来确定a b c d A B C D

[ Last edited by fspdlh on 2009-3-20 at 20:02 ]
4楼2009-03-20 19:59:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hht666

金虫 (正式写手)

谢谢fspdlh,   提供的数据第一组与实际相差太远.第二组不满足A×a+B×b+C×c+D×d=0,
我自己用EXCEL拼凑了几个数据,初算了一下, A=125025, a=1.18,  B=-110170, b=21.1, C=-1113011,c=0.196, D=31961557,d=0.07498,
caculated
29848802
28782013
26780451
24921904
23187854
21568459
20056542
18645797
17330281
16104289
14962330
13899139
12909690
11989199
11133132
A+B+C+D=29792400可以有一定误差, 但A×a+B×b+C×c+D×d=0需要满足,我的数据就是等于1269,有点偏差, 其中数据范围还是符合要求的.我手上没MATLAB软件,请好心人重新再算算.
5楼2009-03-21 11:56:42
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

fspdlh

金虫 (正式写手)

我用你的初值算了一下,计算前r=0.9948,计算后:
>> [abcd ABCD r ]=datafit()

X=

29536194.64
29008722.91
27667836.28
26159305.46
24573391.62
22962735.76
21365229.00
19808246.42
18310992.69
16886312.45
15542139.46
14282656.69
13109221.39
12021099.45
11016044.88



abcd =

    0.0537    2.5824    0.1597    0.1595


ABCD =

  1.0e+010 *

    0.0010   -0.0000   -2.2485    2.2505


r =

    0.9999

>> vpa(abcd,10)

ans =

[ .5373300786e-1,    2.582437202,    .1596700667,    .1595465853]


>> vpa(ABCD,10)

ans =

[     9947271.618,    -381453.9706, -.2248515653e11,  .2250538312e11]


>> vpa(r,10)

ans =

.9999282088


>> sum(ABCD)

ans =

  2.9792e+007

>> vpa(sum(ABCD))

ans =

29792400.000007629394531250000000


>> abcd*ABCD'

ans =

  9.5367e-007

>>
你这结果相关性已经这么高了,再变动A B C D的总和就没必要了,你说我的结果那个式子不满足,这是显示位数的问题。
1.bmp

[ Last edited by fspdlh on 2009-3-21 at 13:44 ]
6楼2009-03-21 12:54:48
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hht666

金虫 (正式写手)

我用残差法,假设在T 足够大时,前面三项都趋近于零,就可以得出那把D 和  d固定下来.D=EXP(17.449)=37844643, d=0.0882.请再算算情况怎么样.

[ Last edited by hht666 on 2009-3-21 at 13:54 ]
7楼2009-03-21 13:49:48
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

fspdlh

金虫 (正式写手)

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
hht666(金币+25,VIP+0):非常感谢 3-21 18:00
29523786.32
29036825.91
27771350.75
26262045.13
24633185.40
22967894.59
21320683.29
19726250.87
18205655.33
16770632.53
15426613.99
14174830.41
13013772.85
11940203.43
10949850.22

>> vpa(abcd,10)

ans =

[    .3505511676,    1211.256494,    1171.246711, .8820000000e-1]


>> vpa(ABCD,10)

ans =

[ -7969456.956,  2418493.915, -2492663.807,    37844643.]


>> vpa(r,10)

ans =

.9998853657
8楼2009-03-21 16:17:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hht666

金虫 (正式写手)

相加后还是不为零,也许是模型有问题吧,要不就是实验数据太少,没有设计好,好多能定下来的值没有搞出来
9楼2009-03-21 17:59:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

fspdlh

金虫 (正式写手)

是我刚才搞错了,这回才对
29523786.31
29036825.90
27771350.75
26262045.14
24633185.41
22967894.59
21320683.29
19726250.87
18205655.33
16770632.53
15426613.99
14174830.41
13013772.85
11940203.44
10949850.23



ans =

.9998853657



abcd =

  1.0e+003 *

    0.0004    2.8830    2.7608    0.0001


ABCD =

  1.0e+007 *

   -0.7969    0.1867   -0.1949    3.7845


r =

    0.9999

>> sum(ABCD)

ans =

    29792400

>> abcd*ABCD'

ans =

  5.4948e-008

这是舍入误差造成的

>> vpa(abcd,10)

ans =

[    .3505511687,    2883.001975,    2760.841200, .8820000000e-1]


>> vpa(ABCD,10)

ans =

[ -7969456.980,  1866514.535, -1949300.555,    37844643.]


>>

[ Last edited by fspdlh on 2009-3-21 at 18:20 ]
10楼2009-03-21 18:14:10
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 hht666 的主题更新
普通表情 高级回复(可上传附件)
信息提示
请填处理意见