Julia数据分析入门

编辑:李东风

日期:2020-04-04

Julia比较适合用作数值计算, 编程既有Python、R、Matlab这样的语言的简洁, 又有C++这样的编译语言的运行效率。 统计数据分析、作图需要用到许多复杂的算法, 有些算法耗时很多, 比如MCMC等。 大量数据的分析、计算、测试都需要易用的编程和高效的运行效率, Julia在这两点都很适合。

Julia用作统计数据分析, 缺点是其问世时间还比较短, 许多统计模型还没有现成的App包。 但是, 基本的数据管理、作图、统计分布、随机模拟、最优化等统计计算必须的功能已经很完善, 所以能够自己编写的统计计算任务很适合用Julia编写。

Julia目前的版本是1.4.0。

本文承接“Julia语言入门”一文, 从统计数据分析角度简单先容如何用Julia读入数据, 进行数据清理,汇总,简单的统计分析。

数据框

数据分析任务的大部分时间并不是消耗在建模、计算上面, 而是对数据进行整理。 这包括从网络、文本文件、关系数据库、NoSQL数据库等位置获取数据, 将其转换成统计过程需要的数据输入形式, 清理无效数据, 将数据格式进行一致化, 对数据进行简单的探索性分析(EDA), 比如各个属性(变量)的类型,分布,属性之间的关系, 观测之间的分类, 异常值(离群值)识别等。

统计分析最常见的输入数据格式是数据框(Data Frame), 这是与矩阵类似的数据结构, 但是可以同时保存不同类型的数据, 同一列的数据类型相同。 R语言中称为数据框, SAS语言中称为数据集(dataset)。 数据框的每一行称为一个观测, 每一列称为一个变量。

Julia语言的DataFrames包(package)提供了数据框数据类型。 如果没有安装此包可用如下命令安装:

import Pkg
Pkg.add("DataFrames")

数据框的有关数据类型

数据框中的列与统计学中随机变量的分类类似, 分为类别(categorical)变量与数值(continuous)变量。 类别变量主要用来对观测分组, 也可以作为模型中的自变量和因变量。 对类别变量不能进行通常的四则运算和其它数学计算。 字符型列只能作为类别变量, 某些离散取值的数值型列在建模时也可以作为类别变量, 比如,用1代表男性,2代表女性, 这样的列就应该当作类别变量使用。 对数值变量则可以进行各种数学运算。

统计数据与一般的数值计算问题有一个重要差别, 就是统计数据中常常含有缺失值, 如记录遗失、无应答、无交易、不适用等。 用来保存统计数据的数据结构必须能表示缺失值, 这是Julia基本的向量(一维数组)不能直接提供的功能。

DataFrames包调用了Missings包提供的一个Missing数据类型, 数据值missing就是Missing类型,表示缺失值。 这个数据类型允许将正常值与missing组合在一起成为一个向量,如

In [1]:
using DataFrames
[1,2,missing,4]
Out[1]:
4-element Array{Union{Missing, Int64},1}:
 1
 2
  missing
 4

数据框生成与访问

DataFrame()可以直接生成小的数据框,如

In [2]:
using DataFrames
da1orig = DataFrame(name=["张三", "李四", "王五", "赵六"], 
    age=[33, 42, missing, 51],
    sex=["M", "F", "M", "M"])
da1 = copy(da1orig)
Out[2]:

4 rows × 3 columns

nameagesex
StringInt64?String
1张三33M
2李四42F
3王五missingM
4赵六51M

其中copy()函数制作一个da1orig数据框的副本。 数据框属于可变数据类型(mutable type), 如果直接将da1orig赋值给一个变量da1, 则这两个变量实际指向同一个数据框, 改变其中一个变量会同时修改另一个变量的值。

df为数据框,和矩阵类似, 可以用size(df,1)求数据框行数(观测个数), 用size(df,2)求数据框列数(变量个数)。如

In [3]:
(size(da1, 1), size(da1, 2))
Out[3]:
(4, 3)

可以用first(df,k)取出df前面k行组成的子集, 用last(df,k)取出df最后k行组成的子集。

可以像对矩阵那样取数据框元素或子集。 比如, df[2,1]表示数据框df第2行第1列元素。 列序号也可以用变量名符号表示, 变量名符号就是将变量名前面添加冒号, 不用字符串表示,如

In [4]:
da1[1, :sex]
Out[4]:
"M"

用单独的叹号或冒号表示所有行或所有列,比如,df[!,2]df[:,2]取出df的第二列:

In [5]:
da1[!,2]
Out[5]:
4-element Array{Union{Missing, Int64},1}:
 33
 42
   missing
 51
In [6]:
da1[:, 2]
Out[6]:
4-element Array{Union{Missing, Int64},1}:
 33
 42
   missing
 51

可以用变量名指定取出一列,行下标位置写冒号或叹号,如

In [7]:
da1[!,:age]
Out[7]:
4-element Array{Union{Missing, Int64},1}:
 33
 42
   missing
 51

对行、列下标可以指定范围, 用下标向量或者变量名符号向量, 如

In [8]:
da1[2:4, [1,3]]
Out[8]:

3 rows × 2 columns

namesex
StringString
1李四F
2王五M
3赵六M
In [9]:
da1[2:4, [:name, :sex]]
Out[9]:

3 rows × 2 columns

namesex
StringString
1李四F
2王五M
3赵六M

CSV文件读写

用CSV包的read()函数可以读入CSV文件到数据框, 用write()函数将数据框保存到CSV文件。

设文件class9.csv内容如下:

name,sex,age,height,weight
Sandy,F,11,130.0,23.0
Karen,F,12,143.0,35.0
Kathy,F,12,152.0,38.0
Alice,F,13,144.0,38.0
Thomas,M,11,146.0,39.0
James,M,12,146.0,38.0
John,M,12,150.0,45.0
Robert,M,12,165.0,58.0
Jeffrey,M,13,159.0,38.0

用如下命令将其读入为数据框:

In [10]:
using CSV
d_class = DataFrame(CSV.File("class9.csv"))
Out[10]:

9 rows × 5 columns

namesexageheightweight
StringStringInt64Float64Float64
1SandyF11130.023.0
2KarenF12143.035.0
3KathyF12152.038.0
4AliceF13144.038.0
5ThomasM11146.039.0
6JamesM12146.038.0
7JohnM12150.045.0
8RobertM12165.058.0
9JeffreyM13159.038.0

如果数据项之间是用一个空格分隔的, 可以加选项用delim=' '的方式读入。 不支撑多个空格作为分隔。 如果数据项之间是用一个制表符分隔的, 可以加选项delim='\t')的方式读入。 如果数据中没有列名, 可以加选项datarow=1

将数据框写入CSV文件class2.csv,命令如

CSV.write("class2.csv", d_class)

简单修改

rename!()函数可以为数据框变量改名, 如

In [11]:
da1 = copy(da1orig)
rename!(da1, [:sex => :性别, :age => :年龄])
Out[11]:

4 rows × 3 columns

name年龄性别
StringInt64?String
1张三33M
2李四42F
3王五missingM
4赵六51M

可以直接修改变量或者添加新变量,如

In [12]:
dcl = copy(d_class)
dcl[!,:height] ./= 100
first(dcl, 3)
Out[12]:

3 rows × 5 columns

namesexageheightweight
StringStringInt64Float64Float64
1SandyF111.323.0
2KarenF121.4335.0
3KathyF121.5238.0
In [13]:
dcl = copy(d_class)
dcl[!,:ratio] .= dcl[!,:weight] ./ dcl[!,:height]
dcl[!,:classno] .= 1
first(dcl,3)
Out[13]:

3 rows × 7 columns

namesexageheightweightratioclassno
StringStringInt64Float64Float64Float64Int64
1SandyF11130.023.00.1769231
2KarenF12143.035.00.2447551
3KathyF12152.038.00.251

信息查询

利用Query包可以实现从各种数据源的信息查询, 这也包括从数据框查询。

注意:Query包目前还在Julia 0.6阶段,不能用在Julia 1.0.0版。

查询以@from宏开始, 由若干个查询命令连接而成, @from宏可以指定行子集条件, @select宏可以指定列子集并为结果列变量改名, 比如, 从d_class数据框中查询所有性别为女、年龄在12岁及以下的人的姓名、身高:

In [14]:
using Query

df_tmp = @from i in d_class begin
    @where i.sex=="F" && i.age <= 12
    @select {i.name, 体重=i.weight}
    @collect DataFrame
end
Out[14]:

3 rows × 2 columns

name体重
StringFloat64
1Sandy23.0
2Karen35.0
3Kathy38.0

查询的结果保存成了一个数据框。 如果省略@collect命令, 结果会是一个Julia迭代器, 可以用在for循环中遍历。

排序

sort!()函数将数据框安装某一列或某几列排序。 这个函数会修改其输入, 输入的数据框会被修改为排序后的值。

用第二个自变量指定按哪一列排序:

In [15]:
typeof(d_class)
Out[15]:
DataFrame
In [16]:
sort!(d_class, :age)
d_class
Out[16]:

9 rows × 5 columns

namesexageheightweight
StringStringInt64Float64Float64
1SandyF11130.023.0
2ThomasM11146.039.0
3KarenF12143.035.0
4KathyF12152.038.0
5JamesM12146.038.0
6JohnM12150.045.0
7RobertM12165.058.0
8AliceF13144.038.0
9JeffreyM13159.038.0

可以指定多列,当前一列相同时按后一列排序:

In [17]:
sort!(d_class, (:sex, :age))
d_class
Out[17]:

9 rows × 5 columns

namesexageheightweight
StringStringInt64Float64Float64
1SandyF11130.023.0
2KarenF12143.035.0
3KathyF12152.038.0
4AliceF13144.038.0
5ThomasM11146.039.0
6JamesM12146.038.0
7JohnM12150.045.0
8RobertM12165.058.0
9JeffreyM13159.038.0

在指定排序变量时, 可以用order()函数指定选项, 如rev=true指定降序, by=uppercase指定按转换为大写后排序。 例如

In [18]:
sort!(d_class, [:sex, order(:age, rev=true)])
d_class
Out[18]:

9 rows × 5 columns

namesexageheightweight
StringStringInt64Float64Float64
1AliceF13144.038.0
2KarenF12143.035.0
3KathyF12152.038.0
4SandyF11130.023.0
5JeffreyM13159.038.0
6JamesM12146.038.0
7JohnM12150.045.0
8RobertM12165.058.0
9ThomasM11146.039.0

横向合并

join()函数可以用来按照关键列对两个数据框进行横向合并。

一对一横向合并

In [19]:
dfj1 = DataFrame(id=["a", "b"], X=[11,12])
Out[19]:

2 rows × 2 columns

idX
StringInt64
1a11
2b12
In [20]:
dfj2 = DataFrame(id=["b", "a"], Y=[21,22])
Out[20]:

2 rows × 2 columns

idY
StringInt64
1b21
2a22
In [21]:
dfj3 = join(dfj1, dfj2, on = :id)
Out[21]:

2 rows × 3 columns

idXY
StringInt64Int64
1a1122
2b1221

一对多横向合并

In [22]:
dfj4 = DataFrame(id=["a", "a", "b"], Y=[41,42,43])
Out[22]:

3 rows × 2 columns

idY
StringInt64
1a41
2a42
3b43
In [23]:
dfj5 = join(dfj1, dfj4, on = :id)
Out[23]:

3 rows × 3 columns

idXY
StringInt64Int64
1a1141
2a1142
3b1243

多对多横向合并

In [24]:
dfj6 = DataFrame(id=["a", "a", "b"], X=[61,62,63])
Out[24]:

3 rows × 2 columns

idX
StringInt64
1a61
2a62
3b63
In [25]:
dfj7 = DataFrame(id=["a", "a", "b"], Y=[71,72,73])
Out[25]:

3 rows × 2 columns

idY
StringInt64
1a71
2a72
3b73
In [26]:
dfj8 = join(dfj6, dfj7, on = :id)
Out[26]:

5 rows × 3 columns

idXY
StringInt64Int64
1a6171
2a6172
3a6271
4a6272
5b6373

横向合并选项

join(a,b, on=...)默认将左右能匹配的观测保留, 不能匹配的删除, 这在数据库术语中称为“内连接”(inner join)。 可以用选项kind=指定类别, 如kind = :left, kind = :right, kind = :outer

长宽表转换

考虑如下数据:

"subject","x_1","x_2","x_3","y_1","y_2","y_3"
1,5,7,8,9,7,6
2,8,2,10,1,1,9
3,7,2,5,10,8,3
4,1,5,6,10,1,1
5,9,7,10,8,8,10

读入为数据框:

In [27]:
d_wide = DataFrame(CSV.File("widetab.csv"))
Out[27]:

5 rows × 7 columns

subjectx_1x_2x_3y_1y_2y_3
Int64Int64Int64Int64Int64Int64Int64
11578976
228210119
337251083
441561011
5597108810

这是5位病人3次去医院的记录数据, 每次去有xy两个测量值。

stack()函数将其中的测量值合并到一列value, 并增加一列表示值所对应的变量名的字符型列。 第二自变量指定要合并的列,可以用列号范围、列号向量、列号名称。 第三自变量可选,指定保持不变的区分病人的变量。 如

In [28]:
d_long4w = stack(d_wide, 2:7, :subject)
Out[28]:

30 rows × 3 columns

variablevaluesubject
SymbolInt64Int64
1x_151
2x_182
3x_173
4x_114
5x_195
6x_271
7x_222
8x_223
9x_254
10x_275
11x_381
12x_3102
13x_353
14x_364
15x_3105
16y_191
17y_112
18y_1103
19y_1104
20y_185
21y_271
22y_212
23y_283
24y_214
25y_285
26y_361
27y_392
28y_333
29y_314
30y_3105

也可以写成

d_long4w = stack(d_wide, [:x_1, :x_2, :x_3, :y_1, :y_2, :y_3], :subject)

这个表格的形式虽然是长表, 但并不是大家最终所求。 大家希翼将变量名与随访时间分为两个变量, 可用如下代码:

In [29]:
d_long4w[!,:time] = [parse(Int, String(x)[3:3]) for x in d_long4w[!, :variable] ];
d_long4w[!, :variable] = [String(x)[1:1] for x in d_long4w[!, :variable] ];

按照病人号、随访时间、测量值名排序,长表现在变成了:

In [30]:
sort!(d_long4w, (:subject, :time, :variable) )
d_long4w
Out[30]:

30 rows × 4 columns

variablevaluesubjecttime
StringInt64Int64Int64
1x511
2y911
3x712
4y712
5x813
6y613
7x821
8y121
9x222
10y122
11x1023
12y923
13x731
14y1031
15x232
16y832
17x533
18y333
19x141
20y1041
21x542
22y142
23x643
24y143
25x951
26y851
27x752
28y852
29x1053
30y1053

unstack()可以将放在同一列的变量值合并到一行中。 第二自变量表示按照那些变量分组, 第三自变量表示保存了变量名的列(类型是Symbol), 第四自变量是保存的值。 如

In [31]:
d_long4w2 = unstack(d_long4w, [:subject, :time], :variable, :value)
Out[31]:

15 rows × 4 columns

subjecttimexy
Int64Int64Int64?Int64?
11159
21277
31386
42181
52221
623109
731710
83228
93353
1041110
114251
124361
135198
145278
15531010

如果使用R的tidyr包, 可以用pivot_longerpivot_wider更容易地完成转换。 SAS的PROC TRANSPOSE也可以比较容易地完成这个数据的转换。

总之, Julia DataFrames的宽表变长表用stack()函数, 长表变宽表用unstack()函数, 拆分变量名与时间可以用字符串功能以及类型转换功能来实现。

变量名在Julia中是Symbol数据类型, 不是字符串。 用String()将Symbol类型转换成字符串类型, 用parse(Int, s)将字符串s转换成整数, 用parse(Float64, s)将字符串s转换成浮点数。

缺失值管理

Julia语言提供了missing表示缺失值, missing的功能与R语言中的NA, 数据库SQL语言中的null类似。

数据框中可以用missing表示缺失元素, 见上面的da1数据框。 但是,如果数据框中某列原来没有缺失值, 则不能将其中的值赋值为missing

ismissing()判断某个值是否缺失值, 加点以后可以判断某个向量的每个元素,如

In [32]:
da1 = copy(da1orig)
ismissing.(da1[!,:age])
Out[32]:
4-element BitArray{1}:
 0
 0
 1
 0

函数any()可以判断布尔值向量元素是否存在真值, all()可以判断布尔值向量元素是否都是真值。如

In [33]:
any(ismissing.(da1[!,:age]))
Out[33]:
true

missing参与的四则运算、比较、数学计算一般返回缺失值。 用skipmissing()可以在计算时删除指定列中的缺失值进行计算, 如

In [34]:
da1[!,:age]
Out[34]:
4-element Array{Union{Missing, Int64},1}:
 33
 42
   missing
 51
In [35]:
sum(da1[!,:age])
Out[35]:
missing
In [36]:
sum(skipmissing(da1[!,:age]))
Out[36]:
126

Missings.replace()可以在计算中用指定值替换缺失值,如

In [37]:
sum(Missings.replace(da1[!,:age], 0))
Out[37]:
126

为了获得整个数据框df中不含任何缺失值的观测的序号向量, 用DataFrames.completecases(df)函数。 将整个数据框中含有缺失值的观测都删去, 用DataFrames.dropmissing!(df)函数。 如

In [38]:
da1 = copy(da1orig)
DataFrames.dropmissing!(da1)
da1
Out[38]:

3 rows × 3 columns

nameagesex
StringInt64String
1张三33M
2李四42F
3赵六51M

为了将去掉缺失值后的内容转换成一个普通向量, 用collect()函数,如

In [39]:
da1 = copy(da1orig)
collect(skipmissing(da1[!,:age]))
Out[39]:
3-element Array{Int64,1}:
 33
 42
 51

含有缺失值的向量实际是Array{Union{Missing, T},1}类型, 其中T是String, Int64, Float64等值类型。如

In [40]:
typeof(da1[!,:age])
Out[40]:
Array{Union{Missing, Int64},1}

Missings.nonmissingtype(eltype(x))函数可以求带有missing值的向量x中非缺失值的数据类型,如

In [41]:
Missings.nonmissingtype(eltype(da1[!,:age]))
Out[41]:
Int64

这样,为了从不允许缺失值的向量生成允许缺失值的向量,可以用如下方法进行转换:

In [42]:
x = copy(da1orig[!,:sex])
x = convert(Array{Union{Missing, eltype(x)}, 1}, x)
Out[42]:
4-element Array{Union{Missing, String},1}:
 "M"
 "F"
 "M"
 "M"

可以将上述操作写成自定义函数,并作用到指定列:

In [43]:
allowmissing(x) = convert(Array{Union{Missing, eltype(x)}, 1}, x)
da1 = copy(da1orig)
da1[!,:sex] = allowmissing(da1[!,:sex])
Out[43]:
4-element Array{Union{Missing, String},1}:
 "M"
 "F"
 "M"
 "M"

反过来,为了将允许有缺失值但是实际上不包含缺失值的向量转换成不允许有缺失值的向量,程序如

In [44]:
denymissing(x) = convert(Array{Missings.nonmissingtype(eltype(x)), 1}, x)
da1 = copy(da1orig)
da1[!,:sex] = denymissing(da1[!,:sex])
Out[44]:
4-element Array{String,1}:
 "M"
 "F"
 "M"
 "M"

事实上,如果df是数据框,则DataFrames.allowmissing!(df)使得数据框中所有列都允许取缺失值。 DataFrames.allowmissing!(df, col)可以使得第col列取缺失值, DataFrames.allowmissing!(df, cols)可以使得下标向量cols指定的那些列允许取缺失值。 DataFrames.disallowmissing!()函数起到反向的作用。 如

In [45]:
da1 = copy(da1orig)
DataFrames.allowmissing!(da1, 2:3)
typeof(da1[!,:sex])
Out[45]:
Array{Union{Missing, String},1}

数据框变量概括

describe()函数对数据框各个变量进行简单概括。如

In [46]:
describe(d_class)
Out[46]:

5 rows × 8 columns

variablemeanminmedianmaxnuniquenmissingeltype
SymbolUnion…AnyUnion…AnyUnion…NothingDataType
1nameAliceThomas9String
2sexFM2String
3age12.01112.013Int64
4height148.333130.0146.0165.0Float64
5weight39.111123.038.058.0Float64

describe()对每个变量给出了数据类型, 缺失值个数和缺失百分比, 对字符型变量给出了不同值个数, 对数值型变量给出了数据个数、均值、最小值、最大值、中位数、四分之一分位数和四分之三分位数。 类似于R语言的summary()函数。

分组汇总

对数据框经常需要按某一个或几个分类变量进行分类汇总。 用by()函数进行分类汇总, 第一自变量是输入数据框, 第二自变量是分组变量, 第三自变量是汇总用的函数。 汇总用函数的输入是一组的子数据框。 如

In [47]:
by(d_class, :sex, size)
Out[47]:

2 rows × 2 columns

sexx1
StringTuple…
1F(4, 5)
2M(5, 5)

可以用无名函数作为第三自变量,如

In [48]:
using StatsBase
by(d_class, :sex, 
   df -> DataFrame(
        n=size(df, 1), 
        mh=mean(df[!, :height]), 
        mw=mean(df[!, :weight])))
Out[48]:

2 rows × 4 columns

sexnmhmw
StringInt64Float64Float64
1F4142.2533.5
2M5153.243.6

按多个分类变量分组,如

In [49]:
by(d_class, [:sex, :age], 
   df -> DataFrame(
        n=size(df, 1), 
        maxh=maximum(df[!, :height]), 
        maxw=maximum(df[!, :weight])))
Out[49]:

6 rows × 5 columns

sexagenmaxhmaxw
StringInt64Int64Float64Float64
1F131144.038.0
2F122152.038.0
3F111130.023.0
4M131159.038.0
5M123165.058.0
6M111146.039.0

分类变量

某些统计模型允许使用分类变量, 比如, 逻辑斯谛回归的因变量是分类变量。 用字符串表示分类变量很直观, 但是不方便用在数学模型计算当中, 在数据量很大时保存许多有重复的字符串往往也比较浪费存储空间。

CategoricalArrays包提供了CategoricalArray(分类数组)类型, 专门用来表示分类变量。 这种类型类似于R的因子(factor)类型, 将分类变量的不同值编码为整数号码, 然后保存这些号码以及号码与变量实际值之间的对应表。 用CategoricalArray()函数将分类变量转换成分类数组类型, 如

In [50]:
using CategoricalArrays
dcl = copy(d_class)
dcl[!,:sex] = CategoricalArray(dcl[!,:sex])
Out[50]:
9-element CategoricalArray{String,1,UInt32}:
 "F"
 "F"
 "F"
 "F"
 "M"
 "M"
 "M"
 "M"
 "M"

levels()求分类数组的各个不同值,如

In [51]:
levels(dcl[!,:sex])
Out[51]:
2-element Array{String,1}:
 "F"
 "M"

实际上, 还可以用categorical!()函数将数据框中指定的某一列或几列转换为因子,如

In [52]:
dcl = copy(d_class)
categorical!(dcl, [:sex, :age])
Out[52]:

9 rows × 5 columns

namesexageheightweight
StringCategorical…Categorical…Float64Float64
1AliceF13144.038.0
2KarenF12143.035.0
3KathyF12152.038.0
4SandyF11130.023.0
5JeffreyM13159.038.0
6JamesM12146.038.0
7JohnM12150.045.0
8RobertM12165.058.0
9ThomasM11146.039.0
In [53]:
levels(dcl[!,:sex])
Out[53]:
2-element Array{String,1}:
 "F"
 "M"
In [54]:
levels(dcl[!,:age])
Out[54]:
3-element Array{Int64,1}:
 11
 12
 13

DataFrames.eltype.(eachcol(df))函数求数据框df每一列的类型。 如

In [55]:
eltype.(eachcol(dcl))
Out[55]:
5-element Array{DataType,1}:
 String
 CategoricalString{UInt32}
 CategoricalValue{Int64,UInt32}
 Float64
 Float64

日期和时间类型

有些数据中包括日期, 或者日期时间。 Julia的Dates包提供了Date类型和DateTime类型。 日期不包含时区, 可以调用TimeZones包的功能来支撑时区。

可以从年月日的整数值用Date()转换成日期, 也可以从日期字符串按照某个模板转换成日期。 如:

In [56]:
import Dates
Dates.Date(2018)
Out[56]:
2018-01-01
In [57]:
Dates.Date(2018, 10)
Out[57]:
2018-10-01
In [58]:
Dates.Date(2018, 10, 31)
Out[58]:
2018-10-31
In [59]:
Dates.Date("2018-10-31", "y-m-d")
Out[59]:
2018-10-31
In [60]:
Dates.Date("20181031", "yyyymmdd")
Out[60]:
2018-10-31
In [61]:
Dates.Date.([2018, 2018], [3, 10], [15, 31])
Out[61]:
2-element Array{Dates.Date,1}:
 2018-03-15
 2018-10-31
In [62]:
Dates.Date.(["20180315", "20181031"], "yyyymmdd")
Out[62]:
2-element Array{Dates.Date,1}:
 2018-03-15
 2018-10-31

可以用DateTime()函数将年、月、日、时、分、秒、毫秒整数值转换成日期时间, 精确到1毫秒。 如

In [63]:
Dates.DateTime(2018, 10, 31)
Out[63]:
2018-10-31T00:00:00
In [64]:
Dates.DateTime(2018, 10, 31, 12, 15, 30)
Out[64]:
2018-10-31T12:15:30
In [65]:
Dates.DateTime(2018, 10, 31, 12, 15, 30, 136)
Out[65]:
2018-10-31T12:15:30.136

两个日期或者时间可以比较大小, 可以相减。结果带有单位, 用Dates.value()转换为表示天数或者毫秒数的整数值。 如

In [66]:
Dates.Date(2018, 10, 31) - Dates.Date(2018)
Out[66]:
303 days
In [67]:
Dates.value(Dates.Date(2018, 10, 31) - Dates.Date(2018))
Out[67]:
303
In [68]:
Dates.DateTime(2018, 10, 31, 12, 15, 30, 136) - Dates.DateTime(2018)
Out[68]:
26223330136 milliseconds
In [69]:
Dates.value(Dates.DateTime(2018, 10, 31, 12, 15, 30, 136) - Dates.DateTime(2018)) / (3600*1000*24)
Out[69]:
303.510765462963

Date和DateTime还有许多伴随的函数和功能, 比如, 从日期和时间中提取其组成部分, 星期几, 一个月中第几个星期几, 所在月的天数, 时间区间计算, 等等。 详见Julia手册中标准库部分的日期时间库说明。

数据框的其它功能

DataFrames包中还有许多函数,详见其文档的函数部分。 可用于DataFrame类型的函数包括(未特别说明都属于DataFrames包):

  • Base中的filter, filter!, join, similar, sort, sort!, unique
  • aggregate
  • allowmissing!, disallowmissing!调整列使得其允许或不允许包含缺失值
  • by将观测分组汇总
  • colwise对每列实行一些函数
  • combine, groupby可以将数据框分组处理后合并处理结果
  • completecases, dropmissing, dropmissing!用于获得完全观测
  • eachrow用来逐行遍历数据框
  • eltype.(eachcol(df))获得每列的数据类型
  • first(df, k), last(df, k)返回数据框头部或者尾部的若干行
  • melt, meltdf, stack, stackdf用于将宽表变成长表, unstack用于将长表变成宽表
  • names!为数据框修改所有变量名
  • unique, unique!, nonunique用于处理重复行或重复的变量值组合
  • rename, rename!指定修改某些变量名
  • describe对每列进行简单描述

描述统计

有一系列与统计有关的扩展包, 包括DataFrames, CSV, StatsBase, Distributions, HypothesisTests, KernelDensity, Loess, GLM, CategoricalArrays, Distances, Clustering, MultivariateStats, TimeSeries等。 StatsBase包提供了一些常见的统计函数。

Julia标准库一个Statistics包, 包含了均值、中位数、分位数、标准差、方差、相关系数等函数。

在一个统计分析项目中, 统计学家对手头的数据往往只有大致的了解, 但是每个变量是什么类型、能取什么样的值, 是分类变量还是连续变量, 分布情况, 变量之间的关系如何, 观测有无明显的分组, 这些都需要了解。 这样的工作称为探索性数据分析(EDA)。 Julia的一些函数可以用来进行这样的分析。

x是数值型向量, 或者用skipmissing()处理过的含有missing值的向量。 Base包和StatsBase包中的函数提供了一些简单概括函数, 如求和、平均、标准差等。

为了查看某个函数有哪些方法以及各个方法来自那些源文件, 可以调用methods()函数,如:

In [70]:
using StatsBase
methods(mean)
Out[70]:
# 8 methods for generic function mean:

sum(x)求和,prod(x)求乘积。

mean(x)求均值,median(x)求中位数。

std(x)求标准差,var(x)求方差,variation(x)求变异系数, mad(x)求平均绝对偏差, iqr(x)返回四分位间距(或称四分位极差), 即四分之三分位数减去四分之一分位数。

skewness(x)计算偏度, kurtosis(x)计算超额峰度。 从源程序看,偏度计算公式为 $$ \text{skewness} = \frac{\frac{1}{n} \sum (x_i - \bar x)^3} {\left( \frac{1}{n} \sum (x_i - \bar x)^2 \right)^{3/2}} $$ 峰度计算公式类似,是标准化的四阶矩减去3。

minimum(x)求最小值,maximum(x)求最大值。 quantile(x)返回5数概括(最小值、四分之一分位数、中位数、 四分之三分位数、最大值), quantile(x, p)返回对应于概率p的分位数。如

In [71]:
quantile(rand(1000), [0.01, 0.99])
Out[71]:
2-element Array{Float64,1}:
 0.00850012747502421
 0.9878773035560613

其中rand()函数用来生成标准均匀分布随机数向量。

zscore(x)返回x的Z分数。

对整数向量xsort(unique(x))获得x的所有不同值, 而counts(x)获得这些值的出现次数, proportions(x)获得这些值的出现比例。 对一般向量x, 可以用countmap(x)函数获得不同值到对应计数的字典。

有几个求秩得分的函数。 ordinalrank()返回同名次时不同秩的结果, competerank()返回同名次时同秩的结果, tiedrank()返回同名次时平均秩的结果。

对向量x, y, cor(x,y)计算相关系数, corspearman(x, y)计算Spearman秩相关系数, corkendall(x, y)计算Kendall τ统计量。

autocor(x)计算看成时间序列的向量x的自相关函数列(ACF)。

置信区间和假设检验

HypothesisTests包提供了基本的置信区间和假设检验计算。

单正态总体方差已知时均值的Z检验

OneSampleZTest()函数实行此检验。 下面模拟生成N($10,2^2)$分布的100个样本点, 在已知$\sigma=2$的条件下对$H_0: \mu=11$作双侧检验Z检验:

In [72]:
using Random, Distributions
using StatsBase
using HypothesisTests
Random.seed!(101)
x = rand(Normal(10, 2), 100);
tres11 = OneSampleZTest(mean(x), 2.0, length(x), 11.0)
Out[72]:
One sample z-test
-----------------
Population details:
    parameter of interest:   Mean
    value under h_0:         11.0
    point estimate:          10.258558563291658
    95% confidence interval: (9.8666, 10.6506)

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           0.0002

Details:
    number of observations:   100
    z-statistic:              -3.707207183541712
    population standard error: 0.2

可以用StatsBase.confint()函数从以上检验结果中抽取置信区间, 置信区间不依赖于检验所用的$\mu_0$:

In [73]:
confint(tres11)
Out[73]:
(9.866565766383646, 10.65055136019967)

可以用pvalue()函数从检验结果中提取p值,如

In [74]:
pvalue(tres11)
Out[74]:
0.0002095574982142723

pvalue()函数中可以用tail = :left指定左侧检验, 用tail = :right指定右侧检验。如

In [75]:
pvalue(tres11, tail = :left)
Out[75]:
0.00010477874910713615

单正态总体均值的t检验

OneSampleZTest()函数实行此检验。 如

In [76]:
Random.seed!(101)
x = rand(Normal(10, 2), 100);
tres12 = OneSampleTTest(x, 11.0)
Out[76]:
One sample t-test
-----------------
Population details:
    parameter of interest:   Mean
    value under h_0:         11.0
    point estimate:          10.258558563291658
    95% confidence interval: (9.8031, 10.7141)

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           0.0017

Details:
    number of observations:   100
    t-statistic:              -3.2297939849657924
    degrees of freedom:       99
    empirical standard error: 0.22956307435076084

单个比例的假设检验

BinomialTest()函数作此检验。

假设随机抽查的1000个人中有750个人支撑某项政策, 设总体支撑率为$p$, 要检验零假设$H_0: p=0.8$。 可用如下程序:

In [77]:
tres13 = BinomialTest(750, 1000, 0.80)
Out[77]:
Binomial test
-------------
Population details:
    parameter of interest:   Probability of success
    value under h_0:         0.8
    point estimate:          0.75
    95% confidence interval: (0.7219, 0.7766)

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           0.0001

Details:
    number of observations: 1000
    number of successes:    750

利用检验结果可以提取置信区间和p值。 比例的推断有各种不同的方法, 所以confint()函数和pvalue()函数可以加选项method=,包括:

  • method = :clopper_pearson: 这是利用二项分布作精确检验,求置信区间比较保守;
  • method = :wald: 利用正态近似,当总体比例接近0或者1时效果会很差;
  • method = :wilson: Wilson置信区间是一种改进的正态近似;
  • method = :jeffrey: 是基于Jeffrey无信息先验的贝叶斯credible区间;
  • method = :agresti_coull: 也是对正态近似的改进;
  • method = :arcsine: 用反正弦变换使得分布方差不依赖于比例的方法。

如:

In [78]:
confint(tres13, tail=:both, method=:wilson)
Out[78]:
(0.7222397197409138, 0.7758469010163087)
XML 地图 | Sitemap 地图