数学研发网设为首页收藏本站
切换到窄版

 找回密码
 欢迎注册
数学研发论坛»论坛 › 【数学研究】 › 数学软件 › NIST非线性拟合-----接力赛
123下一页
返回列表 发新帖
查看: 27612|回复: 27

[讨论] NIST非线性拟合-----接力赛

[复制链接]
wayne
发表于 2011-2-21 21:02:47 | 显示全部楼层 |阅读模式

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要 登录 才可以下载或查看，没有账号？欢迎注册

×
原始数据以及 待拟合的 目标公式都在该链接里，
http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml

我们不妨来 接力攻破一下~~
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复

使用道具 举报

forcal
发表于 2011-2-21 21:08:22 | 显示全部楼层
给楼主一个建议：能否一个一个地给出题目的数据及中文说明？即给出一个，解决一个，再给出一个。

像我这样英文水平的人就也可以进行尝试了。也能方便更多像我这样的人。

评分

参与人数 1鲜花 +2 收起 理由
wayne + 2 好建议！1

查看全部评分

毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

wayne
 楼主| 发表于 2011-2-21 21:26:41 | 显示全部楼层

低难度组第1题： Misra1a

本帖最后由 wayne 于 2011-2-21 21:30 编辑

参考链接：http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml
题目背景意义： 关于牙医研究中单分子的吸收率。

=================================
待拟合的公式： y = b1*(1-exp[-b2*x])
其中b1,b2是参数，自变量x是体积，因变量y是压强。共14组观测数据：
Data:   y               x
      10.07E0      77.6E0
      14.73E0     114.9E0
      17.94E0     141.1E0
      23.93E0     190.8E0
      29.61E0     239.9E0
      35.18E0     289.0E0
      40.02E0     332.8E0
      44.82E0     378.4E0
      50.76E0     434.8E0
      55.05E0     477.3E0
      61.01E0     536.8E0
      66.40E0     593.1E0
      75.47E0     689.1E0
      81.78E0     760.0E0
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

forcal
发表于 2011-2-21 21:28:25 | 显示全部楼层
第一题Forcal求解：
  1. !using["fcopt","math"];
  2. init(::Array,max)=
  3. {
  4.     max=14,
  5.     Array=arrayinitns{max,2 :
  6.         "      10.07E0      77.6E0
  7.       14.73E0     114.9E0
  8.       17.94E0     141.1E0
  9.       23.93E0     190.8E0
  10.       29.61E0     239.9E0
  11.       35.18E0     289.0E0
  12.       40.02E0     332.8E0
  13.       44.82E0     378.4E0
  14.       50.76E0     434.8E0
  15.       55.05E0     477.3E0
  16.       61.01E0     536.8E0
  17.       66.40E0     593.1E0
  18.       75.47E0     689.1E0
  19.       81.78E0     760.0E0"
  20.     }.free()
  21. };
  22. f(b1, b2 : i,s,x,y,z : Array,max)=
  23. {
  24.     s=0,i=0,(i<max).while{
  25.         y=Array[i,0], x=Array[i,1],
  26.         s=s+[ b1*(1-exp[-b2*x])-y]^2,
  27.         i++
  28.     },
  29.     s
  30. };
  31. Opt[HFor("f")];
复制代码
结果（前面是最优参数，最后一个数是误差，Forcal结果基本都是这样的）：
238.9421251876012         5.501564435330968e-004    0.1245513889445084

从Forcal的角度，难度比较小。
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

wayne
 楼主| 发表于 2011-2-22 12:01:14 | 显示全部楼层
3# wayne
用Mathematica解：
1、直接从网络抓取数据，不用拷贝：
  1. data = Import["http://www.itl.nist.gov/div898/strd/nls/data/LINKS/DATA/Misra1a.dat", "Data"];
复制代码
aa = Cases[data, {_Real, _Real}]
2、数据有理化
  1. xy = Cases[Rationalize[aa], {x_, y_} -> {y/100, x/10}]
复制代码
3、拟合参数值：
  1. {10 b1, b2/100} /. FindFit[xy, b1*(1 - Exp[-b2*x]), {b1, b2}, x, WorkingPrecision -> 20]
复制代码
===================================================
{238.94212917913832416, 0.00055015643180517008397}
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

forcal
发表于 2011-2-22 18:22:45 | 显示全部楼层
即便没有用过Mathematica，也开眼界，长见识。Mathematica的优化能力在软件四大家中恐怕要居首位吧？个人浅陋的理解，虽然都没有用过。

该推出第二题了吧？速度稍快些也无妨，呵呵。
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

forcal
发表于 2011-2-22 19:48:27 | 显示全部楼层
英文看起来还是比较费劲的，我来第二题，错了大家给指出哦。

低难度组第2题： Chwirut2
参考链接：http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml
=================================
待拟合的公式： y = exp(-b1*x)/(b2+b3*x)
其中b1,b2,b3是拟合参数。共54组观测数据：
  1. Data:  y             x
  2.       92.9000E0     0.500E0
  3.       57.1000E0     1.000E0
  4.       31.0500E0     1.750E0
  5.       11.5875E0     3.750E0
  6.        8.0250E0     5.750E0
  7.       63.6000E0     0.875E0
  8.       21.4000E0     2.250E0
  9.       14.2500E0     3.250E0
  10.        8.4750E0     5.250E0
  11.       63.8000E0     0.750E0
  12.       26.8000E0     1.750E0
  13.       16.4625E0     2.750E0
  14.        7.1250E0     4.750E0
  15.       67.3000E0     0.625E0
  16.       41.0000E0     1.250E0
  17.       21.1500E0     2.250E0
  18.        8.1750E0     4.250E0
  19.       81.5000E0      .500E0
  20.       13.1200E0     3.000E0
  21.       59.9000E0      .750E0
  22.       14.6200E0     3.000E0
  23.       32.9000E0     1.500E0
  24.        5.4400E0     6.000E0
  25.       12.5600E0     3.000E0
  26.        5.4400E0     6.000E0
  27.       32.0000E0     1.500E0
  28.       13.9500E0     3.000E0
  29.       75.8000E0      .500E0
  30.       20.0000E0     2.000E0
  31.       10.4200E0     4.000E0
  32.       59.5000E0      .750E0
  33.       21.6700E0     2.000E0
  34.        8.5500E0     5.000E0
  35.       62.0000E0      .750E0
  36.       20.2000E0     2.250E0
  37.        7.7600E0     3.750E0
  38.        3.7500E0     5.750E0
  39.       11.8100E0     3.000E0
  40.       54.7000E0      .750E0
  41.       23.7000E0     2.500E0
  42.       11.5500E0     4.000E0
  43.       61.3000E0      .750E0
  44.       17.7000E0     2.500E0
  45.        8.7400E0     4.000E0
  46.       59.2000E0      .750E0
  47.       16.3000E0     2.500E0
  48.        8.6200E0     4.000E0
  49.       81.0000E0      .500E0
  50.        4.8700E0     6.000E0
  51.       14.6200E0     3.000E0
  52.       81.7000E0      .500E0
  53.       17.1700E0     2.750E0
  54.       81.3000E0      .500E0
  55.       28.9000E0     1.750E0
复制代码
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

forcal
发表于 2011-2-22 19:52:13 | 显示全部楼层
第二题Forcal求解：
  1. !using["fcopt","math"];
  2. init(::Array,max)=
  3. {
  4.     max=54,
  5.     Array=arrayinitns{max,2 :
  6.         "      92.9000E0     0.500E0
  7.       57.1000E0     1.000E0
  8.       31.0500E0     1.750E0
  9.       11.5875E0     3.750E0
  10.        8.0250E0     5.750E0
  11.       63.6000E0     0.875E0
  12.       21.4000E0     2.250E0
  13.       14.2500E0     3.250E0
  14.        8.4750E0     5.250E0
  15.       63.8000E0     0.750E0
  16.       26.8000E0     1.750E0
  17.       16.4625E0     2.750E0
  18.        7.1250E0     4.750E0
  19.       67.3000E0     0.625E0
  20.       41.0000E0     1.250E0
  21.       21.1500E0     2.250E0
  22.        8.1750E0     4.250E0
  23.       81.5000E0      .500E0
  24.       13.1200E0     3.000E0
  25.       59.9000E0      .750E0
  26.       14.6200E0     3.000E0
  27.       32.9000E0     1.500E0
  28.        5.4400E0     6.000E0
  29.       12.5600E0     3.000E0
  30.        5.4400E0     6.000E0
  31.       32.0000E0     1.500E0
  32.       13.9500E0     3.000E0
  33.       75.8000E0      .500E0
  34.       20.0000E0     2.000E0
  35.       10.4200E0     4.000E0
  36.       59.5000E0      .750E0
  37.       21.6700E0     2.000E0
  38.        8.5500E0     5.000E0
  39.       62.0000E0      .750E0
  40.       20.2000E0     2.250E0
  41.        7.7600E0     3.750E0
  42.        3.7500E0     5.750E0
  43.       11.8100E0     3.000E0
  44.       54.7000E0      .750E0
  45.       23.7000E0     2.500E0
  46.       11.5500E0     4.000E0
  47.       61.3000E0      .750E0
  48.       17.7000E0     2.500E0
  49.        8.7400E0     4.000E0
  50.       59.2000E0      .750E0
  51.       16.3000E0     2.500E0
  52.        8.6200E0     4.000E0
  53.       81.0000E0      .500E0
  54.        4.8700E0     6.000E0
  55.       14.6200E0     3.000E0
  56.       81.7000E0      .500E0
  57.       17.1700E0     2.750E0
  58.       81.3000E0      .500E0
  59.       28.9000E0     1.750E0"
  60.     }.free()
  61. };
  62. f(b1, b2, b3 : i,s,x,y : Array,max)=
  63. {
  64.     s=0,i=0,(i<max).while{
  65.         y=Array[i,0], x=Array[i,1],
  66.         s=s+[ exp(-b1*x)/(b2+b3*x)-y]^2,
  67.         i++
  68.     },
  69.     s
  70. };
  71. Opt[HFor("f")];
复制代码
结果：
0.1665774524075732        5.165354127091879e-003    1.214996964994243e-002    513.0480294504323
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

forcal
发表于 2011-2-22 20:02:14 | 显示全部楼层
低难度组第3题： Chwirut1
参考链接：http://www.itl.nist.gov/div898/strd/nls/data/chwirut1.shtml
=================================
待拟合的公式： exp[-b1*x]/(b2+b3*x)
其中b1,b2,b3是拟合参数。共214组观测数据：
  1. Data:  y            x
  2.      92.9000E0     0.5000E0
  3.      78.7000E0     0.6250E0
  4.      64.2000E0     0.7500E0
  5.      64.9000E0     0.8750E0
  6.      57.1000E0     1.0000E0
  7.      43.3000E0     1.2500E0
  8.      31.1000E0     1.7500E0
  9.      23.6000E0     2.2500E0
  10.      31.0500E0     1.7500E0
  11.      23.7750E0     2.2500E0
  12.      17.7375E0     2.7500E0
  13.      13.8000E0     3.2500E0
  14.      11.5875E0     3.7500E0
  15.       9.4125E0     4.2500E0
  16.       7.7250E0     4.7500E0
  17.       7.3500E0     5.2500E0
  18.       8.0250E0     5.7500E0
  19.      90.6000E0     0.5000E0
  20.      76.9000E0     0.6250E0
  21.      71.6000E0     0.7500E0
  22.      63.6000E0     0.8750E0
  23.      54.0000E0     1.0000E0
  24.      39.2000E0     1.2500E0
  25.      29.3000E0     1.7500E0
  26.      21.4000E0     2.2500E0
  27.      29.1750E0     1.7500E0
  28.      22.1250E0     2.2500E0
  29.      17.5125E0     2.7500E0
  30.      14.2500E0     3.2500E0
  31.       9.4500E0     3.7500E0
  32.       9.1500E0     4.2500E0
  33.       7.9125E0     4.7500E0
  34.       8.4750E0     5.2500E0
  35.       6.1125E0     5.7500E0
  36.      80.0000E0     0.5000E0
  37.      79.0000E0     0.6250E0
  38.      63.8000E0     0.7500E0
  39.      57.2000E0     0.8750E0
  40.      53.2000E0     1.0000E0
  41.      42.5000E0     1.2500E0
  42.      26.8000E0     1.7500E0
  43.      20.4000E0     2.2500E0
  44.      26.8500E0     1.7500E0
  45.      21.0000E0     2.2500E0
  46.      16.4625E0     2.7500E0
  47.      12.5250E0     3.2500E0
  48.      10.5375E0     3.7500E0
  49.       8.5875E0     4.2500E0
  50.       7.1250E0     4.7500E0
  51.       6.1125E0     5.2500E0
  52.       5.9625E0     5.7500E0
  53.      74.1000E0     0.5000E0
  54.      67.3000E0     0.6250E0
  55.      60.8000E0     0.7500E0
  56.      55.5000E0     0.8750E0
  57.      50.3000E0     1.0000E0
  58.      41.0000E0     1.2500E0
  59.      29.4000E0     1.7500E0
  60.      20.4000E0     2.2500E0
  61.      29.3625E0     1.7500E0
  62.      21.1500E0     2.2500E0
  63.      16.7625E0     2.7500E0
  64.      13.2000E0     3.2500E0
  65.      10.8750E0     3.7500E0
  66.       8.1750E0     4.2500E0
  67.       7.3500E0     4.7500E0
  68.       5.9625E0     5.2500E0
  69.       5.6250E0     5.7500E0
  70.      81.5000E0      .5000E0
  71.      62.4000E0      .7500E0
  72.      32.5000E0     1.5000E0
  73.      12.4100E0     3.0000E0
  74.      13.1200E0     3.0000E0
  75.      15.5600E0     3.0000E0
  76.       5.6300E0     6.0000E0
  77.      78.0000E0      .5000E0
  78.      59.9000E0      .7500E0
  79.      33.2000E0     1.5000E0
  80.      13.8400E0     3.0000E0
  81.      12.7500E0     3.0000E0
  82.      14.6200E0     3.0000E0
  83.       3.9400E0     6.0000E0
  84.      76.8000E0      .5000E0
  85.      61.0000E0      .7500E0
  86.      32.9000E0     1.5000E0
  87.      13.8700E0     3.0000E0
  88.      11.8100E0     3.0000E0
  89.      13.3100E0     3.0000E0
  90.       5.4400E0     6.0000E0
  91.      78.0000E0      .5000E0
  92.      63.5000E0      .7500E0
  93.      33.8000E0     1.5000E0
  94.      12.5600E0     3.0000E0
  95.       5.6300E0     6.0000E0
  96.      12.7500E0     3.0000E0
  97.      13.1200E0     3.0000E0
  98.       5.4400E0     6.0000E0
  99.      76.8000E0      .5000E0
  100.      60.0000E0      .7500E0
  101.      47.8000E0     1.0000E0
  102.      32.0000E0     1.5000E0
  103.      22.2000E0     2.0000E0
  104.      22.5700E0     2.0000E0
  105.      18.8200E0     2.5000E0
  106.      13.9500E0     3.0000E0
  107.      11.2500E0     4.0000E0
  108.       9.0000E0     5.0000E0
  109.       6.6700E0     6.0000E0
  110.      75.8000E0      .5000E0
  111.      62.0000E0      .7500E0
  112.      48.8000E0     1.0000E0
  113.      35.2000E0     1.5000E0
  114.      20.0000E0     2.0000E0
  115.      20.3200E0     2.0000E0
  116.      19.3100E0     2.5000E0
  117.      12.7500E0     3.0000E0
  118.      10.4200E0     4.0000E0
  119.       7.3100E0     5.0000E0
  120.       7.4200E0     6.0000E0
  121.      70.5000E0      .5000E0
  122.      59.5000E0      .7500E0
  123.      48.5000E0     1.0000E0
  124.      35.8000E0     1.5000E0
  125.      21.0000E0     2.0000E0
  126.      21.6700E0     2.0000E0
  127.      21.0000E0     2.5000E0
  128.      15.6400E0     3.0000E0
  129.       8.1700E0     4.0000E0
  130.       8.5500E0     5.0000E0
  131.      10.1200E0     6.0000E0
  132.      78.0000E0      .5000E0
  133.      66.0000E0      .6250E0
  134.      62.0000E0      .7500E0
  135.      58.0000E0      .8750E0
  136.      47.7000E0     1.0000E0
  137.      37.8000E0     1.2500E0
  138.      20.2000E0     2.2500E0
  139.      21.0700E0     2.2500E0
  140.      13.8700E0     2.7500E0
  141.       9.6700E0     3.2500E0
  142.       7.7600E0     3.7500E0
  143.       5.4400E0     4.2500E0
  144.       4.8700E0     4.7500E0
  145.       4.0100E0     5.2500E0
  146.       3.7500E0     5.7500E0
  147.      24.1900E0     3.0000E0
  148.      25.7600E0     3.0000E0
  149.      18.0700E0     3.0000E0
  150.      11.8100E0     3.0000E0
  151.      12.0700E0     3.0000E0
  152.      16.1200E0     3.0000E0
  153.      70.8000E0      .5000E0
  154.      54.7000E0      .7500E0
  155.      48.0000E0     1.0000E0
  156.      39.8000E0     1.5000E0
  157.      29.8000E0     2.0000E0
  158.      23.7000E0     2.5000E0
  159.      29.6200E0     2.0000E0
  160.      23.8100E0     2.5000E0
  161.      17.7000E0     3.0000E0
  162.      11.5500E0     4.0000E0
  163.      12.0700E0     5.0000E0
  164.       8.7400E0     6.0000E0
  165.      80.7000E0      .5000E0
  166.      61.3000E0      .7500E0
  167.      47.5000E0     1.0000E0
  168.      29.0000E0     1.5000E0
  169.      24.0000E0     2.0000E0
  170.      17.7000E0     2.5000E0
  171.      24.5600E0     2.0000E0
  172.      18.6700E0     2.5000E0
  173.      16.2400E0     3.0000E0
  174.       8.7400E0     4.0000E0
  175.       7.8700E0     5.0000E0
  176.       8.5100E0     6.0000E0
  177.      66.7000E0      .5000E0
  178.      59.2000E0      .7500E0
  179.      40.8000E0     1.0000E0
  180.      30.7000E0     1.5000E0
  181.      25.7000E0     2.0000E0
  182.      16.3000E0     2.5000E0
  183.      25.9900E0     2.0000E0
  184.      16.9500E0     2.5000E0
  185.      13.3500E0     3.0000E0
  186.       8.6200E0     4.0000E0
  187.       7.2000E0     5.0000E0
  188.       6.6400E0     6.0000E0
  189.      13.6900E0     3.0000E0
  190.      81.0000E0      .5000E0
  191.      64.5000E0      .7500E0
  192.      35.5000E0     1.5000E0
  193.      13.3100E0     3.0000E0
  194.       4.8700E0     6.0000E0
  195.      12.9400E0     3.0000E0
  196.       5.0600E0     6.0000E0
  197.      15.1900E0     3.0000E0
  198.      14.6200E0     3.0000E0
  199.      15.6400E0     3.0000E0
  200.      25.5000E0     1.7500E0
  201.      25.9500E0     1.7500E0
  202.      81.7000E0      .5000E0
  203.      61.6000E0      .7500E0
  204.      29.8000E0     1.7500E0
  205.      29.8100E0     1.7500E0
  206.      17.1700E0     2.7500E0
  207.      10.3900E0     3.7500E0
  208.      28.4000E0     1.7500E0
  209.      28.6900E0     1.7500E0
  210.      81.3000E0      .5000E0
  211.      60.9000E0      .7500E0
  212.      16.6500E0     2.7500E0
  213.      10.0500E0     3.7500E0
  214.      28.9000E0     1.7500E0
  215.      28.9500E0     1.7500E0
复制代码
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

forcal
发表于 2011-2-22 20:05:49 | 显示全部楼层
第三题Forcal求解：
  1. !using["fcopt","math"];
  2. init(::Array,max)=
  3. {
  4.     max=214,
  5.     Array=arrayinitns{max,2 :
  6.         "     92.9000E0     0.5000E0
  7.      //（此处略去10000字）
  8.      28.9500E0     1.7500E0"
  9.     }.free()
  10. };
  11. f(b1, b2, b3 : i,s,x,y : Array,max)=
  12. {
  13.     s=0,i=0,(i<max).while{
  14.         y=Array[i,0], x=Array[i,1],
  15.         s=s+[ exp[-b1*x]/(b2+b3*x)-y]^2,
  16.         i++
  17.     },
  18.     s
  19. };
  20. Opt[HFor("f")];
复制代码
结果：
0.1902757520340506        6.131347827474484e-003    1.053102614699184e-002    2384.477139692282
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

123下一页
返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-4-25 19:17 , Processed in 0.132647 second(s), 18 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表