Hatena::ブログ(Diary)

iAnalysis 〜おとうさんの解析日記〜

2009 | 11 |
2010 | 02 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2011 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2012 | 01 | 03 | 04 | 05 | 10 | 12 |
2013 | 01 | 02 | 04 |
2014 | 03 | 12 |
2016 | 03 |

2016-03-26

異常検知について

この本を読みました。

いろんな場合でどんな手法があるかが書いてあるので、とても整理されています。

まとめると、

  • 1次元で山がひとつの異常検知:ガンマ分布で異常度を計算する
  • 多次元または山が複数ある異常検知:クラスタリング
  • 不要な次元を含むデータからの異常検知:主成分分析による異常検知
  • 目的変数のある異常検知:リッジ回帰
  • 周期のわかる時系列データの異常検知:自己回帰モデル
  • 周期のわからない時系列データの異常検知:部分状態空間モデル

って感じです。

多くの場合は、クラスタリングと主成分分析異常検知でまかなえます。

時系列異常検知のできる部分状態空間モデルは魅力的ですが、本によるとRのパッケージは充実してないみたい。要調査。

2016-03-19

ggplot2でロジスティック回帰の近似曲線を描く

バージョンアップで昔の指定方法ができなくなってた。

調べたら下記のようにやれば描ける。ggplot2は綺麗なんだけど仕様変更が多いから困る。

binomial_smooth <- function(...) {
  geom_smooth(method = "glm", 
              method.args = list(family = "binomial"), ...)
}

ggplot(data, aes(x, y)) + geom_point() + binomial_smooth()

2016-03-12

メモ:ROC曲線の最適カットオフを計算する関数

epiライブラリROC関数ではカットオフまで計算してくれなかったので、ROC関数を計算するように書き換え。

library(epi)

ROC1 <- function (test = NULL, stat = NULL, form = NULL, plot = c("sp", 
    "ROC"), PS = is.null(test), PV = TRUE, MX = TRUE, MI = TRUE, 
    AUC = TRUE, grid = seq(0, 100, 10), col.grid = gray(0.9), 
    cuts = NULL, lwd = 2, data = parent.frame(), ...) 
{
    rnam <- if (!missing(test)) 
        deparse(substitute(test))
    else "lr.eta"
    if (is.null(form)) {
        if (is.null(stat) | is.null(test)) 
            stop("Either 'test' AND 'stat' OR 'formula' must be supplied!")
        lr <- glm(stat ~ test, family = binomial)
        resp <- stat
        Model.inf <- paste("Model: ", deparse(substitute(stat)), 
            "~", deparse(substitute(test)))
    }
    else {
        lr <- glm(form, family = binomial, data = data)
        resp <- eval(parse(text = deparse(form2)), envir = lr$model)
        Model.inf <- paste("Model: ", paste(paste(form)[c(2, 
            1, 3)], collapse = " "))
    }
    m <- as.matrix(base::table(switch(PS + 1, test, lr$fit), 
        resp))
    m <- addmargins(rbind(0, m), 2)
    fv <- c(-Inf, sort(unique(switch(PS + 1, test, lr$fit))))
    nr <- nrow(m)
    m <- apply(m, 2, cumsum)
    sns <- (m[nr, 2] - m[, 2])/m[nr, 2]
    spc <- m[, 1]/m[nr, 1]
    pvp <- m[, 2]/m[, 3]
    pvn <- (m[nr, 1] - m[, 1])/(m[nr, 3] - m[, 3])
    res <- data.frame(cbind(sns, spc, pvp, pvn, fv))
    ddaattaa <- data.frame(lr[21], lr$fit)
    names(res) <- c("sens", "spec", "pvp", "pvn", rnam)
    auc <- sum((res[-1, "sens"] + res[-nr, "sens"])/2 * abs(diff(1 - 
        res[, "spec"])))
    
    mx <- max(res[, 1] + res[, 2])
    mhv <- which((res[, 1] + res[, 2]) == mx)
    mxf <- fv[mhv]
    
    cutoff <- ddaattaa[ddaattaa$lr.fit == mxf, ]
    
    invisible(list(res = res, AUC = auc, lr = lr, cutoff = cutoff, optimal = res[mhv, ], 
                   n = nrow(ddaattaa), freq = table(ddaattaa[, 1])))
}

b <- ROC(form = a[, 1] ~ a[, 2], plot="ROC")
c <- ROC1(form = a[, 1] ~ a[, 2], plot="ROC")
c$cutoff
c$AUC
c$optimal
c$n
c$freq

2014-12-12

【厳選】ビジネスを成功に導くための14のデータ分析メソッド日経ムック

日経ビッグデータから出版されました『実践!”超”分析の教科書』の監修と執筆をしました。この本は、「分析はスゴいって騒がれてるけど、自分には難しそう。でも興味はある」と思っているような初心者のために、分析って何なのか分かりやすく伝えたいという思いで、監修・執筆しました。


自分が担当したのはケーススタディキーワード解説の部分です。分量でいうと本の約半分くらいです。あと、最初のインタビューにも出ています。

ケーススタディはまず、私がiAnalysis(アイアナリシス)社で3年間コンサルをしてきて、ビジネスに効果が出た分析方法だけをピックアップしました。毎日ビッグデータとかアナリティクスとか騒がれてますが、「結局どんなことをするの?」的なことを疑問に思われてる方が多いかと思います。そういう方に、「データを活用するとこんなことができますよ」と説明しています。それも、ビジネスに役立つものだけを厳選しているので、これだけ押さえていれば、まずは充分だと思います。いろんなところで報道されている事例は、だいたいこれらのどれかに入るんじゃないでしょうか。

  1. モニタリング
  2. KPI厳選
  3. マーケティング調査
  4. 顧客分析
  5. テキストマイニング
  6. 仮説検証
  7. 広告効果測定
  8. 施策効果検証試験
  9. 反応スコアリング
  10. ライフタイム分析
  11. 需要予測
  12. 故障予測
  13. 個別最適化
  14. レコメンデーション

そしてその手法にあった事例を、たくさんの企業の中から選ばせてもらいました。分析の初心者にも興味を持ってもらうため、誰でも知っている大企業と、業界で注目を集めている企業を織りまぜて選びました。実に多くの企業が、データを活用して事業を拡大しています。

  1. イーグルバス
  2. セリア
  3. りそな銀行
  4. ローソン
  5. 明治安田生命/富国生命保険
  6. セブン&アイ・ホールディングス
  7. DeNA
  8. 日本エスリード
  9. 日本航空
  10. オルビス
  11. 大林組
  12. 大阪ガス
  13. アサヒビール
  14. 良印計画

これからより一層データ活用を促進していこうと思われている方には、とても参考になる内容だと思います。お付き合いのある大手システム開発会社様では、数十冊購入し、社内教育に利用するとおっしゃっていました。


また、iAnalysis社では”オリオン”という、ビジネスツールにアナリティクス機能を付加させたウェブサービスを開発していますオリオンは、日常業務のビジネスツールとして利用しつつ、高度なアナリティクスを簡単に業務に活かすことを目指したサービスです。あらゆる現場で使ってもらうために、スマホタブレット・パソコンのどれからでも利用できます。現在(2014/12/12)、オリオンはクローズドα版なので、URLは一般公開していませんが、”超”分析の教科書のP37に限定でデモ版を公開しています。ユーザー登録できますし、登録したアカウントは正式サービス後も使えますので、興味がある方は教科書を買って登録してみて下さい。


オリオンの一部機能は、日経ビッグデータ様で取り上げて頂きました。

「有力サイエンティストが独自開発ツールで事業拡大、iAnalysisがIoTベンダーと連携」

オリオン自体は、あらゆる業種業態のビジネスシーンで活用することを考えています。その中で、IoTやM2Mで発生するデータも扱えるようにしていきたいと考えています。


オリオンの詳細は、URLも含めて年明けあたりに公開しようと思っています。そのときに、オリオンの機能も順次紹介していこうと思います。

2014-03-12

マレーシア航空旅客機の衛星写真を画像解析

マレーシア航空のいろんな話題が飛び交っている中、Twitterで「【助けて】全世界のインターネットユーザーに協力を呼びかけ! この写真から「消息を絶ったマレーシア航空の旅客機」を見つけてください」という記事を見つけました。


衛星写真の画像が大量にあるので、人海戦術で破片を見つけよう、という趣旨です。私は分析屋なので、データサイエンスを使って手助けできないか?と思い、少し分析してみました。何かの一助になればと思い、ブログで公開します。


データサイエンスを利用するには、まずデータが必要です。衛星画像はこのサイトにあります。


そのサイトから、まずは特徴的な画像を拾ってきて、分析してみます。サイトを見てみると、画像はだいたい3つのパターンに分かれているようでした。

  1. 雲の画像
  2. 海の画像
  3. 何かの物体の画像

3つ目の”何かの物体”を画像から判定できれば、それが旅客機かもしれません。分析のロジックを作るために、ある船が含まれる画像を使って分析してみました。その画像が次のものです。サイト上ではこのページ(23975番)の画像です。

f:id:isseing333:20140312211139p:image


衛星画像をいくつか適当に眺めて、たまたま見つけたものです。まずは特徴的なこの画像を使って、分析ロジックを検討します。船とヘリポートが写っており、これを画像解析によって見つけられるか試します(画像データのダウンロードはできなかったので、Web画面のスクリーンショットを保存しました)。


分析には、いつものように統計解析ソフトであるRを利用します。Rでの画像処理に関しては、@wdkzさんの「Rでウォーリーを探してみた」の記事が大変参考になりました。この場を借りてお礼申し上げます。


まず、飛行機が写っているとしたら画像では白っぽくなるだろうと思い、画像から”白っぽい物体”を抽出するロジックを考えました。いろいろ試行錯誤しましたが、結果的に、次のような順で行うことで”白っぽい物体”を抽出できました。


1.白っぽい部分以外を1色に塗りつぶす

f:id:isseing333:20140312211140p:image


2.白っぽい部分も1色にする

f:id:isseing333:20140312212139p:image


3.色を少しぼやけさせる

f:id:isseing333:20140312212316p:image


4.ぼやけた部分をくっきりさせる

f:id:isseing333:20140312212356p:image


5.ノイズを除去する

f:id:isseing333:20140312212427p:image



最終的に出来た画像が、5の処理で行ったものです。赤く協調されているものが”白っぽいもの”です。もともとあった船とヘリポートがしっかりと認識されています。ただ、他にもいくつか赤い点があります。もとの画像をみると、恐らく海の波か何かなんじゃないかと思いますが、ちょっと分かりません。ただ、1つ目の画像処理の時点で、何か浮き出てますね。。何だろ、、、もしかしたら目的のものの一部?


これで”白っぽいものを認識する”ロジックはできたはずなので、確認のために他の画像にも当てはめてみます。先ほどの海域の左にある、雲の画像に当てはめてみたものが、次の結果です。


もとの画像

f:id:isseing333:20140312212754p:image


解析後の画像

f:id:isseing333:20140312212819p:image



2枚目が解析後の画像ですが、雲が1つの物体と認識されてないので明らかにおかしいです。多分、判定ロジックのどこかがおかしく、まだ改良しなくてはならなさそうです。これがうまくいったら次は、「雲かどうか」も何かしらのロジックで判別しなくてはならないですね。


とりあえず今思い付く方法はこれくらいです。これらをまとめると、次の3点です。

  1. 画像の白色判定はうまくいったので、飛行機を見つけるときの補助ができる(人が目で見る分には、処理1か2だけで良いかもしれない)
  2. 【課題】雲を見分ける必要がある(雲が無い衛星写真があればベスト)
  3. 【課題】サイトから画像が全部ダウンロードできないので、自動で取得してくる必要がある

今できるのはこれくらいですが、改良すればもっと実用的になるかもしれません。何かコメントがあればTwitterの方でどうぞ。全画像が手に入れば、ビッグデータ解析になります。

なるべく早く事案が解決することを願っています。



以下、Rのコードです。ご参考下さい。

# ImageMagickなどのインストールImageMagickは結局使ってない)
library(installr)
install.ImageMagick()

install.packages(c("abline", "biOps"))
source("http://www.bioconductor.org/biocLite.R")
biocLite("EBImage")
source("http://rimagebook.googlecode.com/svn/installRImageBook.R")
installRImageBook()


require(EBImage)
require(biOps)
require(RImageBook)

# setwd直下に、画像番号をファイル名にしたキャプチャpng画像を置いておく
# setwd直下に画像の番号を付けたフォルダを作る必要がある

WhiteExtract <- function(number){
  MalayAir <- readImage(file=paste("MalayAir", number, ".png", sep=""))
  # display(MalayAir)
  # gtkImageMagickの設定が難しい。writeImageで代用する。
  
  
  # 白色抽出
  writeImage(MalayAir, paste(number, "/MalayAirTest0.png", sep=""))
  
  MalayAir.white <- MalayAir
  MalayAir.white[MalayAir.white > 0.1] <- 1
  writeImage(MalayAir.white, paste(number, "/MalayAirTest1.png", sep=""))
  
  MalayAir.white2 <- thresh(MalayAir.white, w=50, h=50, offset=0.1)
  writeImage(MalayAir.white2, paste(number, "/MalayAirTest2.png", sep=""))
  
  MalayAir.white3 <- gblur(MalayAir.white2, sigma=1)
  writeImage(MalayAir.white3, paste(number, "/MalayAirTest3.png", sep=""))
  
  MalayAir.white4 <- (MalayAir.white3 > 0.1)
  writeImage(MalayAir.white4, paste(number, "/MalayAirTest4.png", sep=""))
  
  kern <- makeBrush(7, shape="diamond")
  MalayAir.white5 <- closing(opening(MalayAir.white4, kern), kern)
  writeImage(MalayAir.white5, paste(number, "/MalayAirTest5.png", sep=""))
  writeImage(MalayAir.white5, paste("MalayAir", number, "_判定後.png", sep=""))
  
  
  MalayAir.white5.label <- bwlabel(MalayAir.white5)
  print(max(MalayAir.white5.label))
}

WhiteExtract(21387)
WhiteExtract(24889)