应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 专业学科区 » 数学 » 概率论与统计 » 关于用Maple或Mathematica实现符号多项式拟合
232/3返回列表上一页123下一页
查看: 2436  |  回复: 22
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

赵梦92

木虫 (著名写手)

【答案】应助回帖

maple擅长符号运算，但不擅长矩阵处理，而matlab新版融合maple符号运算以后，并两者优势兼而有之，用matlab编码更简捷一些

发自小木虫Android客户端
赞一下 回复此楼
11楼2016-01-08 13:20:35
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xzczd

木虫 (小有名气)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
loujing: 金币+40, ★★★很有帮助 2016-01-08 17:04:29
看看楼上的某些发言就觉得Mathematica推广工作任重而道远……楼主你这个问题是可以用Mathematica推导的。不过在给出解法前得先和你讨论清楚一件事：你给出的
的表达式是不是错了。我手算和用软件推导（当然编程里面混入了我的分析，或许是我的分析本身不对）的结果都是：
赞一下(1人) 回复此楼
小木虫Mathematica版块已毁（当然原本也不咋的），建议大家前往百度贴吧或Stackexchange。
12楼2016-01-08 14:38:05
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xzczd

木虫 (小有名气)

【答案】应助回帖

引用回帖:
12楼: Originally posted by xzczd at 2016-01-08 14:38:05
看看楼上的某些发言就觉得Mathematica推广工作任重而道远……楼主你这个问题是可以用Mathematica推导的。不过在给出解法前得先和你讨论清楚一件事：你给出的
w_j的表达式是不是错了。我手算和用软件推导（当然编程 ...

刚不放心所以又带具体的数算了一遍，如果各个式子的定义真的如你所写，那么分母各项的指数毫无疑问就是2j而不是j+1。
赞一下 回复此楼
小木虫Mathematica版块已毁（当然原本也不咋的），建议大家前往百度贴吧或Stackexchange。
13楼2016-01-08 14:44:20
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

loujing

铁杆木虫 (正式写手)
引用回帖:
13楼: Originally posted by xzczd at 2016-01-08 14:44:20
刚不放心所以又带具体的数算了一遍，如果各个式子的定义真的如你所写，那么分母各项的指数毫无疑问就是2j而不是j+1。...

您好，感谢您的回复。w_j最终的闭解，是原书上给出的。所以我想利用软件验证一下。

我在网上搜了一下，李航老师在他新浪博客上给出了勘误表，http://blog.sina.com.cn/s/blog_7ad48fee01017dpi.html，如截图。
勘误表上，原错误w_j的结果分子上是x_i，而我的书上是x_i^j。看来，我这本书是已经修订过的了（购于2015年6月，2015年3月印刷的版本）。
关于用Maple或Mathematica实现符号多项式拟合

您能否将您的代码给出，或将.nb以附件形式上传，让我研究一下，非常感谢。
赞一下 回复此楼
14楼2016-01-08 17:03:58
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

loujing

铁杆木虫 (正式写手)
引用回帖:
13楼: Originally posted by xzczd at 2016-01-08 14:44:20
刚不放心所以又带具体的数算了一遍，如果各个式子的定义真的如你所写，那么分母各项的指数毫无疑问就是2j而不是j+1。...

如4楼所贴的书截图一样，一般遇到这种，利用输入x构造增广矩阵X；然后直接得到w的闭解 inv(X'X)X'y。所以还没去细研究过每个w_j到底是什么。

既然书上给出了w_j的闭解形式，又老听说Mathematica可以用来辅助公式推导，我就想用这个练练手的。

还望您能给出源码，谢谢了。
赞一下 回复此楼
15楼2016-01-08 17:17:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xzczd

木虫 (小有名气)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
loujing: 金币+30, ★★★很有帮助 2016-01-08 18:15:58
引用回帖:
14楼: Originally posted by loujing at 2016-01-08 17:03:58
您好，感谢您的回复。w_j最终的闭解，是原书上给出的。所以我想利用软件验证一下。

我在网上搜了一下，李航老师在他新浪博客上给出了勘误表，http://blog.sina.com.cn/s/blog_7ad48fee01017dpi.html，如截图。
...

我觉得你对这段勘误的理解是不是有误……这种勘误格式应该是相当于把wj的表达式整个都从原书里去掉了。总之这里给出我推导出12楼结果的代码及解说。

这个问题总的来说还是很有意思的。楼主的计算之所以会失败大致有两个原因。其一是原推导中不动声色地对抽象求和使用了一些相对特殊的处理，Sum[...w[j], {j, 1, M}]里的w[j]的j其实是个局部变量，仅用于由1遍历至M，而“对w[j]求导”这一运算中的j是一个“可以取值为1到M的任意整数的变量”，也就是说Sum内外的w[j]的所指其实是不同的；其二，Sum在未计算为封闭形式时并不会将内部的常数项自动向外提取（并且似乎也没有哪个化简函数会进行这种操作），一个简单的例子是
CODE:
Sum[1/2 a[b], {b, c, e}]

其中1/2并不会自动计算到外面。这是可以接受的，毕竟这种计算并没有使表达式变得更简单，而对于更一般的情形而言，这种“不变换”甚至可以说是必要的，比如：
CODE:
Sum[ff a[b], {b, c, e}]

这里的ff在表面上虽然是一个常数项，但是考虑到它在后续编程中可能被赋为任何值（包括与b有关的值）所以冒然将它提出显然会引发很多麻烦——于是我们就不难理解Solve为什么会求解失败了：对楼主的问题而言，若要求解方程的话，将符号常量提取到Sum外面无疑是必须的。

好的，说了这么多，我们到底应该怎么在Mathematica中完成这个推导呢？前已述及，直接基于内置函数的推导之所以会失败，是因为这一推导中用到了一些相对特殊的规则，那么，我们自己把这些规则编写出来就可以了，Mathematica核心语言中的模式匹配就是干这个的：
CODE:
Clear[m, n]

f[i_, m_] := Sum[w[j]*x[i]^j, {j, 0, m}]
l = 1/2 Sum[(f[i, m] - y[i])^2, {i, n}];

(*以下规则对应前面提及的第一点*)
rule[w_] := Sum[expr_, {i_, downup__}] /; ! FreeQ[w, i] :> expr
(*以下三条规则对应第二点*)
ruleExpand = Sum[(a_ + b_) c_, i_] :> Sum[(a + b) c // Expand, i];
ruleDistribute = Sum[expr1_ + expr2_, i_] :> Sum[expr1, i] + Sum[expr2, i];
ruleExtract = Sum[a_ b_., {i_, downup__}] /; FreeQ[a, i] :> a Sum[b, {i, downup}];
(*把相应的变换规则代入：*)
D[l /. rule[w@j], w[j]] //. {ruleExpand, ruleDistribute, ruleExtract}
Solve[% == 0, w[j]]

上面这个算是比较中规中矩的编法，其实，还可以更疯狂一点的。（注意，以下内容仅为展示一种可能性，对解决楼主的问题并非必要。）因为在这里D和Sum的固有规则对我们的计算造成了麻烦，那我们何不自己写一套D和Sum呢：
CODE:
rule2[w_] := sum[expr_, {i_, downup__}] /; ! FreeQ[w, i] :> expr

(*d[a_ +b_,x_]:=d[a,x]+d[b,x]*)
d[a_ b_, x_] /; FreeQ[a, x] := a d[b, x]
d[a_^b_, x_] /; FreeQ[b, x] := b a^(b - 1) d[a, x]
d[sum[expr_, i_], x_] := If[FreeQ[i, x], sum[d[expr, x], i], d[expr, x]]
d[expr_, i_] := D[expr, i]

sum[(a_ + b_) c_, i_] := sum[(a + b) c // Expand, i]
sum[expr1_ + expr2_, i_] := sum[expr1, i] + sum[expr2, i]
sum[a_ b_., {i_, downup__}] /; FreeQ[a, i] := a sum[b, {i, downup}]

f2[i_, m_] := sum[w[j]*x[i]^j, {j, 0, m}]
l2 = 1/2 sum[(f2[i, m] - y[i])^2, {i, n}];

d[l2 /. rule2[w@j], w[j]]

Solve[% == 0, w[j]] /. sum -> Sum
赞一下 回复此楼
小木虫Mathematica版块已毁（当然原本也不咋的），建议大家前往百度贴吧或Stackexchange。
16楼2016-01-08 17:24:49
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

loujing

铁杆木虫 (正式写手)
引用回帖:
12楼: Originally posted by xzczd at 2016-01-08 14:38:05
看看楼上的某些发言就觉得Mathematica推广工作任重而道远……楼主你这个问题是可以用Mathematica推导的。不过在给出解法前得先和你讨论清楚一件事：你给出的
w_j的表达式是不是错了。我手算和用软件推导（当然编程 ...

原书截图
关于用Maple或Mathematica实现符号多项式拟合-1
关于用Maple或Mathematica实现符号多项式拟合-2
赞一下 回复此楼
17楼2016-01-08 17:28:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xzczd

木虫 (小有名气)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
loujing: 金币+30, ★★★很有帮助 2016-01-08 23:56:20
引用回帖:
17楼: Originally posted by loujing at 2016-01-08 17:28:16
原书截图

...

……刚发现我为了凑出和你的顶楼相似的结果犯了个错误。盯着四楼的图看了半天（真伤眼），现在我基本确定，正确的推导是：
CODE:
f[i_, m_] := Sum[w[j]*x[i]^j, {j, 0, m}]
l = 1/2 Sum[(f[i, m] - y[i])^2, {i, n}];

rulej1[w_] := Sum[expr_, {i_, downup__}] /; ! FreeQ[w, i] :> sum[{i, downup}][expr]
rulej2 = sum[_]'[_] :> 1;
rulej3 = sum[i_][expr_] :> Sum[expr, i];

ruleExpand = Sum[(a_ + b_) c_, i_] :> Sum[(a + b) c // Expand, i];
ruleDistribute = Sum[expr1_ + expr2_, i_] :> Sum[expr1, i] + Sum[expr2, i];
ruleExtract = Sum[a_ b_., {i_, downup__}] /; FreeQ[a, i] :> a Sum[b, {i, downup}];

(D[l /. rulej1[w@j], w[j]] /. {rulej2, rulej3} //. {ruleExpand, ruleDistribute,
     ruleExtract} // Expand) == 0



这和4楼结果是一致的，显然，光靠一个方程是定不了w[j]的，必须联立方程组。w[j]的简单的封闭解是不存在的。
赞一下 回复此楼
小木虫Mathematica版块已毁（当然原本也不咋的），建议大家前往百度贴吧或Stackexchange。
18楼2016-01-08 17:59:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xzczd

木虫 (小有名气)

【答案】应助回帖

引用回帖:
17楼: Originally posted by loujing at 2016-01-08 17:28:16
原书截图

...

16楼已经说了，这种勘误格式等于是把这一整段从书中删掉了……
赞一下 回复此楼
小木虫Mathematica版块已毁（当然原本也不咋的），建议大家前往百度贴吧或Stackexchange。
19楼2016-01-08 18:01:14
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

loujing

铁杆木虫 (正式写手)
引用回帖:
18楼: Originally posted by xzczd at 2016-01-08 17:59:09
……刚发现我为了凑出和你的顶楼相似的结果犯了个错误。盯着四楼的图看了半天（真伤眼），现在我基本确定，正确的推导是：

f := Sum
l = 1/2 Sum;

rulej1 := Sum /; ! FreeQ :> sum
rulej2 = sum' :> ...

追问一下，在得到这里的联立方程组后，符号推导的工作是不是就可以算结束了？

如果接下来，我给出具体的一组(x_i,y_i)，如何再利用这里的推导，计算出具体的w_j？
赞一下 回复此楼
20楼2016-01-08 18:15:38
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 loujing 的主题更新
232/3返回列表上一页123下一页
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭