✈ – ページ 2

numo-gnuplot playground

投稿日: 2017年8月6日2017年8月6日投稿者: k2

require 'daru'
require 'numo/gnuplot'

df = Daru::DataFrame.from_csv "https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/datasets/AirPassengers.csv"

Numo.gnuplot do
  set :title, "Air Passengers", :font, "Arial, 20"
  set :xlabel, "Year"
  set :ylabel, "Passengers"
  set :grid
  set :key, :left
  plot df.time, df.AirPassengers, w: :linespoints, pt:6, lw:2, lc:'rgb "blue"', t:"Passengers"
end

require 'daru'

require 'numo/gnuplot'

df = Daru::DataFrame.from_csv "https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/datasets/AirPassengers.csv"

Numo.gnuplot do

set :title, "Air Passengers", :font, "Arial, 20"

set :xlabel, "Year"

set :ylabel, "Passengers"

set :grid

set :key, :left

plot df.time, df.AirPassengers, w: :linespoints, pt:6, lw:2, lc:'rgb "blue"', t:"Passengers"

end

df = Daru::DataFrame.from_csv "https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/datasets/BJsales.csv"

Numo.gnuplot do
  set :title, "BJsales", :font, "Time, 20"
  set :xlabel, "Time"
  set :ylabel, "BJsales"
  unset :key
  plot df.time, df.BJsales, w: :boxes
end

df = Daru::DataFrame.from_csv "https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/datasets/BJsales.csv"

Numo.gnuplot do

set :title, "BJsales", :font, "Time, 20"

set :xlabel, "Time"

set :ylabel, "BJsales"

unset :key

plot df.time, df.BJsales, w: :boxes

end

df = Daru::DataFrame.from_csv "https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/datasets/CO2.csv"

Numo.gnuplot do
  set title: "Carbon Dioxide Uptake in Grass Plants ", font: "Times, 20"
  set xlabel: "Day"
  set ylabel: "CO2"
  set xrange: [91..334]
  unset :key
  set yrange: [350..400]
  plot df.Day, df.CO2, w: "filledcurve x1", lc: 'rgb "navy"'
end

df = Daru::DataFrame.from_csv "https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/datasets/CO2.csv"

Numo.gnuplot do

set title: "Carbon Dioxide Uptake in Grass Plants ", font: "Times, 20"

set xlabel: "Day"

set ylabel: "CO2"

set xrange: [91..334]

unset :key

set yrange: [350..400]

plot df.Day, df.CO2, w: "filledcurve x1", lc: 'rgb "navy"'

end

今日の落書き red-arrowで遊ぶ

投稿日: 2017年8月5日2017年8月6日投稿者: k2

難しそうなので今日まで敬遠していたred-arrowを使ってみた。
公式の通りインストール。手順通りコピペでインストールしているだけなので何もわかっていないけど、そういうアホなユーザーも世の中一定数いるから仕方ないよね。Pythonの方もcondaでPyArrowをインストールした。

難しそうなので公式のexampleを試すだけ

パイソン

import pyarrow as pa
import numpy as np

ndarray = np.random.randn(10,6)
print(ndarray)

tensor = pa.Tensor.from_numpy(ndarray)
with pa.OSFile("/dev/shm/tensor.arrow", "wb") as sink:
    pa.write_tensor(tensor, sink)

import pyarrow as pa

import numpy as np

ndarray = np.random.randn(10,6)

print(ndarray)

tensor = pa.Tensor.from_numpy(ndarray)

with pa.OSFile("/dev/shm/tensor.arrow", "wb") as sink:

pa.write_tensor(tensor, sink)

ルビー

require 'arrow'
require 'arrow-numo-narray'

Input = Arrow::MemoryMappedInputStream

tensor = nil

Input.open("/dev/shm/tensor.arrow") do |input| 
  tensor = input.read_tensor(0)
end

p tensor.to_narray

require 'arrow'

require 'arrow-numo-narray'

Input = Arrow::MemoryMappedInputStream

tensor = nil

Input.open("/dev/shm/tensor.arrow") do |input|

tensor = input.read_tensor(0)

end

p tensor.to_narray

おお、普通に動くじゃないか…。
/dev/shmというのが何か気になったのでインターネッツでググるとRAMディスクのことらしい（正確には違うらしい）。

ということは…なにこれすごい∑( ﾟдﾟ )

公式スライド

RubyもApache Arrowでデータ処理言語の仲間入り

正直Webサイト製作に関するプロの道具であることを求めなければ、Rubyでデータ解析できる状況は十分に整ってきていると感じる。

機械学習に関してはscikit-learnのようなライブラリの力を借りる必要があり、pycallを呼び出して適宜red-arrowでデータを変換していくような使い方が望ましいのかもしれない。

RubyやPythonといった便利な道具があることを知らない or キャリアをまともにこなしていると学習コストを支払えない制度設計になっているために、プログラムでやった方が早いことも含めてプログラムと縁のない方法で仕事をこなさざるをえない非エンジニアの世界の方が、エンジニアの世界よりも広いのである。。

データを用意する

https://raw.githubusercontent.com/pandas-dev/pandas/master/pandas/tests/data/iris.csv

wget https://raw.githubusercontent.com/pandas-dev/pandas/master/pandas/tests/data/iris.csv

1	wget https://raw.githubusercontent.com/pandas-dev/pandas/master/pandas/tests/data/iris.csv

Jupyter Notebookを起動する

source activate fugaenv
jupyter-notebook

1 2	source activate fugaenv jupyter-notebook

※ホビーでRubyを手軽に使うのにJupyter Notebookが必ずしもベストとは思わない…

ライブラリを読み込む

require 'numo/narray'
require 'numo/gnuplot'
require 'statsample'
require 'nyaplot'

N = Numo
D = Daru
S = Statsample

require 'numo/narray'

require 'numo/gnuplot'

require 'statsample'

require 'nyaplot'

N = Numo

D = Daru

S = Statsample

CSVファイルを読み込む

df = D::DataFrame.from_csv "iris.csv"
df.head

1 2	df = D::DataFrame.from_csv "iris.csv" df.head

	SepalLength	SepalWidth	PetalLength	PetalWidth	Name
0	5.1	3.5	1.4	0.2	Iris-setosa
1	4.9	3.0	1.4	0.2	Iris-setosa
2	4.7	3.2	1.3	0.2	Iris-setosa
3	4.6	3.1	1.5	0.2	Iris-setosa
4	5.0	3.6	1.4	0.2	Iris-setosa
5	5.4	3.9	1.7	0.4	Iris-setosa
6	4.6	3.4	1.4	0.3	Iris-setosa
7	5.0	3.4	1.5	0.2	Iris-setosa
8	4.4	2.9	1.4	0.2	Iris-setosa
9	4.9	3.1	1.5	0.1	Iris-setosa

アヤメの品種を確認する

df.Name.uniq

1	df.Name.uniq

	Name
0	Iris-setosa
50	Iris-versicolor
100	Iris-virginica

基礎統計量を表示してみる

df.describe.round 3

1	df.describe.round 3

	PetalLength	PetalWidth	SepalLength	SepalWidth
count	150.0	150.0	150.0	150.0
mean	3.759	1.199	5.843	3.054
std	1.764	0.763	0.828	0.434
min	1.0	0.1	4.3	2.0
max	6.9	2.5	7.9	4.4

グループにしてみる

gpb = df.group_by("Name")
gpb.mean.round 3

1 2	gpb = df.group_by("Name") gpb.mean.round 3

	SepalLength	SepalWidth	PetalLength	PetalWidth
Iris-setosa	5.006	3.418	1.464	0.244
Iris-versicolor	5.936	2.77	4.26	1.326
Iris-virginica	6.588	2.974	5.552	2.026

アヤメの品種ごとにデータをわける

setosa     = df.where(df.Name.eq "Iris-setosa")
versicolor = df.where(df.Name.eq "Iris-versicolor")
virginica  = df.where(df.Name.eq "Iris-virginica")

setosa = df.where(df.Name.eq "Iris-setosa")

versicolor = df.where(df.Name.eq "Iris-versicolor")

virginica = df.where(df.Name.eq "Iris-virginica")

SepalLengthのヒストグラム DataFrameから

Nyaplotを使う場合

df.SepalLength.plot type: :histogram do |plot, diagram|
  plot.xrange [3.5, 8.5]
  plot.yrange [0, 25]
  plot.y_label "Frequency"
  plot.x_label "SepalLength"
end

df.SepalLength.plot type: :histogram do |plot, diagram|

plot.xrange [3.5, 8.5]

plot.yrange [0, 25]

plot.y_label "Frequency"

plot.x_label "SepalLength"

end

※ Daruから NyaplotやGnuplotrb、Plotlyを通してグラフを描けるらしい. 残念ながらNyaplotの開発は止まっている…

SepalLengthのヒストグラム　Statsampleから

# SVG IRubyに表示
d = S::Graph::Histogram.new(setosa.SepalLength, bins:13)
# ターミナルに表示
puts S::Graph::Histogram.new(setosa.SepalLength, bins:13).summary

# SVG IRubyに表示

d = S::Graph::Histogram.new(setosa.SepalLength, bins:13)

# ターミナルに表示

puts S::Graph::Histogram.new(setosa.SepalLength, bins:13).summary

StatsampleがDaruと統合する前に独自に持っていたグラフ機能が残っている

setosaの散布図を描く

Nyaplotを使うとき

setosa.plot type: :scatter, x: 'SepalLength', y: 'SepalWidth' do |plot, diagram|
  plot.x_label "SepalLength"
  plot.y_label "SepalWidth "
end

setosa.plot type: :scatter, x: 'SepalLength', y: 'SepalWidth' do |plot, diagram|

plot.x_label "SepalLength"

plot.y_label "SepalWidth "

end

Numo::Gnuplotを使うとき

N.noteplot do
  set title: "Setosa 散布図"
  set xlabel: "SepalLength"
  set ylabel: "SepalWidth"
  set key: :left
  plot [setosa.SepalLength.to_a, setosa.SepalWidth.to_a, w: :points, pt:7, t:"setosa"]
end

N.noteplot do

set title: "Setosa 散布図"

set xlabel: "SepalLength"

set ylabel: "SepalWidth"

set key: :left

plot [setosa.SepalLength.to_a, setosa.SepalWidth.to_a, w: :points, pt:7, t:"setosa"]

end

Daru::Vector→RubyArrayへの変換が必要。

相関係数を求める

setosa.corr.round(3)

1	setosa.corr.round(3)

	SepalLength	SepalWidth	PetalLength	PetalWidth
SepalLength	1.0	0.747	0.264	0.279
SepalWidth	0.747	1.0	0.177	0.28
PetalLength	0.264	0.177	1.0	0.306
PetalWidth	0.279	0.28	0.306	1.0

散布図行列

現状では簡単に描く方法は存在していないとおもわれる。

単回帰分析を行う

lr = S::Regression.simple(setosa.SepalLength, setosa.SepalWidth)
puts lr.summary

1 2	lr = S::Regression.simple(setosa.SepalLength, setosa.SepalWidth) puts lr.summary

= Regression of SepalLength over SepalWidth
  Table 1
+----------+--------+
| Variable | Value  |
+----------+--------+
| r        | 0.747  |
| r^2      | 0.558  |
| a        | -0.623 |
| b        | 0.807  |
| s.e      | 0.256  |
+----------+--------+

= Regression of SepalLength over SepalWidth

Table 1

+----------+--------+

| Variable | Value |

+----------+--------+

| r | 0.747 |

| r^2 | 0.558 |

| a | -0.623 |

| b | 0.807 |

| s.e | 0.256 |

+----------+--------+

Numo::Gnuplotで回帰直線を描く

N.noteplot do
  set title: "散布図と単回帰分析"
  set xlabel: "SepalLength"
  set ylabel: "SepalWidth"
  set key: :left
  plot [setosa.SepalLength.to_a, setosa.SepalWidth.to_a, w: :points, pt:7, t:"setosa"],
  ["#{lr.b}*x+#{lr.a}", lw: 3, t:"単回帰"]
end

N.noteplot do

set title: "散布図と単回帰分析"

set xlabel: "SepalLength"

set ylabel: "SepalWidth"

set key: :left

plot [setosa.SepalLength.to_a, setosa.SepalWidth.to_a, w: :points, pt:7, t:"setosa"],

["#{lr.b}*x+#{lr.a}", lw: 3, t:"単回帰"]

end

多分Nyaplotでも描ける

重回帰分析

# df[0..-2] Nameを取り除く
lr = S::Regression.multiple(df[0..-2], "SepalWidth")
puts lr.summary

# df[0..-2] Nameを取り除く

lr = S::Regression.multiple(df[0..-2], "SepalWidth")

puts lr.summary

= Multiple reggresion of SepalLength,PetalLength,PetalWidth on SepalWidth
  Engine: Statsample::Regression::Multiple::RubyEngine
  Cases(listwise)=150(150)
  R=0.724
  R^2=0.524
  R^2 Adj=0.514
  Std.Error R=0.302
  Equation=1.035 + 0.608SepalLength + -0.588PetalLength + 0.568PetalWidth
  == ANOVA
    ANOVA Table
+------------+--------+-----+-------+--------+-------+
|   source   |   ss   | df  |  ms   |   f    |   p   |
+------------+--------+-----+-------+--------+-------+
| Regression | 14.667 | 3   | 4.889 | 53.487 | 0.000 |
| Error      | 13.345 | 146 | 0.091 |        |       |
| Total      | 28.013 | 149 | 4.980 |        |       |
+------------+--------+-----+-------+--------+-------+

  Beta coefficients
+-------------+--------+--------+-------+--------+
|    coeff    |   b    |  beta  |  se   |   t    |
+-------------+--------+--------+-------+--------+
| Constant    | 1.035  | -      | 0.269 | 3.842  |
| SepalLength | 0.608  | 1.160  | 0.062 | 9.823  |
| PetalLength | -0.588 | -2.394 | 0.062 | -9.526 |
| PetalWidth  | 0.568  | 1.000  | 0.122 | 4.669  |
+-------------+--------+--------+-------+--------+

= Multiple reggresion of SepalLength,PetalLength,PetalWidth on SepalWidth

Engine: Statsample::Regression::Multiple::RubyEngine

Cases(listwise)=150(150)

R=0.724

R^2=0.524

R^2 Adj=0.514

Std.Error R=0.302

Equation=1.035 + 0.608SepalLength + -0.588PetalLength + 0.568PetalWidth

== ANOVA

ANOVA Table

+------------+--------+-----+-------+--------+-------+

| source | ss | df | ms | f | p |

+------------+--------+-----+-------+--------+-------+

| Regression | 14.667 | 3 | 4.889 | 53.487 | 0.000 |

| Error | 13.345 | 146 | 0.091 | | |

| Total | 28.013 | 149 | 4.980 | | |

+------------+--------+-----+-------+--------+-------+

Beta coefficients

+-------------+--------+--------+-------+--------+

| coeff | b | beta | se | t |

+-------------+--------+--------+-------+--------+

| Constant | 1.035 | - | 0.269 | 3.842 |

| SepalLength | 0.608 | 1.160 | 0.062 | 9.823 |

| PetalLength | -0.588 | -2.394 | 0.062 | -9.526 |

| PetalWidth | 0.568 | 1.000 | 0.122 | 4.669 |

+-------------+--------+--------+-------+--------+

【Ruby】 Numo::Gnuplot で箱ひげ図 (boxplot)

投稿日: 2017年7月9日2017年7月9日投稿者: k2

忘れないうちにメモ

1.配列をさっと比較したいとき

require 'numo/gnuplot'

a = Array.new(100) { rand }
b = Array.new(100) { rand * 2 }
c = Array.new(100) { rand * 3 }

Numo.gnuplot do
  plot [a, u: '(0):1', w: :boxplot],
       [b, u: '(1):1', w: :boxplot],
       [c, u: '(2):1', w: :boxplot]
end

gets

require 'numo/gnuplot'

a = Array.new(100) { rand }

b = Array.new(100) { rand * 2 }

c = Array.new(100) { rand * 3 }

Numo.gnuplot do

plot [a, u: '(0):1', w: :boxplot],

[b, u: '(1):1', w: :boxplot],

[c, u: '(2):1', w: :boxplot]

end

gets

2.ラベル付きのデータのとき

require 'numo/gnuplot'

label = [*'A'..'J'] * 5
data  = Array.new(50) { |i| i }

Numo.gnuplot do
  plot label, data, u: '(0):2:(0):1', w: :boxplot
end

gets

require 'numo/gnuplot'

label = [*'A'..'J'] * 5

data = Array.new(50) { |i| i }

Numo.gnuplot do

plot label, data, u: '(0):2:(0):1', w: :boxplot

end

gets

using の設定がややこしいけど、丸暗記するしかないっぽいね

医療とAIの未来について再び妄想する

投稿日: 2017年6月4日2017年8月7日投稿者: k2

頭の整理をするために、また医療とAIの未来について考えている。

医療データサイエンス 1.0

これは、近い将来にきっと起きる変化だ。AIの医療応用を考える多くの人が思い描いている未来だ。
具体的には
・CTやMRI画像などが、Deep Learning技術で自動的に診断できるようになる
・医療用のデータセットを作成する患者団体や、病院、学会が現れる
・APIを用いて医療データを取得する手段が整備されはじめる
・データセンタを持つ企業が医療データを低価格で診断するサービスを始めてデータを収集・蓄積しはじめる　
・対抗して、よりオープンな医療データセットやオープンなモデルも公開される
・大病院が工学技師としてデータサイエンティストを雇いはじめる
・もしくは医者がデータサイエンスをかじりはじめる
・少しずつ人工知能が普及しはじめる
・画像診断の精度、血液検査の精度、内視鏡の精度、遺伝子の解析、心電図モニタリングetc
　すべてでコンピュータが人間を上回る
・個別化された医療が行われるようになる
・診断だけではなく、治療にロボットの活用がはじまる
・確実にどんどん成果を生まれ、治療成績も激的に向上していく
・にもかかわらず、不思議なほど医療そのものには何らの変化がない　なんでだろう？
・保険点数が厳しく減額され、利益が減少したりする
・IT企業に対する患者の医療訴訟が多発するようになる。高額の賠償が請求される事例も発生する
・IT企業も、次第に状況を理解して（最初は不合理だと思っていた）製薬会社のやり方を真似するようになる
・もしくはもともと無料でサービスを提供できる体力のあるところだけが他の会社を焼き払う形で生き残る
・最初からこうなることを見越していた老舗企業もしぶとく生き残る
・患者とエンジニアの間にゆっくりと失望が広がっていく
・業界の熱気が冷めて、少しずつ落ち着いてくる
・なんかいい感じになる
といった感じだろうか。とってもワクワクしないよね。あまりにも想像するのが容易だ。
妄想にしても陳腐だ。
こういったことに無邪気に夢を感じるには、少し歳をとりすぎた気がする。

問題は、ここからだと思っている。上記の筋書きはちゃんと考えれば誰もが見通せる素直な未来だ。
未来は、だいたい思い描いた通り素直には歩いてくれない。

医療データサイエンス 2.0

上で描いた1.0は人工知能は、あまり医療の本質に触れていない。
つまり1.0は来るが、そのままの形ではこない。

僕が思いつくだけで大きな焦点は２つあると思う。
１つ目の焦点は、割と簡単だ。簡単だということは、乗り越えうるということである。
乗り越えた先は、たとえばアルファ碁のあとの囲碁界みたいな事になっていると思う。ある意味楽しい世界である。

２つ目の焦点も簡単だが、もう少し医療の本質に関わる問題である。
この問題が乗り越えられるかどうかはわからない。乗り越えるとしたら、人工知能は俺が考えているよりもやばい。
少なくとも、今の段階では人工知能がこういった問題を乗り越えることは想像できないな。多分、人工知能はこの問題を乗り越えられずに迂回すると思う。
少なくとも人間は、こういう問題を直接乗り越えるよりは迂回することを好む。

この問題について、ちょうどいいツイートがあったので引用しておく

人工知能の医療応用の壁は、技術を越えたあとの「命と命の比較」、「命と経済の比較」という答えのない、そして文化や思想信条に関わる質問に、答えを出して可視化しなければならないと言うことだと思う。
今はカンファレンスの中で医師の常識や成り行き等で処理されてきたが、そこが難敵になるはず。

— 人工知能放射線科応用研 (@71LP6ETOwR8U9tg) 2017年7月19日

そしてその根本にある問題や、その他の問題もおぼろげに見える気がする。

いつものことだが妄想日記にしても考えがよくまとまらない。

考えれば考えるほど、やっぱり人工知能も今までの技術革新と同じで、結局は医療は変わらないという説の方を採用したくなるね。

Rubyを使ってDICOM画像で遊んでみる

投稿日: 2017年5月25日2017年5月30日投稿者: k2

RubyでDICOM画像を扱う情報がほとんどないので少しずつメモしてみる。

ruby-dicomhttps://github.com/dicom/ruby-dicom

インストール

gem install dicom

1	gem install dicom

DICOM画像を読み込む

require 'dicom'
dcm = DICOM::DObject.read("test.dcm")

1 2	require 'dicom' dcm = DICOM::DObject.read("test.dcm")

メタデータを取り出す

dcm.modality　# {"Modality"=>"CT"}

1	dcm.modality　# {"Modality"=>"CT"}

画像データを取り出す

pixels = dcm.pixels
pixels_narray = dcm.narray # 旧NArray ：画像によっては正しく動かないこともあるので注意

1 2	pixels = dcm.pixels pixels_narray = dcm.narray # 旧NArray ：画像によっては正しく動かないこともあるので注意

IRuby Notebook上にDICOM画像を表示する

imageメソッドで、ImageMagickオブジェクトに変換できる。
内部ではrmagickもしくはmini_magickが動いている。
rmagickは更新が滞っているので、mini_magickがおすすめ。

image = dcm.image # Magick::Image
# normalizeするともっと見栄えがよくなる
image = image.normalize

image = dcm.image # Magick::Image

# normalizeするともっと見栄えがよくなる

image = image.normalize

※画像はDSB2017に使用されたsampleの1つです

ヒストグラムを表示する

Guido ZuidhofFull Preprocessing Tutorial
の序盤をRubyに書き換えてみた。とりあえず動いてくれればOKという発想なので、例によって間違いあるかも。

require 'numo/narray'
require 'numo/gnuplot'
require 'dicom'
require 'histogram/array'

# Some contents
INPUT_FOLDER = 'stage1'
patients = Dir.entries(INPUT_FOLDER)
patients = patients[2..-1] # remove './', '../'

require 'numo/narray'

require 'numo/gnuplot'

require 'dicom'

require 'histogram/array'

# Some contents

INPUT_FOLDER = 'stage1'

patients = Dir.entries(INPUT_FOLDER)

patients = patients[2..-1] # remove './', '../'

class DICOM::Element
  def values
    value.split("\\").map(&:to_f)
  end
end

def load_scan(path)
  file_paths = Dir.glob(File.join(INPUT_FOLDER, path, '/*'))
  slices = file_paths.map{ |fp| DICOM::DObject.read(fp)}
  slices.sort_by!{ |s| s.image_position_patient.values[2]}
  slice_thickness = (
    slices[0].image_position_patient.values[2] -
      slices[1].image_position_patient.values[2] ||
      slices[0].slice_location.values[2] -
      slices[1].slice_location.values[2]
    ).abs
  slices.each do |s|
    s.slice_thickness = slice_thickness
  end
end

class DICOM::Element

def values

value.split("\\").map(&:to_f)

end

def load_scan(path)

file_paths = Dir.glob(File.join(INPUT_FOLDER, path, '/*'))

slices = file_paths.map{ |fp| DICOM::DObject.read(fp)}

slices.sort_by!{ |s| s.image_position_patient.values[2]}

slice_thickness = (

slices[0].image_position_patient.values[2] -

slices[1].image_position_patient.values[2] ||

slices[0].slice_location.values[2] -

slices[1].slice_location.values[2]

).abs

slices.each do |s|

s.slice_thickness = slice_thickness

end

def get_pixels_hu(slices)
  image = Numo::Int16[*slices.map(&:pixels)]
  slices.each_with_index do |slice, i|
    intercept = slice.rescale_intercept.value.to_f
    slope = slice.rescale_slope.value.to_f
    if slope != 1.0
      image[i,true] = slope * image[i,true] # 自動でfloatにキャストされ、なおかつInt16に戻るみたい
      image[i,true] += intercept
    end
  end
  image
end

def get_pixels_hu(slices)

image = Numo::Int16[*slices.map(&:pixels)]

slices.each_with_index do |slice, i|

intercept = slice.rescale_intercept.value.to_f

slope = slice.rescale_slope.value.to_f

if slope != 1.0

image[i,true] = slope * image[i,true] # 自動でfloatにキャストされ、なおかつInt16に戻るみたい

image[i,true] += intercept

end

image

end

first_10_patients = Array.new(10) { |pt| load_scan(patients[pt]) }
first_10_patients_pixels = first_10_patients.map! { |pt| get_pixels_hu(pt) }

1 2	first_10_patients = Array.new(10) { \|pt\| load_scan(patients[pt]) } first_10_patients_pixels = first_10_patients.map! { \|pt\| get_pixels_hu(pt) }

# histogram 旧NArray対応のみ...なのでRuby配列に一旦変換
first_10_patients_pixels.each do |ptp|
  first_patient_pixels_array = ptp.flatten.to_a
  bins,freq = first_patient_pixels_array.histogram(50)
  bins.map!{|b| b.round.to_s}

  noteplot = Numo.noteplot do
    set title:"ヒストグラム"
    set style:[:fill, :solid, {border:-1}]
    set xtics:["rotate by": -45]
    plot bins, freq, u:"2:xtic(1)", w: :boxes, lc:{rgb:"orange"}, t:"分布"
  end
  
  IRuby.display noteplot
end

# histogram 旧NArray対応のみ...なのでRuby配列に一旦変換

first_10_patients_pixels.each do |ptp|

first_patient_pixels_array = ptp.flatten.to_a

bins,freq = first_patient_pixels_array.histogram(50)

bins.map!{|b| b.round.to_s}

noteplot = Numo.noteplot do

set title:"ヒストグラム"

set style:[:fill, :solid, {border:-1}]

set xtics:["rotate by": -45]

plot bins, freq, u:"2:xtic(1)", w: :boxes, lc:{rgb:"orange"}, t:"分布"

end

IRuby.display noteplot

end

今日の落書き Ruby の PyCall で Keras の Mnist 画像部類

投稿日: 2017年5月3日2017年5月3日投稿者: k2

RubyでPythonっぽいコードを書くと、RubyでPythonが実行できるというPyCallというGemがあります。
このPyCallでDeep Learning用のライブラリであるKerasが動かせるということなので、やってみました。

こちらのサイトを参考にしました。
PyCallを使えばRubyでもKerasでDeep Learningができる（洋食の日記）

require 'pycall/import'
require 'optparse'
# Mnist データ取得用
require 'ruby_brain'
require 'ruby_brain/dataset/mnist/data'
# Numo::NArray
require 'numo/narray'
require 'tempfile'

# NArrayをNumPyに変換するモンキーパッチ
class Numo::SFloat
  def to_numpy(type, xp)
    shape = self.shape
    p shape
    string = self.to_string
    file = Tempfile.new('tmp')
    begin
      File.binwrite(file.path, string)
      ndarray = xp.fromfile.(file.path, dtype:type).reshape.(shape)
    ensure
      file.close
      file.unlink
    end
    ndarray
  end
end

batch_size = 128
epochs = 20

opt = OptionParser.new
opt.on('--batchsize=VAL', '-b VAL', 'Number of images in each mini-batch'){ |v|
  batch_size = v.to_i
}
opt.on('--epoch=VAL', '-e VAL', 'Number of sweeps over the dataset to train'){|v|
  epochs = v.to_i
}
opt.parse!(ARGV)

training_dataset, test_dataset = RubyBrain::DataSet::Mnist::data
x_train = Numo::SFloat[*training_dataset[:input]] / 255
x_test  = Numo::SFloat[*test_dataset[:input]] / 255
y_train = Numo::SFloat[*training_dataset[:output]]
y_test  = Numo::SFloat[*test_dataset[:output]]

include PyCall::Import
pyimport :numpy, as: :np
pyimport :keras
pyfrom 'keras.models', import: 'Sequential'
pyfrom 'keras.layers', import: ['Dense', 'Dropout']
pyfrom 'keras.optimizers', import: 'RMSprop'


x_train = x_train.to_numpy(:float32, np)
x_test  = x_test.to_numpy(:float32, np)
y_train = y_train.to_numpy(:float32, np)
y_test  = y_test.to_numpy(:float32, np)

model = Sequential.()
model.add.(Dense.(512, activation: :relu, input_shape:[784]))
model.add.(Dropout.(0.2))
model.add.(Dense.(512, activation: :relu))
model.add.(Dense.(10, activation: :softmax))

model.summary.()

model.compile.(loss: :categorical_crossentropy,
               optimizer: RMSprop.(),
               metrics:['accuracy'])

model.fit.(x_train, y_train,
           batch_size: batch_size,
           epochs: epochs,
           verbose: 1,
           validation_data: [x_test, y_test])

score = model.evaluate.(x_test, y_test, verbose:0)
print(sprintf("Test loss: %.6f\n", score[0]))
print(sprintf("Test accuracy: %.6f\n", score[1]))

require 'pycall/import'

require 'optparse'

# Mnist データ取得用

require 'ruby_brain'

require 'ruby_brain/dataset/mnist/data'

# Numo::NArray

require 'numo/narray'

require 'tempfile'

# NArrayをNumPyに変換するモンキーパッチ

class Numo::SFloat

def to_numpy(type, xp)

shape = self.shape

p shape

string = self.to_string

file = Tempfile.new('tmp')

begin

File.binwrite(file.path, string)

ndarray = xp.fromfile.(file.path, dtype:type).reshape.(shape)

ensure

file.close

file.unlink

end

ndarray

end

batch_size = 128

epochs = 20

opt = OptionParser.new

opt.on('--batchsize=VAL', '-b VAL', 'Number of images in each mini-batch'){ |v|

batch_size = v.to_i

}

opt.on('--epoch=VAL', '-e VAL', 'Number of sweeps over the dataset to train'){|v|

epochs = v.to_i

}

opt.parse!(ARGV)

training_dataset, test_dataset = RubyBrain::DataSet::Mnist::data

x_train = Numo::SFloat[*training_dataset[:input]] / 255

x_test = Numo::SFloat[*test_dataset[:input]] / 255

y_train = Numo::SFloat[*training_dataset[:output]]

y_test = Numo::SFloat[*test_dataset[:output]]

include PyCall::Import

pyimport :numpy, as: :np

pyimport :keras

pyfrom 'keras.models', import: 'Sequential'

pyfrom 'keras.layers', import: ['Dense', 'Dropout']

pyfrom 'keras.optimizers', import: 'RMSprop'

x_train = x_train.to_numpy(:float32, np)

x_test = x_test.to_numpy(:float32, np)

y_train = y_train.to_numpy(:float32, np)

y_test = y_test.to_numpy(:float32, np)

model = Sequential.()

model.add.(Dense.(512, activation: :relu, input_shape:[784]))

model.add.(Dropout.(0.2))

model.add.(Dense.(512, activation: :relu))

model.add.(Dense.(10, activation: :softmax))

model.summary.()

model.compile.(loss: :categorical_crossentropy,

optimizer: RMSprop.(),

metrics:['accuracy'])

model.fit.(x_train, y_train,

batch_size: batch_size,

epochs: epochs,

verbose: 1,

validation_data: [x_test, y_test])

score = model.evaluate.(x_test, y_test, verbose:0)

print(sprintf("Test loss: %.6f\n", score[0]))

print(sprintf("Test accuracy: %.6f\n", score[1]))

短いコードですが、半年ぐらいプログラミングに触れていなかったのでかなり時間がかかりました。

工夫したところ

・RubyからPython呼んでいる時点で無駄が多いことをしている。だから何も気にせずに、ライブラリを贅沢に使おうと思った。
・MnsitをRubyの配列で呼ぶためだけにruby_brainを使用した。
・255で割るためだけにNumo/NArrayを使用した。
・結果として、mnistの元ファイル ⇢ Ruby配列 ⇢ NArray配列 ⇢ テキストファイル ⇢ numpy配列という情報の変換、受け渡しにとても時間がかかっている。
　（・しかし、これは、PCを使わない人工知能にも広く見られる現象だ）
・ハマったのは、Kerasでmnist読み込む時と、ruby_brainでnnist呼びこむ時のデータ形式の違い。

よかったところ

・うごく！
・Rubyでディープラーニングする手段がある。0と1では大違い。

不満が残るところ

・PyCallでKerasを実行している最中に Ctrl+C でターミナルを停止できない
・エラーメッセージがもっと改行してくれると嬉しい
・やはり Ruby ←→ Pyrhon の情報のやりとりに不自由を感じる
　・Rubyで前処理して一度ファイルに保存、Pythonで処理するのが現実的かもしれない

ASUS ZenBook に Keras を入れるまで

投稿日: 2017年4月29日2017年5月4日投稿者: k2

常用していた中古のノートパソコンをMacbookに買い換えました。周囲のMac使用率が50%に近づいている現在、Macを購入する最大のメリットは、周囲とのコミュニケーションの心理的な促進作用だと思います。Macには人と人とのつながり、コミュニケーションを誘発する作用があります。そこにはスペックだけ見ているとわからない、お金にできない、大きな価値があります。それに、Macを使っていると気分がいいです。Macについて語るとき、みんな笑顔になります。洗練されたGUIとシンプルなソフトウェアに囲まれて作業していると、落ち着いた気持ちになります。そう、間違いなくMacを選択するメリットは大きいのです。そうだMacbook Proを購入しよう。2年ぶりにUbuntuからMacに戻ろう。鼓動の高なりを感じながら、僕は家を出た。

そして、ASUSとともに帰宅。

Intelのシールの右に、緑色のNVIDIAのシールが張られている。GeForce 940MX搭載である。つまりCudaが使える。それからMacは高い。そりゃあなた、お金は大事ですよ。（しかし現行MacにGeForce搭載されていたら本当にMacbookを買っていたと思う）

以下困った時のためにKeras導入までの自分用メモ。

① まず購入して真っ先に、Intel等のシールをはがす。Macでは不要な作業。

② ブートUSBの作成
UNetbootinではなくこちらを使用した。とくに不具合なし。Ubuntu 16.04
UbuntuのブートUSB作成 (usb-creator を使用)

③ 画面の解像度がおかしい
→ソフトウェアアップデートして再起動

④ Cuda8.0のインストール
第456回　Ubuntu 16.04 LTSでCUDA 8.0を使用する 2017年2月1日柴田充也
→ こちらの情報が適確で初心者にもわかりやすいのでその通りにした
→ deb(network)を選択した
→ 自動的にnvidia-375がインストールされる。NVIDIAのサイトでsupported productにGeForce 940MXが含まれていることを確認
→ PATHの設定、LD_LIBRARY_PATHの設定
→ nvcc –version や echo $LD_LIBRARY_PATH で確認

⑤ nvidia-smiコマンドが見つからない
→ F2でBIOSの設定に入り secure boot をdisableに変更して起動

⑥ cuDNN のインストール
CUDA 8.0とcuDNN 5.1をUbuntu 16.04LTSにインストールする
→ こちらの通りにした。
→ cuDNN のバージョンは要注意

⑦ Anacondaの導入
公式サイトの通り

⑧ 新しい環境の作成:
conda create -n keras .anaconda python=3.5
→ 最新過ぎるとはまりやすいという感覚で何となく3.5に

⑧ tensorflowのインストール
公式サイトの通り(Ubuntu/Anaconda使用版)

⑨ pythonでtensorflowを読み込んだ時のエラー
ImportError: libcudnn.5: cannot open shared object file: No such file or directory
→ まずはLD_LIBRARY_PATHの設定を見直し → 改善せず
→ stackoverflowの中段の解決策が有効
→ そもそもディレクトリのなかにlibcudnn.5が存在しないので、シンボリックリンクを作成
→ ln -s libcudnn.so.6* libcudnn.so.5
→ import tensorflow as tf が動作することを確認

⑩ Keras のインストール
公式サイトの通り
pip install keras
→ 公式exampleのmnist(MLP)が動作することを確認
→ nvidia-smi で確かにCPUではなくGPUで動作していることの確認

細かい不具合あるかもしれないが、これで手元でKerasを勉強できる環境が整った。

UbuntuでRuby2.4.1にTkをインストールする

投稿日: 2017年4月5日2017年4月5日投稿者: k2

TkがRuby本体から分離されましたが、オプションは今までの通りでも大丈夫なようです。

gem install tk — –with-tcltkversion=8.6 \
–with-tcl-lib=/usr/lib/x86_64-linux-gnu \
–with-tk-lib=/usr/lib/x86_64-linux-gnu \
–with-tcl-include=/usr/include/tcl8.6 \
–with-tk-include=/usr/include/tcl8.6 \
–enable-pthread

ゼロから作るDeep Learning ―Pythonで学ぶディープラーニングの理論と実装を買った

投稿日: 2016年10月20日2016年10月23日投稿者: k2

書店をぶらついていると、真っ赤なタイトルの本が目に飛び込んできた。

ページをめくってみた。すごい。勉強になる。数式読めない人けど、なんとなくDeep Learningが理解できるかも知れないきざしのきざしの兆候を感じる。すごい。うじゃひょぇえええええ！！！！！期待にドキドキ胸が鼓動を早める。そう、これは…買っただけで勉強しない時の予感。

ふと、ページをめくる指が震える。
おかしい。何かがおかしい。かすかに感じるこの違和感はなんだろう。もしかして、この素晴らしい本には、ひょっとして一点だけ欠陥と言わざるをえないところがあるのではないか。すべて整っているのに肝心のところで、どこか画竜点睛を欠くのではないか。

その原因は、最初はぼんやりとしていたが、次第に誰の目にもわかるほど明瞭になっていった。

どういうわけか、本文中のコードに Ruby + NArray が一行も使われてないのである。そのかわり Python + numpy という一風変わった言事とライブラリが使用されている。Python? パイソン？聞いたことがない新鮮な響きだ。Rubyはどこへいったのだろう？公園だと思ったら行ってみたら墓地だった。そんな衝撃を想像してほしい。

これでは作者の美的センスを疑わざるをえない状況である。それだけではない。近代文明社会への冒涜、民主主義への挑戦というものではなかろうか。久々に怒りで額に汗が流れた。グルグルと景色が回りはじめ、書棚に並ぶ本たちの笑い声がドッと湧いてこだました。

※ 本文はフィクションです。謹んで勉強させていただきます。

無知をさらす
・損失関数についてあまり意識してなかった
　・交差エントロピー誤差というものを知らなかった
　・なるほどsoftmaxと交差エントロピー誤差はセットなのね
・数値微分という言葉を知らなかった
・勾配によるニューラルネットワークの更新がよく理解できない
・chainerって連鎖率（chain rule）から名前がつけられたのだろうか
　・合成関数の微分は合成関数を構成するそれぞれの関数の微分の積によって表すことができる
・ハイパーパラメータという言葉はよく聞くが、なぜハイパーというのか知らなかった
・ミニバッチはランダムにサンプルを取得するものだと知らなかった（順番にサンプル取得するものだとばかり思っていた）
・加算、乗算の計算グラフによる逆伝播を理解した
・Adamが理解できないぞ？
・Xavierの初期値は、sigmoidやtanhに有効だが、ReLuには「Heの初期化」がある。へー。
・Batch Normalization 今度やってみよう
・ハイパーパラメータの最適化は結局総当り式？
・特徴マップという言葉を知らなかった
・im2colねえ…。これRubyだけで実装するとやっぱりfor(each)使わざるを得ない気がするが。

自信をもったところ
なんだかんだですでに実装したことがあるような内容も多い。①パーセプトロン、②2層ニューラルネットワーク、③計算グラフ　④バッチ処理　⑤CNNみたいな順番で勉強してきてよかったなと。経験的に好きな事はつまみ食い的に勉強した方が、かえって効率的に勉強できる気がする。むしろこの本だけでDeep Learning勉強しようとしていたらどこまで理解出来てたか疑問だ。

この本のおかげでchain ruleと計算グラフについて少し理解した。
加算ノードの逆伝播と、乗算ノードの逆伝播理解したから。Numo::NArrayでLSTMにもう一回挑戦してみたいな。

データを用意する

Jupyter Notebookを起動する

ライブラリを読み込む

CSVファイルを読み込む

アヤメの品種を確認する

基礎統計量を表示してみる

グループにしてみる

アヤメの品種ごとにデータをわける

SepalLengthのヒストグラム DataFrameから

SepalLengthのヒストグラム Statsampleから

setosaの散布図を描く

相関係数を求める

散布図行列

単回帰分析を行う

重回帰分析

1.配列をさっと比較したいとき

2.ラベル付きのデータのとき

医療データサイエンス 1.0

医療データサイエンス 2.0

IRuby Notebook上にDICOM画像を表示する

ヒストグラムを表示する

工夫したところ

よかったところ

不満が残るところ

SepalLengthのヒストグラム　Statsampleから