Rで散布図とヒストグラム

irisのデータから散布図とヒストグラムを作る
前回も使ったけどirisのデータセットって？

irisの正体を調べた方がいた
irisの正体 (R Advent Calendar 2012 6日目) - どんな鳥も
"アヤメ属の3種 (Iris setosa，Iris versicolor，Iris virginica) の"、
"sepal (がく片) とpetal (花弁) の長さ(length)・幅(width)"の測定結果
ということだそう。

前回は種ごとに長さと幅の平均値を計算してグラフにしたが
今回は散布図とヒストグラムを書いて分布を調べる

結果はこんな感じ

setosa(緑)は他の2種にくらべてがく片の幅が大きいが
それ以外は小さく分布も離れている
versicolor(赤)とvirginica(青)はがく片のサイズは似ているが
花弁のサイズで区別がつきそう
virginicaは外れ値なのかもう1グループありそうにも見える

コードはこちら
MakingSDGraph()からGraphSD()とGraphHist()を呼び出している
irisのデータセットの読み込み方は前回と同じ

Data.Iris <- iris
# setosa, versicolor, virginica
Species.List <- as.character(unique(Data.Iris$Species))
# Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, Species
Col.Data <- colnames(Data.Iris)
# グラフの作成
MakingSDGraph(Data.Iris, Species.List, Col.Data[5], Col.Data)

MakingSDGraph()はpngファイルの書き込む用意をして
GraphSD()とGraphHsit()に渡すデータセットを整形する

MakingSDGraph <- function(data.all, sp.list, col.name, col.list){

  # データごとの色の定義
  data.color <- data.all[col.name]

  # グラフを書き込むファイルの用意, 5x5の領域を確保する
  png.name = sprintf("SD-Hist.png")
  png(png.name, width = 1200, height = 1200)
  par(mfcol = c(4, 4), lwd = 2, mar = c(5, 5, 1, 1), ps = 12)

  n <- c(); m <- c()
  for (i in 1:4){
    for (j in 1:4){
      par1 = as.character(col.list[i]); par2 = as.character(col.list[j])

      # 各列を変数に代入
      xT <- data.all[[par1]]; yT <- data.all[[par2]]
      xyT <- as.data.frame(cbind(xT, yT, data.color))

      # 最小値、最大値を求める
      minimum <- c(min(xT), min(yT)); maximum <- c(max(xT), max(yT))
      minimum.sign <- sign(minimum); maximum.sign <- sign(maximum)
      minimum <- floor(minimum * minimum.sign) * minimum.sign - 1
      maximum <- ceiling(maximum * maximum.sign) * maximum.sign + 1

      if(par1 != par2){
        GraphSD(par1, par2, minimum, maximum, xyT, sp.list)
      }
      else{
        GraphHist(par1, minimum, maximum, xyT, sp.list)
      }
    }
  }
  dev.off()
}

GraphSD()で散布図を作成する

GraphSD <- function(p1, p2, min, max, d, s){
    plot(0, 0, type = "n",
    xlim = c(min[1], max[1]),
    ylim = c(min[2], max[2]),
    xlab = p1,
    ylab = p2,
    cex.lab = 2,
    cex.axis = 2
  )
  # 点を描く
  d[,3] <- gsub(s[1], "#9bbb5980", d[,3])
  d[,3] <- gsub(s[2], "#c0504d80", d[,3])
  d[,3] <- gsub(s[3], "#4f81bd80", d[,3])
  points(d[,1], d[,2], col = d[,3], pch = 16, cex = 2.5)
}

GraphHist()でヒストグラムを作成する
ビン幅はkで定義しているがちょうどよいところは何度か試す
スタージェスの公式で求めた値そのままよりも3分の1にしたほうが良さげ

GraphHist <- function(p1, min, max, x, s){
  Lower <- min[1] - 1 ; Higher <- max[1] + 1
  k <- seq(Lower, Higher, (Higher - Lower) / (1 + log2(nrow(x))) / 3)

  # リガンドごとにデータを分ける
  sub_xr <- subset(x, (x[,3] == s[2]))
  sub_xb <- subset(x, (x[,3] == s[3]))
  sub_xg <- subset(x, (x[,3] == s[1]))

  # 分布の最大値を求める
  hr <- hist(as.numeric(sub_xr[,1]), breaks = k, plot = FALSE)
  hb <- hist(as.numeric(sub_xb[,1]), breaks = k, plot = FALSE)
  hg <- hist(as.numeric(sub_xg[,1]), breaks = k, plot = FALSE)
  max.hist <- max(hr$count, hb$count, hg$count)

  # ヒストグラム赤
  hist(as.numeric(sub_xr[,1]), col = "#c0504d60",
    breaks = k, axes = F,
    xlim = c(Lower, Higher),
    ylim = c(0, max.hist),
    main = "", xlab = p1,
    cex.lab = 2,
  )

  # ヒストグラム青
  hist(as.numeric(sub_xb[,1]), col = "#4f81bd60",
    breaks = k, axes = F,
    xlim = c(Lower, Higher),
    main = "", xlab = p1,
    add = TRUE
  )

  #  ヒストグラム緑
  hist(as.numeric(sub_xg[,1]), col = "#9bbb5960",
    breaks = k, axes = F,
    xlim = c(Lower, Higher),
    main = "", xlab = p1,
    add = TRUE
  )

  Axis(side = 2, cex.axis = 2)
  Axis(side = 1, at = (Lower + 1):(Higher - 1),
    cex.axis = 2,
    lty = NULL)
}

追加でライブラリを読み込ませずにグラフを作っているけど
ggplotだとどうだろう

12 March 2019追記
こういう書き方もある
散布図行列を描くには （corrplot, pairs, GGally） - StatModeling Memorandum
追記ここまで

23 June 2020 追記
python + pandas + seabornだとこんな感じ
【python】iris(アヤメ)のデータセットをpandasとseabornを使って可視化する
追記ここまで

Rで棒グラフにエラーバーをつける

何かしらデータを取ったらグラフ作ってエラーバーを加えて
統計的有意差があるかを確認をするのだが
データ量が増えるとエクセルだとめんどい

Rで棒グラフを書いてエラーバーを作れるようにした

irisのデータでグラフを書く
結果はこんな感じ

元のデータこれ

  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa
~省略~

Rのコンソールで
> iris
で確認できる
アヤメ属に含まれる3種の花弁とがく片の測定データで
setosa, versicolor, virginicaそれぞれ50点ずつデータがある

まずデータを読み込み、種の名前を取り出す
種名はsetosa, versicolor, virginicaの3種
次にデータの項目
Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, Species
を取り出してグラフを書く関数に渡す

Data.Iris <- iris
Species.List <- as.character(unique(Data.Iris$Species))
Col.Data <- colnames(Data.Iris)
MakingBarGraph(Data.Iris, Species.List, Col.Data[5], Species.List[1])

関数MakingBarGraph()は引数に
irisのデータ、
種名のリスト、
種名の入ったカラムの名前、
統計的な比較をする種名、
を使う

今回はsetosaと他の2種の測定結果とに
統計的有意差があればスターをつける

エラーバー、スターの書き方は以下のサイトを参考にした
棒グラフ | カテゴリーで分けられるデータの視覚化（棒グラフおよびエラーバー）
出力はpng形式でファイル名はbargraph.png

グラフを書く関数がこちら

MakingBarGraph <- function(data.all, sp.list, col.name, sp.name){
  tmp <- subset(data.all, data.all[[col.name]] == sp.name)
  tmp <- tmp[, colnames(data.all) != col.name]

  data.mean <- c()
  data.error <- c()

  for(i in sp.list){
    tmp <- subset(data.all, data.all[col.name] == i)
    tmp <- tmp[, colnames(data.all) != col.name]
    # 統計量を求める、列方向に計算
    # StatCalは別に定義した関数
    data.all.stat <- (apply(tmp, 2, StatCal))
    data.all.stat <- data.frame(t(data.all.stat))

    data.mean <- rbind(data.mean, data.all.stat$Mean)
    data.error <- rbind(data.error, data.all.stat$SEM)
  }
  colnames(data.mean) <- colnames(tmp)

  rownames(data.mean) <- sp.list
  colnames(data.error) <- paste(colnames(tmp), "error", sep = ".")

  print(data.mean)
  print(data.error)

  png.name = sprintf("bargraph.png")
  png(png.name, width = 800, height = 600)

  bar.mean <- barplot(data.mean,
    col=c("#9bbb5960", "#c0504d60","#4f81bd60"), # bar color
    #main = File.Name.First, # graph title
    ylim = c(0, max(data.mean + data.error) + 1),
    beside = T
  )

  legend("topleft", legend = sp.list,
    fill =c("#9bbb5960", "#c0504d60","#4f81bd60"),
    col = c("#9bbb5960", "#c0504d60","#4f81bd60"),
    border = c("#9bbb5960", "#c0504d60","#4f81bd60"),
    box.lwd = 0, box.lty = 0
  )

  arrows(bar.mean, data.mean - data.error, bar.mean, data.mean + data.error,
    code = 3, lwd = 1, angle = 90, length = 0.1
  )

  # 有意差の検定
  T12 <- c()
  T13 <- c()
  for(i in 1:4){
    t1 <- data.all[data.all[[col.name]] == sp.list[1],]
    t2 <- data.all[data.all[[col.name]] == sp.list[2],]
    t3 <- data.all[data.all[[col.name]] == sp.list[3],]

    # まずはF検定 var.test()
    if(var.test(t1[,i], t2[,i])$p.value < 0.05){
      F12 = F
    }else{
      F12 = T
    }
    if(var.test(t1[,i], t3[,i])$p.value < 0.05){
      F13 = F
    }else{
      F13 = T
    }

    # T検定, t.test()
    T12 <- append(T12, t.test(t1[,i], t2[,i], var.equal = F12)$p.value)
    T13 <- append(T13, t.test(t1[,i], t3[,i], var.equal = F13)$p.value)
  }

  # **, *をグラフにつける

  for (i in 1:length(T12)) {
    xi <- bar.mean[2, i] # x axis
    yi <- data.mean[2, i] + data.error[2, i] + 0.5
    if (T12[i] < 0.05) {
      if (T12[i] < 0.01) {
        text(xi, yi, "**")
      } else if (T12[i] < 0.05) {
        text(xi, yi, "*")
      }
    }
  }
  # **, *をグラフにつける
  for (i in 1:length(T13)) {
    xi <- bar.mean[3, i] # x axis
    yi <- data.mean[3, i] + data.error[3, i] + 0.5
    if (T13[i] < 0.05) {
      if (T13[i] < 0.01) {
        text(xi, yi, "**")
      } else if (T13[i] < 0.05) {
        text(xi, yi, "*")
      }
    }
  }
  dev.off()
}

基本統計量はsummary()で求められるが
SEMは出てこないので基本統計量を求める関数をつかう
この回答を参考にした
How can I summarizing data statistics using R - Stack Overflow

StatCal <- function(x, Na.Remove=TRUE){
  result <- c(Mean = mean(x, na.rm=Na.Remove),
              SD = sd(x, na.rm=Na.Remove),
              Median = median(x, na.rm=Na.Remove),
              Min = min(x, na.rm=Na.Remove),
              Max = max(x, na.rm=Na.Remove),
              N = length(x),
              SEM = sd(x, na.rm=Na.Remove)/sqrt(length(x)))
}

もし規格化するなら次の関数を使う
呼び出し方はこんな感じ
Data.Iris <- NormData(Data.Iris, Col.Data[5], Species.List[1])
この場合はsetosaの平均値で規格化する
19 Septembe 2020 追記
平均値で割っているだけなので規格化とは異なるかも
normalizationと呼んでいる場合もあるようだけど
Normalizations: dividing by mean - Cross Validated
How, When, and Why Should You Normalize / Standardize / Rescale… – Towards AI

追記ここまで

NormData <- function(data.all, col.name, sp.name){
  # 列名の一時保存
  tmp.col <- data.all[, colnames(data.all) == col.name]
  tmp.col.name <- colnames(data.all)
  # 基準になるデータセットの取得と平均値の計算
  tmp <- subset(data.all, data.all[[col.name]] == sp.name)
  data.all.stat <- apply(tmp[, colnames(tmp) != col.name], 2, StatCal)
  data.all.stat <- data.frame(t(data.all.stat))
  data.all.mean <- data.all.stat$Mean
  # 規格化
  data.all <- t(apply(data.all[, colnames(data.all) != col.name], 1, function(x){x/data.all.mean}))
  data.all <- data.frame(data.all)
  data.all <- cbind(data.all, tmp.col)
  colnames(data.all) <- tmp.col.name

  return(data.all)
}

結果はこうなる

ggplotとか使えばもっと楽できれいにできるのかもだけど

次はRで

前回pythonで書いたことを今度はRに移植

p <- c("PEN", "APPLE", "PINEAPPLE")
word <- ""

while(word != "PEN, PINEAPPLE, APPLE, PEN"){
  n <- as.integer(runif(4, min = 1, max = 4))

  word <- c(p[n[1]])
  for(j in 2:4){
    word <- paste(word, p[n[j]], sep = ", ")
  }

  print(word)
}

print("PPAP!")

ベクトル間の比較でTRUE/FALSEは返すが
その数を調べる方法がわからなかった
乱数で単語を選んで文字列に結合
そして正解と比較した
結果はこちら

[1] "PINEAPPLE, APPLE, APPLE, APPLE"
[1] "APPLE, PEN, PEN, PEN"
[1] "APPLE, PINEAPPLE, APPLE, PEN"
~中略~
[1] "PINEAPPLE, APPLE, PEN, APPLE"
[1] "PEN, PEN, APPLE, PINEAPPLE"
[1] "PEN, PINEAPPLE, APPLE, PEN"
[1] "PPAP!"

やっぱりRで書いている人もいるが
1行で書いてたりする
Rでズンドコ - Qiita

python手始め

学生の時触っていたのはCだった
pythonは初めてだったので、教科書を探してみた　

pythonを勉強するにあたっての教科書はこれを使った
第3版を使っていたけど、第4版が出ましたね
O'Reilly Japan - Pythonチュートリアル 第4版
(O'Reillyのサイトに飛びます)

オライリーからいろいろ本が出てるがとりあえず薄いのでプレッシャーにならなかった
細かくは書かれてないのであとはネットで調べる

このときの環境はiPadで動くPythonista3
‎「Pythonista 3」をApp Storeで

結構いろいろできる
iOS上で動作する革命的ものづくり環境「Pythonista 3」の魅力をとくと語る

11 March 2019追記
顔認識もできたのは感動した
コードもそんなに長くない
togetterまとめ
Pythonista(iOSで動くPython)で顔認識 - Togetter
ソースコード
FaceDetectTest2.py · GitHub

現在はwindows10 + ATOM + scriptを使用中

で、一通りやってみたら自分でなにか書いてみたくなる
こういうのがあるわけですが
ズンドコキヨシまとめ - Qiita
ズンドコキヨシ with python / ruby / Lua - Qiita
pythonでズンドコキヨシ - Qiita
うまく流行りに乗れないというか、天の邪鬼の自分はこんな感じで

import random

def choicewords(words):
	# from random import choice
	# choice(['Pen', 'Apple', 'Pineapple'])
	choose = []
	for i in range(4):
		w = random.choice(words)
		choose.append(w)
	return choose

def ppap():
	words = ['PEN', 'APPLE', 'PINEAPPLE']
	choose = choicewords(words)
	while choose != ['PEN', 'PINEAPPLE', 'APPLE', 'PEN']:
		choose = choicewords(words)
		print(choose)		

	print('PPAP!')

if __name__ == "__main__":
	ppap()

実行結果は

['APPLE', 'PEN', 'PEN', 'PINEAPPLE']
['PEN', 'PEN', 'APPLE', 'PEN']
['PINEAPPLE', 'PINEAPPLE', 'APPLE', 'PINEAPPLE']
~中略, 83行~
['PINEAPPLE', 'APPLE', 'APPLE', 'PINEAPPLE']
['APPLE', 'PINEAPPLE', 'APPLE', 'APPLE']
['PEN', 'PINEAPPLE', 'APPLE', 'PEN']
PPAP!

わざわざランダムに言葉を選ぶ関数を書いたが

random.choice()

を使うともっと短くなる

初心者プログラマ

いろいろあって職種を変えて転職予定。
今度はプログラミングをします。
今までも触り程度にコーディングしていたけれど仕事になるほどできるんだろうか・・・。

今後のために書いたコード、これから調べたことを記録していきます。

まずはコードの埋め込みってどうやるんだろうと思って調べてみたら
はてな記法ってのがあるんですね。

はてな記法一覧 - はてなダイアリーのヘルプ
記事中へのソースコードの埋め込み - 駆け出しプログラマの備忘録
はてなブログにソースコードを埋め込む - どこまでいってもプログラマ
はてなブログにソースコードを載せる3つの方法 - つばさのーと

">||"と"||<"でコードを挟むといいらしい。

例えば
print("Hello World!")
を">||"と"||<"で挟むと

print("Hello World!")

こうなる。
">|python|"とすると色分けしてくれる。

print("Hello World!")

Rも">|r|"でいけた。

print("Hello World!")