24小时热门版块排行榜    

查看: 3932  |  回复: 8

cheng1378653

新虫 (初入文坛)

[求助] 用matlab里的nlinfit如何拟合分段函数

Sample Text
t=[0 0.0833 0.25 0.5 1 2];
c=[0 0.0833 1.2491 1.7025 0.2623 0.0757];
函数形式为c=A*(1-exp(-K1*t))       0<=t<=T
                c=A*(1-exp(-k1*T))*exp(-k2*(t-T))     t>=T
如何用脚本式或是函数式M文件拟合,求出A,k1,k2,(如果可以请划一下图形)
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖


感谢参与,应助指数 +1
fegg7502: 金币+1, 鼓励交流 2013-04-02 09:24:49
1stOpt更简单:
CODE:
Function c=if(t Data;
t=[0 0.0833 0.25 0.5 1 2];
c=[0 0.0833 1.2491 1.7025 0.2623 0.0757];

均方差(RMSE): 0.0282872879644411
残差平方和(SSE): 0.00480102396229929
相关系数(R): 0.999311354215736
相关系数之平方(R^2): 0.998623182664489
决定系数(DC): 0.998184285936899
卡方系数(Chi-Square): 0.0580621164226486
F统计(F-Statistic): 367.165006987895

参数        最佳估算
----------        -------------
t1        0.317859946002189
a        -0.0366738585259589
k1        -14.2282073650289
k2        3.70145325521581

jv1.jpg

2楼2013-03-29 08:50:59
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lgycjpcqu

金虫 (正式写手)

【答案】应助回帖


感谢参与,应助指数 +1
fegg7502: 金币+1, 鼓励交流 2013-04-02 09:24:59
引用回帖:
2楼: Originally posted by dingd at 2013-03-29 08:50:59
1stOpt更简单:

Function c=if(t<T1,A*(1-exp(-K1*t)),A*(1-exp(-k1*T1))*exp(-k2*(t-T1)));
Data;
t=;
c=;


均方差(RMSE): 0.0282872879644411
残差平方和(SSE): 0.00480102396229929
相关系数(R): ...

为啥我的1sopt运算结果不一致呢?
参数        最佳估算
----------        -------------
t1                 0.500000049987681
a                  -0.0568673833778535
k1                 -6.86274581292278
k2                 3.70079073587775
不是说1sopt得到的是唯一解吗?
3楼2013-03-29 10:38:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lgycjpcqu

金虫 (正式写手)

【答案】应助回帖

★ ★ ★ ★
cheng1378653: 回帖置顶 2013-03-29 18:09:29
cheng1378653: 金币+2, 有帮助 2013-03-30 22:02:11
fegg7502: 金币+2, 应助指数+1, 鼓励交流 2013-04-02 09:24:40
还是2楼的结果好,用刚才我用1sopt运算的结果作为初始值,用MATLAB进行拟合得到结果和1楼一致,看来非线性最优化还是1sopt强大啊
A =
-0.0367
k1 =
-14.2282
k2 =
3.7015
T =
0.3179
你需要MATLAB 代码的话如下
clear
tic
fun1=inline('beta(1).*(1-exp(-beta(2).*t))','beta','t');
fun2=inline('beta(1).*(1-exp(-beta(2).*beta(4)))*exp(-beta(3)*(t-beta(4)))','beta','t');
fun=@(beta,t)((t<0.5).*fun1(beta,t)+(t>=0.5).*fun2(beta,t));
t=[0 0.0833 0.25 0.5 1 2];
c=[0 0.0833 1.2491 1.7025 0.2623 0.0757];
beta=nlinfit(t,c,fun,[-0.05  -6.8 3.7 0.5]);
A=beta(1);
k1=beta(2);
k2=beta(3);
T=beta(4);
warning off all
c1=beta(1).*(1-exp(-beta(2).*t(1:4)));
c2=beta(1).*(1-exp(-beta(2).*beta(4))).*exp(-beta(3).*(t(4:6)-beta(4)));
cfit(1:4)=c1;
cfit(4:6)=c2;
plot(t,c,'*')
hold on;
plot(t,cfit,'r-')
toc
4楼2013-03-29 10:47:34
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖


fegg7502: 金币+1, 应助指数+1, 3ks 2013-04-02 09:25:13
引用回帖:
3楼: Originally posted by lgycjpcqu at 2013-03-29 10:38:11
为啥我的1sopt运算结果不一致呢?
参数        最佳估算
----------        -------------
t1                 0.500000049987681
a                  -0.0568673833778535
k1                 -6.86274581292278
k2                 3.70079073587775
...

如果你用的1stOpt老版本(如1.5),可是如下代码:
CODE:
Variable t,c;
Function c=if(t Data;
0        0
0.0833        0.0833
0.25        1.2491
0.5        1.7025
1        0.2623
2        0.0757

5楼2013-03-29 11:15:58
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lgycjpcqu

金虫 (正式写手)

引用回帖:
5楼: Originally posted by dingd at 2013-03-29 11:15:58
如果你用的1stOpt老版本(如1.5),可是如下代码:

Variable t,c;
Function c=if(t<T1,A*(1-exp(-K1*t)),A*(1-exp(-k1*T1))*exp(-k2*(t-T1)));
Data;
0        0
0.0833        0.0833
0.25        1.2491
0.5        1.7025
1        0. ...

这次对了,看来是我的版本太老了
6楼2013-03-29 13:53:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cheng1378653

新虫 (初入文坛)

引用回帖:
4楼: Originally posted by lgycjpcqu at 2013-03-29 10:47:34
还是2楼的结果好,用刚才我用1sopt运算的结果作为初始值,用MATLAB进行拟合得到结果和1楼一致,看来非线性最优化还是1sopt强大啊
A =
-0.0367
k1 =
-14.2282
k2 =
3.7015
T =
0.3179
你需要MATLAB 代码的 ...

你好,你的这个代码是默认t=0.5时取得最大值,而我是想要根据拟合的结果,求出取得最大值得T值,谢谢你的帮助
7楼2013-03-30 22:01:05
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lgycjpcqu

金虫 (正式写手)

【答案】应助回帖


fegg7502: 金币+1, 应助指数+1, 鼓励交流 2013-04-02 09:25:26
引用回帖:
7楼: Originally posted by cheng1378653 at 2013-03-30 22:01:05
你好,你的这个代码是默认t=0.5时取得最大值,而我是想要根据拟合的结果,求出取得最大值得T值,谢谢你的帮助...

fun=@(beta,t)((t<0.5).*fun1(beta,t)+(t>=0.5).*fun2(beta,t)); 这一句中的0.5应该换成beta(4),这个是一个参数,我忘改过来了,不好意思,不过和其它参数一样也必须有一个合理的初始值才能得到好的结果。
8楼2013-03-31 10:34:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cheng1378653

新虫 (初入文坛)

引用回帖:
8楼: Originally posted by lgycjpcqu at 2013-03-31 10:34:12
fun=@(beta,t)((t<0.5).*fun1(beta,t)+(t>=0.5).*fun2(beta,t)); 这一句中的0.5应该换成beta(4),这个是一个参数,我忘改过来了,不好意思,不过和其它参数一样也必须有一个合理的初始值才能得到好的结果。...

谢谢
9楼2013-03-31 12:25:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 cheng1378653 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 面上意见出来了 +4 黄鸟于飞Chao 2026-08-29 8/400 2026-08-29 14:50 by 黄鸟于飞Chao
[博后之家] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +4 ASdOkHsho7FD 2026-08-28 5/250 2026-08-29 14:24 by jCd0dEvKHShX
[考博] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +4 ASdOkHsho7FD 2026-08-28 7/350 2026-08-29 14:15 by jCd0dEvKHShX
[硕博家园] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +5 ASdOkHsho7FD 2026-08-28 9/450 2026-08-29 14:03 by jCd0dEvKHShX
[公派出国] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 ASdOkHsho7FD 2026-08-28 4/200 2026-08-29 14:02 by jCd0dEvKHShX
[找工作] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +3 ASdOkHsho7FD 2026-08-28 6/300 2026-08-29 11:51 by jCd0dEvKHShX
[教师之家] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +4 ASdOkHsho7FD 2026-08-28 5/250 2026-08-29 11:49 by jCd0dEvKHShX
[考研] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +4 ASdOkHsho7FD 2026-08-28 6/300 2026-08-29 11:40 by jCd0dEvKHShX
[找工作] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +4 ASdOkHsho7FD 2026-08-28 6/300 2026-08-29 11:29 by jCd0dEvKHShX
[考研] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +6 ASdOkHsho7FD 2026-08-28 9/450 2026-08-29 11:18 by jCd0dEvKHShX
[基金申请] 中青基了要发朋友圈吗? +3 349506619 2026-08-28 3/150 2026-08-29 10:49 by 木风9012
[基金申请] 国自然评审意见 +13 wangmingqi 2026-08-28 19/950 2026-08-29 10:22 by Poppy1104
[基金申请] 我就是申请一个面上项目而已,这评审意见是按照杰青的条件评的吧? +6 gouxfjh 2026-08-28 9/450 2026-08-29 09:59 by jklily
[教师之家] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +4 ASdOkHsho7FD 2026-08-28 7/350 2026-08-29 09:26 by G6APbkg8SA6w
[基金申请] 基金未中,这种答复是模板吗? +5 zhaosm1982 2026-08-27 6/300 2026-08-27 16:00 by lfy8008
[基金申请] 能否退出参与的面上项目解除限项 +23 koalala 2026-08-24 26/1300 2026-08-26 14:29 by 宝贝虫子
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 国合可查了 +3 paperzjh 2026-08-26 3/150 2026-08-26 10:41 by LemmonTr
[基金申请] 牛来!米来!面来! +8 beefly 2026-08-26 8/400 2026-08-26 08:37 by xuzhipiao
[基金申请] 没有任何消息-是不是就凉了 +9 图啦图啦 2026-08-24 10/500 2026-08-25 11:59 by 南海小哥
信息提示
请填处理意见