裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

Julia でデータ分析

2020年12月30日 | ブログラミング

統計解析

julia> describe(df)
1500×7 DataFrame
  Row │ variable  mean     min      median   max      nmissing  eltype   
      │ Symbol    Float64  Float64  Float64  Float64  Int64     DataType 
──────┼──────────────────────────────────────────────────────────────────
    1 │ X1        50.0424     6.38   50.01     92.58         0  Float64
    2 │ X2        49.9949     3.97   50.02     90.16         0  Float64
    3 │ X3        50.0425     6.8    50.03     92.79         0  Float64
   :

julia> combine(df, names(df) .=> sum)
1×1500 DataFrame
 Row │ X1_sum     X2_sum     X3_sum     X4_sum     ⋯
     │ Float64    Float64    Float64    Float64    ⋯
 ────┼───────────────────────
   1 │ 5.00424e6  4.99949e6  5.00425e6  5.00398e6  ⋯

julia> sum(df.X2) # 合計
4.99949053e6

julia> var(df.X3) # 不偏分散
99.8737108092948

julia> var(df.X3, corrected=false) # 分散
99.87271207218672

julia> std(df.X4) # 標準偏差
10.012003767127581

julia> cor(df.X1, df.X3) # ピアソンの積率相関係数
-0.0027621173217232155

julia> using StatsBase
[ Info: Precompiling StatsBase [2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91]

julia> corspearman(df.X1, df.X3) # スピアマンの順位相関係数
-0.0033731070923354164

julia> corkendall(df.X1, df.X3) # ケンドールの順位相関係数
-0.002256851537930442

散布図

julia> xy = df[:, [:X1, :X2]]
100000×2 DataFrame
    Row │ X1       X2      
        │ Float64  Float64 
────────┼──────────────────
      1 │   51.26    44.35
      2 │   41.03    59.43
      3 │   61.89    29.32
      4 │   56.88    59.77
   :

julia> plot(df.X1, df.X2, st=:scatter, xlabel="X1", ylabel="X2", label="", title="scattergram")

ヒートマップ

julia> using StatsBase

julia> using Plots

julia> fit(Histogram, (df.X1, df.X2)) |> Plots.plot

単回帰分析

julia> using GLM

julia> ols = lm(@formula(X2 ~ X1), df)
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

X2 ~ 1 + X1

Coefficients:
─────────────────────────────────────────────────────────────────────────────────
                    Coef.  Std. Error       t  Pr(>|t|)    Lower 95%    Upper 95%
─────────────────────────────────────────────────────────────────────────────────
(Intercept)  50.0038       0.161703    309.23    <1e-99  49.6868      50.3207
X1           -0.000177375  0.00316896   -0.06    0.9554  -0.00638851   0.00603376
─────────────────────────────────────────────────────────────────────────────────

重回帰分析

julia> ols = lm(@formula(X2 ~ X1 + X3 + X10), df)
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

X2 ~ 1 + X1 + X3 + X10

Coefficients:
─────────────────────────────────────────────────────────────────────────────────
                    Coef.  Std. Error       t  Pr(>|t|)     Lower 95%    Upper 95%
─────────────────────────────────────────────────────────────────────────────────
(Intercept)  49.63         0.275956    179.85    <1e-99  49.0892       50.1709
X1           -0.000161144  0.00316895   -0.05    0.9594  -0.00637226    0.00604997
X3            0.00151638   0.00316354    0.48    0.6317  -0.00468412    0.00771688
X10           0.00594262   0.00316205    1.88    0.0602  -0.000254952   0.0121402
─────────────────────────────────────────────────────────────────────────────────

回帰分析のまとめ

julia> response(ols) # 従属変数
100000-element Array{Float64,1}:
 44.35
 59.43
 29.32
 59.77
 56.75
 33.05
  :
 
julia> coeftable(ols)
──────────────────────────────────────────────────────────────────────────────────
                    Coef.  Std. Error       t  Pr(>|t|)     Lower 95%    Upper 95%
─────────────────────────────────────────────────────────────────────────────────
(Intercept)  49.63         0.275956    179.85    <1e-99  49.0892       50.1709
X1           -0.000161144  0.00316895   -0.05    0.9594  -0.00637226    0.00604997
X3            0.00151638   0.00316354    0.48    0.6317  -0.00468412    0.00771688
X10           0.00594262   0.00316205    1.88    0.0602  -0.000254952   0.0121402
─────────────────────────────────────────────────────────────────────────────────

julia> coefnames(ols) # 独立変数の名前
4-element Array{String,1}:
 "(Intercept)"
 "X1"
 "X3"
 "X10"

julia> coef(ols) # 回帰係数
4-element Array{Float64,1}:
 49.63004506573552
 -0.00016114387378327575
  0.0015163795173625358
  0.005942618126945406

julia> stderror(ols) # 偏回帰係数の標準誤差
4-element Array{Float64,1}:
 0.2759555477601921
 0.00316895479304485
 0.003163542394883129
 0.00316204541408164

julia> round.(stderror(ols), digits=5)
4-element Array{Float64,1}:
 0.27596
 0.00317
 0.00316
 0.00316

julia> vcov(ols) # 偏回帰係数の分散・共分散行列
4×4 Array{Float64,2}:
  0.0761515    -0.000504924  -0.000499421  -0.000497993
 -0.000504924   1.00423e-5    2.75781e-8    2.00964e-8
 -0.000499421   2.75781e-8    1.0008e-5    -5.57047e-8
 -0.000497993   2.00964e-8   -5.57047e-8    9.99853e-6

julia> confint(ols) # 回帰係数の信頼区間
4×2 Array{Float64,2}:
 49.0892       50.1709
 -0.00637226    0.00604997
 -0.00468412    0.00771688
 -0.000254952   0.0121402

julia> predict(ols) # 予測値
100000-element Array{Float64,1}:
 49.97502667684598
 49.998951264807296
 50.010448604441656
 49.95823034383087
 50.04644378530134
 50.04986504302499
 50.0524343170104
  :

julia> round.(predict(ols), digits=5)
100000-element Array{Float64,1}:
 49.97503
 49.99895
 50.01045
 49.95823
 50.04644
 50.04987
  :

julia> residuals(ols) # 残差
100000-element Array{Float64,1}:
  -5.625026676845977
   9.431048735192704
 -20.690448604441656
   9.811769656169133
   6.703556214698658
 -16.99986504302499
   :

julia> r2(ols) # 決定係数
3.775037972164608e-5

julia> adjr2(ols) # 自由度調整済み R2  R とちょっと値が異なる
7.75031223032574e-6

julia> aic(ols) # AIC
744259.4650342457

julia> aicc(ols) # 調整済み AIC
744259.4656342817

julia> bic(ols) # BIC
744307.0296615706

julia> loglikelihood(ols) # 対数尤度
-372124.7325171229

julia> deviance(ols) # デビアンス
9.994475506014155e6

julia> nulldeviance(ols) # ヌル・デビアンス
9.994852815503202e6

julia> nobs(ols) # 分析に用いたサンプルサイズ
100000.0

julia> dof(ols) # 独立変数の個数 + 2
5

julia> dof_residual(ols)  # 残差の自由度
99996.0


ヒストグラム(度数分布)

julia> fit(Histogram, df.X1)
Histogram{Int64,1,Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}
edges:
  5.0:5.0:95.0
weights: [2, 16, 93, 465, 1616, 4403, 9067, 15026, 19263, 19100, 14993, 9225, 4447, 1603, 540, 121, 18, 2]
closed: left
isdensity: false

julia> fit(Histogram, df.X1, nbins=10)
Histogram{Int64,1,Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}
edges:
  0.0:10.0:100.0
weights: [2, 109, 2081, 13470, 34289, 34093, 13672, 2143, 139, 2]
closed: left
isdensity: false

julia> using Plots

julia> plot(df.X1, st=:histogram, label="", xlabel="observed value", ylabel="frequencies", title="histogram of X1")

コメント   この記事についてブログを書く
« Julia のベンチマークテスト | トップ | Julia で統計学分布関数と仮... »
最新の画像もっと見る

コメントを投稿

ブログ作成者から承認されるまでコメントは反映されません。

ブログラミング」カテゴリの最新記事