【幾何問題】二次元平面で2つの固定長の線の軌跡を求める

python ニュートンラフソン法 幾何学

ふと思い出した

DLのプログラムを書いていて、 ふと大学時代にmatlabで書いていた幾何問題を解く問題（※）に似ているな、と思い、 改めてpythonで簡単な問題を解いてみた。

※お車の話。最後に余談として記載。

最新MATLABハンドブック第六版

最新MATLABハンドブック第六版

二次元平面で2つの固定長の線の軌跡

一体なんなのか。

こんな線ACと線BCがあったとする。

f:id:kazuhitogo:20190211115918p:plain

線は剛体で固定長だとする。

そのとき、Bを動かしたとき、 Cはどこに行くの？

f:id:kazuhitogo:20190211115901p:plain

という問題。

ちなみに数学的にこの問題を解こうとすると、 Aを中心とする長さACの円と、 B'を中心とする長さBCの円の交点を求めればよいので、 下記方程式を解けばよい…のだが、それなりにめんどくさい。

\displaystyle ( x - A_x ) ^ 2 + ( y - A_y ) ^ 2 = (\vec{AC}) ^ 2

\displaystyle ( x - B ' _ x ) ^ 2 + ( y - B ' _ y ) ^ 2 = (\vec{BC}) ^ 2

というかこのような非線形な方程式は私には解けない。 で、それを現実的にどう解くか…というのが今回のエントリ。

ちなみにさきほどの方程式の解は基本的に２つある（２つの円の交点は基本２つ）が、 現実世界においてはBをB'に移動したらCはCに近い方の解にしか移動しないため、 そちらの座標を計算する。

どうやるのか

ある程度の誤差を許容しつつ（現実世界では0.0001mmの誤差くらいは許容、というより製造限界)、 コンピュータは簡単な問題を繰り返し解くのが得意なので、その方法で解く。

ニュートンラフソン法というのがある。 ja.wikipedia.org

  1. 簡単にいうと複雑な問題（今回で言えば二元二次方程式）を、線形の連立方程式として解く。

  2. 複雑な問題を簡易に問いているので当然結果には誤差が出る。

  3. 解き終わった後誤差を評価し、また同じ方法で解く…を繰り返し、ある一定の誤差より低くなったらやめる

というものである。

1の連立方程式として解くところがミソであるのだが、 Cの座標を微小に動かしたときに、それぞれの線の長さがどう変わるか、を算出し、 連立方程式化する。

で、DLやってて思い出した理由がここにあって、 誤差を評価して微小（＝つまり微分）に動かして修正して、 を繰り返して行くのが、まさにDLの発想だな、と思い出した次第。

改めて図示すると以下の通り。

初期位置

f:id:kazuhitogo:20190211115918p:plain

Bを移動

f:id:kazuhitogo:20190211185659p:plain

Cはどこに行く？

f:id:kazuhitogo:20190211185732p:plain

Bを動かすと赤に誤差が生じている

f:id:kazuhitogo:20190211185815p:plain

CをX方向に微小に動かして線分AC,BCの長さがどうなるかを求める

f:id:kazuhitogo:20190211185853p:plain

CをY方向に微小に動かして線分AC,BCの長さがどうなるかを求める

f:id:kazuhitogo:20190211185917p:plain

そこからは連立方程式を立てて解きましょう…。

コード

from matplotlib import pyplot as plt
import math
import numpy as np

# point A,B,C
pa = np.array([1,1])
pb =  np.array([3,-1])
pc = np.array([6,3])

# 元のデータ(描画用)
pa_org = np.array([1,1])
pb_org =  np.array([3,-1])
pc_org = np.array([6,3])


# 2点間の距離を求める
# a = np.array([x1,y1])
# b = np.array([x2,y2])
def len_2p(a,b):
    return np.linalg.norm(b-a)

# 2つの誤差の絶対値の合計を返す（誤差関数とする）
def e_2(a,b):
    return np.abs(a) + np.abs(b)

# 描画
def plt3p(a,b,c):
    x = np.array([a[0],b[0],c[0]])
    y = np.array([a[1],b[1],c[1]])
    plt.gca().set_aspect('equal', adjustable='box')
    plt.plot(x,y,"-*")

# 初期の長さ(i=init)
lab_i = len_2p(pa,pb)
lbc_i = len_2p(pb,pc)

# bをx方向に1動かす(動かしたときにcがどこに行くかを計算)
pb = pb + np.array([1,0])

# 移動後長さ(a=after)
lac_a = len_2p(pa,pc)
lbc_a = len_2p(pb,pc)

# 誤差(e=error)
lac_e = lac_i - lac_a
lbc_e = lbc_i - lbc_a

# 誤差計算
e = e_2(lac_e,lbc_e)

# 微小に動かす単位(今回はx,yともに1動かす)
unit_d = 1

# cが繰り返しの演算の中でどう動いていったかを残す配列
pc_hist = np.empty((2,0))
pc_hist = np.append(pc_hist,pc.reshape(2,-1),axis=1)

# 誤差が一定の数値以下になるまで回し続ける
# 解けない問題（先の連立方程式だと解が無いくらい離した場合など）は無限ループになるので気をつけるべし
while(e > 0.000000001): 
    # x方向に微小に動かしたpc
    pc_dx = pc + np.array([unit_d,0])
    # y方向に微小に動かしたpc
    pc_dy = pc + np.array([0,unit_d])
    
    # x方向に微小に動かしたpcの長さ評価
    lac_dx = len_2p(pa,pc_dx)
    lbc_dx = len_2p(pc_dx,pb)

    # y方向に微小に動かしたpcの長さ評価
    lac_dy = len_2p(pa,pc_dy)
    lbc_dy = len_2p(pc_dy,pb)

    # 連立方程式を解く
    # lac_dx * x + lac_dy * y  = lac_e * unit_d
    # lbc_dx * x + lbc_dy * y = lbc_e * unit_d
    # を解く
    move = np.linalg.solve(np.array([[lac_dx,lac_dy],[lbc_dx,lbc_dy]]),np.array([lac_e,lbc_e])) * unit_d
    pc = pc + move
    pc_hist = np.append(pc_hist,pc.reshape(2,-1),axis=1)
    
    # 移動後長さ(a=after)
    lac_a = len_2p(pa,pc)
    lbc_a = len_2p(pb,pc)
    
    # それぞれの線の長さの誤差(e=error)
    lac_e = lac_i - lac_a
    lbc_e = lbc_i - lbc_a
    
    # 先の長さ２つを誤差関数に入れる
    e = e_2(lac_e,lbc_e)

# 描画
plt3p(pa_org,pc_org,pb_org)
plt3p(pa,pc,pb)
plt.plot(pc_hist[0],pc_hist[1])

するとこうなる。 f:id:kazuhitogo:20190211221809p:plain 青が元の位置、オレンジが動かした位置、緑が連立方程式を解いたCの推移である。 かなり速く収束し、 164回で0.000000001未満の誤差に落ち着く。 参考までに0.001未満の精度だと53回である。

余談

大学時代はこれを６次元でやっており、 X,Y,Zの空間に、さらに回転のロール・ヨー・ピッチを入れて計算していた。

車のサスペンションのジオメトリ計算で、 マルチリンク式サスペンションを、 バンプした時、ステアを切った時、それぞれどう動くかを計算していたのである。

前処理18本ノック

前処理 bash shell sql

いろいろ前処理で泣きそうになりながらやってきた証をここに記す（お墓みたいだな・・・） 前処理X本ノックとあるがこのXの数は随時増える予定。

前処理大全[データ分析のためのSQL/R/Python実践テクニック]

前処理大全[データ分析のためのSQL/R/Python実践テクニック]

bash

1. csvの連結

日毎にCSVが出力されるが、DBにロードする際、ETL等で１件ずつロードがめんどくさいよね・・・ってときに、 CSVを一気に連結する方法。 すべて同じ列構成で、ヘッダ行がある前提で、ヘッダ行は残す。

$ # カレントディレクトリ内の拡張子がcsvのものをすべて連結する。
$ awk 'NR==1 || FNR!=1' *.csv >> concat.csv

2. テキストファイルの分割

行数で分割したい場合（リモートのDBにロード等で、分割しないとタイムアウトしたりする場合） ただしヘッダ行は考慮されない。

$ # 10000行毎に分割。分割あとはhoge-xx.csvになる。
$ split -l 10000 hoge.csv -d  --additional-suffix=.csv  hoge-

3. csvの不要列を削除

データ容量が大きく、不要列を削除して容量を削減したい場合に使用。 必要な列数目を指定して、残す処理。

$ # 4-6列目,8-19列目を残す
$ cat hoge.csv | cut -d "," -f 4-6,8-19 >> extract.csv

4. csvの不要行を削除

データ容量が大きく、不要列を削除して容量を削減したい場合に使用。 必要な行の条件を指定して、残す処理

$ #2列目がfugaのレコードだけを残す
$ cat hoge.csv | awk -F, '$2 == fuga {print $0}' > done.csv

5. tsvをcsvに変換

水平タブをカンマに置換

$ sed -E 's/\t/,/g' hoge.tsv >> hoge.csv

6. csvをtsvに変換

カンマを水平タブに置換

$ sed -E 's/,/\t/g' hoge.csv >> hoge.tsv

7. 固定長で改行を入れる

レガシーシステムでしか見ないけれども、

aaabbbccc

を

aaa
bbb
ccc

のように固定長で改行をいれたい場合。(csv化のための前処理とか)

$ # 3文字毎の折返し
$ fold -w 3 test.txt

バイト指定の場合はオプションが-b

$ # 3byte折返し
# fold -b 3 test.txt

8. 空白埋めされたcsvの空白除去

こんなのを

aaa   ,bbb    ,ccc    ,ddd  
eee   ,ff     ,g      ,h

こうしたい場合

aaa,bbb,ccc,ddd
eee,ff,g,h
$ tr -d " " < hoge.csv >> hoge_del_space.csv

ただしこのままだとyyyy/mm/dd hh:mm:ssみたいな形の値が入っていた場合も、 yyyy/mm/ddhh:mm:ssとdとhの間も埋められてしまう。

そのときはsedでの置換をする必要がある。 ダブルクオーテーションでくくってあれば下記でできる。

sed -E 's/\s+"/"/g' hoge.csv >> hoge_del_space.csv

9. 文字コード変換

sjisをutf8に変換

$ iconv -f SJIS -t UTF-8

[改訂第3版]シェルスクリプト基本リファレンス ──#!/bin/shで、ここまでできる (WEB+DB PRESS plus)

[改訂第3版]シェルスクリプト基本リファレンス ──#!/bin/shで、ここまでできる (WEB+DB PRESS plus)

SQL

1. 履歴型のマスタで最新値だけを取ってくる

履歴型のマスタ（という言葉がすでに崩壊しているが）等でよくあるパターン

登録日 商品コード 価格
2018/01/01 A 100
2018/10/01 A 120
2018/12/01 A 90
2018/11/01 B 120

みたいなのがあったとき、

登録日 商品コード 価格
2018/12/01 A 90
2018/11/01 B 120

こんな感じで最新のデータだけ取り出したい場合の処理。

-- データ作成
drop table if exists 商品マスタ;
create table 商品マスタ as
with t1 ( 
select to_date('20180101','YYYYMMDD') "登録日",'A' "商品コード",100 "価格" union all
select to_date('20181001','YYYYMMDD') "登録日",'A' "商品コード",120 "価格" union all
select to_date('20181201','YYYYMMDD') "登録日",'A' "商品コード",90 "価格" union all
select to_date('20181101','YYYYMMDD') "登録日",'B' "商品コード",120 "価格"
)
select * from t1
;

-- クエリ
select distinct
          first_value(登録日) over (partition by 商品コード order by 登録日 desc nulls last) "登録日"
         ,商品コード "商品コード"
         ,first_value(価格) over (partition by 商品コード order by 登録日 desc nulls last) "価格"
from
          商品マスタ
;

商品コード毎に新しい値で全部塗り替えてdistinctで重複を除外している。

2. 重複IDがあるレコード/単一IDのレコードだけを抽出する

ID 名前
1 hoge
1 fuga
2 chome
3 hogehoge

みたいなテーブルがあったとき、 IDがダブる1のレコードだけを取りたいとき。 with句を使うと楽（？）

-- データ作成
drop table if exists name_master;
create table name_master as 
with t1 (
select 1 "ID", 'hoge' "名前" union all
select 1 "ID", 'fuga' "名前" union all
select 2 "ID", 'chome' "名前" union all
select 3 "ID", 'hogehoge' "名前"
)
select * from t1
;
-- クエリ
with t1 as (
select 
    ID
from
    name_master
group by ID
having count(1) > 1
)
select
    a.*
from
     name_master a
     inner join t1 b on a.ID = b.ID
;

これで

ID 名前
1 hoge
1 fuga

が取り出せる。

逆にhaving count=1にすると単一レコードだけを取り出せる。

drop table if exists name_master;
create table name_master as 
with t1 (
select 1 "ID", 'hoge' "名前" union all
select 1 "ID", 'fuga' "名前" union all
select 2 "ID", 'chome' "名前" union all
select 3 "ID", 'hogehoge' "名前"
)
select * from t1;

with t1 as (
select distinct
    ID
from
    name_master
group by ID
having count(1) = 1
)
select
    a.*
from
     name_master a
     inner join t1 b on a.ID = b.ID
;

結果は下記になる。

ID 名前
2 chome
3 hogehoge

3. from-toで持っているテーブルの時間を展開

よくあるログデータの持ち方で、

user_id from_date to_date
A 2018/12/01 2018/12/03
B 2018/12/02 2018/12/05

のようなデータがあったとき、 （あるツールをAさんは12/1-3の間使いました、Bさんは2018/12/2-5の間使いました、のような意味）

user_id from_date
A 2018/12/01
A 2018/12/02
A 2018/12/03
B 2018/12/02
B 2018/12/03
B 2018/12/04
B 2018/12/05

のように加工したい場合は下記のようなSQLで対応可能。 (postgresqlの場合)

-- データ作成
drop table if exists user_log;
create table user_log as
select 'A' "user_id",to_date('20181201','YYYYMMDD') "from_date" ,to_date('20181203','YYYYMMDD') "to_date" union all
select 'B' "user_id",to_date('20181202','YYYYMMDD') "from_date" ,to_date('20181205','YYYYMMDD') "to_date"
;
-- クエリ
with RECURSIVE rt(user_id,from_date,to_date) as(
select user_id,from_date,to_date
  from user_log
union all
select user_id,from_date+1,to_date
  from rt
 where from_date+1 <= to_date)
select user_id,from_date "valid_date" from rt
order by user_id,from_date;

試してはいないが、MSSQL(sqlserver)の場合は上記のSQLからrecuresive(~~)を除外すればできる（はず）。 簡易的な解説をすると、 再帰SQLを使い、条件を満たしている間（from_dateがto_dateを超えるまで)は、 インクリメントさせ続ける。

ただし、この処理は重いので、大きいテーブルに対してはあまりやるべきではない。

↓のようにpythonでもできるがこちらも重い。 kazuhitogo.hateblo.jp

毎回型の確認等が発生するため遅い。 ので、コンパイル言語でやったほうがよい。 あるいはjuliaで関数作るか。

julia,コンパイル言語 >> python > SQL である。

4. pivot（縦→横）

よくあるあれです。

今回はこんな例を用意

店 商品 値段
A りんご 100
A みかん 50
A バナナ 200
B りんご 110
B みかん 40
B バナナ 100

これを

店 りんご みかん バナナ
B 110 40 100
A 100 50 200

こうしたい。

postgresqlだとcrosstab使えだとか、 oracle,mssqlだとpivot使えだとか、 いろいろ書き方はありますが、 どれもこれも商品が増えるとめんどくさいので、最初から動的SQLで。 （今回の例はpostgresqlだが、oracleもMSSQLも似たような書き方はできる。

-- データ作成
drop table if exists 値段マスタ;
create table 値段マスタ as
with t1 as (
select 'A' "店",'りんご' "商品", 100 "値段" union all
select 'A' "店",'みかん' "商品", 50 "値段" union all
select 'A' "店",'バナナ' "商品", 200 "値段" union all
select 'B' "店",'りんご' "商品", 110 "値段" union all
select 'B' "店",'みかん' "商品", 40 "値段" union all
select 'B' "店",'バナナ' "商品", 100 "値段"
)
select * from t1
;

-- 横持ちデータ作成
drop table if exists 値段マスタ横持ち;
DO $$
DECLARE
  vSQL text;
  CR RECORD;
BEGIN
  vSQL := 'create table 値段マスタ横持ち as select 店 ';
  FOR CR IN (SELECT distinct "商品" from "値段マスタ") LOOP
    vSQL := vSQL || ',max(case when "商品" = ''';
    vSQL := vSQL || CR.商品;
    vSQL := vSQL || '''';
    vSQL := vSQL || 'then 値段 else null end) ';
    vSQL := vSQL || CR.商品;
  END LOOP;
  vSQL := vSQL || ' from 値段マスタ group by 店';
  RAISE NOTICE '%', vSQL;
  EXECUTE vSQL;
END$$;

-- selectする
select * from 値段マスタ横持ち;

解説すると、 発行したいSQLは下記の通り。

create table 値段マスタ横持ち as 
select 
          店 
         ,max(case when "商品" = 'りんご'then 値段 else null end) りんご
         ,max(case when "商品" = 'みかん'then 値段 else null end) みかん
         ,max(case when "商品" = 'バナナ'then 値段 else null end) バナナ 
from 
          値段マスタ 
group by 
          店

このSQLをvSQL変数で作り、最後にEXECUTE vSQLで実行している。 商品が増えたときにも対応できるよう、max(~~)のところを、

SELECT distinct "商品" from "値段マスタ"

でカーソルに収め、自動で発行している。

ただし、今回はりんご、みかん、バナナくらいなので、 手書きやキーマクロ等でSQLを生成したほうが速いのは言うまでも無い。

あくまで何が入るか予測できなかったり、 列の数が巨大（といっても列数の上限はoracleが1000,MSSQLが1024,postgresが1600しか無いが)のときに、 このような動的SQLを発行するとよい。

5. unpivot(横→縦)

pivotの逆。

これを

店 りんご みかん バナナ
B 110 40 100
A 100 50 200

こうしたい。

店 商品 値段
A りんご 100
A みかん 50
A バナナ 200
B りんご 110
B みかん 40
B バナナ 100 
-- データ作成
drop table if exists 値段マスタ横持ち;
create table 値段マスタ横持ち as
with t1 as (
select 'A' "店",100 "りんご",50 "みかん", 200 "バナナ" union all
select 'B' "店",110 "りんご",40 "みかん", 100 "バナナ"
)
select * from t1

-- テーブル作成
drop table if exists 値段マスタ縦持ち;
DO $$
DECLARE
  vSQL text;
  CR RECORD;
BEGIN
  vSQL := 'create table 値段マスタ縦持ち as ';
  FOR CR IN (SELECT attname FROM pg_attribute WHERE attrelid = '値段マスタ横持ち'::regclass and attstattarget = -1 and attname != '店') LOOP
      vSQL := vSQL || 'select 店, ''';
      vSQL := vSQL || CR.attname;
      vSQL := vSQL || ''' 商品 ,';
      vSQL := vSQL || CR.attname;
      vSQL := vSQL || ' from 値段マスタ横持ち union all ';
  END LOOP;
  vSQL := left(vSQL,-11);
  RAISE NOTICE '%', vSQL;
  EXECUTE vSQL;
END $$;

-- クエリ

select * from 値段マスタ縦持ち order by 店,商品;

これもpivotと同様、発行したいSQLがあり、それをカーソルループで作る。

発行したいSQLは下記の通り。

create table 値段マスタ縦持ち as 
select 
          店
         , 'りんご' 商品 
         ,りんご 
from 
           値段マスタ横持ち 
union all 
select 
          店
         , 'みかん' 商品 
         ,みかん 
from 
           値段マスタ横持ち 
union all 
select 
          店
         , 'バナナ' 商品 
         ,バナナ 
from 
          値段マスタ横持ち

これを作成するために、動的SQLで、列リストが必要。 列リストは下記SQLで取れる。（カーソル）

SELECT 
          attname 
FROM 
          pg_attribute 
WHERE 
          attrelid = '値段マスタ横持ち'::regclass 
      and attstattarget = -1 
      and attname != '店'

これはpostgresにしかできないが、 pg_attributeテーブルのatterelidにテーブルを指定し、 attstattarget = -1の条件をつけることで、 該当テーブルの列リストを作ることができる。

ただし、今回のところで言う、 "店"カラムは動的SQLのパラメータには不要なため、 これもまたwhere句で除外する。

postgresqlにしかできない、と書いたが、 OracleはALL_TAB_COLUMNSから、MSSQLはsys.columnsから列一覧が取れるので、 似たようなSQLは記述できる。

ただし、今回はりんご、みかん、バナナくらいなので、 手書きやキーマクロ等でSQLを生成したほうが速いのは言うまでも無い。

あくまで何が入るか予測できなかったり、 列の数が巨大（といっても列数の上限はoracleが1000,MSSQLが1024,postgresが1600しか無いが)のときに、 このような動的SQLを発行するとよい。

6. 過去データと新データを新データ優先してマージ

下記のような２つのデータがあったとする。

過去データ

登録日 売上
2018/12/01 1000
2018/12/02 1001
2018/12/03 1002
2018/12/04 1003
2018/12/05 1004

新データ

登録日 売上
2018/12/04 1010
2018/12/05 1011
2018/12/06 2000

ここから新しいデータを優先して(重複する場合は新データを採用)、 縦にくっつけたい(union)場合。

求める結果は下記。

登録日 売上
2018/12/01 1000
2018/12/02 1001
2018/12/03 1002
2018/12/04 1010
2018/12/05 1011
2018/12/06 2000

sqlは下記で行ける。

with 過去データ as (
select to_date('20181201','YYYYMMDD') "登録日",1000 売上 union all
select to_date('20181202','YYYYMMDD') "登録日",1001 売上 union all
select to_date('20181203','YYYYMMDD') "登録日",1002 売上 union all
select to_date('20181204','YYYYMMDD') "登録日",1003 売上 union all
select to_date('20181205','YYYYMMDD') "登録日",1004 売上
),新データ as (
select to_date('20181204','YYYYMMDD') "登録日",1010 売上 union all
select to_date('20181205','YYYYMMDD') "登録日",1011 売上 union all
select to_date('20181206','YYYYMMDD') "登録日",2000 売上
)
select   
         *
from 
         新データ
union all
select
         *
from
         過去データ a
where
         not exists(select 1 from 新データ b where a.登録日 = b.登録日 )
order by 登録日

まずは新データを全部持ってきて、新データと過去データを比較して、無いものだけを過去データから持ってくる。 (not existsの部分)

ただし、not existsの処理は重いので注意。（まぁ仕方なしの場合が多いけれども）

7. 多数決でデータを作る

なんのこっちゃ、かと思いますが、例えば下記。

商品 原価
A 100
A 100
A 120
B 50
B 50
B 70

こんなデータがトランザクションで入ってきたとして、 でも原価って一律決まってて誰かの入力ミスだよね…。

みたいなときに、多数決で決めれば多分合ってんだろ、という場合、 下記のようなデータをどうやって抜き出すか、という問題。

商品 原価
A 100
B 50

実装はこれでいける。

with 原価データ as (
select 'A' "商品",100 "原価" union all
select 'A' "商品",100 "原価" union all
select 'A' "商品",120 "原価" union all
select 'B' "商品",50 "原価" union all
select 'B' "商品",50 "原価" union all
select 'B' "商品",70 "原価"
),t1 as (
select
          "商品"
         ,"原価"
         ,count(1) "件数"
from
         原価データ
group by 
          "商品"
         ,"原価"
)
select distinct
         "商品"
         ,first_value("原価") over (partition by "商品" order by "件数" desc) "原価" 
from
         t1

その原価の件数を取って多いものソートし、それをfirst_valueで塗り替える。

(本来こんなデータがあっちゃいけないんですがそれはさておき)

8. データが無い日を補完する

こんな情報があったとする。

日付 売上
2018/01/01 1000
2018/01/02 2000
2018/01/03 3000
2018/01/05 5000
2018/01/07 7000

4日や6日は売上がなくてトランザクションが存在しない場合NULLを入れたくなる。 まともなシステムならカレンダマスタ的なのが存在するが、 無い場合どうしよ、ってエントリ。

無いものは作ったれ…とカレンダマスタを作る。 ただし、不要なカレンダは作りたくないので、最小と最大部分のマスタを作る実装。

前述の再帰withでもできるけれども、 postgresqlにはgenerate_seriesがあるので、それを使ったもの。

generate_seriesはざっくりいうと、 日付に限らず連続するものを最小と最大とインターバルを指定すると、 その値の集合を返してくれる。

https://www.postgresql.jp/document/9.1/html/functions-srf.html

ただしそのままだと固定値を都度入れることになるので、 pl/pgsqlで可変で使えるようにする。

-- データ作成

drop table if exists 売上情報;
create table 売上情報 as 
with t1 as (
select to_date('20180101','YYYYMMDD') "日付",1000 "売上" union all
select to_date('20180102','YYYYMMDD') "日付",2000 "売上" union all
select to_date('20180103','YYYYMMDD') "日付",3000 "売上" union all
select to_date('20180105','YYYYMMDD') "日付",5000 "売上" union all
select to_date('20180107','YYYYMMDD') "日付",7000 "売上"
)
select * from t1
;

-- カレンダ作成

DO $$
DECLARE
  min_date date;
  max_date date;
BEGIN
  select into min_date min("日付") from 売上情報;
  select into max_date max("日付") from 売上情報;
  RAISE NOTICE '%', min_date;
  RAISE NOTICE '%', max_date;
  drop table if exists カレンダマスタ;
  create table カレンダマスタ as
  SELECT
            generate_series "日付"
  FROM 
            generate_series(min_date,max_date, '1 days')
  ;
END
$$
;

-- カレンダと売上の結合

select
         a.日付
         ,b.売上
from
         カレンダマスタ a
         left join 売上情報 b on a.日付 = b.日付

結果は下記が返る。

日付 売上
2018/01/01 1000
2018/01/02 2000
2018/01/03 3000
2018/01/04 NULL
2018/01/05 5000
2018/01/06 NULL
2018/01/07 7000

oracleもsqlserverも似たような方法はある。

oracle doruby.jp

sqlserver 夏椰の庵 - Secluded Spot of Kaya -

スッキリわかる SQL 入門 ドリル215問付き! (スッキリシリーズ)

スッキリわかる SQL 入門 ドリル215問付き! (スッキリシリーズ)

python

1. リストと辞書の組み合わせからcsvに変換

下記のような辞書のキーに予測がつかないようなリストがあったとする。

x = [{"hoge":1,"fuga":2},{"hoge":3,"fuga":4},{"hoge":5,"hogehoge":6}]

これをキー一覧から表形式に変換する。

fuga hoge hogehoge 　
2 1 NaN
4 3 NaN
NaN 5 6 
# pandasを使う
import pandas as pd
# データ作成
x = [{"hoge":1,"fuga":2},{"hoge":3,"fuga":4},{"hoge":5,"hogehoge":6}]
# キー（カラム）一覧を作成
key_set = set([ii for i in x for ii in i.keys()])
# indexのみのdataframe作成
csv_df = pd.DataFrame(index=[i for i in range(len(x))], columns=[])
# 各カラムのデータを入れていく
for ii in key_set:
    csv_df[ii] = pd.DataFrame([i.get(ii) for i in x])
# csv出力
csv_df.to_csv("hoge.csv",index = False)

リスト内包表記を連発しているが、簡単に解説すると下記の通り。

  1. キーをまずは重複ありでもいいのでとにかくとってきて、それをset関数で重複を排除する
  2. 空のデータフレームを作成する
  3. キー順にデータを入れていく。

もう少し速い処理もかけそうな気もするけれども一旦。

入門 Python 3

入門 Python 3

コンサル寄りのデータサイエンティスト（自称）によるmac環境整備

mac macbookair データサイエンティスト

調教…！調教…！

データサイエンティストならずとも、 プログラムを書く人にとってPCは神聖なるツールで、 買ったばかりのものは調教し無くてはならないかと思います。 (神聖なものは調教しないといけないという価値観）

私はvimやemacsの宗教戦争に巻き込まれたくはないのですが、 普段winで使っている相当の環境を準備くらいはしないと使い物にならないので、 導入するもの、手順を残しました。

普段の開発環境がwinで、久しぶりにmacを使った人の備忘録として…。 （もっとこんな手順あるよ！とか教えてください）

設定

  • macの３本指スワイプを有効にする これで有効化しました。 pc-karuma.net

  • ダークモード mojaveからダークモードが使えるようになったそうで、消費電力的にも目にも優しいそうなので。 初回起動時に設定もできますが、再設定方法はこちら。 http://ascii.jp/elem/000/001/747/1747977/ascii.jp

入れたもの

  • google chrome

    safariなんぞクソくらえ…！ インストール方法も特に難しくないので割愛

  • google日本語入力

    基本的にグーグル先生に身も心も売ってます。

  • visual studio code

    特にエディタにこだわりはないですし、 普段Winではサクラエディタ使ったり、IDE使ったりしてますけど、 最近のvscodeはよいので…。

  • terminal起動ショートカット

    ターミナルを使うことは多いので、下記リンクを参考にコマンド+option+Tで起動できるようにしました。 ry0.github.io

  • xcode-select

    入れておかないとhome-brewが使えないとのことなので一応。 下記コマンドでインストール shell /usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"

  • git

    自分のバージョン管理はもちろんのこと、アプリもgithubで配布されることも多いので入れます。

    shell brew install git

  • PlantUML

    テキストでUMLは管理したいよね。表示に難はあるけど。

qiita.com

  • python絡み anacondaを普段使っているので、macでも…と思うのですが、 condaコマンドを有効にするのにスマートなやり方がなかったので、下記手順で実施。 pyenvいれて、anaconda入れて、pathとおして…的な。

    shell brew install pyenv pyenv install -l | grep anaconda pyenv install anaconda3-5.3.0 echo 'export PYENV_ROOT="${HOME}/.pyenv"' >> ~/.bash_profile echo 'export PATH="${PYENV_ROOT}/bin:$PATH"' >> ~/.bash_profile echo 'eval "$(pyenv init -)"' >> ~/.bash_profile echo ". ~/.pyenv/versions/anaconda3-5.3.0/etc/profile.d/conda.sh" >> ~/.bash_profile source ~/.bash_profile pyenv global anaconda3-5.3.0

  • postgres 最近postgresがお気に入り（oracleとかmssqlとか,mysqlhはライセンスが逐一めんどくて）。 mariadbはなぜか遅くて…。 brewで入れてセットアップ。

$ # postgresqlをインストール
$ brew install postgresql
$ # utf-8でDBを初期化(日本語圏ではロケールの意味がないらしいので--no-localeオプションを入れておく
$ initdb /usr/local/var/postgres -E utf8 --no-locale
$ # インストールの確認。DB一覧が表示されればOK。
$ psql -l
$ # サービスを起動
$ brew services start postgresql
$ #データディレクトリの環境変数を設定(.bash_profileに記載）
$ echo "export PGDATA=/usr/local/var/postgres" >> ~/.bash_profile

あとは何かと便利なのでGUIツール(pgcommander)も導入しておく

PG Commander, a PostgreSQL client for Mac

あとは起動して、nicknameは適当に、hostをlocalhost,portは5432,userはOSのユーザ名, passwordなし、databaseはpostgresで取りあえず使える。

随時更新中…

０から作るOS開発　Vol.0　環境準備編　OSを作る環境の準備と設定

０から作るOS開発　Vol.0　環境準備編　OSを作る環境の準備と設定

SBI証券口座の資産残高を自動で取得する（python & selenium)〜selenium基礎編〜

python selenium SBI証券

資産がいくらあるのかを毎回確認するのはめんどくさい

みなさん、自分の口座においくら入ってますか？？ …と問われたときにパッと答えられる人はなかなかいませんし、私も５桁円以下は全く覚えてません。

というか外貨もあるので、レートとか考えなきゃいけません。

めんどいので自動化しましょう。 というのが、今回のエントリ。

Selenium

Webサイトの回帰テストでよく用いられる、Seleniumを使いました。

Selenium実践入門 ―― 自動化による継続的なブラウザテスト (WEB+DB PRESS plus)

Selenium実践入門 ―― 自動化による継続的なブラウザテスト (WEB+DB PRESS plus)

実践 Selenium WebDriver

実践 Selenium WebDriver

Seleniumとはなんぞや？ってのはググればいくらでも出てくるので割愛ですが、 簡単に言うと人間の操作を自動でやってくれるもので、 スクショとかも勝手に撮らせることもできます。

環境構築

以前はWinマシンを使ってたのですが、 ついこないだmacbook airを買ったばっかりなので、 → kazuhitogo.hateblo.jp

環境を整えるところから…。

pythonとanacondaを使える前提から。

下記をmacのterminalから実行してseleniumのライブラリとchromedriver(chromeをseleniumで動かすためのアプリ)をインストールする。 （一応仮想環境も別途作っておく

# selenium用仮想環境作成
conda create -n selenium python=3.7 anaconda
# selenium用仮想環境起動
conda activate selenium
# seleniumインストール
conda install selenium
# chromedriverインストール
brew cask install chromedriver

…簡単ですね。

さて、ここからはpythonのコード。

# seleniumのWebdriverインポート
from selenium import webdriver

# chrome driverインスタンス生成
driver = webdriver.Chrome()

# URL遷移
driver.get('https://site2.sbisec.co.jp/ETGate/')

これだけでブラウザが起動し、 SBI証券のWebサイトが開きます。

f:id:kazuhitogo:20181125182535p:plain

ここから ユーザーネームとパスワードを入力して、 ログインボタンを押下したいわけです。

そのためには、人間の操作と同じ通り、下記ステップで行います。

  1. ユーザーネームテキストボックスを選択
  2. ユーザーネームを入力
  3. パスワードテキストボックスを選択
  4. パスワードを入力
  5. ログインボタンを押下

幸いSBI証券はユーザーネーム、パスワード、ログインボタン全てにHTMLでIDを振ってくれているので、 IDをたどって選択して、入力＆押下しましょう。

IDを知らない人はここらへんを確認ください。

HTMLタグ/共通属性/要素にIDを付ける - TAG index

# ユーザーネームとパスワードを設定
username = "hoge"
password = "hogepass"

# ユーザーネームテキストボックスを選択
user_id = driver.find_element_by_id("user_input")
# ユーザーネームを入力
user_id.find_element_by_name("user_id").send_keys(username)
# パスワードテキストボックスを選択
user_pass = driver.find_element_by_id("password_input")
# パスワードを入力
user_pass.find_element_by_name("user_password").send_keys(password)
# ログインを押下
login = driver.find_element_by_class_name("ov").click()

これでユーザーネームにhoge,パスワードにhogepassを入力してくれます。 ソースコードの中にユーザーネームやパスワードを直書きするのに抵抗があるかたは、 例えばこのソースコードをSBI.pyとしたときに、

python SBI.py hoge hogepass

のようにpython実行時の引数として渡すことも可能です。 そのときはこのように修正するといけます。

# コマンドライン引数を受け取る
import sys
args = sys.argv

# username直打ちをやめる
# username = "hoge"
username = args[1]
# password直打ちをやめる
# password = "hogepass"
password = args[2]

さて、これでログインできたはずです。 私は円建てとドル建ての資産があるので、 それぞれチェックの必要があります。

まずは右上にある「口座管理」のリンクに飛びましょう。

こちらはIDが無いので、xpathでアクセスしましょう xpathを知らない人はこちらを確認ください。

XPathのまとめ、要素の参照方法いろいろ │ Web備忘録

ちなみにxpathを調べるの、ブラウザの開発者ツールとかを使うのもよいのですが、 chromeの拡張機能である、「xPath Finder」を使うと便利でした。

xPath Finder - Chrome Web Store

さて、「口座管理」にxpathを使って飛びます。

# 口座管理に遷移
driver.find_element_by_xpath("/html/body/div[1]/div[1]/div[2]/div/ul/li[3]/a/img").click()

これを打つことで飛ぶことができます。

さて、口座管理画面を開くと、もう円建ての資産が表示されています。 （スクショは個人情報満載なので割愛）

ここからxpathを使って、総資産額をとってきましょう。

# 円建ての資産総額を取得
yen_all_asset = float(driver.find_element_by_xpath("/html/body/div[1]/table/tbody/tr/td[1]/table/tbody/tr[2]/td/table[1]/tbody/tr/td/form/table[2]/tbody/tr[1]/td[2]/table[5]/tbody/tr/td[1]/table[3]/tbody/tr[8]/td[2]/div/b").text.replace(',',''))
# 円建ての総額を表示
print(yen_all_asset)

これで円建ての資産がコンソールに表示されます。

xpathで総額が表示されているパスまで行くことができ、 textという変数に要素の中身が格納されているので、それを取得します。

ただし、こちらは文字列形式で格納されているうえ、 カンマがあってそのままじゃ数値型にキャストできないため、 replaceメソッドでカンマを除去した上、 float型にキャストしています。

さて、同様にドル建ての資産に行きましょう。

「口座（外貨建）」と書かれたリンクに飛べばよいので、 同様にxpathで遷移し、そこに書かれた総額を取得します。

# 口座（外貨建）に遷移
driver.find_element_by_xpath("/html/body/div[1]/div[3]/div/div/ul/li[2]/div").click()

# ドル建ての資産総額を取得
doller_all_asset = float(driver.find_element_by_xpath("/html/body/div[1]/table[2]/tbody/tr/td[1]/table/tbody/tr/td[2]/table[2]/tbody/tr[1]/td[2]/table[5]/tbody/tr/td[1]/table[3]/tbody/tr[2]/td[3]/table/tbody/tr[2]/td[2]/b").text.replace(',',''))
# ドル建ての総額を表示
print(doller_all_asset)

これでドル建ての資産が表示されます。

最後にchromedriverを閉じましょう。 （閉じないとログインしっぱなしの画面が残ることに） また、合計額も併せて表示しておきましょう。

# 閉じる
driver.close()

# summary
print("円建資産： " + str(yen_all_asset) + "円")
print("ドル資産建： " + str(doller_all_asset) + "円")
print("総額：" + str(yen_all_asset + doller_all_asset)+ "円")

これで完了です。

一応最後にまとめておきます。

from selenium import webdriver
import sys

args = sys.argv

#ログインID、パスワードはコマンドライン引数から
username = args[1]
password = args[2]

# chrome driverインスタンス生成
driver = webdriver.Chrome()

# URL遷移
driver.get('https://site2.sbisec.co.jp/ETGate/')

# user_idを入力する場所探索
user_id = driver.find_element_by_id("user_input")
user_id.find_element_by_name("user_id").send_keys(username)
user_pass = driver.find_element_by_id("password_input")
user_pass.find_element_by_name("user_password").send_keys(password)
login = driver.find_element_by_class_name("ov").click()

# 口座管理に遷移
driver.find_element_by_xpath("/html/body/div[1]/div[1]/div[2]/div/ul/li[3]/a/img").click()

# 円
yen_all_asset = float(driver.find_element_by_xpath("/html/body/div[1]/table/tbody/tr/td[1]/table/tbody/tr[2]/td/table[1]/tbody/tr/td/form/table[2]/tbody/tr[1]/td[2]/table[5]/tbody/tr/td[1]/table[3]/tbody/tr[8]/td[2]/div/b").text.replace(',',''))

# 口座（外貨建）に遷移
driver.find_element_by_xpath("/html/body/div[1]/div[3]/div/div/ul/li[2]/div").click()

# ドル建ての資産総額を取得
doller_all_asset = float(driver.find_element_by_xpath("/html/body/div[1]/table[2]/tbody/tr/td[1]/table/tbody/tr/td[2]/table[2]/tbody/tr[1]/td[2]/table[5]/tbody/tr/td[1]/table[3]/tbody/tr[2]/td[3]/table/tbody/tr[2]/td[2]/b").text.replace(',',''))

driver.close()

print("円建資産： " + str(yen_all_asset) + "円")
print("ドル資産建： " + str(doller_all_asset) + "円")
print("総額：" + str(yen_all_asset + doller_all_asset)+ "円")

あとは、この結果をメールとかで飛ばしてもいいですね。 ↓のプログラムを多少改変すればできます。

kazuhitogo.hateblo.jp

はじめての株1年生 新・儲かるしくみ損する理由がわかる本 (アスカビジネス)

はじめての株1年生 新・儲かるしくみ損する理由がわかる本 (アスカビジネス)

  • 作者: 竹内弘樹
  • 出版社/メーカー: 明日香出版社
  • 発売日: 2009/12/22
  • メディア: 単行本（ソフトカバー）
  • 購入: 2人 クリック: 18回
  • この商品を含むブログを見る

ざっつおーるせんきゅー

macbook air(2018)買いました〜17年ぶりのmac開封の儀でおふざけ〜

…！

f:id:kazuhitogo:20181121083050j:plain

ズズズ…！

f:id:kazuhitogo:20181121083102j:plain

バーンッ！

f:id:kazuhitogo:20181121083112j:plain

このマックではありません。 いや、買ったものですが。 （美味しくいただきました。）

買ったのはマックのホットアップルパイでもなく、 ラズベリーパイでもなく、

Raspberry Pi3 Model B ボード＆ケースセット 3ple Decker対応 (Clear)-Physical Computing Lab

Raspberry Pi3 Model B ボード＆ケースセット 3ple Decker対応 (Clear)-Physical Computing Lab

apple社製のmacbook airを買ってしまいました。

総額193104円…！ ストレージ256GB,メモリ16GB…！

アップル 13．3インチ MacBook Air(1．8GHz Dual Core i5 / 8GB / 128GB) MQD32J/A

アップル 13．3インチ MacBook Air(1．8GHz Dual Core i5 / 8GB / 128GB) MQD32J/A

アップルの公式サイトから買うと定価でしか買えないので、 ポイントが5%つくビックカメラで買うかー、なんて思ってました。

だが…だが…

「カスタマイズ（今回で言うメモリ増設）をする場合はポイント付かない上、 納期がアップル公式サイトより遅れます^^本当にメモリ16GBいるの？ お前パソコンわかってんの？どうせネットサーフィンしかしないでしょ？」

と店員に遠回しに煽り口調で言われたので 腹たち紛れにビックカメラのアップルブースで、 アップルの公式サイトで速攻ポチりました。

メモリ8GBなんて、 ちょっとDBインスタンス建てて、 DBからデータをロードして、 ちょっと処理しただけで超えるわ！

あるいは、BMPとか扱ってopencvとかで、加工した瞬間すぐ超えるわ！ と心の中でイキリながら…。

注文してから待つこと一週間、 深センから送られてきたので開封の儀の様子をお届けします。

まず箱。

f:id:kazuhitogo:20181121083123j:plain

写真からは伝わらないのですが、なんか臭いっす。 １週間位風呂入っていない人の匂いがした気がします。

梱包用の箱はとっととお払い箱。 開けます。

f:id:kazuhitogo:20181121083136j:plain

お、こころ踊る。

f:id:kazuhitogo:20181121083147j:plain

f:id:kazuhitogo:20181121083157j:plain

あ・・・れ・・・？ ゴールド買ったつもりがスペースグレー…？

伝票見るとたしかにスペースグレー…。 きっと色選択間違いだな…。

まぁ変に尖ったゴールドよりスペースグレーの落ち着いた色のほうがよい、と、 脳内心理で思ってスペースグレーを無意識に買ったのでしょう。

f:id:kazuhitogo:20181121083211j:plain うーん、リサイクルアルミらしいけど、違いはわからない。 強いて言えば、気温が落ちてきてアルミはちょっと冷たいな、と小学生並みの感想を。

f:id:kazuhitogo:20181121083222j:plain お、これがtype-cの充電ケーブルか。 スマホ以外では初めて。 私はandroiderで比較的しょっちゅう買い換えるので、type-cケーブルはたくさん持ってるんですけどね。（自慢） なら、なぜmac買ったか、ってのはさておき。

f:id:kazuhitogo:20181121083232j:plain

充電器も出てきました。これでandroid充電したら早いのかなぁ？（駄目、ゼッタイ）

さて、開封＆接続。

f:id:kazuhitogo:20181121083243j:plain

おいぃぃぃ！ 接続した瞬間起動しやがった！！！ｗｗｗ

電源ボタン押させてやー！

まぁいいや。

macbook air(2018)はtype-cを二本させるだけで、 拡張性が皆無という事前情報を知っていたので（特にすごくない）、 拡張性をもたせるべくドックを買いました。

この記事によると、 www.gizmodo.jp macbook pro用のドックが挿さる（ただし使えるとは言っていない）とのことで、 macbook pro用のを確認もせずに買いました。

あ、ポチっとな。

Wavlink USB C Hub Type C adapter USB Cハブ 13

Wavlink USB C Hub Type C adapter USB Cハブ 13"/15" MacBook Pro 2016/2017 thunderbolt 3 充電ポート PD機能付きサンダーボルト3 ハブ ドッキングステーション 4K HDMIポート 2* USB 3.0 ポートMicro SD/SDカードリーダー アルミニウム 7 in1 （シルバー） …

f:id:kazuhitogo:20181121083254j:plain

じゃん。

f:id:kazuhitogo:20181121083307j:plain

ささった…！

使えた…！

よく口コミで聞いていたような、 Wi-Fiが使えなくなる、とかそういった症状もなく、 HDMI端子もtype-cも使えました。

microSDやtype-aはこれから試します。

さて、このままではビックカメラの店員の想定していた使い方以下しかしていないので、 macbook air調教編をそのうち。

よくわからずmodel.compile(loss='binary_crossentropy',xxx...)している方むけのcross entropy入門

ニューラルネットワーク python

私です

すみません、出落ちです。

私のようなコピペプログラマにとって、 他人が書いたソースコードを適当に自分用に変数を変えて切った貼ったした上、 さも「私がやりました(ｷﾘｯ」と言う人は、 model.compile(以下略)を、 C言語を始めた人ががおまじないのように

#include <stdio.h>

int main(int argc, char *args[])
{
    printf("Hello, world!\n");
    return 0;
}

のように書いて、 includeで何？stdio.hって何？intって何？mainって何？argcって何？（以下略 わからないけどとりあえずHello,World!って出るね！

みたいな感じで、分類問題なら

model.compile(optimizer=tf.train.AdamOptimizer(),loss='binary_crossentropy',metrics=['accuracy')

って書いてmodel.fit(train_x,train_y…)すれば学習して、model.predict(x,y,…)で予測できるんでしょ！

ってやってる人は多いかと思います。（いや、私だけ？） Adam等はさておき、cross entropyについて、 理解したつもりになったのでこちらにて。

cross entropyって何？

さて、cross entropyって何？という話で、

  • 負の対数尤度関数 (Negative Loglikelihood Function）
  • 交差エントロピー誤差関数ともいう

らしいです。 なんのこっちゃ。

数式は以下の通りです。(binary cross entropyの場合)

\displaystyle E = -\sum_{i=1}^N \left( t_i \log y_i + (1 - t_i ) \log ( 1 - y_i )\right)

なんのこっちゃ。

cross entropyは分類問題を0～1の連続値で返してほしいときに使うので、二値問題で考えてみましょう。

二値分類問題なので、あるデータを与えられたとき、それが0なのか、1なのかを判別します。 (多値分類問題もone-hot化＆softmaxするので、各値が1(True)なのか0(False)なのかの程度を返すので同じですが)

さてもう一度、先程の数式。

\displaystyle E = -\sum_{i=1}^N \left( t_i \log y_i + (1 - t_i ) \log ( 1 - y_i )\right)

yが予測結果でtが真の値です。 予測する対象が一つ(オンライン処理)であればΣが不要なので、こうなります。

\displaystyle E = -\left( t \log y + (1 - t) \log ( 1 - y )\right)

さて、ここで二値分類なのでyは0～1を、tは0か1が入ってきます。 tが0の場合と、tが1の場合で、それぞれyを0～1へ推移させた場合、Eはどのような推移をしていくのかを見てみましょう。

import numpy as np
from matplotlib import pyplot as plt

# yの値を0(.00001)～1に推移
y = np.arange(0.00001,1,0.00001).reshape(-1,1)

# T=0をyの配列分用意
t_0 = np.zeros(y.shape[0]).reshape(-1,1)
# T=1をyの配列分用意
t_1 = np.ones(y.shape[0]).reshape(-1,1)

# np.logをそのまま使うとnp.log(0)のとき発散してしまうので、最小値を1e-100に固定
log = lambda x:np.log(np.clip(a=x, a_min=1e-100, a_max=x))

# 交差エントロピー誤差関数(公式通りsumしておく)
cross_entropy = lambda t,y: np.sum(-t*log(y)-(1-t)*log(1-y),axis=1)

# エントロピーの値をそれぞれ出す
E_0 = cross_entropy(t_0,y)
E_1 = cross_entropy(t_1,y)

# 描画
plt.figure()
plt.plot(y,E_0,c='red',label = "T=0")
plt.plot(y,E_1,c='blue',label = "T=1")
plt.xlabel("y")
plt.ylabel("E")

f:id:kazuhitogo:20181101150608p:plain

これを見ると、値が遠ければ遠いほどものすごく大きい値を返す(T-y=|1|のときは発散する)ことがわかります。 すなわち間違いが大きい場合は誤差Eも跳ね上がるし、近ければ誤差は少ないよ、ということです。

学習して学習して、そしてこの誤差を0に出来たとき人は幸せになれるわけです。

で、0にするにはwとbをいじって0になる値を探すわけですが….。

それを数学的に解くのは難しいので、

  1. 微分してその刹那wやbをどっちにどれくらいずらせば良い方向に行くのかを計算
  2. 実際にずらす

を繰り返して、 最適なwとbを探すのがニューラルネットの学習方法なわけですね。

なのでEをwとbでそれぞれ偏微分してどれくらいwとbをずらせばよいのか計算しましょう。 Eの式にwは含まれないので一旦yをはさみます。

\displaystyle \frac{dE}{dw} = \frac{dE}{dy}\frac{dy}{dw}

さて、\displaystyle E = -\left( t \log y + (1 - t) \log ( 1 - y )\right) でしたので、 対数関数の微分は対数の関数を分数に持っていけばよいので（高校数学ですね！）、

\displaystyle \frac{dE}{dy} = \frac{t}{y} - \frac{1-t}{1-y}

ですね。

それを先程の式に当てはめると、

\displaystyle \frac{dE}{dy}\frac{dy}{dw} = -(\frac{t}{y} - \frac{1-t}{1-y})\frac{dy}{dw}

となります。

さて、yをwで微分ですが、二値分類で確率の出力をしている前提の話ですので、 出力層の活性化関数はsigmoid関数を使っているはずです。

sigmoid関数は

\displaystyle y = \frac{1}{1+e^{-z}}

の通りなので、 zにwx + bを入れ、

\displaystyle y = \frac{1}{1+e^{-wx-b}}

となります。 …

…

…

微分の仕方忘れたので、ここからパクリます。 これをwで微分すると、

\displaystyle \frac{dy}{dw} = y(1-y)x となります。

さて、合わせるとこうなります。

\displaystyle \frac{dE}{dy}\frac{dy}{dw} = -(\frac{t}{y} - \frac{1-t}{1-y})(y(1-y)x)

texを書くの疲れてきたので、整理を省略すると、

\displaystyle \frac{dE}{dy}\frac{dy}{dw} = -(t-y)x

となります。

エントロピーを小さくするためにどれくらいずらしたのか、を計算した結果が、 - (真の値 - 予測の値) * xになってしまうんですね。 これをwに加算して、再度学習すればよいわけです。

同様のことをbについても行うと、

\displaystyle \frac{dE}{db}=-(t-y)

となり、こちらもほぼ同様な(真の値 - 予測の値)になります。 これをbに加算して（以下略

ちなみにこのような美しい数式になるのは、 誰かが楽に計算したいがために、このような数式を見つけ出したのでしょう。 …たぶん。

さて、これを持って学習を可視化してみましょう。

データが多すぎると大変なので、ANDゲートをロジスティック回帰で解く問題としてみましょう。 (((1,1)が1,(0.0),(0,1),(1,0)が0と返すようになればOK)

%matplotlib nbagg
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.animation as animation
rs = np.random.RandomState(4545)

# tは正解,yは予測結果
cross_entropy = lambda t,y: np.sum(-t*log(y)-(1-t)*log(1-y))

# np.logをそのまま使うとnp.log(0)のとき発散してしまうので、最小値を1e-100に固定
# cross_entropy関数で使う
log = lambda x:np.log(np.clip(a=x, a_min=1e-100, a_max=x))

# 活性化関数sigmoid
sigmoid = lambda x: 1.0 / (1.0 + np.exp(-x))

# 予測結果
predict = lambda w,b,x:sigmoid(np.matmul(w,x)+b)

# データ
train_x = np.array([[0,0],[0,1],[1,0],[1,1]],dtype=float)
train_y = np.array([0,0,0,1])

# w,bの初期化
w = np.array([1.0,1.0])
b = 1.0
# 学習率
e = 0.05
# 可視化用履歴変数
w_hist = []
b_hist = []
i_hist = []
y_hist = []
cost_hist = []
w_hist.append(w)
b_hist.append(b)
dw = np.zeros(2,dtype=float)
db = 0.0
# epoch数
epoch = 40
for j in range(epoch):
    # 学習順序ガチャ
    learn_order = rs.choice(np.arange(train_x.shape[0]),train_x.shape[0],replace=False)
    for i in learn_order:
        t = train_y[i]
        x = train_x[i]
        y = predict(w,b,x)
        cost = cross_entropy(t,y)
        delta = -(t-y)
        dw = e * delta * x
        db = e * delta
        w = w - dw
        b = b - db
        i_hist.append(i)
        w_hist.append(w)
        b_hist.append(b)
        y_hist.append(y)
        cost_hist.append(cost)

# wx + bの線履歴を計算
x0_hist = [-3,3]
x1_hist = np.zeros(len(w_hist)*2).reshape(len(w_hist),2)
for i in range(len(w_hist)):
    x1_hist[i,0] = -(w_hist[i][0]/w_hist[i][1])*x0_hist[0] - b_hist[i]/w_hist[i][1]
    x1_hist[i,1] = -(w_hist[i][0]/w_hist[i][1])*x0_hist[1] - b_hist[i]/w_hist[i][1]

def update_fig(i):
    print(i)
    plt.cla()
    plt.xlim(-1,2)
    plt.ylim(-1,2)
    plt.scatter(train_x[:,0],train_x[:,1])
    # 線を引く
    plt.plot(x0_hist,x1_hist[i],c='black')
    plt.scatter(train_x[i_hist[i],0],train_x[i_hist[i],1],c='red')
    # cost表示
    plt.text(-0.5, -0.5, str(cost_hist[i]))
    
    if i < 4000:
        # 修正後の線を引く
        plt.plot(x0_hist,x1_hist[i+1],c='red')
    for i,(x,y) in enumerate(zip(train_x[:,0],train_x[:,1])):
        plt.annotate(str(train_y[i]),(x,y))
    
    
fig = plt.figure(figsize = (6, 6))
anim = animation.FuncAnimation(fig,update_fig,interval=10,frames=len(w_hist))

すごい勢いでコードを貼りましたが、基本的にはコメントの通り。 出力結果はこんな感じ。

https://lh3.googleusercontent.com/alwL14dCOR1HiGcjJ7q6A_0IcDUuij8mhvS6Jj_2iOmpy6s9g2OsKC1S1ame6uGsrsUmSgWeKRUFvmpVy9U1_63tzvj2pfIKKxh3Q7oP01kLFgsQixMPx9A8gFBUTFImNLwYxQYWOqI6ZLo-9pg4SNMSF0VeMMyIDJhO3BAxgfuuUOX8EUmMDj-XuC0bzNnvNrt3MLpnDk5aMksI05tYVFP_YQ8KFcFr2bE5hBkxl9whbzgVUls0qIpEgFJRqxAHBNbBH-zhyksRQD0KsqfpZnvo6a1rSDicjwlPi_WyfZIWtf7pMuVEkm4eGU6y_IfDcABarsXo9RTr63JjPUEd_crLxerxK3f1VmSV_jlyfY1xfIqUrUhs2kDRueyrEo5lMuPuzScBp_UavOEoFtvJHcDbijgSi00IkoFYm1HRHYUYEi7A9oahkWLevWbPWCO8UxT0hn3hg0HbW_YIOFxPlXMm8TVF9nIEfsiEKvBausG4kE1aFDpTu7M-xwWlanmhmmK2cRY-tOHhPtpcGID-iCE_6ZPqlCEkwJWyqVev27a675tGO5eERCA2qYWc1wAoduMT6zXVqdQJUJ1PFsK3oi730oo-pp7XlJZwqhbj-U76M7idiiG7a_L4afkUd_fbvqhYde3-N3ZTQk73psT237nIXUhhb-4u5ji8RIUbjVNRuybXFU7NpqxnTWFTLz7DWExcJK8QSPa4t_1r=s600-no

黒線が学習前、赤線が学習後、数値が変わっているのがcross_entropy、赤点が学習している点、青がそれ以外の点、で、 1,1だけが分離される線がひければOK、です。

なんとなく感覚的につかめたのではないでしょうか、クロスエントロピー。

ざっつおーるせんきゅー。

詳解 ディープラーニング ~TensorFlow・Kerasによる時系列データ処理~

詳解 ディープラーニング ~TensorFlow・Kerasによる時系列データ処理~

juliaで単純パーセプトロンをライブラリを使わずに実装（を多少きれいに）

python julia 単純パーセプトロン ニューラルネットワーク

juliaで単純パーセプトロンを実装した例が僕の日本語ググり力では見つからなかったので、 こちらに多少きれいにまとめておく。

julia初心者過ぎてclassとかの作り方が不明(moduleをどうやら使うらしい)なので、 関数のみ！関数、バンザイ…！only function! no more module!

simple perceptron on julia 

#必要なライブラリ(CSVとStatsBaseの干渉には参った…)
using CSV
using StatsBase
import Base: * , >

# ステップ関数
function fn_step(x)
  return 1 * (x > 0)
end

# 予測
function fn_predict(x,w,b)
  return fn_step(w' * x +b)
end

#学習
function fn_fit(x::Array,y::Array,learning_rate = 0.01,w=true,b=true)
  local dim = size(x)[2]
  local N = size(y)[1]
  local epoch_counter = 0
  # weigth initialize
  if w
    w = ((rand(dim) - ones(dim)/2) * 2)
  end
  # bias initialize
  if b
    b = (rand(1)[1]-0.5) * 2
  end
  while true
    epoch_counter += 1
    local fin_flag = true
    local learning_order = sample(1:N,N,replace=false)
    for i in learning_order
      delta_w = learning_rate * (y[i] - fn_predict(w,x[i,:],b)) * x[i,:]
      delta_b = learning_rate * (y[i] - fn_predict(w,x[i,:],b))
      w += delta_w
      b += delta_b
      fin_flag *= (delta_w == zeros(dim)) * (delta_b == 0)
    end
    if fin_flag
      return w,b,epoch_counter
      break
    end
  end
end

以外とシンプル（？）に書けました。

あとは、こんな感じで書けば使えます。

# train_xは二次元の表で横方向にデータ種、縦方向にデータ量
train_x = convert(Array,CSV.read("train_x.csv", header=false, delim=','))
# train_yは1,0ラベルのみの一次元配列
train_y = convert(Array,CSV.read("train_y.csv", header=false, delim=','))

# 学習
# weight,bias,epoch回数を返す
w,b,e = fn_fit(train_x,train_y)

# 予測
# 予測結果を返す
predict = fn_predict(train_x[1,:],x,b)

毎回関数にw,bを乗っけるのはアホやな…。 module使ってちゃんと作りましょう。そのうち。