公告:网站程序已升级到1.8.3,修复了提问时可能报错的问题,请清除浏览器缓存
2019-11-10

欢迎来到 Mathematica 问答社区

提问时请贴上文本代码

语法高亮:在编辑器中点击

被禁止的话题:广告破解

请阅读:《提问的智慧》

备用域名:mma.ooo

支持LaTex数学公式
行内公式标识符:\$ 或“$\backslash ($”+“$\backslash )$”,
行间公式标识符:\$\$ 或 “$\backslash [$”+“$\backslash ]$”

社区建议QQ群:365716997

分类

+2 投票
752 浏览

背景:

实验室得到的一组原子比,已经归一化(和为1),但由于实验会有误差,所以比值不会完全准确,我希望把比例处理成一组整数比,使得数值尽量小的前提下,与原始数值的误差尽量小(我并不知道误差的精确定义),比如如果为了精确,我们得到了200:200:299,那么我更愿意误差再大一点而得到2:2:3。

下面模拟一组2:2:3的实数比例
 

SeedRandom[1]
list={2,2,3}/.t_?NumberQ:>t+RandomReal[{-.1,.1}]//#/Total[#]&
{0.292955,0.27291,0.434135}

然后有什么好办法可以让它还原成2:2:3不?

用户: Lozmlve*永 (1.1k 分)
修改于 用户:Lozmlve*永

3 个回答

+5 投票
 
已采纳

先发代码:

xyz = {0.292955, 0.27291, 0.434135};
data = {#, Total[(# - Round@#)^2]/3 &[# xyz]} & /@ Range[1, 1000, 0.1];
len = Length@data;
sigma = 1;
ans = If[Last@# < sigma,
         Sow@{#, Round[xyz First@#], 
           k = xyz.Round[xyz First@#]/(#.# &@
               xyz), #.# &@(# - Round@# & /@ (k xyz))/3}; 
         sigma = Last@#] & /@ data // Reap // Last // Last // 
   DeleteDuplicates[#, Last@#1 == Last@#2 &] &;
TableForm[Flatten /@ ans, 
 TableHeadings -> {None, {"k-approx", "\[Sigma]-approx", "a", "b", 
    "c", "k", "\[Sigma]"}},
 TableAlignments -> Center]

运行结果为

====================以下是我的思路的分割线(太长不看)====================

首先,原始的数据应为{a,b,c},然而由于实验误差,测得的值为{a+e,b+e,c+e}(e为服从NormalDistribution[0,sigma]的正态分布变量),再放缩至{a+e,b+e,c+e}/k(注意这里k可以是任意实数而不限于正整数),最后表现为{x,y,z}的形式。

那么,如何从{x,y,z}还原至{a,b,c}呢?先需要放大至k{x,y,z},用Round[k{x,y,z}]来拟合{a,b,c}的值,其偏离值#-Round@#&/@(k{x,y,z}可视为e的三次取样,不防设为d1,d2,d3。由概率统计中的极大似然思想,可知三次取样为d1,d2,d3服从的正态分布的参数sigma^2的估计值为(d1^2+d2^2+d3^2)/3。显然sigma越小,拟合效果越好。

这就有点思路了,可以让k逐渐增大,对每个k可以对应一个sigma的估计值。而我们的目的是要找最接近{x,y,z}比例的{a,b,c},也即任给{a1,b1,c1}<{a,b,c},对应的sigma1>sigma。换句话说,在sigma=f[k]这个函数中找到某个特定的k0,对应的sigma0是k取[1,k0]间的值的最小值,即k<k0时必有f[k]<f[k0]。此时对应的{a,b,c}即是比之前拟合效果更好的拟合。

然而,sigma=f[k]=#.#&@(#-Round@#&/@(k{x,y,z}))/3,这样涉及数论函数Round的一个式子无法用导函数为0的方法求出k值。果断耍个流氓,先对k值离散化,用k=Range[1,1000,0.1]来求出一系列的sigma值,然后找出所有sigma0使得sigma0比前面的sigma值都小,即为一个局部极小值的近似值,此时对应的{a,b,c}即为所求。以此k值(即图中的k-approx)为初值求出真正使得sigma0值最小的k和sigma即可。根据sigma=(kx-a)^2+(ky-b)^2+(kz-c)^2,k作微小改变时{a,b,c}不会改变(否则k{x,y,z}必有一个接近n+0.5,此时sigma值必然很大,矛盾),可将a,b,c视作常数并对k求导,于是k=(ax+by+cz)/(x^2+y^2+z^2)时即为k-approx附近的极小值点,并以此计算出sigma即可。

上述代码就是基于这个思路来写的。先找出所有sigma近似最小值及对应的k近似值,以此算出{a,b,c},去重(可能会有好几个k对应同一组{a,b,c}),最后算出真正的k和sigma。

做出表格后,可根据需要的sigma值为基准进行选择,比如只要sigma^2>0.004即可,那{2,2,3}即为满足需要的整数比;若要求sigma^2取0.001,则需要的整数比应为{29,27,43},以此类推。

 

用户: 瓦屋青衣 (321 分)
采纳于 用户:Lozmlve*永
你这统计学,用得杆杆的。。
目测是有限范围内暴力枚举并选优?
对,本质上仍是暴力法,实在不会对含Round的函数进行最优化分析……从结果来猜测似乎与连分数的最佳逼近也有所关联,如果能用上其结论就可以直接求值了,然而数论学得实在太渣,完全不知从何下手。苹果有什么好的想法没?
+2 投票
findApproximateRatio[list_, err_] := 
 Rationalize[Divide[#, First@#] &@list, err] // #/GCD @@ # &

误差的话调整err的值就可以了

V2:多一点选择,分别以列表中所有的数进行参考

findApproximateRatio[list_, 
  err_] := (Rationalize[#, err] // #/GCD @@ # &) & /@ 
  Transpose@Outer[Divide, list, list]

 

用户: happyfish (1.8k 分)
修改于 用户:happyfish
比如实际数值有一组是{1.08, 0.95, 97.965},你觉得这个比例是多少呢?我觉得这里用Rationalize来处理可能会不太恰当,如果某个值与原数据的误差为0.04达到最合题意时,我们的err取为0.1的话,这个方法是不会取到这个0.04的,我感觉这应该是个优化方面的问题,拙见,望指点
这个比例应该是多少不应该看起来像是几就是几,至少应该有一套严谨的理论,然而相关统计学知识我并不了解。现在如果不能拿出一个标准或者算法,我也不知道要怎么改
@Lozmlve*永 您这个例子放到上面的方法需要更大的k,如果相差太大,k也需要更大,方法还是有问题的
0 投票
xyz = {0.292955, 0.27291, 0.434135};
Manipulate[xy = Rationalize[xyz[[1]]/xyz[[2]], d];
 yz = Rationalize[xyz[[2]]/xyz[[3]], d];
 lcm = LCM[Denominator[xy], Numerator[yz]]; {xy*lcm, lcm, lcm/yz}, {d,
   0.01, 0.1}]
用户: wengdeping (41 分)
修改于 用户:wengdeping
...