So-net無料ブログ作成
検索選択
確率・統計、データマイニング ブログトップ
前の30件 | -

統計フォーラム2017 [確率・統計、データマイニング]

イベントの紹介です。

統計フォーラム2017
http://adv-analytics.com/statisticsday/

10月18日は、統計の日ということらしい。
大枠は決まっているんだけど、内容はこれから作っていきますw

nice!(8)  コメント(0) 
共通テーマ:学問

ネスティッドロジットモデル [確率・統計、データマイニング]

ネスティッドロジットモデル(Nested Logit)をoptimを使って解きました。
階層構造(入れ子)になっている選択を通常のロジットで解くとIIA問題があり、上手く推定できません。

対処方法としては、プロビットモデルで解くか、ネスティッドロジットモデル(Nested Logit)を使って解くかのアプローチになるわけですが、ネスティッドロジットモデルに対応しているパッケージなどが少なく、自分で尤度関数を作成して解くのが一番分かりやすそうです。
Rだとmlogitというパッケージがありますね。

【ポイント 1】
log sum で上がってくる部分のλは0≦λ≦1である必要があるため、
exp(bx) / (1 + exp(bx))
という形式で与えることで、0≦λ≦1を保証しています。

【ポイント 2】
optim 関数はデフォルトで最小化を行うため、control 引数に fnscale = -1 を代入して、最大化を行うようにする

【Nested LogitのRコード】
hh <- nrow(dat_ss) # レコード数
b0 <- c(0, 0, 0, 0, 0) # パラメータの初期化

# ロジットモデルの尤度関数
fr <- function(x)
{
b1 <- x[1] # 階層1の説明変数 1
b2 <- x[2] # 階層1の説明変数 2
b3 <- x[3] # 階層2の説明変数 1
b4 <- x[4] # 階層2の説明変数 2
bx <- x[7] # 0 <= bx <= 1

LL = 0
for(i in 1:hh) {
# 効用の計算
  U_K2 <- b3 * dat_ss[i,3] + b4*dat_ss[i,4]
U_K1 <- b1 * dat_ss[i,1] + b2*dat_ss[i,2] + (exp(bx)/(1+exp(bx)))*log(1+exp(U_K2))

# 選択確率の計算
PP1 <- exp(U_K1)/(exp(U_K1)+1)
PP2 <- exp(U_K2)/(exp(U_K2)+1)

# 対数尤度の計算
LLL <- dat_ss[i,5]*log(PP1)+(1-dat_ss[i,5])*log(1-PP1)+dat_ss[i,5]*(dat_ss[i,6]*log(PP2)+(1-dat_ss[i,6])*log(1-PP2))
LL <- LL + LLL
}
return(LL)
}

# 対数尤度関数の最大化
res_nl <- optim(b0, fr, method = "BFGS", hessian = TRUE, control = list(fnscale = -1))

nice!(44)  コメント(0)  トラックバック(0) 
共通テーマ:学問

Rを使って重回帰分析を色々な方法で解く その2 [確率・統計、データマイニング]

Rを使って重回帰分析を色々な方法で解く その1
http://skellington.blog.so-net.ne.jp/2017-05-30

こちらの続き。
汎用最適化関数 optim() を使って解く方法。

### optim を使ったパラメータ推定 ###
b0 <- c(0, 0, 0, 0) # 推定するパラメータの初期化
num_row <- nrow(dat) # レコード数

y_ols = dat$Sepal.Length
x_ols = cbind(1, dat$Sepal.Width, dat$Petal.Length, dat$Petal.Width)

ols <- function (y_ols, x_ols) {
min_rss <- function (param) {
sum((y_ols - x_ols %*% param)^2)
}
k <- ncol(x_ols)
return(optim(par = rep(0, k), fn = min_rss, method = "BFGS", hessian = TRUE))
}

result <- ols(y_ols, x_ols)

result
$par
[1] 2.3518922 0.6548347 0.2375590 0.2521268

$value # 対数尤度値
[1] 2.586648

$counts
function gradient
24 8

$convergence # 収束判定
[1] 0

$message
NULL

$hessian
[,1] [,2] [,3] [,4]
[1,] 100.0 342.80 146.20 24.60
[2,] 342.8 1189.20 502.32 85.24
[3,] 146.2 502.32 216.70 36.56
[4,] 24.6 85.24 36.56 7.14

nice!(1)  コメント(0)  トラックバック(0) 
共通テーマ:学問

Rを使って重回帰分析を色々な方法で解く その1 [確率・統計、データマイニング]

まずは、パッケージ lm() を使って解く方法と行列を使って解く方法から。

# irisのデータを使用
head(iris)

# Species == "setosa" 50レコードを使用
dat <- subset(iris, Species == "setosa", c(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width ))

# 相関行列
round(cor(dat), 3)

Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 1.000 0.743 0.267 0.278
Sepal.Width 0.743 1.000 0.178 0.233
Petal.Length 0.267 0.178 1.000 0.332
Petal.Width 0.278 0.233 0.332 1.000


### lm を使った重回帰式 ###
# Length = β1 * Width + β2 * Length + β3 *Width + β0
fit.lm <- lm(Sepal.Length ~ . , data = dat)
summary(fit.lm)

Call:
lm(formula = Sepal.Length ~ ., data = dat)

Residuals:
Min 1Q Median 3Q Max
-0.40662 -0.17721 0.01222 0.13388 0.49693

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.35189 0.39287 5.986 3.03e-07 ***
Sepal.Width 0.65483 0.09245 7.083 6.83e-09 ***
Petal.Length 0.23756 0.20802 1.142 0.259
Petal.Width 0.25213 0.34686 0.727 0.471
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2371 on 46 degrees of freedom
Multiple R-squared: 0.5751, Adjusted R-squared: 0.5474
F-statistic: 20.76 on 3 and 46 DF, p-value: 1.192e-08


### 行列を使って解く方法 ###
solve(t(x_ols) %*% x_ols) %*% t(x_ols) %*% y_ols
[,1]
[1,] 2.3518898
[2,] 0.6548350
[3,] 0.2375602
[4,] 0.2521257

nice!(42)  コメント(0)  トラックバック(0) 
共通テーマ:学問

コンジョイント分析 [確率・統計、データマイニング]

8章 コンジョイント分析について


現代マーケティング・リサーチ -- 市場を読み解くデータ分析

現代マーケティング・リサーチ -- 市場を読み解くデータ分析

  • 作者: 照井 伸彦
  • 出版社/メーカー: 有斐閣
  • 発売日: 2013/11/22
  • メディア: 単行本(ソフトカバー)



よくやる間違い
属性水準(説明変数)として
・バッテリ時間
・保証期間
・色
・価格
を考えたときに、「価格」が入っているのは間違い

ダメな理由は、コンジョイント分析は、属性間の独立性を担保する必要があり、
価格は他の属性との相関があるため使えない。

nice!(48)  コメント(0)  トラックバック(0) 
共通テーマ:日記・雑感

重回帰で回帰係数の符号が逆転している場合の意味 [確率・統計、データマイニング]

個々の変数の相関が正(負)なのに、重回帰式において、その係数が負(正)になることがあります。

そんな時に、どう考えるか?

試しにネットに落ちている中古車データを拾ってきました。

20170517_img01.jpg

相関係数を計算すると

20170517_img02.jpg

となります。

走行距離や年数が古くなれば、中古車価格は下落するという事が分かります。

ここで、中古車価格を目的変数に、走行距離と年数を説明変数にすると
 中古車価格 = -10.0 × 走行距離 - 10.7 × 年数 + 252.0
という重回帰式が得られます。

走行距離が1万km増えると、中古車価格は10.0万円下落し、
同様に年数が1年増えると、中古車価格は10.7万円下落することが分かります。

さて、次に、年数を目的変数に、走行距離と中古車価格を説明変数にすると
 年数 = -0.0499 × 価格 - 0.154× 走行距離 + 13.86
という結果になりました。

年数と走行距離の相関係数は、0.562と正の値なのですが、
回帰式のパラメータは、マイナスとなっています。

つまり、「走行距離が1万キロ増えると、年数が0.154年減少する」となっていて感覚と値が違う結果になってしまいます。

これをどう考えるか?です。

重回帰の場合、解釈としては、注意が必要で、「他の変数が同じ値だった場合(他の変数を影響を統制した場合)の効果」をみていることになります。

つまり、「走行距離が同じであったとしたら、走行距離が1万キロ増えた場合、年数はどうなるか?」を表していることになります。

例として、ある車種が、
 【車種 A】走行距離:10万キロ、価格:100万円、年数:10年
だったとした場合、
 【車種 B】走行距離:20万キロ、価格:100万円、年数:???年
車種Bは、20万キロも走っているけど、車種Aと同じです。
つまり、年数が10年ではなく、もっと小さい数字であることが期待されます。

この構造を表しているのが重回帰のパラメータということになります。

一方、一般的には、走行距離が大きくなれば、年数も同時に大きくなる(正の相関がある)という事を表しているのが、相関係数となります。

nice!(38)  コメント(0)  トラックバック(0) 
共通テーマ:学問

KDD2017の申し込みが開始! [確率・統計、データマイニング]

KDD2017の申し込みが開始!
http://www.kdd.org/kdd2017/registration

今年はカナダのHalifaxです。

6/28までに申込すると、早期割引になりますね。
学割もかなりの額なので、学割で行くべきか悩みます。
(一応、大学院です。。。)

nice!(2)  コメント(0)  トラックバック(0) 
共通テーマ:学問

統計手法から統計モデリングへ その3 モデルの違いをパス図で理解 [確率・統計、データマイニング]

# パッケージの読み込み
library(lavaan)
 
# データの準備
lower <- '
1.000
0.774 1.000
0.383 0.496 1.000
'
full <- getCov(lower, names=c("Y", "X1", "X2"))




モデル 1
# モデル1
model_1 <- '
X1 ~ Y
X2 ~ X1
X2 ~ Y
X1 ~~ X1; X2 ~~ X2
'
 
fit_1 <- lavaan(model_1, sample.cov = full, sample.nobs = 20)
summary(fit_1, standardized=T)


ここで Y から X2のパス係数は有意でないので次の様なモデルを考える。



モデル 2
# モデル2
model_2 <- '
X1 ~ Y
X2 ~ X1
X1 ~~ X1; X2 ~~ X2
'
 
fit_2 <- lavaan(model_2, sample.cov = full, sample.nobs = 20)
summary(fit_2, standardized=T)






ここで面白いのが、モデル1とモデル2のYからX1へのパス係数の値は同じである。
そして、X1からX2のパス係数の値は、0.498から0.496に変化をしている。

これは、0.496 = 0.498 - 0.002ということだが、X2へのYの影響は、全てX1を経由しての「間接効果」を考えていることになる。

nice!(6)  コメント(0)  トラックバック(0) 
共通テーマ:学問

統計手法から統計モデリングへ その2 熱伝搬モデル [確率・統計、データマイニング]

y = a + b * x
という単回帰分析を考えた場合、パラメータa, bは分散共分散行列から計算することができます。



a = (yの平均) - σxy/σxx * (xの平均)
b = σxy/σxx

下記のような問題を考えた場合、
式1:y = a1 + b1 * x

式2:x = a2 + b2 * y
⇒ y = -a2/b2 + 1/b2 * x

回帰直線の変数xとyを入替えた場合の傾きは、一致しません。

つまり、因果の方向性を知っている場合に、
式1:(原因)= a + b * (結果)
式2:(結果)= a + b * (原因)

式2を考えた方が自然な発想になります。

ここで炉内温度、炉外温度の話に戻して、、、

炉内温度:Y [通常は、測定できない]
↓(熱伝搬)
炉外温度:X1 [測定可能]
↓(熱伝搬)
炉外温度:X2 [測定可能]

Yは炉内の真の温度で測定できない ⇒ 潜在変数
X1, X2は、観測変数となります。

因子分析や、パス解析、共分散構造分析を使ってモデリングをするのが次のステップとなります。

nice!(1)  コメント(0)  トラックバック(0) 
共通テーマ:学問

統計手法から統計モデリングへ その1 因果推論のモデリング [確率・統計、データマイニング]

統計手法から統計モデリングへ
https://www.jstage.jst.go.jp/article/rssj1981/18/1/18_1_57/_pdf

椿先生がよく講義で使われる元ネタ

こちらを自分なりに色々試して理解しました。

考え方はこんな感じ。

炉内温度:Y [通常は、測定できない]
↓(熱伝搬)
炉外温度:X1 [測定可能]
↓(熱伝搬)
炉外温度:X2 [測定可能]



このようなデータがあった時に、Yの温度をX1とX2を使って、どう推定するか?という問題になります。

そこで、
Y = a1 * X1 + a2 * X2 + a0
という重回帰分析を使って、モデル化しようというのが最初の第一歩。

結果は、
Y = 1.214 * X1 - 0.001 * X2 + 10.8
となるのですが、X2の係数は有意でなく、再度
Y = a1 * X1 + a0
というモデルを考え、パラメータを推定します。

その結果は、
Y = 1.213 * X1 + 10.8
となります。

そこで、Yの予測結果(炉内温度)は、上の式で与えられると考えるのですが、このモデルのおかしいところは、X1やX2が決まれば、Yが決まるという構造を表していることです。

本来は、Yが決まってX1やX2が決まるので、原因と結果が逆転していることになります。

nice!(1)  コメント(0)  トラックバック(0) 
共通テーマ:学問

IBM Watson Summit 2017 [確率・統計、データマイニング]

IBM Watson Summit 2017 に参加してきました。
実際に、自分でWatsonを呼び出すハンズオンセッションがあったりと、なかなか良かったです。

自分は、機械学習よりも統計モデル/ビジネス側の人間なので、そのあたりのセッションが減ってきているのは残念ではありますが、どちらか一方というよりかは使い分けな気がします。

片一方だけ知っていてはダメで、時間はかかるものの両方それなりに深いレベルまで理解しておくことが重要だと思っていましたし、今回のセッションを受けてやはりそうだなって思いました。

nice!(1)  コメント(0)  トラックバック(0) 
共通テーマ:学問

SPSS Modeler、シミュレーションノードを用いたモンテカルロ法の話 その2 [確率・統計、データマイニング]

SPSS Modeler、シミュレーションノードを用いたモンテカルロ法の話
http://skellington.blog.so-net.ne.jp/2017-04-21

前回は、正規分布を使ったシミュレーションでした。

今回は、ベータ分布を使ったシミュレーションの例になります。

ベータ分布は、ベルヌーイ分布や二項分布の事前分としての相性がとても良い分布です。
共役分布(きょうやくぶんぷ)と読みます。

3月の視聴率を30%と「点」で与えるのではなく、「30%付近に分布している」と分布で考えます。



シミュレーションノードを設定するとイメージが湧くかと思います。



例えば、形状1と2を1/10倍(サンプルサイズが1/10倍)にすると、分布は広がってきます。
つまり、30%よりも離れた値を取る割合が増えます。

形状の計算は、60人の30%ということで、18人、60人-18人=42人と計算します。
18+1=19
42+1=43



一方、形状1と2を10倍(サンプルサイズが10倍)にすると、分布の幅は狭くなります。

形状の計算は、6000人の30%ということで、1800人、6000人-1800人=4200人と計算します。
1800+1=1801
4200+1=4201



後は、「 theta_3 > theta_2 and theta_2 > theta_1 」となる確率を計算すれば、OKです。



朝野先生の教科書では、10万回シミュレーションされいて、0.57213という値でした。
SPSS Modelerを使ったシミュレーションでは、0.57299となり、非常に近い値を得ることができました。

nice!(4)  コメント(0)  トラックバック(0) 
共通テーマ:学問

SPSS Modeler、シミュレーションノードを用いたモンテカルロ法の話 [確率・統計、データマイニング]

朝野先生が書かれた『ベイズ統計学』をSPSS Modelerに実装しました。


ビジネスマンがはじめて学ぶ ベイズ統計学 ―ExcelからRへステップアップ―

ビジネスマンがはじめて学ぶ ベイズ統計学 ―ExcelからRへステップアップ―

  • 作者: 朝野 煕彦
  • 出版社/メーカー: 朝倉書店
  • 発売日: 2017/02/20
  • メディア: 単行本(ソフトカバー)



まずは、マルコフ連鎖モンテカルロ法ではない、普通のモンテカルロ法の話。

6章(6.1 モンテカルロ法より)
店舗A, B, C, Dの売上データあり、それを集計すると
 店舗A(平均100, 分散20)
 店舗B(平均120, 分散20)
 店舗C(平均110, 分散15)
 店舗D(平均100, 分散10)
となっていました。

なんとなく店舗Bの平均が高いのですが、1000日間の中で、店舗Bの売上が一番高い日はどれくらいあるでしょうか?といった問題をモンテカルロ法を使って計算します。

実際のストリームは、いたってシンプルです。



「シミュレーション生成」ノードの中身
分布で、正規分布を選んで平均と分散のパラメータを設定



売上なので、マイナスの値が出てくるとおかしいので、最小値を0.0としておきました。
実際は、平均100、分散10~20程度だと0以下になることは極めて稀なのですが、念のためです。w

相関のところは特に設定していませんが、こちらも設定することは可能です。
例えば、週末は売上が高くなるけど、平日は売上が落ちる、みたいな店舗の売上に相関があると仮定するならば、設定しても良いかもしれません。

拡張オプションでは、どれくらいの乱数を発生させるか、つまり、シミュレーションを行うか設定することができます。
デフォルトは10万。

「_max(フィールド作成)」ノードの中身
続いて、どの店舗の売上が最も高いかを計算しています。



@FIELD = max_n([Field_1, Field_2, Field_3, Field_4]


それぞれのフィールドで最大だったら1、それ以外は0という設定です。



「レコード集計とソート」ノードの中身
後は、各店舗別に1の数を合計しているだけなので、省略。

結果は、このようになります。



10万レコードあるので、各行のレコード数を10万で割ると、確率が出てきます。
教科書では、Bの店舗が最大になる確率は、54%となっていました。

今回、SPSS Modelerで実装した例では、54.75%ということで、ほぼ同じ結果となりました。

nice!(42)  コメント(0)  トラックバック(0) 
共通テーマ:学問

平成29年度統計数理研究所公開講座 [確率・統計、データマイニング]

平成29年度統計数理研究所公開講座が発表されました。

http://www.ism.ac.jp/lectures/kouza.html

B. ベイズ統計の理論・モデリング・評価について
D. 統計モデルと赤池情報量規準 AIC 1

気になっているのはこの2つ。

BもDも需要が高そうなので、抽選になりそうです。(^^;
抽選で当たると良いのですが。

nice!(4)  コメント(0)  トラックバック(0) 
共通テーマ:学問

All of Statistics: A Concise Course in Statistical Inference [確率・統計、データマイニング]

いくつか統計まわりの教科書を買いました。

All of Statistics: A Concise Course in Statistical Inference (Springer Texts
in Statistics)
https://www.amazon.co.jp/dp/0387402721/

All of Statistics: A Concise Course in Statistical Inference (Springer Texts in Statistics)

All of Statistics: A Concise Course in Statistical Inference (Springer Texts in Statistics)

  • 作者: Larry Wasserman
  • 出版社/メーカー: Springer
  • 発売日: 2004/10/21
  • メディア: ハードカバー




統計的学習の基礎 ―データマイニング・推論・予測―
https://www.amazon.co.jp/dp/432012362X/

統計的学習の基礎 ―データマイニング・推論・予測―

統計的学習の基礎 ―データマイニング・推論・予測―

  • 作者: Trevor Hastie
  • 出版社/メーカー: 共立出版
  • 発売日: 2014/06/25
  • メディア: 単行本




マーケティングの数理モデル (経営科学のニューフロンティア)
https://www.amazon.co.jp/dp/4254275161/

マーケティングの数理モデル (経営科学のニューフロンティア)

マーケティングの数理モデル (経営科学のニューフロンティア)

  • 作者:
  • 出版社/メーカー: 朝倉書店
  • 発売日: 2001/06
  • メディア: 単行本




実践 ベイズモデリング -解析技法と認知モデル
https://www.amazon.co.jp/dp/4254122209/

実践 ベイズモデリング -解析技法と認知モデル-

実践 ベイズモデリング -解析技法と認知モデル-

  • 作者: 豊田 秀樹
  • 出版社/メーカー: 朝倉書店
  • 発売日: 2017/01/25
  • メディア: 単行本(ソフトカバー)




岩波データサイエンス Vol.5
https://www.amazon.co.jp/dp/4000298550/

岩波データサイエンス Vol.5

岩波データサイエンス Vol.5

  • 作者:
  • 出版社/メーカー: 岩波書店
  • 発売日: 2017/02/16
  • メディア: 単行本(ソフトカバー)




ビジネスマンがはじめて学ぶ ベイズ統計学 ―ExcelからRへステップアップ―
https://www.amazon.co.jp/dp/4254122217/

ビジネスマンがはじめて学ぶ ベイズ統計学 ―ExcelからRへステップアップ―

ビジネスマンがはじめて学ぶ ベイズ統計学 ―ExcelからRへステップアップ―

  • 作者: 朝野 煕彦
  • 出版社/メーカー: 朝倉書店
  • 発売日: 2017/02/20
  • メディア: 単行本(ソフトカバー)



nice!(54)  コメント(0)  トラックバック(0) 
共通テーマ:学問

実践 ベイズモデリング [確率・統計、データマイニング]

豊田秀樹先生の新しい本が発売されていますね。

実践 ベイズモデリング -解析技法と認知モデル-2017/1/25


実践 ベイズモデリング -解析技法と認知モデル-

実践 ベイズモデリング -解析技法と認知モデル-

  • 作者: 豊田 秀樹
  • 出版社/メーカー: 朝倉書店
  • 発売日: 2017/01/25
  • メディア: 単行本(ソフトカバー)



こちらの基礎からのベイズ統計学の続編かと思われます。


基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

  • 作者: 豊田 秀樹
  • 出版社/メーカー: 朝倉書店
  • 発売日: 2015/06/25
  • メディア: 単行本



基礎とありますが、途中は歯ごたえがあります。

ただ、原理をきちんと押さえることなく、ツールを使ってなんとなくそれっぽい解析をするのはとても危険なことなので、色々な本を読んで基礎を固めたうえで、実践をしていきたいものです。

nice!(50)  コメント(0)  トラックバック(0) 
共通テーマ:学問

放送大学「心理統計法」 [確率・統計、データマイニング]

豊田秀樹先生が4月からの放送大学で「心理統計法」という授業があるようです。

http://www.ouj.ac.jp/hp/kamoku/H29/kyouyou/C/sinri/1529196.html

心理学はデータに基づいて心のメカニズムを研究する学問です。この目的のためのデータ分析法について講義します。従来の心理統計法の初年度の講義は、有意性検定の利用を前提としていました。しかし本講義には有意性検定が登場しません。ベイズ流のアプローチで学習系列が展開されます。もちろんt分布・F分布・カイ2乗分布は登場しません。その点で本講義はとてもユニークですから、はじめてデータ分析に入門する方ばかりでなく、長くデータ分析をしてきた方の統計学再入門のための授業としても利用していただけます。


2017年4月1日から毎週土曜日 14時30分~15時15分
全15回、各回45分の講義

p値を使わないってことで、こちらの本が参考になりそうです。

はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―2016/6/2


はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―

はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―

  • 作者: 豊田 秀樹
  • 出版社/メーカー: 朝倉書店
  • 発売日: 2016/06/02
  • メディア: 単行本(ソフトカバー)



nice!(5)  コメント(0)  トラックバック(0) 
共通テーマ:学問

ビッグデータ時代のマーケティング ベイジアンモデリングの活用 [確率・統計、データマイニング]

樋口先生と佐藤先生の本。

研究で使う手法について、この本に書かれているってことで、頂きました。
内容は難しそうですが、少しずつ読み進めていこうと思います。


ビッグデータ時代のマーケティング―ベイジアンモデリングの活用 (KS理工学専門書)

ビッグデータ時代のマーケティング―ベイジアンモデリングの活用 (KS理工学専門書)

  • 作者: 佐藤 忠彦
  • 出版社/メーカー: 講談社
  • 発売日: 2013/01/22
  • メディア: 単行本(ソフトカバー)



Kindle版




nice!(5)  コメント(0)  トラックバック(0) 
共通テーマ:学問

SPSS Modelerで逆関数法を使った乱数発生 [確率・統計、データマイニング]

RやSPSS Modelerには乱数生成の機能があります。
しかし、用意されていない乱数を発生したい場合、どう乱数を発生させるのか?

分布関数F(x)に従う乱数Yを生成する手順
1. (0, 1)の一様乱数Xを生成する

2. F-1(x)を求める(必要がある)。

ここで、三角分布の乱数(min=0, max=4, mode=1)を作成しました。

逆関数法を使って求めた三角分布の乱数





SPSS Modelerには乱数発生のシミュレーションノードがあるのですが、
調べたところ、なんと!三角分布の乱数発生がありました。w

直接、シミュレーションノードを使って求めた三角分布の乱数です。



逆関数法のシミュレーション結果と同じですね。

逆関数法を使うことで、ラプラス分布なども発生可能となります。
(ラプラス分布は、シミュレーションノードに入っていませんでした。)

nice!(54)  コメント(0)  トラックバック(0) 
共通テーマ:学問

正規分布かどうかの検定 [確率・統計、データマイニング]

色々な検定方法がありますが、ややこしいのは、微妙に似ている分布だと
同じデータであるにも関わらず、検定結果が採択になる場合があったり、
棄却される場合があったりしてしまいます。

適合度の検定--正規分布への適合度の検定
http://aoki2.si.gunma-u.ac.jp/lecture/GoodnessOfFitness/normaldist.html

上記以外の方法としては、Rのパッケージ tseries に含まれているJarque-Bera Test(ジャック・ベラ検定)などもあります。

library(tseries)
dat <- read.table("data.csv", header=T)

jarque.bera.test(dat)

## Jarque Bera Test
##
## data: data
## X-squared = 30.881, df = 2, p-value = 1.969e-07

上記の場合は、p値が 1.969e-07 なので、有意水準 5% だったとすると、
帰無仮説は棄却されることになります。

nice!(3)  コメント(0)  トラックバック(0) 
共通テーマ:学問

欠測データの統計学 [確率・統計、データマイニング]

統計数理研究所の公開講座『欠測データの統計学』に行ってきました。
http://www.ism.ac.jp/lectures/28m.html

高井先生、星野先生、野間先生によって書かれたこちらの本がベースになっているようです。


欠測データの統計科学――医学と社会科学への応用 (調査観察データ解析の実際 第1巻)

欠測データの統計科学――医学と社会科学への応用 (調査観察データ解析の実際 第1巻)

  • 作者:
  • 出版社/メーカー: 岩波書店
  • 発売日: 2016/04/20
  • メディア: 単行本(ソフトカバー)



今年もいろいろ公開講座を受講しましたが、今年度の受講はこちらが最後。
来年度の公開講座が楽しみです。

nice!(2)  コメント(0)  トラックバック(0) 
共通テーマ:学問

欠損値を含むレコードの処理 その6 オススメの補完方法 [確率・統計、データマイニング]

★ 過去の記事 ★
欠損値を含むレコードの処理 その1 欠測データの生成方法
http://skellington.blog.so-net.ne.jp/2016-12-19

欠損値を含むレコードの処理 その2 C&RT Treeを使った欠測処理
http://skellington.blog.so-net.ne.jp/2016-12-20

欠損値を含むレコードの処理 その3 リストワイズ削除
http://skellington.blog.so-net.ne.jp/2016-12-21

欠損値を含むレコードの処理 その4 平均値代入
http://skellington.blog.so-net.ne.jp/2016-12-22

欠損値を含むレコードの処理 その5 回帰代入
http://skellington.blog.so-net.ne.jp/2016-12-23


こちらの方法ではなんらかの不具合が生じてしまいました。

優れている手法としては、完全情報最尤推定(Full-Information Maximum Likelihood: FIML)や多重代入法(Multiple Imputation: MI)を使うのが良いようです。

多重代入法は、SPSS Modelerにはなく、SPSS Statisticsのオプションである Missing Values に入っているようです。

M.欠測データの統計科学:基礎理論と実践的な方法論
http://www.ism.ac.jp/lectures/28m.html


欠測データの統計科学――医学と社会科学への応用 (調査観察データ解析の実際 第1巻)

欠測データの統計科学――医学と社会科学への応用 (調査観察データ解析の実際 第1巻)

  • 作者:
  • 出版社/メーカー: 岩波書店
  • 発売日: 2016/04/20
  • メディア: 単行本(ソフトカバー)



nice!(1)  コメント(0)  トラックバック(0) 
共通テーマ:学問

欠損値を含むレコードの処理 その5 回帰代入 [確率・統計、データマイニング]

その2 C&RT Treeを使った欠測処理 で行った決定木を使う補完方法の場合、同じ枝の値は全部同じ値になってしまうため、上手く補完ができませんでした。

だったら、決定木ではなく回帰分析を行ってその値で補正しましょうというのは自然な流れです。



[真値]
yの平均:124
yの標準偏差:25.0
xとyの相関:0.592

[MCAR]
yの平均:124 ○
yの標準偏差:19.2 ×(小さい)
xとyの相関:0.793 ×(高い)



[MAR]
yの平均:122 ○
yの標準偏差:21.1 ×(小さい)
xとyの相関:0.816 ×(高い)



[MNAR]
yの平均:148 ×(高すぎる)
yの標準偏差:9.86 ×(小すぎる)
xとyの相関:0.671 ×(高い)



平均値に関して言えば、MCARとMARは大丈夫そうです。
yの標準偏差やxとyの相関に関しては上手く行っていません。

これは、CRT代入や平均値代入と同じく、一つの値で補正をしてしまっているため、母数にバイアスがないのが原因です。

★ 過去の記事 ★
欠損値を含むレコードの処理 その1 欠測データの生成方法
http://skellington.blog.so-net.ne.jp/2016-12-19

欠損値を含むレコードの処理 その2 C&RT Treeを使った欠測処理
http://skellington.blog.so-net.ne.jp/2016-12-20

欠損値を含むレコードの処理 その3 リストワイズ削除
http://skellington.blog.so-net.ne.jp/2016-12-21

欠損値を含むレコードの処理 その4 平均値代入
http://skellington.blog.so-net.ne.jp/2016-12-22

nice!(1)  コメント(0)  トラックバック(0) 
共通テーマ:学問

欠損値を含むレコードの処理 その4 平均値代入 [確率・統計、データマイニング]

リストワイズ削除は、データを削除するのでもったいない感じがします。
そこで、平均値で代入してみてはどうか?という発想です。

まずは、結果から。

[真値]
yの平均:124
yの標準偏差:25.0
xとyの相関:0.592

[MCAR]
yの平均:127 ○
yの標準偏差:14.5 ×(低い)
xとyの相関:0.322 ×(低い)

[MAR]
yの平均:141 ×(高い)
yの標準偏差:13.2 ×(低い)
xとyの相関:0.111 ×(低い)

[MNAR]
yの平均:152 ×(高い)
yの標準偏差:7.99 ×(低い)
xとyの相関:0.194 ×(低い)

なぜ、このような事が起こっているのか散布図を書いてみます。

[MCAR]の場合
平均値は真値に近いです。
これは、xとyが完全にランダムであるから。
しかし、yの値を一定の値にしてしまっているために
yの標準偏差が低くなったり、xとyの相関も低くなったりします。

CRT Treeを使った欠測処理で、CRTが分岐せずに全部同じ値になっていることと同様の結果ですね。



[MAR]の場合(x <= 135を欠損させる)
yの平均が高くなる理由ですが、xとyが相関があります。
欠測している個所は、欠損していない箇所よりも平均値が低くなっているはずです。
しかし、欠損していない個所の平均値で欠損箇所を埋めてしまっているために平均値が低くなってしまいます。



[MNAR]の場合(y <= 135を欠損させる)
MARよりもさらに平均値が高いです。
今回は、低いyの値を欠損させているにもかかわらず、高いyの値で補完していることが原因です。



★ 過去の記事 ★
欠損値を含むレコードの処理 その1 欠測データの生成方法
http://skellington.blog.so-net.ne.jp/2016-12-19

欠損値を含むレコードの処理 その2 C&RT Treeを使った欠測処理
http://skellington.blog.so-net.ne.jp/2016-12-20

欠損値を含むレコードの処理 その3 リストワイズ削除
http://skellington.blog.so-net.ne.jp/2016-12-21

nice!(64)  コメント(0)  トラックバック(0) 
共通テーマ:学問

欠損値を含むレコードの処理 その3 リストワイズ削除 [確率・統計、データマイニング]

リストワイズ削除とは、一つでも欠損があったら全部のレコードを削除するという方法になります。
なんだかもったいない気もします。。。

リストワイズ削除をすると、どのような影響になるのか、SPSS Modelerで検証します。



[真値]
yの平均:124
yの標準偏差:25.0
xとyの相関:0.592

[MCAR]
yの平均:126 ○
yの標準偏差:25.2 ○
xとyの相関:0.557 ○

[MAR]
yの平均:141 ×(高すぎる)
yの標準偏差:22.8 ×(低い)
xとyの相関:0.380 ×(低い)

[MNAR]
yの平均:152 ×(高すぎる)
yの標準偏差:14.2 ×(低すぎる)
xとyの相関:0.401 ×(低い)

xとyが完全ランダムの場合のみyの平均もyの標準偏差も同じになっていますが、
MARやMNARの場合は上手く復元できません。

特にxとyの相関係数は、切断効果により小さくなってしまいます。

リストワイズは手っ取り早いのですが、単純にデータを削除すると、平均や標準偏差など得られた結果がおかしい場合があるので注意が必要ですね。

★ 過去の記事 ★
欠損値を含むレコードの処理 その1 欠測データの生成方法
http://skellington.blog.so-net.ne.jp/2016-12-19

欠損値を含むレコードの処理 その2 C&RT Treeを使った欠測処理
http://skellington.blog.so-net.ne.jp/2016-12-20

nice!(7)  コメント(0)  トラックバック(0) 
共通テーマ:学問

欠損値を含むレコードの処理 その2 C&RT Treeを使った欠測処理 [確率・統計、データマイニング]

IBM SPSS Modelerには欠損値検査が付いています。
(手法としては、全部、微妙なのであまり使わない方が良さそうです。)

[データ検査]ノードというものがあります。
実行すると、このようなアウトプットが出てきます。



欠損値検査タブをクリックすると、どのようなアルゴリズムで欠損値を補完するかを選択することができます。

用意されているのは、下記の4種類。
1. 固定値
2. 無作為
3. 式(回帰式)
4. アルゴリズム(CRT)

1.~3.は、また、別の機会に調べるとして、まずは、4.のアルゴリズムから。
元々、データマイニングツールということもあり、補完するアルゴリズムは決定木のCRT(CA&RT)になっています。

なんで、CARTかといえば、classification and regression treesの名前の通り、
分類と回帰の両方を扱えます。

つまり、名義変数の補完の場合はclassification、連続値の補完としてregressionを使うことができるかです。

一見、上手く行きそうですが、細かく見ていくとおかしなことが起こっています。

今回、xとyは、相関0.600の関係になるようにサンプリングしています。
欠測を与える前のxとy値は0.592となっています。

欠測データを補完した際にxとyの相関も0.600(0.592)付近になることが期待されるのですが、0.315という結果になっていました。

出来上がったモデルを見ると、下図のようになっています。


つまり、欠損データ部分にすべて同じ値を埋めましょう、ということを意味しています。

次に、散布図を書いてみると、このようになります。


緑色の部分が欠損値のデータを補完した部分になります。
このようにすべて同じ値で埋めてしまうため、xとyの相関が0.315と低い値になってしまいました。

次にランダムシードを変更して実行すると、今度は、木が分岐し、散布図は下図のようになっています。


先ほど違ってすべて同じ値ではないのですが、やはり多くの個所を一つの値で置換していることがわかります。

また、別の問題として、C&RTを使って分岐するということは、
ちょっとした決定木の分岐の違いによって補完される値が変わってしまいます。
その結果、相関係数の値が実行するたびに大きく変化することも問題です。

ということで、一見、もっともらしいアルゴリズム(CRT)で置換しているように見えますが、細かく見ていくと、変なことが起こっているので、この置換方法はやめた方が良いように思えます。

★ 過去の記事 ★
欠損値を含むレコードの処理 その1 欠測データの生成方法
http://skellington.blog.so-net.ne.jp/2016-12-19

nice!(1)  コメント(0)  トラックバック(0) 
共通テーマ:学問

欠損値を含むレコードの処理 その1 欠測データの生成方法 [確率・統計、データマイニング]

レコードの中に欠損値があることは、よくあります。

IBM SPSS Modelerでも欠損値がある場合、それを補正するノードは用意されているのですが、今回、欠測データ解析の講義を受けて、なるほど!と思ったことがあったので、まとめたいと思います。

基本、SPSS Modelerでシミュレーションしていきます。

欠測データ解析と言えば、Rubin(1976)があまりにも有名です。

欠損値を3つのタイプに分けることができます。
今、(X, Y)の2変量データ、Yに欠測があるとします。

・MCAR(Missing Completely At Random)
YにもXにも依存しない。

・MAR(Missing At Random)
Yには依存しないが、Xに依存する。

・MNAR(Missing Not At Random)
Yに依存する。

いろいろな方法がありますが、上手く行くのはMCARとMARの場合です。

まずは、欠測データ解析用のデータを作ります。


[入力]タブにあるシミュレーション生成ノードを使います。
2変数(x, y)はともに平均が125, 標準偏差が25, 相関は0.6とします。

項目の選択
シミュレーションしたフィールド
ここで平均が125, 標準偏差が25を設定


相関
ここで相関0.6を設定


拡張オプション
発生させる乱数として1000個のデータを生成
ランダムシードを12345と設定


・MCAR(Missing Completely At Random)の作り方
randというフィールドを作成します。
random(100) と書くと、0 < x <= 100となります。
実際は、1から100までの100通りの乱数が生成されます。



random0(100) と書くと、0 <= x <= 100となります。
実際は、0から100までの101通りの乱数が生成されます。
通常は、randomを使うことが多いかと思われます。

YにもXにも依存しない。
つまり、randの値が66以下の場合にデータを欠損させます。


・MAR(Missing At Random)の作り方
Yには依存しないが、Xに依存する。
つまり、x <= 135 の場合に、yの値を欠損させます。


・MNAR(Missing Not At Random)
Yに依存する。
つまり、y <= 135 の場合に、yの値を欠損させます。


以下、欠測データを補完した場合にどうなるかをシミュレーションしていきます。

ちなみに、このデータは、相関は0.6で作っているので、y = a x + b という単回帰を考えた場合、回帰係数は0.6となるはずです。

完全データで線形回帰を行った場合、
y = 0.5962 * x + 50.83
となりました。



確かに、回帰係数は 0.6 となっています。
欠測データじゃないので、当然の結果と言えば当然の結果ですが。。。(^^;

欠測データがある場合、
・そのレコードを削除する
・平均値や最頻値などの固定値で置換する
・回帰などのアルゴリズムを使う
が考えられます。

SPSS Staticsでは、Missing Valuesを購入すれば、EMアルゴリズムで推定したり
多重代入を使って置換することができるようです。

SPSS Modelerに搭載されているアルゴリズムは、C&RT Treeが使えます。

nice!(50)  コメント(0)  トラックバック(0) 
共通テーマ:学問

統計・機械学習における確率的最適化 [確率・統計、データマイニング]

統計数理研究所の公開講座『統計・機械学習における確率的最適化』に行ってきました。
上期に行われた『スパース推定』の続編的な講座です。

ただ、かぶっている部分も多くあって、それぞれの先生によっていろんな角度から説明されるので、非常に参考になりました。

下期は、あと『欠測データの統計科学:基礎理論と実践的な方法論』を受講する予定です。

nice!(1)  コメント(0)  トラックバック(0) 
共通テーマ:学問

コーシー分布:乱数を使ったシミュレーション [確率・統計、データマイニング]

Rを使って正規分布とコーシー分布の乱数を1,000個発生させます。

# 正規分布の乱数を発生
x <- rnorm(1000)

# コーシー分布の乱数を発生
y <- rcauchy(1000, location = 0, scale = 1)

グラフを書いてみるとコーシー分布が正規分布に比べて、裾が重い分布ということが良くわかります。

まずは、正規分布のグラフ(青)。


これにコーシー分布(赤)を重ねてみると、、、

確かに0を中心とした乱数が発生されていますが、正規分布と比較してより大きな値もたくさん出ています。

さらに分かりやすいのが、グラフのy軸をデータに合わせてみると・・・

時々、ものすっごく大きな値が出てきています。

このグラフは正規分布のグラフと一緒に書いているのですが、正規分布の方は0付近でおとなしく値が出ています。

グラフを書いてみると、コーシー分布は少しやんちゃな分布のように見えます。

nice!(3)  コメント(0)  トラックバック(0) 
共通テーマ:学問
前の30件 | - 確率・統計、データマイニング ブログトップ