英语原文共 14 页，剩余内容已隐藏，支付完成后下载完整资料
利用Python进行科学计算：灵活高效的求取输入为多元随机值的统计特征方程的值
绪论：本文通过用于评估随机变量的统计特性的多维集成示例来检查基于Python的高级应用程序对数值密集型应用程序的可行性。我们讨论使用高级脚本代码和低级编译代码来实现数学公式增量表达式的方法。由于Python语言的动态类型，算法的组件可以很容易地以通用的方式编码为算法模板。使用Enthought开发套件，它们可以有效地被组装到一个灵活的计算框架中，该框架可以配置为执行用于任意组合集成方案和实例化代码版本的代码。本文使用包含不连续性的随机双参数函数的平均值的简单运行示例来描述开发周期。此示例还用于比较可用算法和执行功能的性能。包含更多示例和性能研究结果的已实施的软件包已通过免费存储库[1]和CPCP库提供。
正文：
1、介绍
用于科学计算的高级语言为应用程序开发人员提供了一种方便高效的数学模型制定和实施工具。用于科学应用程序快速原型开发的广泛使用的高级软件示例包括Maple，Matlab，Octave，R和S 。这些工具提供了丰富的文档，可视化工具，符号运算符和大量数值方法。在[2]中，已经讨论了Matlab工具包对数字应用原型和教学课程的适用性。与像FORTRAN和C或C 这样的低级编程语言相比，这些高级工具使开发更加高效，特别是在应用程序开发的早期阶段。另一方面，这样的高级工具很少会与以编译语言编写的应用程序的性能相竞争。灵活性和效率之间的折衷是应用程序员在科学计算领域的日常问题。理想的开发环境应该在开发的早期阶段提供高生产力，同时在应用程序达到成熟状态时提供快速代码执行的简单方法。
在过去的十年中，开源软件领域出现了一个特别吸引人的开发环境。 脚本语言Python已经建立了开发和集成两个丰富的数字库numpy和scipy的平台。这些库体现了过去几十年在编译语言FORTRAN和C （http://www.netlib.org）中开发的算法和方法。教科书已经充分地介绍了应用于科学应用的Python脚本编程的灵活性[3]。正如[4,5]中进一步记载的那样，脚本语言的高度灵活性并不一定伴随着编译代码的低效率。事实上，当使用紧凑阵列表示[6]以矢量化形式实现数学表达式时，灵活性和效率之间的冲突被降低到可接受的水平[7]。
除了与编译的低级语言的比较之外，Python开发环境的生产力和效率已经使用几个基准示例与Matlab和Octave进行了比较[8,9]。这些研究得出的结论是，基于Python环境的功能与商业工具的功能相当。然而，作者认为Python环境在一个关键领域具有决定性优势：它的语言背后有一个开放和可扩展的对象模型，以便它可以进一步发展。使用元类扩展语言中的类定义的可能性允许将设计模式结合到语言中，从而进一步提高开发人员的生产力。尤其是，Enthought Tool Suite（http://www.enthought.com）引入了一个使用预定义特征的类的精确定义[10]。类特征与面向对象语言中的常用属性相对应，但是它还包含有关初始化、验证、变更通知和可视化的信息。此扩展属性规范只需最少的编程工作就能将模型-视图-控制器设计模式合并到语言并提供自动对象可视化和状态变化的动态控制。这种语义丰富的语言支持允许开发人员直接从类视图规范自动生成用户界面，并对对象的持久性、状态依赖关系的声明以及对数据可视化提供有效支持。
本文主要阐述了基于Python的环境在科学应用开发中的可行性。这里提出的概念、观察和结论与对三个学科感兴趣的读者有关：统计矩的估计，科学计算应用程序的编程以及纤维束和/或脆性基质复合材料的模拟。在整篇文章中，展示了应用实施概念在估计一组平行纤维的应力 - 应变响应的统计矩的作用，该任务对应于使用随机参数估计“基本”函数的统计矩。图1中显示了这样一个函数的例子，该函数描述了一个拉伸载荷下纤维的应力 - 应变响应，其中包含两个独立的、相同分布的刚度和强度的随机参数。单光纤响应的随机实现显示在图1中。目的是有效地估计光纤束内光纤的平均响应，如图1右侧黑色实线所示。
这个具体例子的实现已经被推广，因此它适用于具有独立随机参数的任意函数的统计表征。与pychem [11]一样，实施的框架属于多变量分析软件包的范畴。我们将证明，当使用上述实施环境时，数学公式的抽象特性可以反映在代码中，而不会有任何额外的性能损失。本文的特别目的是
- 以简单的Python脚本开始显示应用程序的快速原型，
- 演示如何将脚本概括并放大到可配置的算法对象，该算法对象可以包含随机采样和选项的多种变体以加快代码的速度，
- 为实施的所有可用功能提供结果验证框架和执行效率研究框架。
本文的结构安排如下。首先，在第二节中介绍了所实现算法的符号规范。然后在第3节构造纯Python中算法的简单实现，另一方面使用numpy和scipy包。接下来，算法的实现被推广到包含其他类型的抽样和语言（第4节）。然后使用Enthought作为中间件，将该算法的引入变体合并到第5节中用于多维集成的通用交互式框架中。可配置程序包对第8节中描述的几种版本采样类型的效率进行彻底比较。比较执行效率的三种代码生成和编译当前响应函数及其在第9部分的随机化代码段。整个文件中提供的代码段构成一个可执行的脚本，包括spirrid包中可供下载在线CPCP库和github.com存储库[1]。
2. 用独立随机变量估计函数的统计矩
目标是估计作为控制变量ε和构成随机域Theta;的随机变量theta;1...theta;m的函数给出的随机问题Q的统计矩：
其中q（ε;theta;）进一步被称为响应函数。一个随机问题的第k个原始统计矩是给定的
由于这里只考虑独立的随机变量，所以随机向量theta;的联合概率密度函数g（theta;）（PDF）可以用单变量边际密度表示
积分将作为分布在离散化随机域上的离散值的总和进行数值计算
其中Delta;Gi（theta;i）表示取决于下面指定的特定采样类型的加权因子。随机域内积分点Theta;m的分布可以表示为
其中是第i个变量的离散化点的数量，j是对一个变量的离散化点进行计数的数值。有很多方法通过采样点来覆盖随机域。为了演示语言的表达能力并讨论可能的实现技术的计算效率，我们将以四种不同的方式实现积分（3）：
- 随机变量的等距网格（TGrid）。第i个随机变量被规则间隔值覆盖。与积分单元相关的概率给定为
表示离散点之间的间隔
- 随机变量的非等距网格。采样点是通过采样概率（PGrid）的等距网格生成的。将采样概率转换为等式（1）中的随机域Theta;。（5）针对使用逆累积分布（也称为百分点函数）列举为i = 1··m的每个随机变量执行：
其中的值可以得到
通过这种取样呈现的整合单元对于每个变量具有相同的概率（权重）：
这种类型的采样在此处作为网格类型采样与具有恒定积分单元概率的蒙特卡罗类型采样的概念
- 粗蒙特卡罗采样（MCS）。（8）中给出的的结构化网格，而不是采用随机选择采样概率，根据式（7）选择相应的采样点。 当所有随机变量的采样点数时，Theta;i数组可以变平坦，并且可以在单个循环中完成积分，如
- 拉丁超立方体抽样（LHS）。蒙特卡洛类型抽样的增强版本，它使用输入随机变量的分布函数的分层来确保抽样概率的均匀覆盖[12-14]。一个LHS样本倾向于通过 单位超立方体的抽样概率大于蒙特卡罗样本。因此，与MCS相比，方程（11）中的统计量估计量的方差通常会减小。实际上，斯坦因[15]显示LHS 与简单随机抽样（粗蒙特卡罗）相比，方差减小量增加。方差减少量随函数q所依赖的随机量中的可加性程度和单调性增加。
类似于MCS的情况，基础数据结构可以被平坦化并且可以按照方程（11）中的规定在单个循环中执行积分。方程（7）中的采样点通过以下计算的采样概率来获得：
表示包含序列1，...的随机置换的矩阵。通过使用这种方案，每个层的中位数被选为采样点。
该方法的进一步改进可以通过选择每个阶层的概率手段来实现[16]。 我们注意到抽样概率的随机排序不能保证模拟数据的精确单位相关矩阵，特别是对于小样本量，参见文献[17]。在[16]中已经提出了一种可能的补救措施，即采用一种算法来减少不希望的相关性例子。
为了说明实现的说明，让我们介绍一个简单的双参数响应函数，如图1（左）所示
该函数定义了线性弹性脆性纤维的给定正控制应变ε的应力，刚度参数lambda;和断裂应变xi;。符号H（x）表示当x lt;0时产生零的Heaviside函数，xge;0时产生一个。我们人为选择一个包含不连续性的函数来研究积分算法重现平均响应的平滑性的能力由非光滑函数控制的组成响应。如果许多纤维平行地起作用，就像复合材料中的裂纹桥表现出缺陷。这意味着材料参数lambda;和xi;必须被认为是随机的。使用方程（4）可以得到一束光纤中的平均应力-应变行为
这个表达式的结果绘制在图1中（右边，黑色曲线）。它表明长丝的平均响应是非线性的。所描述的积分展示了描述纱线行为的基于应变纤维束模型的结构和复合材料[18]。作者已经使用该程序对复丝纱线拉伸试验进行建模[19]。
3.脚本实现等距集成方案
为了说明实现数值模型的可能方法，让我们构造一个简短的脚本来传递式（14）的结果。在第一步中，我们将参数lambda;和xi;的概率分布定义为正常，并生成 具有20个值的随机区域Theta;lambda;和Theta;xi;的等距离离散以及边缘密度glambda;和gxi;的值。
from scipy.stats.distributions import norm 1
import numpy as np 2
# set the mean and standard deviation of la and xi 3
m_la, std_la = 10.0, 1.0 4
m_xi, std_xi = 1.0, 0.1 5
# construct objects representing normal distributions 6
pdistrib_la = norm(loc=m_la, scale=std_la) 7
pdistrib_xi = norm(loc=m_xi, scale=std_xi) 8
# get operators for probability density functions 9
g_la = pdistrib_la.pdf 10
g_xi = pdistrib_xi.pdf 11
# number of integration points set equal for both variables 12
n_i = 10 13
# generate midpoints of n_i intervals in the range (-1,1) 14
theta_arr = np.linspace(-(1.0 - 1.0 / n_i), 15
1.0 - 1.0 / n_i, n_i) 16
# scale up theta_arr to cover the random domains 17
theta_la = m_la 4 * std_la * theta_arr 18
theta_xi = m_xi 4 * std_xi * theta_arr 19
# get the size of the integration cells 20
d_la = (8 * std_la) / n_i 21
d_xi = (8 * std_xi) / n_i 22
注意，随机域已经被等距集合点集合覆盖，该集合点分两步产生：首先，范围lang;-1,1rang;在第15行被集成区间的中点阵列覆盖。第二，范围lang;-1,1rang;是以平均值为中心并在第18-19行放大四个标准偏差，以充分覆盖正态分布的非零域。响应函数本身与方程（13）一致地实现，
def Heaviside(x): 23
'''Heaviside function.''' 24
return x gt;= 0.0 25
26
def q_eq13(eps, la, xi): 27
'''Response function of a single fiber.''' 28
return la * eps * Heaviside(xi - eps) 29
请注意，应该避免使用控制语句if和else。相反，使用返回布尔值的Heaviside函数来允许数组作为响应函数的参数。这对于允许对函数进行矢量化评估是必要的，稍后会解释。
公式（14）中总和表达式将由三个循环组成：一个循环用于收集沿控制变量ε的响应值和两个沿着随机变量域Theta;lambda;和Theta;xi;的循环：
def mu_q_eq14_loops(eps_arr): 30
'''Loop-based calculation of mean values.''' 31
mu_q_arr = np.zeros_like(eps_arr) 32
for i, eps in enumerate(eps_arr): 33
mu_q = 0.0 34
for la in theta_la: 35
for xi in theta_xi: 36
dG = g_la(la) * g_xi(xi) * d_la * d_xi 37
mu_
全文共19572字，剩余内容已隐藏，支付完成后下载完整资料
资料编号：[12711]，资料为PDF文档或Word文档，PDF文档可免费转换为Word
以上是毕业论文外文翻译，课题毕业论文、任务书、文献综述、开题报告、程序设计、图纸设计等资料可联系客服协助查找。
