Sage统计快速入门

这个 Sage 快速入门教程是为MAA准备工作坊“Sage：对本科生使用开源数学软件”开发的（由NSF提供资金，到期日：0817071）。它是根据Creative Commons Attribution -ShareLiked 3.0许可证授权的 (CC BY-SA ）

尽管Sage最初是作为一个代数和几何学的项目开始的，但它在统计学和金融学方面有许多功能。尤其是由于 R project 作为Sage的一个组成部分，我们拥有非常强大的统计技术。

基本描述性统计

一些基本的统计函数是内置的。

sage: mean([1,2,3,5])
11/4

sage: std([1,2,2,4,5,6,8]) # The standard deviation
sqrt(19/3)

一旦我们超越了这些东西，通常有几种方法来完成事情——这些方法可能很复杂，但也很强大。

分布

让我们使用本机Python随机生成器从给定的连续分布类型生成一个随机样本。

  • 我们使用最简单的方法从（正态）平均值为2的对数正态分布中生成随机元素 \(\sigma=3\) .

  • 注意，真的没有办法绕过某种循环。

sage: my_data=[lognormvariate(2,3) for i in range(10)]
sage: my_data # random
[13.189347821530054, 151.28229284782799, 0.071974845847761343, 202.62181449742425, 1.9677158880100207, 71.959830176932542, 21.054742855786007, 3.9235315623286406, 4129.9880239483346, 16.41063858663054]

我们可以检查数据的对数平均值是否接近2。

sage: mean([log(item) for item in my_data]) # random
3.0769518857697618

下面是一个使用 Gnu scientific library 在引擎盖下面。

  • 让 dist 是分配给标准偏差为3的连续高斯/正态分布的变量。

  • 然后我们使用 .get_random_element() 方法10次，每次加2次，平均值等于2。

sage: dist=RealDistribution('gaussian',3)
sage: my_data=[dist.get_random_element()+2 for _ in range(10)]
sage: my_data # random
[3.18196848067, -2.70878671264, 0.445500746768, 0.579932075555, -1.32546445128, 0.985799587162, 4.96649083229, -1.78785287243, -3.05866866979, 5.90786474822]

现在，直接获取这些东西的直方图有点烦人。在这里，我们得到了这个分布的一个更大的抽样，并用10个箱子绘制了一个柱状图。

sage: my_data2 = [dist.get_random_element()+2 for _ in range(1000)]
sage: T = stats.TimeSeries(my_data)
sage: T.plot_histogram(normalize=False,bins=10)
Graphics object consisting of 10 graphics primitives

要访问离散分布，我们访问Sage的另一部分，它内置了统计信息： Scipy .

  • 我们必须 import 这个模块。

  • 我们使用 binom_dist 表示20次试验和5%预期失败率的二项分布。

  • 这个 .pmf(x) 方法给出了 \(x\) 失败，然后在条形图中标出 \(x\) 从0到20。（别忘了 range(21) 表示所有整数 零到二十 ）

sage: import scipy.stats
sage: binom_dist = scipy.stats.binom(20,.05)
sage: bar_chart([binom_dist.pmf(x) for x in range(21)])
Graphics object consisting of 1 graphics primitive

这个 bar_chart 函数执行直方图的一些职责。

Scipy的统计数据也可以做其他事情。在这里，我们找到一个早期数据集的中值（作为第50个百分位）。（我们用 Python int 在Numpy中解决一个bug。）

sage: scipy.stats.scoreatpercentile(my_data, int(50)) # random
0.51271641116183286

这里要记住的关键是查看文档！

  • 特别是对于Scipy来说，Sage中的所有内容都不是用简单的命令“包装”的，所以您可能需要做一些实验。

  • 改进这些文档将是让学生参与其中的一个很好的方法。

从Sage内部使用R

Sage还有其他几个部分具有统计功能，但到目前为止，最重要的是 R project ，这是统计分析的行业和学术标准 all 种类。

有几种方法可以访问R。

  • 最简单的方法之一就是 r() 把你想做成统计对象的东西，然后。。。

  • 通过使用R命令 r.method() 把它们交给Sage做进一步的处理。

以下Kruskal-Wallis测试的示例直接来自 r.kruskal_test? 在笔记本里。

sage: x=r([2.9, 3.0, 2.5, 2.6, 3.2]) # normal subjects
sage: y=r([3.8, 2.7, 4.0, 2.4])      # with obstructive airway disease
sage: z=r([2.8, 3.4, 3.7, 2.2, 2.0]) # with asbestosis
sage: a = r([x,y,z]) # make a long R vector of all the data
sage: b = r.factor(5*[1]+4*[2]+5*[3]) # create something for R to tell which subjects are which
sage: a; b # show them
 [1] 2.9 3.0 2.5 2.6 3.2 3.8 2.7 4.0 2.4 2.8 3.4 3.7 2.2 2.0
 [1] 1 1 1 1 1 2 2 2 2 3 3 3 3 3
Levels: 1 2 3

sage: r.kruskal_test(a,b) # do the KW test!
    Kruskal-Wallis rank sum test

data:  sage17 and sage33
Kruskal-Wallis chi-squared = 0.7714, df = 2, p-value = 0.68

看来我们不能拒绝无效假设。

认真使用R的最佳方法是使用所谓的“百分比指令”，简单地要求每个单元在R中完全求值。这是一个来自johnverzani的线性回归示例 simpleR 文本。注意R也使用 # 表示注释的符号。

sage: %r
....: x = c(18,23,25,35,65,54,34,56,72,19,23,42,18,39,37) # ages of individuals
....: y = c(202,186,187,180,156,169,174,172,153,199,193,174,198,183,178) # maximum heart rate of each one
....: png() # turn on plotting
....: plot(x,y) # make a plot
....: lm(y ~ x) # do the linear regression
....: abline(lm(y ~ x)) # plot the regression line
....: dev.off()     # turn off the device so it plots
Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x
   210.0485      -0.7977

null device
          1
prep/Quickstarts/../media/Rplot001.png

在R中得到一个完整的工作表（并且能够忽略 % )，您也可以 r 菜单中靠近顶部的选项 sage 在里面。

（还有另一个Python接口，称为 rpy2 接口，但我们目前不建议将其与Sage一起使用。）