65.9K
CodeProject 正在变化。 阅读更多。
Home

高级数学的交互式环境

dawright
starIconstarIconstarIconstarIcon
emptyStarIcon
starIcon

4.68/5 (11投票s)

2009 年 7 月 14 日

Ms-PL

8分钟阅读

viewsIcon

42306

结合 .NET 数学库和 .NET 动态编程语言，创建类似 Matlab/Mathematica 的数学和数据分析交互环境。

引言

像 Mathematica 和 Matlab 这样的系统允许用户与数学对象进行交互式操作：计算表达式、创建和使用变量，并在一个立即执行输入命令并打印结果的界面中对它们进行操作。

用于科学编程的 Meta.Numerics 库提供了一些昂贵专有系统的高级数学和统计功能，但以 .NET 库 API 的形式存在，您可以从自己的程序中调用。没有“前端”允许您以 Mathematica 和 Matlab 的方式交互式地访问 Meta.Numerics 功能。构建一个“前端”用户界面是一项昂贵且耗时的任务。有没有办法让我们现有的一些 .NET 交互式解析器能够交互式地暴露 Meta.Numerics 功能？

IronPython 提供了一个现成的答案。IronPython 是使用 .NET 动态语言运行时 (DLR) 的 Python 编程语言的实现。它的作者编写了一个所谓的读-求值-打印-循环（REPL）界面，允许用户交互式地编写 Python 程序。由于 IronPython 是一种 .NET 语言，它可以调用 .NET API。因此，即使我们对用 Python 编程不感兴趣，我们也可以使用 IronPython 的 REPL 来与 Meta.Numerics 对象进行交互。本教程将展示如何做到这一点。

入门

安装 IronPython（可在 http://ironpython.codeplex.com/Release/ProjectReleases.aspx 获取）和 Meta.Numerics（可在 http://metanumerics.codeplex.com/Release/ProjectReleases.aspx 获取）。

不幸的是，IronPython 安装程序不会将 IronPython 添加到您的程序菜单中，但在安装后，您将在“程序文件”目录下的 IronPython 文件夹中找到可执行文件（例如，C:\Program Files\IronPython 2.0.1\ipy.exe）。启动 IronPython 可执行文件；您将看到一个如下所示的交互式界面。

要能够在此交互式界面中访问 Meta.Numerics 对象，您需要添加对 Meta.Numerics 程序集的引用并导入您想要的这些对象。如果您已将 Meta.Numerics 程序集安装到计算机的 GAC（如果使用 MSI 安装程序，则会自动发生），您可以直接输入

>>> import clr
>>> clr.AddReference("Meta.Numerics")

（为求清晰：您无需输入“>>>”——此处仅为区分您在 IronPython 提示符下输入的内容与 IronPython 响应的内容。）

另一方面，如果 Meta.Numerics 程序集文件仅位于您的硬盘驱动器上但不在您的 GAC 中（这在您解压 ZIP 文件时会发生），您可以改为输入

>>> import clr
>>> clr.AddReferenceToFileAndPath(r"C:\Program Files\Meta.Numerics\Meta.Numerics.dll")

当然，如果 Meta.Numerics 程序集在您计算机上的路径不同，您应将其替换为实际路径。

现在 IronPython 已成功链接到 Meta.Numerics 程序集，我们只需告诉它将哪些 API 导入到默认命名空间。这基本上类似于发出 C# 的 using 或 VB 的 Imports 语句。我们现在就将这样做，以导入我们计划使用的所有 API。

>>> from Meta.Numerics import *
>>> from Meta.Numerics.Functions.AdvancedMath import *
>>> from Meta.Numerics.Functions.AdvancedIntegerMath import *
>>> from Meta.Numerics.Functions.OrthogonalPolynomials import *
>>> from Meta.Numerics.Functions.FunctionMath import *
>>> from Meta.Numerics.Statistics import *
>>> from Meta.Numerics.Matrices import *

语句 from Meta.Numerics.Statistics import * 导入 Meta.Numerics.Statistics 命名空间中的所有类。语句 from Meta.Numerics.Functions.AdvancedMath import * 导入 FunctionMath 类中的所有静态 *方法*，因此可以通过直接输入它们的名字（无需类前缀）来访问它们。例如，您将在下面看到，我们只需输入 Gamma(0.5)，而无需输入 AdvancedMath.Gamma(0.5)。

为了避免反复输入所有这些内容，只需创建一个包含您到目前为止输入的所有 Python 命令的文本文件。（我无法为您完成此操作，因为路径字符串将取决于您的安装。）将文件命名为“MetaNumerics.py”并将其放在（比如）Meta.Numerics 目录中。下次您想在 IronPython 中交互式地使用 Meta.Numerics 时，只需启动 IronPython，然后输入...

>>> execfile(r"C:\Program Files\Meta.Numerics\MetaNumerics.py")

...然后您就可以开始工作了。

基础与高级数学

由于简单的算术表达式就是 Python 表达式，您只需直接输入即可。

>>> 2 + 3 * 4
14

请注意，Python 尊重运算顺序。您还可以赋值变量并计算包含它们的表达式。

>>> a = 6 / 2
>>> a + 2
5

当然，这根本没有用到 Meta.Numerics。为了使用它，让我们调用一个高级函数。

>>> Gamma(0.5)
1.77245385091
>>> Gamma(5) 
24.0

请注意，Γ(1/2) 是 π1/2，而 Γ(5) 是 4!。您还可以计算不完全 Gamma 函数、Beta 函数和 Psi 函数。那么黎曼 Zeta 函数呢？

>>> RiemannZeta(3.0)
1.20205690316

没问题。您还可以计算狄利克雷 Eta 函数。有关特殊函数可用列表的完整信息，请参阅项目网站上的 Meta.Numerics 文档。

现在让我们做得更花哨一些。让我们找到零阶贝塞尔函数 J0 的第一个根。Meta.Numerics 的 FindZero 函数期望接收一个函数作为参数，该函数接受一个参数：一个双精度数。但 Meta.Numerics 的 BesselJ 函数接受两个参数：一个整数（阶数）和一个双精度数（参数）。为了克服这个小障碍，我们只需定义一个返回 BesselJ(0,x) 的 Python 函数。

>>> def J0(x): return(BesselJ(0,x))
>>> x0 = FindZero(J0, 3.0)
>>> x0
2.40482555770

（FindZero 函数的第二个参数只是一个根搜索的起点。）现在，让我们计算 J0 从原点到第一个根的积分。

>>> Integrate(J0, Interval.FromEndpoints(0.0,x0))
1.47030004338

很酷。

统计

现在，让我们转向一些数据分析。Meta.Numerics 为不同类型实验数据定义了不同的类。每个类都提供适当的描述性统计、统计检验和拟合程序。我们将举几个例子。

示例

假设我们有一些组件的寿命测量值。让我们将它们放入合适的数据容器中。

>>> s = Sample()
>>> s.Add(1.0)
>>> s.Add(1.3)
>>> s.Add(1.5)
>>> s.Add(1.7)
>>> s.Add(2.9)

哎呀，最后一个数据点应该是 1.9。幸好我们在交互式框架中工作。

>>> s.Remove(2.9)
True
>>> s.Add(1.9)

那么，平均寿命是多少？

>>> s.PopulationMean
1.48 ± 0.156204993518133

请注意，Meta.Numerics 提供了对其总体参数估计的不确定性。它还将提供任意高阶矩的估计以及不确定性。

假设我们想保证平均寿命高于 1.25。我们的数据能否让我们有 95% 的信心做到这一点？为此，进行 t 检验。

>>> t = s.StudentTTest(1.25)
>>> t.Statistic
1.47242411923
>>> t.LeftProbability
0.892554887728

还不够。我们的数据让我们有 89% 的信心满足要求，但不是 95% 的信心。

DataSet

接下来，我们处理带有误差线的数据。假设我们已经进行了以下测量，例如，水库水位随时间的变化。

时间 (h) 水位 (m)
2.0 2.5 ± 0.5
3.0 2.0 ± 1.0
5.0 -1.0 ± 0.5
8.0 -4.0 ± 0.5
10.0 -1.5 ± 0.5
12.0 0.0 ± 1.0

创建一个 DataSet 对象来处理这些数据。

>>> D = DataSet()
>>> D.Add(2.0,2.5,0.5)
>>> D.Add(3.0,2.0,1.0)
>>> D.Add(5.0,-1.0,0.5)
>>> D.Add(8.0,-4.0,0.5)
>>> D.Add(10.0,-1.5,0.5)
>>> D.Add(12.0,0.0,1.0)
>>> D.Count6

让我们尝试用稳定的变化率来模拟数据。

>>> f1 = D.FitToLine()
>>> f1.GoodnessOfFit.LeftProbability
0.999999999654

该模型基本上被排除了：在原假设（模型正确）下，获得如此糟糕拟合的概率仅为 0.0000000001。那么正弦振荡呢？

>>> def ff(a,x): return(a[0] * Math.Sin(2.0 * PI * x / a[1] + a[2]))
>>> f2 = D.FitToFunction(ff,(2,10,0))
>>> f2.GoodnessOfFit.LeftProbability
0.690194092616

这拟合得相当好。（请注意，(2,10,0) 只是振幅、周期和相位参数的初始猜测。只要您的初始猜测不太离谱，其精确值就不重要，Meta.Numerics 会找到最佳拟合点。）那么，最佳拟合振荡周期是多少？

>>> f2.Parameter(1)
14.5075811109964 ± 0.989537307149617

大约 14 小时。

BinaryContingencyTable

让我们再处理一种实验数据。假设我们随机分配患者接受或不接受实验性治疗，并获得以下结果。

康复 死亡 总计
治疗 12 6 18
未治疗 8 15 23
总计 20 21 41

这是一个 2x2 列联表，Meta.Numerics 为这种类型的实验提供了一个定制类。

>>> B = BinaryContingencyTable()
>>> B[0,0] = 12
>>> B[0,1] = 6
>>> B[1,0] = 8
>>> B[1,1] = 15

治疗似乎有帮助，但样本量很小：我们能有多确定这不是统计上的巧合？某些条目太小，无法满足 Pearson 检验的假设。幸运的是，Meta.Numerics 实现 Fisher 精确检验，即使对于小样本也是合适的。

>>> B.FisherExactTest().LeftProbability
0.061603593338

在原假设（治疗无效）下，获得如此显著结果的几率仅为 6%。因此，我们可以相当有信心，但还不够达到通常医学接受惯例要求的 95% 的信心。

那么，接受治疗比不接受治疗的机会大约提高了多少？

>>> B.OddsRatio
3.75 ± 2.49217525467211

大约是 4 倍，但存在相当大的不确定性。

请参阅 Meta.Numerics 文档，了解这些数据分析类的其他功能，并了解更多数据分析类。

矩阵运算

Meta.Numerics 为不同类型的矩阵定义了特定的类：例如，SquareMatrix、SymmetricMatrix、TridiagonalMatrix。每个类都具有适合该类型矩阵的方法，并针对利用矩阵结构进行了优化实现。

例如，让我们创建一个 3x3 的对称矩阵。

>>> M = SymmetricMatrix(3)
>>> M[0,0] = 6
>>> M[0,1] = 4
>>> M[0,2] = 1
>>> M[1,1] = 5
>>> M[1,2] = 2
>>> M[2,2] = 3
>>> M
{  6  4  1  }
{  4  5  2  }
{  1  2  3  }

请注意，由于 Meta.Numerics 知道我们的矩阵是对称的，我们只需告诉它一半的矩阵元素。它根据对称性填充了其余元素，并且在底层，它只使用了 SquareMatrix 存储相同条目所需的存储空间的一半。

Cholesky 分解是将矩阵 M 分解为 M = S ST，其中 S 是下三角矩阵。（S 通常被称为矩阵 M 的“平方根”）。它仅适用于正定对称矩阵，但如果您有这样的矩阵，Cholesky 分解比 LU 分解更快、更稳定。SymmetricMatrix 类提供了此分解。

>>> CD = M.CholeskyDecomposition()
>>> S = CD.SquareRootMatrix()
>>> S
{   2.44948974278               0              0  }
{   1.63299316186   1.52752523165              0  }
{  0.408248290464  0.872871560944  1.43924583426  }

让我们验证分解是否正确。

>>> S * S.Transpose()
{  6  4  1  }
{  4  5  2  }
{  1  2  3  }

请注意，表达式的计算需要矩阵乘法。您可以使用 Cholesky 分解来求解方程组...

>>> CD.Solve((1.0,2.0,3.0))
{  0  }
{  0  }
{  1  }

...或求矩阵的逆。

>>> MI = CD.Inverse()

Meta.Numerics 还可以找到特征值和特征向量。有关使用 Markov 矩阵对角化分析 Monopoly 这个棋盘游戏的一个有趣应用，请参阅 CodeProject 文章“Markov Monopoly”。

除了让我们能够动态交互我们的库之外，IronPython 界面还允许我们以一些有用的方式反思我们的库。想知道上一个结果的类型吗？

>>> type(MI)
<type 'SymmetricMatrix'>

想知道该类型提供了哪些属性和方法吗？

>>> dir(MI)
['CholeskyDecomposition', 'Clone', 'ColumnCount', 'Dimension', 'Eigensystem', 'E
igenvalues', 'Equals', 'GetHashCode', 'GetType', 'Inverse', 'Item', 'MemberwiseC
lone', 'ReferenceEquals', 'RowCount', 'ToString', 'Trace', 'Transpose', '__add__
', '__class__', '__delattr__', '__div__', '__doc__', '__eq__', '__getattribute__
', '__getitem__', '__hash__', '__init__', '__mul__', '__ne__', '__new__', '__rad
d__', '__reduce__', '__reduce_ex__', '__repr__', '__rmul__', '__rsub__', '__seta
ttr__', '__setitem__', '__str__', '__sub__']

结论

您已经看到了如何在 IronPython 界面中使用 Meta.Numerics 库进行交互。这项技术以极少的努力增加了大量价值。这个想法可以轻松应用于允许任何 .NET 库的交互式使用。

© . All rights reserved.