ME差
ふとipadに入れているアプリの『血圧ノート』ってのを見たら、随分進化してた。オイラーの 自作版ソフトと同じように、平均とか偏差を出している。この結果を医者に見せてくださいって、自慢気に紹介してたぞ。
それはいいんだけど、新しい概念が発明されてた。ME差ですって。何処かの医療コンサルタントが、捻り出したんだな。薬が沢山売れて、製薬会社が儲かるようにね。
Mってのはモーニングの略、Eってのはイブニングの略。心臓がギューと縮まって血液を送り出す時、血管に圧力がかかる。これが収縮期血圧とか、最高血圧とか言われるものだ。
朝と夜でこの差が大きいと、血管プッツンして危ないですよってのが、ME差らしい。大体、この差が15から20ぐらいより大きいと、問題ですよって脅かす訳だ。幅が有るのは、まだ提唱中だから、世間一般に認知されていないからだろうね。
これと、朝晩の血圧の平均値の135を組み合わせる。135ってのは、高血圧学会なりWHOなりがお墨付きを与えた、高血圧の指標だ。
| 平均135以下 | 平均135越え
-----------+-----------------+----------------
ME差15越え | 早朝血圧上昇型 | 早朝優位高血圧
ME差15以下 | 正常血圧 | 持続性高血圧
起きる時はゆっくりと、気温差を無くす、早朝高血圧に適した薬の服用（これを言いたいのが見え見え）。某タレントさんが罹患した、悪性高血圧は、、、と、脅かすのに余念がないな。
そういう発明に乗って、オイラーも調べてみるか。幸いたっぷりとデータが有るからね。 時系列にため込んだCSVファイルになってる。テキストなんで、awkと相性が良いでしょう。 ってんで、awkを思い出しながら、スクリプトを書いてみた。
旅行に行ってるとかして、朝晩のデータが欠落してる場合が有るので、それは除外するようにしたよ。ってか、そうしないといけない。まあ、朝晩のサンプリング。かっこよく言うと、離散的だな。
腹筋とかやると血圧が上昇する事を、体で感じるんだけど、こういうのはどうなのよ。血圧上昇するから、そんな無茶するなって事？
#! /usr/bin/awk -f
BEGIN { FS = ","; }
{
hh = substr($1,7,2);
ymd= substr($1,1,6);
if ( hh < 12 ){ # Case AM
mb = $2;
today = ymd;
next;
}
if (today == ymd){ # Case PM
eb = $2;
av = (mb + eb) / 2;
dl = mb - eb;
if (av > 135 && dl > 15) r1 +=1;
if (av <= 135 && dl > 15) r2 +=1;
if (av > 135 && dl<= 15) r4 +=1;
if (av <= 135 && dl<= 15) r3 +=1;
ttl++;
}
}
END {
printf "%8d%8d\n%8d%8d\n", r2,r1,r3,r4;
printf " Total count = %d\n", ttl;
}
実行時間を測りつつ、走らせてみた。
[ham@cent tmp]$ time ./me.awk current.csv
856 4
1773 75
Total count = 2708
real 0m0.012s
user 0m0.010s
sys 0m0.001s
リナ陣営のやつは便利な機能が搭載されてるんだけど、知らず知らずのうちにそれを使ってしまって動かない事があるので、shだとかawkのコード書きはBSDで行うようにしてる。
by octave
丁度、octaveにハマっているんで、同じ事をoctaveでやってみる。
% me blood
tic();
bld = csvread('current.csv');
am = bld( rem(bld(:,1),100)<12, [1 2] ); % [ [ymdh hi lo ps] ... ]
pm = bld( rem(bld(:,1),100)>12, [1 2] );
am = [fix(am(:,1)/100) am(:,2)]; % [ [ymd hi] ... ]
pm = [fix(pm(:,1)/100) pm(:,2)];
tbl = pm(:,1); % for lookup (must be sorted)
toc()
res = zeros(2,2);
for i=1:length(am)
if ( x = lookup(tbl, am(i,1), 'm') )
av = (am(i,2) + pm(x,2)) / 2.0;
dl = am(i,2) - pm(x,2);
if (av > 135 && dl > 15) res(1,2) +=1; endif
if (av <= 135 && dl > 15) res(1,1) +=1; endif
if (av > 135 && dl<= 15) res(2,2) +=1; endif
if (av <= 135 && dl<= 15) res(2,1) +=1; endif
endif
endfor
res
sum(sum(res))
toc()
tic,tocを忍ばせて、経過時間をざっくりとモニターするようにした。
[ham@cent tmp]$ octave me.m
warning: suggest parenthesis around assignment used as truth value near line 12, column 11 in file '/tmp/me.m'
Elapsed time is 0.00973105 seconds.
res =
856 4
1775 75
ans = 2710
Elapsed time is 0.127146 seconds.
結果がawkのそれと異なっているけど何故？ 気になるぞ！ こういう動くけど、結果が違うってバグを見つけるのは至難の技が必要そう。 別の方法を試してみるかな。
その前に、今回の作の意図をば。なるべく、らしく書くのが目標だった。 前半は、らしく書けたけど、後半はありきたりだなあ。そして、圧倒的に後半が遅い！ forを使ってるから？ いいや、無理してlookupなんてのを使ったからかな。
octave:2> size(am) ans = 2723 2 octave:3> size(pm) ans = 2718 2
朝と晩の測定個数が違う。両方のデータが揃っているのを対象にしたかったんで、晩の日付のテーブルを作り、その中に朝の日付が含まれているかってのを一々lookupしてる。見つかれば、そのindexが得られる。見つからないと０が返る。こういう手法で、同時に存在するデータを得てる。超無駄な事だ。
前半部分も、手軽に何個もマトリックスを作っているけど、場所の確保やら値の設定で無駄を してるな。手軽に出来るからこそ、戒めないとな。
csvread()にどれぐらいの時間がかかっているのだろう？csvread()の後ろにtoc()を置いてみた。
octave> me Elapsed time is 0.00999403 seconds. Elapsed time is 0.01033 seconds.
前半部分で、ほとんどの時間はcsvread()に費やされているとな。IOは高価なオペレーションなんだな。マトリックスを対象にした演算は気にする程でないとな。
まてまて、遅いのはIOがからんでいるって、そんな大雑把な事で委員会？ csvreadって、読み取った数字文字列をバイナリー表現の数値に変換してからマトリックスに格納するんだよな。 何行あるか分からんファイルを相手にするんだから、一行分づつをappendするんだかLispっぽくconsするんだかは知らないけど、高負荷な事をやってるはずだな。
csvread()
何か高速にする方法が無いかと、help csvread する。dlmreadへのラッパーって事が判明。
m/io/csvread.m
function x = csvread (filename, varargin)
x = dlmread (filename, ",", varargin{:});
endfunction
続いて、help dlmread する。指定した所だけ読み込む方法が解説されてた。いわゆる範囲を指定出来る。それもあのEXCELみたいにね。そんな事で、スピード比べ。
clear all
tic()
% bld = csvread('current.csv');
bld = dlmread('current.csv', ',' , 'A1..B999999');
toc()
A列の1行目からB列の99999行分読み取ってくれ。99999ってのは十分と言う意味ね。これで、使いもしない列はオミットされるから、速くなる期待。
octave> Elapsed time is 0.031836 seconds. % by csvread octave> Elapsed time is 0.032347 seconds. octave> Elapsed time is 0.030863 seconds. octave> Elapsed time is 0.0301418 seconds. % by dlmread octave> Elapsed time is 0.02987 seconds. octave> Elapsed time is 0.0301368 seconds.
3回づつ、測定してみた。
期待したんだけど、誤差の範囲内だな。データはint型で格納されてると思ったんだけど、double型だった。デフォルトの数値型はdoubleになってる。今まで気がつかなかったぞ。
octave:40> s="19070121,123\n888,56789\n"
s = 19070121,123
888,56789
octave:41> q = textscan(s, "%d,%d");
octave:42> q
q =
{
[1,1] =
19070121
888
[1,2] =
123
56789
}
こういう機能が有るから使え？ 大カッコで括られていて、cellクラスとか言うらしい。深入りはしないでおこう。
検証
で、一番の問題は、正しい答えが出てこない事。どんなにスピードが速くても、間違っていたんじゃ、目もあてられない。
全てを疑って、元データから検証してみる。朝晩に測定値の有る日数の確認。cutコマンドのバイト位置による抽出でymdを取り出し、uniq して、その個数を数える。
[ham@cent tmp]$ cut -b1-6 current.csv | sort | uniq -c | grep ' 2 ' | wc
2706 5412 40590
[ham@cent tmp]$ cut -b1-6 current.csv | sort | uniq -c | sort
:
2 190630
3 131013
3 131021
自信のある（？）awkの実行結果とも２日程ズレていたので、もしやと思って眺めてみると、だぶりを発見。こういうのに関しては、自作のアプリでは検出出来ないからなあ。
[ham@cent tmp]$ grep 131013 current.csv 13101303,116,63,56 13101304,106,64,59 13101321,107,56,63 [ham@cent tmp]$ grep 131021 current.csv 13102103,116,72,55 13102104,116,73,58 13102121,126,70,60
入力を間違ったのかな？ 修正しておこう。
この修正を施したら、awkが出す結果とoctaveの出す結果が一致した。第三者の眼って大事だなあ。大事に至らなくてよかったわい。
profile
プロファイラーが用意されてるので、そいつのお世話にもなってみる。
octave:1> profile on octave:2> me octave:3> profile off octave:4> profshow # Function Attr Time (s) Time (%) Calls -------------------------------------------------------- 1 me 0.112 83.99 1 3 dlmread 0.010 7.16 1 11 lookup 0.006 4.67 2721 6 binary > 0.002 1.27 8125 14 binary <= 0.001 1.09 8124 12 binary + 0.001 0.64 2708 13 binary - 0.001 0.49 2708 7 binary / 0.001 0.40 2710 16 disp 0.000 0.11 2 4 rem 0.000 0.07 2 2 csvread 0.000 0.04 1 18 profile 0.000 0.03 1 15 display 0.000 0.01 2 8 fix 0.000 0.01 2 5 binary < 0.000 0.01 1 17 sum 0.000 0.01 2 9 zeros 0.000 0.00 1 19 nargin 0.000 0.00 1 21 false 0.000 0.00 1 20 binary != 0.000 0.00 1
profile('clear')で、ログデータがクリア出来る。一定時間毎に割り込んで、何をしてたか 調べているんだろうから、数度実行してみると良いかも知れない。
octave:6> profexplore
Top
1) me: 1 calls, 0.134 total, 0.112 self
2) profile: 1 calls, 0.000 total, 0.000 self
profexplore> 1
Top
me: 1 calls, 0.134 total, 0.112 self
1) csvread: 1 calls, 0.010 total, 0.000 self
2) lookup: 2721 calls, 0.006 total, 0.006 self
3) binary >: 8125 calls, 0.002 total, 0.002 self
4) binary <=: 8124 calls, 0.001 total, 0.001 self
5) binary +: 2708 calls, 0.001 total, 0.001 self
6) binary -: 2708 calls, 0.001 total, 0.001 self
7) binary /: 2710 calls, 0.001 total, 0.001 self
8) display: 2 calls, 0.000 total, 0.000 self
9) rem: 2 calls, 0.000 total, 0.000 self
10) fix: 2 calls, 0.000 total, 0.000 self
11) binary <: 1 calls, 0.000 total, 0.000 self
12) sum: 2 calls, 0.000 total, 0.000 self
13) zeros: 1 calls, 0.000 total, 0.000 self
14) length: 1 calls, 0.000 total, 0.000 self
こちらは、ターゲットを絞った内容確認になる。
ここには出ていないけど、forで回すのが決定的に遅い。forを何とか駆逐出来ないか？
for less
強引にforを表舞台から消してみた。 その仕掛けは、arrayfunと言う関数。アレー（複数指定可）と、各エレメントに対して実行したい関数のハンドルを指定する。
gaucheあたりのforeachだな。関数のハンドルってoctaveの方言が出てきたけど、scheme風に言うと、lambdaだ。
下記で定義した関数の引数xに対して処理を書いていく。関数なので、基本的には関数外は見えない。それじゃ困ると言う事で、関数の内外で、グローバル宣言と言う怪しげ事をする。超不格好で、見るからにこんな事をするんじゃ無いってのがプンプンする。
朝と晩で共通に表れる日を、集合演算 intersectで求めてみた。こんな事も出来るよって見本かな。
clear -a
bld = csvread('current.csv');
am = bld( rem(bld(:,1),100)<12, [1 2] ); % [ ymdh hi lo ps; ... ]
pm = bld( rem(bld(:,1),100)>12, [1 2] );
global am2 = [fix(am(:,1)/100) am(:,2)]; % [ ymd hi; ... ]
global pm2 = [fix(pm(:,1)/100) pm(:,2)];
global tba = am2(:,1); % [ ymd ... ]
global tbp = pm2(:,1);
both= intersect(tba,tbp);
global res = zeros(2,2);
function diside (x)
global am2 pm2 tba tbp res
i = lookup(tba, x, 'm'); mb = am2(i, 2);
j = lookup(tbp, x, 'm'); eb = pm2(j, 2);
av = (mb + eb) / 2;
dl = mb - eb;
if (av > 135 && dl > 15) res(1,2) +=1; endif
if (av <= 135 && dl > 15) res(1,1) +=1; endif
if (av > 135 && dl<= 15) res(2,2) +=1; endif
if (av <= 135 && dl<= 15) res(2,1) +=1; endif
endfunction
arrayfun(@diside, both)
res
sum(sum(res))
気になる実行時間だけど、最初に書いたものより遅くなった。Elapsed time is 0.147579 seconds. ですってさ。多分、lookupが足を引っ張っているんだろうね。
octave:2> whos
Variables in the current scope:
Attr Name Size Bytes Class
==== ==== ==== ===== =====
am 2721x2 43536 double
g am2 2721x2 43536 double
ans 1x1 8 double
bld 5439x4 174048 double
both 2708x1 21664 double
pm 2718x2 43488 double
g pm2 2718x2 43488 double
g res 2x2 32 double
g tba 2721x1 21768 double
g tbp 2718x1 21744 double
Total is 51664 elements using 413312 bytes
スクリプトの冒頭に入れている clear -a は、グローバル宣言した変数も含めて、クリアする為だ。こいつを入れておかないと、結果が累積されて、ホェーとなる。
res =
5136 24
10638 450
ans = 16248
global res = zeros(2,2);
global today mb
function diside (x,y)
global today mb res
ymd = int32(x) / 100;
if (rem(x,100) < 12)
today = ymd;
mb = y;
elseif (today == ymd)
av =( mb + y) /2;
dl = mb - y;
if (av > 135 && dl > 15) res(1,2) +=1; endif
if (av <= 135 && dl > 15) res(1,1) +=1; endif
if (av > 135 && dl<= 15) res(2,2) +=1; endif
if (av <= 135 && dl<= 15) res(2,1) +=1; endif
endif
endfunction
bld = csvread('current.csv');
arrayfun(@diside, bld(:,1), bld(:,2))
res
sum(sum(res))
またawkのそれに似せて書いたスクリプトだと、 Elapsed time is 0.165908 seconds. こんな具合に更に遅くなった。隠れforは元凶だな。
forが遅いのか？ forループの中でインタプリットしてプリミティブ関数を呼び出すのが遅いのか？ それが問題だ。
roots
ME差を肴に、色々やったけど、ちょっと趣向を変えて、４つのエリアに分かれるその原点を満たす血圧はどんな値になるんだろう？ mbを朝血圧、ebを夜血圧とすると
mb + eb = 270 mb - eb = 15
説明では、平均って事で、２で割ってたけど、それを払ってある。こうすると、連立方程式になるな。
左辺の式において、係数だけを抜き出した行列aを作る。右辺の値を列ベクトルで表す。そうしておけは、いろいろな方法で、連立方程式を解けるぞ。(一つ条件がある。det(a) != 0 の場合だけだよ）
octave:1> a = [1 1; 1 -1]
a =
1 1
1 -1
octave:2> b = [270; 15]
b =
270
15
octave:3> inv(a) * b
ans =
142.50
127.50
octave:6> linsolve (a,b)
ans =
142.50
127.50
octave:7> a \ b
ans =
142.50
127.50
N次式の解法も出来る。helpよりの抜粋。
p(x) = x^2 - 5.
c = [1, 0, -5];
roots (c)
=> 2.2361
=> -2.2361
別な例
octave:10> roots([2 3 -4 5]) ans = -2.62482 + 0.00000i 0.56241 + 0.79759i 0.56241 - 0.79759i
2*x^3 + 3*x^2 - 4*x + 5 = 0 の根を求めるんだな。便利、って言っても、身の回りで こんなのの出番は有るのか？ （複素数なんて普通の人は無用の長物だ）