Statistik

Kombinatorik

Funktionen zur Berechnung der Zahl von Permutationen, Kombinationen und Variationen finden sich im Basismodul, bzw. im Modul Random von Julia:

Funktion	Name
Fakultät	`factorial`
Binomialkoeffizient	`binomial`
Entfernung von Duplikaten	`unique`
Zufällige Permutation	`randperm`

Es lassen sich leicht Teilmengen innerhalb von Arrays bilden, da Julia Arrays und Bereiche als Indexargumente erlaubt:

julia> psq = [1; 2; 3; 5; 7; 11; 13];

julia> psq[2:2:end]
3-element Array{Int64,1}:
  2
  5
 11

julia> psq[[2, 5, 7]]
3-element Array{Int64,1}:
  2
  7
 13

Wie sich alle Permutationen einer Menge erzeugen lassen, beschreibt ein Blogbeitrag.

Eine Stichprobe aus einem Array a zieht man mit Random.shuffle. a braucht dabei kein Array zu sein, sondern ein beliebiges iterierbares Objekt, wie zum Beispiel ein Bereich.

import Random: shuffle
a = 5:5:100 # Folge: 5, 10, ... , 95, 100
umfang = 5
shuffle(a)[1:umfang]

Effizienter erhält man eien Stichprobe durch Bernoulli sampling mit Wahrscheinlichkeit p (wobei hier der Stichprobenumfang schwankt)

import Random: randsubseq
gesamt = 1:100_000
p = 0.01
randsubseq(gesamt, p)

Ränge bildet man durch zweimalige Anwendung von sortperm

import Random: sortperm
julia> bign = 7;

julia> seq = rand(100:200, bign)
7-element Array{Int64,1}:
 152
 176
 174
 109
 190
 156
 198

julia> rgseq = sortperm(sortperm(seq))
7-element Array{Int64,1}:
 2
 5
 4
 1
 6
 3
 7

Deskriptive Statistik

Grundfunktionen zur deskriptiven Statistik in Julia finden sich in Base und Statistics.

Funktion	Name
Minimum	`minimum`
Maximum	`maximum`
Quantil zu gegebenem `p`	`quantile`
Mittelwert	`mean`
Median	`median`
empirische Standardabweichung	`std`
empirische Varianz	`var`

Harmonisches Mittel und geometrische Mittel lassen sich leicht definieren als:

harcmean(v) = length(v)/sum(1 ./ v)

geommean(v) = prod(v .^ (1 / length(v)))

oder effizienter

geommean(v) = exp(1 / length(v) * sum(log.(v)))

Siehe auch harmmean und genmean im Paket StatsBase.

Das Paket Dataframes liefert statistische Kenndaten für jede Variable eines Datensatzes mit der Funktion describe.

import DataFrames; Dm = DataFrames;
df = Dm.DataFrame()
df[:A] = collect(1:20);
df[:Squares] = df[:A] .^2
df[:Rand] = rand(20);

Dm.describe(df)

3×8 DataFrames.DataFrame
│ Row │ variable │ mean     │ min      │ median   │ max      │ nunique │ nmissing │ eltype  │
├─────┼──────────┼──────────┼──────────┼──────────┼──────────┼─────────┼──────────┼─────────┤
│ 1   │ A        │ 10.5     │ 1        │ 10.5     │ 20       │         │          │ Int64   │
│ 2   │ Squares  │ 143.5    │ 1        │ 110.5    │ 400      │         │          │ Int64   │
│ 3   │ Rand     │ 0.417612 │ 0.056507 │ 0.374668 │ 0.981625 │         │          │ Float64 │

Häufigkeiten

Die Funktion count liefert die Anzahl von Werten, die eine angegebene Bedingung erfüllen. Ferner können Elemente mit bestimmten Werten mithilfe der Funktion filter extrahiert werden. Absolute und relative Häufigkeiten erhält man mit den Funktionen counts und proportions im Paket StatsBase.

v = rand(2000); count(x -> x > 0.6, v)

filter(x -> x > 0.6, v)

829-element Array{Float64,1}:
 0.813869
 0.780748
 0.965371
 0.661237
 ⋮       
 0.942066
 0.872244
 0.979791
 0.711869
 0.877063
 0.947413

graphische Darstellungen

Ein Blatt-Stiel-Diagramm ist sinnvoll für kleine Datensätze. Es verschafft einen Überblick über die Daten und enthält dennoch die gesamte Information. (Funktionsdefinition stemleaf)

Beispiel: Stundenmittelwerte von Temperaturdaten im Monat Mai.

import FileIO: load, save; using DataFrames
tmai = DataFrame(load("T-2010-05.csv", skiplines_begin = 2))
dattmai = collect(skipmissing(tmai[:T]))

stemleaf(dattmai .* 10)

 1 | 78
 2 | 1
 3 | 15588
 4 | 223333334456788
 5 | 0122233445666677788999
 6 | 0000111111111222222234444455556666666677777777788888889999
 7 | 000000000001111122222222233333333444444444455555556666666666777777888889999
 8 | 00000111111111222222223333333333444444444445555555566666678888888888999999
 9 | 000000000011111111111222222222222334444444444445555555555666666666666667777778888888999
10 | 00000011111122222333333344444444455555666666666777778888888889999
11 | 0000011122333445555556666666677778888888899999
12 | 000000011122222233333344445555556777777788888899999
13 | 00001111111222222222333333444455577777788889999
14 | 0000111222234444445556666677888999999
15 | 001111222344455557888999
16 | 01112223445556777888899
17 | 0001222222334446666678889
18 | 001223444466778899
19 | 01112234477899
20 | 0012346788999
21 | 012233344445669
22 | 566799
23 | 0257
24 | 112234569
25 | 22

Bei größeren Datensätzen werden die Daten in Intervallen zusammengefasst und als Balkendiagramm dargestellt.

import StatsBase: fit, Histogram; import Plots: bar
hgm = fit(Histogram, dattmai, closed = :right, nbins = 25)
Plots.bar(hgm)

Das Graphikpaket stellt auch direkt die Ausgabe eines Histograms aus dem Datensatz bereit.

import Plots: histogram
histogram(dattmai, nbins = 25)

histtemperature

Eine empirische Verteilungsfunktion erhält man mit als die kumulierte Summe:

ecf = cumsum(hgm.weights)
bar(ecf ./ ecf[end], grid = false, title = "Empirische Verteilungsfunktion",
label = "ecf")

empF

Normal-Quantil-Plot

Das Paket Distributions stellt unter anderem die Quantile der Normalverteilung bereit. Wenn x die zu untersuchenden Daten enthält, erhält man mit sort(x) eine aufsteigend sortierte Liste. Diese Werte sind aufzutragen gegen die Quantile der Normalverteilung:

import Distributions: Normal, quantile
distn = Normal(mean(dattmai), std(datmai))

bign = length(dattmai)
zq = [Distributions.quantile(n01, (i - 0.5)/bign) for i in 1:bign];

import Plots: plot, scatter!
plot(zq, zq, title="Normal-Quantil-Plot", label = "")
scatter!(zq, sort(dattmai), markersize = 1.2, label = "Temperatur")

normalquantil

Regressionsanalyse

Funktionen in Julia zur linearen Regression sollen anhand zufällig erzeugter Erhebungsdaten illustriert werden:

vx = [range(-7, stop=7, length=25);]
vy = [18.3, 24.8, 23.1, 19.7, 22.0, 27.4, 31.3, 36.5, 40.3, 37.2, 41.5, 39.7, 35.2, 34.8, 35.4, 34.9, 43.2, 46.1, 44.8, 43.7, 39.6, 43.1, 45.1, 54.8, 51.0]
scatter(vx, vy)

Beispiel-Erhebungsdaten

Gesucht werden Parameter $\alpha$ und $\beta$ , die den Abstand von $\hat{y}$ zu $y$ minimieren für den Ansatz: $\hat{y} = \alpha + \beta * y$ .

Für dieses eindimensionale Beispiel erweitern wir vx mit einer konstanten Spalte:

julia> vxm = hcat(vx, ones(size(vx, 1))

julia> import MultivariateStats: llsq
julia> (beta, alpha) = llsq(vxm, vy, bias = false)
2-element Array{Float64,1}:
  2.0103296703296705
 36.54000000000001
 
julia> scatter(vx, vy, xlim = (-8, 11), label = "y")
julia> plot!(vx, beta .* vx .+ alpha, title = "Daten und Ausgleichsgerade")

Daten und Ausgleichsgerade

Unter der Annahme, dass die Abweichungen der vy normalverteilt sind, sind auch die Schätzer $\hat{\alpha}$ für $\alpha$ und $\hat{\beta}$ für $\beta$ normalverteilt. Es gilt¹

$\frac{\hat{\alpha} - \alpha}{\hat{\sigma}_{\hat{\alpha}}} \sim t(n-2), \quad \frac{\hat{\beta} - \beta}{\hat{\sigma}_{\hat{\beta}}} \sim t(n-2), \quad$

Zu $\hat{\alpha}$ und $\hat{\beta}$ lässt sich jeweils die Standardabweichung und für gegebene Irrtumswahrscheinlichkeit perr ein Vertrauensintervall angeben (Code).

julia> include("DStLinReg")
julia> Lm = DStLinReg;

julia> stdav = sqrt(Lm.varalpha(alpha, beta, vx, vy))
0.8854284906417536

julia> stdbv = sqrt(Lm.varbeta(alpha, beta, vx, vy))
0.21049172464388108

julia> (alpha, beta)
(36.54000000000001, 2.0103296703296705)

julia> perr = 0.05
0.05

julia> Lm.cialpha(alpha, beta, vx, vy, perr)
(34.7083516143521, 38.37164838564793)

julia> Lm.cibeta(alpha, beta, vx, vy, perr)
(1.5748943622148752, 2.445764978444466)

Eine an R angelehnte Implementierung findet sich im Paket GLM.

Funktionen für multiple lineare Regression stellt das Paket MultivariateStats bereit.

Julia Version 1.0 (Funktionsumfang gleich 0.7)

Fahrmeir et al.: Statistik. Berlin 2001. Kapitel 12.1. ↩

Programming Julia