こんにちは。Ｉｓｏｄａです。
このエッセイは主にパソコン関連、
数もの（数学、数式など)、
ファイナンス、科学？等々の話題をメインに
取り上げていきます。
乞うご期待！
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
２０１４年１２月１６日
『ラズベリーパイで遊ぶ』
今回はラズベリーパイと呼ばれる手のひらサイズの
パソコンをご紹介しよう。
ITエンジニア、特にサーバ系のエンジニアはLinuxをいじれないといけない。
Linux自体の導入は比較的楽だが、入れるパソコンが無い場合が多い。
そんな方にもお勧めできるのが、超小、安価なパソコンのラズベリーパイだ。
今回使用するモデルは、
ラズベリーパイB＋
こいつを利用するメリットはいくつかある。
安価である。
OSにLinuxベースのRaspbianを入れられる。
GPIO制御ができる。
以上の３点だ。
特にLinuxやプログラミングを学びたい方には最適。
もともとイギリスのEben Uptonさんという技術者の方が、大学でITを教えていた時、
PCやプログラミングの基本を知らない学生が多すぎるということで、
安価で入手できるこのラズベリーパイの開発を始めたらしい。
なので、プログラミングの基本C言語はもちろん、
Python、Ruby、Javaなんかも使える。
そして、GPIO制御ができるのが筆者が選んだ一番の利点。
写真を見てもらえばわかるように、ラズベリーパイ（以下、ラズパイ）には
電子工作でつかえるピンヘッダーが出ているので、コードを接続していろんな遊びができるのだ。
では、ラズパイでいろいろ遊んでみよう。
まず、ラズパイは基盤むき出しでマウスもキーボードも画面もない。
もちろん、これらのデバイスをUSBなんかで接続すれば使えるのだが、
それじゃあ面白くない。
ここは、WIFI経由でパソコンやタブレットから直接いじれるようにしよう。
と、すぐにいきたいところだが、まずはこいつに色々インストールしないと使えない。
インストするにはラズパイを直接いじらないといけないので、渋々、マウス、キーボード、画面を直接接続する。
で、ラズパイにOSであるRaspbianをインストール。
WIFI等の各種設定をして使えるようにする。
（この辺はネットで検索を）
さて、各種設定済みとして、いよいよラズパイにログインする。
接続するパソコンはwindowsならTeraTarm。
アンドロイドタブレットならConnectBotをインストールしておこう。
では接続します。
TeraTarmなら
ユーザ名 pi
パスワード ○○○
タブレットのConnectBotなら
pi@192.168.0.～
（ちなみに、上記の192.～の数値は各自のネットワークの設定で確認してほしい。
ラズパイが画面に接続してあれば、ifconfigコマンドでWIFIならwlan0のinetアドレス欄に載っている）
パスワードを入力してログイン完了。
これで、ラズパイをマウスやキーボード無しで、winやタブレットから直接操作できる。
基本コマンドで確認。
$ pwd
/home/pi
$ sudo apt-get update
$ sudo apt-get upgrade
この冒頭のsudoはスーパーユーザーでの操作を必要とする時に入れるコマンドです。
ラズパイはrootも取れるので、root取ってsuになればsudoは要らないです。
では、各種ソフトをインストールしてみましょう。
一応、何をやるにも必要なテキストエディタをインスト。
$ sudo apt-get install vim
vim起動
$ vim test.txt
iを押して挿入モードへ。
適当に文字を入力。
koreha test desu
escキーを押して、コマンドモードへ。
:wq
これで保存完了。
ちなみに、アンドロイドタブレットのデフォルトキーボードではescキーが無いので、
hacker’s keyboard等のアプリを入れておきましょう。
あと、もうひとつ。
NEC製のlife touch noteでConnectBotを使うとEscキー、Crtlキーが使えません。
この場合は、画面をタッチするとEscキー、Crtlキーが画面上に出るので、それを押すとよいでしょう。
：マークも通常の位置では出ないので、Shiftキーや Altキー押しながら各マークの場所を確認してください。
さて、本題に戻ります。
$ ls で確認。
test.txtがありましたか。
あったら、cat で確認。
$ cat test.txt
koreha test desu.
どうでしょうか。
ここまで来たら、あなたも初級Linuxユーザーです。
さて、ここからはハード好きな方のために、ラズパイで簡単な電子工作をひとつ。Lチカをしてみます。
ブレッドボードにLEDランプ、３００Ωの抵抗を付けます。
（抵抗値はオームの法則で決まります。無ければ１００Ω～でも可です）
LEDは電流の流れる方向が決まっているので、足の長い方から短い方へ流れるように接続。
ラズパイの GPIO 18→抵抗→LED→ラズパイのGNDへ電流コードを接続します。
プログラミング言語はpythonです。
下記のソースをvim等で書き込んで、testled.pyで保存。
import RPi.GPIO as GPIO
from time import sleep
GPIO.setmode(GPIO.BCM)
GPIO.setup(18,GPIO.OUT)
while True:
GPIO.output(18,GPIO.HIGH)
sleep(2)
GPIO.output(18,GPIO.LOW)
sleep(2)
ちなみにpythonはタブを認識するので、
while True:以降のコードはタブを入れてから入力してください。
で、実行。
$ sudo python testled.py
Ctrlを押しながらCキーで終了。
どうでしょう。Lチカ成功しましたか。
これでラズパイから電流を取り出して、ハードをコントロールする基礎が出来上がりました。
あとは、この電流の流れを回路で制御し、DCモータやサーボモータを動かしてください。
ロボットの頭脳として大活躍してくれるでしょう。
どうでしょう。ラズベリーパイ。
Linuxも使えるし、ハードの制御もできる。
おまけに小さくて持ち運びも楽。
ぜひあなたのデスク横に。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
２０１４年９月２３日
Ｒ言語で、コンデンサの充電電圧を返す関数
を作りました。
ご自由に改変して使ってください。
kondensa=function(v,t,r,c){
ans=v*(1-exp(-(t/(r*c))))
return(ans)
}
使い方
kondensa(v,t,r,c)に代入する。
vは電源電圧
tは充電時間（秒）
rは抵抗値
cはコンデンサの容量
Ｒでの使用例
kondensa(3,2,1000,0.001)
[1] 2.593994
電源３V
2秒後
1000Ωの抵抗値
1000μFのコンデンサ
コンデンサの電圧は、約2.59Vになります。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
Ｒ言語での、
コンデンサの充電時間の時間（秒）を返す関数
を作りました。
tkondensa=function(v,r,c,vt){
#コンデンサの充電式
y=function(t){
v*(1-exp(-(t/(r*c))))-vt
}
#近似解を求める。応用が利くようにニュートン法で。
#c(,)内の引数変更可（初期値-1～10)
ans=uniroot(y, c(-1,10))
#結果を返す
return(ans)
}
使い方
tkondensa(v,r,c,vt)にそれぞれ代入する
vは電源電圧
rは抵抗値
cはコンデンサの容量
vtはコンデンサの電圧
Ｒでの表示
tkondensa(5,1000,0.001,3.16)
$root
[1] 0.9996738
$f.root
[1] 2.762246e-06
$iter
[1] 10
$estim.prec
[1] 6.103516e-05
最初の0.99が近似解です。
コンデンサが3.16Vになるのは、約１秒後です。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
２０１４年６月６日
『phpで遊ぶ』
さて、今回はphpと呼ばれるプログラミング言語を使って、
ちょいと遊んでみましょう。
phpはperlと同じようにサーバ側で動かすプログラムで、
ホームページからの送信情報なんかを受けて→
サーバで処理→結果をホームページに返す、
なんてことを担当しています。
こういったサーバ系のプログラムはテキストエディタさえあれば、
だれでも手軽に作れる上、ネット上で動かすという、
インストール系のソフト＆アプリとは違う性質が
特長なんですね。
さて前置きはともかく、早速作って動かしてみましょう。
まずはお手持ちのテキストエディタに下記を打ち込んでください。
コピペでも可です。
//ここから
<?php
$kekka = $_GET[“koitu”] + $_GET[“aitu”];
echo “計算結果は”.$kekka.”です。”;
?>
//ここまで
上記は足し算プログラムで、
ユーザーの入力値を足し算計算して返してくれます。
（//以降の１行は「コメント」と言って、
プログラム上では読み込みません。
実務上は、読む人にわかりやすいように
コメントを付け加えながらコーディングしていきます。
ちなみに上記画像にはコメント文は入れていません。）
作ったら「test001.php」として保存。
ご自分のホームページスペース上に置いてください。
（ちなみにphpが使えるスペースか確認してください）
で、アップしたphpファイルにアクセスします。
例 http://hogehoge.net/hoge/test001.php
(上記「hoge」部分は各自の環境に置き換えてください)
画面上に「計算結果は0です。」と表示されれば成功です。
次に、
http://hogehogenet/hoge/test001.php?koitu=3&aitu=2
と、先ほどのＵＲＬの末尾に「?koitu=3&aitu=2」を加えます。
すると、「計算結果は5です。」と返してくれます。
ＵＲＬの末尾に?をつけると、
それ以降の文をサーバ側でGET受信してくれます。
よくgoogleなんかの検索サイトで検索語を入力して
enterキーを押すと、
ＵＲＬ欄にhttps～google.co.jp/search?q=%E3%82%A2%E…と
表示されているはずです。
検索ワードを「?～」に付加して、サーバ側に送って
いるんですね。
このphp上では$_GET[“koitu”] と $_GET[“aitu”]に
数字が代入されるので、それらを足し算し、
変数$kekkaに代入して表示させます。
こういったＵＲＬ欄に値を渡してhttp通信をするってのも、
プログラミングのひとつの醍醐味ですね。
さて、これだけじゃ単純なので、
もうちょいと複雑な計算もしてみましょう。
下記を記述して「test002.php」として保存、アップしてください。
//ここから
<?php
$suuji = $_GET[“kaijyou”];
function kj($suuji){
if($suuji==1){
return 1;
} else {
$kekka = $suuji * kj($suuji – 1);
return $kekka;
}
}
echo $suuji.”の階乗の答えは “.kj($suuji).”です。”
?>
//ここまで
どうでしょう。プログラミングっぽくなってきましたね。
上記は階乗計算するプログラムです。
階乗ってのは、例えば、「5!＝5*4*3*2*1」 ってことです。
５以下の整数を乗算していくってやつですね。
では、先ほどと同じく使ってみましょう。
例 http://hogehoge.net/hoge/test002.php?kaijyou=10
「10の階乗の答えは 3628800です。 」と
表示されれば成功です。
ＵＲＬの「?kaijyou=」の数字を適当に変えてみましょう。
http://hogehoge.net/hoge/test002.php?kaijyou=80
８０の階乗の答えが表示されたでしょうか。
（当方で確認したところ１７０の階乗までが限界でした。
それ以上だと無限大エラーが出ます。
ちなみに、手持ちのカシオ製の関数電卓ですと
７０の階乗でエラーが出ました。
階乗計算ってのは限界があるんですな）
でもって、上記の表示に成功したら、各自アップした
phpファイルは削除しておいてください。
使わない動的なファイルをネット上に残しておくのは、
セキュリティ上、好ましくないので。
さて、どうでしょう。
こういったCGI系のサーバサイドプログラミングも
結構面白いものがあります。
ヤフーやグーグル等の検索サイト、
アマゾン、楽天等のショッピングサイト、
その他もろもろ…。
インターネットで何かしようと思ったら必ず必要となる
必須技術ですね。
またこの手のプログラムは、サーバ側のデータベースと
連動させれば、スマホのアプリなんかの通信処理も可能です。
よくスマホゲームで「通信中…」と待たされる時がありますね。
あれなんかは、こういったサーバ側にデータを送ったり、
受信したりしてゲーム情報を更新しているんですね。
最新の情報テクノロジーもまずはここから。
こういった簡単なプログラムから始まります。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
２０１４年５月２日
『androidアプリを作ろう』
さて、今回はアンドロイドのアプリ作成法を
かんたんにご紹介しよう。
使用ソフトはeclipse。
JavaのＧＵＩ開発環境として有名なソフトです。
ちなみにアンドロイドアプリの開発環境の構築は、
専門のサイトがたくさんあるので、そちらを参考にしてください。
さて、まずこのeclipseを立ち上げます。
すると、プロジェクト名を入れろだのなんだのといろいろ
言ってきますが、言われた通りに入れていきましょう。
で、出来上がったのがこちら。
これだけで、すでに画面上にHello world の文字が
表示されています。
プログラマーにとって、このHello worldは
意外と意味のある言葉で、
ソフトを動かすのに基本となるソースの処理が
いろいろ書かれているのです。
例えば開発なんかで使う言語が決まったら、
まずHello worldをどうやって表示させるのか。
大抵のプログラマーはここから入っていきます。
さて、そんなHello world。
起動して一通りの手順を踏めば、自動で入力されている
のですから楽なもんです。
あとはソースコード周辺を見て、フムフムと納得すれば
いいんです。
さて、それじゃソース見てみましょう。
むむっ。どこにmainのソースがあるんじゃい。
CやJavaをＣＵＩベースで使用してきた人たちが最初に
とまどうのが、アンドロイドアプリのmainのソースが
どこにあるのかわからんて事でしょう。
プログラムってのは大抵mainと書かれている箇所があり、
そこからスタートして下の行に渡っていき、
最終的に最後の行で終了という流れになっています。
つまり、上からソースを見てmainを見つければ、
あらかたのソフトの動きがわかるんですね。
ところが、アンドロイドアプリの開発では
このmainの場所がわからない。
これはＶＢなんかのＧＵＩ系のソフト開発でも同じなんですが、
ＶＢの方は例えばボタンをクリックすると
自動でソースを書き込むエディタが出て、書き込み場所まで
指定してくれます。
ところがアンドロイドアプリはそうはいきません。
書き込み場所を自分で見つけないといけないのです。
これが最初の難関でしょう。
まあ、慣れてしまえばonCreateの中とか、終った所とかに
書いていくってのがわかるんですが、
そこまで辿り着くのに結構な時間を費やします。
というより、いざアンドロイドアプリを作ろうと決めて、
開発環境の構築～このHello worldを表示させるポジションに
もっていくまで、相当な時間を要します。
パソコンのスペックが低いと、ただでさえ重いeclipseの動きが
余計に重くなり、インストールしたけど使っていないっていう
ユーザーも多いのではないでしょうか。
しかし、重いソフトならメモリを増やして対処しましょう。
また、ＧＵＩのアプリやソフトを製作するのなら、
このmain無しにも慣れましょう。
イベントドリブンというユーザの操作に応じた処理を
記述していくのが、アプリ系プログラマーの仕事なんです。
さて、話はそれましたが、早速このHello worldを
ちょいといじってみましょう。
まず、valuesフォルダのstrings.xmlファイルを開きます。
string name=”hello_world”を下記のように変更します。
<string name=”hello_world”>こんな感じで表示させます</string>
※上記コードは、実際には１行で記述してあります。
これで一旦保存（保管）してください。
「プロジェクト」→「クリーン」でＯＫです。
グラフィカルレイアウト画面を見ると
下記のように変更されているはずです。
お次はボタンを取り付けて、そのボタンを押すと、
「こんなもんで2」と表示される処理を実装しましょう。
非常にシンプルな処理ですが、全てはここから始まります。
基本の基本をやらないと先には行けません。
これらの処理の詳しい解説は省きますが、
実機にインストールして確認してみます。
ボタンを押すと…。
「こんなもんで2」と表示されます。
どうでしょう。
結構簡単にプログラミングできますね。
アプリ・ソフト開発の基本は、まず既存のソースを
改造すること。
で、ビルド～実行してみて感触をつかんでいく。
慣れてきたらいろんな機能を積み上げて実装していく。
これがコツではないでしょうか。
皆さんもこんな感じで、アプリ製作を楽しんでください。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
２０１４年３月３０日
『Maximaで遊ぶ』
さて今回は、数式処理ソフト Maximaをご紹介します。
世に数式処理ソフトは数あれど、無料でクオリティの
高いソフトとなると少ないです。
このMaximaはその数少ないフリーソフト。
ダウンロード先はこちら。
http://sourceforge.jp/projects/sfnet_maxima/releases/
さて、早速使ってみましょう。
まずは簡単な計算から。
3*4*15
shift+enterキー
180
で、ちょっと複雑な計算もしてみましょう。
入力がちょい面倒ですが。
5*x^3+2*x^2+6*x-171=0
上の式をタブの「高次方程式の解を求める」に入力。
上の方程式欄に入れたら、ＯＫキー。
すると、
[x=3],
[x=(sqrt(851)*%i-17)/10],
[x=-(sqrt(851)*%i+17)/10]
解が３つ出ました。
今回は実数解のみでいいので、
x=3ですね。
どうでしょう。複雑な高次方程式も一発で解けます。
こんな面倒なことしなくても市販の電卓で十分という
意見も聞こえてきそうですが、これじゃなきゃ解けない
問題も多々あるんです。
例えば、ちょいちょい出てくるファイナンス系の問題。
例題
年間６億の賃貸収入が得られる不動産物件が
１２０億円で売りに出ている。
５年間この物件を所持し、５年後に１００億で売りに出す。
この物件のIRR（内部収益率）を求めよ。
また、資本コストを年１％とすると、この不動産物件に
投資すべきか否か判断してください。
（不動産物件なので、かなり話が大きいですがご容赦ください。）
さて、こんなビジネス問題が出たらどうでしょう。
目の前に年６億円の賃貸収入が得られる物件がある。
銀行から融資を受けてでも、この物件に投資する価値は
あるのでしょうか。
こんな問題なら、Maximaで一発です。
まず、数式を出しましょう。
6/(1+x)+6/(1+x)^2+6/(1+x)^3+6/(1+x)^4+106/(1+x)^5-120
=0
これでOKです。
このxを求めれば、IRRが出ます。
しかし、この数式を手計算で解くのは無理ですね。
関数・金融電卓なら解けますが、電卓を持っていなければ、
どうしましょうか。
そうです。パソコンのMaximaを使えばいいのです。
Maximaに入力してみます。
先ほどの方程式欄に、
6/(1+x)+6/(1+x)^2+6/(1+x)^3+6/(1+x)^4+106/(1+x)^5-120
で、ＯＫ。
解が出ました。
[x=0.01783462657614],
IRRが0.01783と出ました。
この不動産に１２０億円投資して、年間1.78％の
収益率があるってことです。
今の銀行預金の金利を考えれば、いけそうな気がしますが、
どうでしょう。
例題では、ここで資本コストが１％かかると書いてあります。
1.78-1=0.78
１２０億円の投資に対して、年間0.78％の実益です。
実質年間９３６０万円の利益が出ると見込んでいます。
どうでしょうか。検討の余地ありと見ますか。
例題は、年に６億円の賃貸収入を見込んでいますが、
賃貸収入はあくまでも予想値です。
テナントが埋まらないリスクもあります。
また、５年後に１００億円でこの物件を売却する予定ですが、
果たして１００億で買い手がつくかどうか。
保証はありません。
もろもろのリスクを踏まえた上で、この投資物件に
ＧＯサインを出せるのか。まさに決断ですね。
ここから先の判断はMaximaでは出せません。
最後は人間の判断になります。
さて、あなたならどうしますか。
（難しい～。そんな決断は、私ら一般人には到底無理ですな）
そんな訳で、今回はMaximaをさらっと紹介してみました。
複雑な数式を解きたい時に、このソフトは大活躍して
くれるんですね。
おっと、もう一つ。
このMaximaをインストールするとグラフ描画ソフトの
gnuplotも自動でインストールされます。
このgnuplotもあると結構便利なので、利用価値大です。
数式だけだと理解しづらい問題も、グラフに表すと簡単
且つわかりやすくなる場合も多々あります。
数式とグラフってのは、ワンセットが一番なんですね。
以上、Maximaの紹介でした。
興味のある方はぜひ使ってみてください。
************************
２０１４年４月１２日追記
統計学上、必要な標本数(サンプル数)を求める関数
ご自由に改変してください。
sample=function(n,e,z,p){
kekka=n/((e/z)^2*((n-1)/(p*(1-p)))+1)
return(kekka)
}
使い方
sample(n,e,z,p）に４つのパラメータを入れます。
n 母集団の数
e 最大誤差 おおむね0.05を代入
z 信頼係数 0.95で出す場合 1.96を代入
p 予想される母平均の比率 通常0.5を代入
pが0.5の時、解が最大になるので、
0.5の代入で無難に最大値を出しておく
例
母集団が2000人いる場合、標本数は何人必要か。
sample(2000,0.05,1.96,0.5)
でenterキー
Ｒでの表示
> sample(2000,0.05,1.96,0.5)
[1] 322.3955
約322人と出ました。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
２０１４年２月９日
『Ｒで解析 陶芸編』
前回は統計ソフトＲをいろいろいじってみましたが、
今回も引き続いてＲで分析してみましょう。
まず下の表をご覧ください。
これは、陶芸の素焼きの焼成表です。
左欄が温度、右欄が単位時間あたりのガス圧の総和。
それぞれ、「ondo」 「neturyou」と勝手に列名を付けています。
※ちなみに「neturyou」＝熱量ではありませんのであしからず。
また、各列の最初の行の数値０は、わざと削除しています。
０値は後で出てくる「ゴンペルツ曲線」に不向きのため。
さて、それをグラフ化したものがこれ。
ゆるやかな曲線を描いて温度が上昇していくのがわかると思います。
今回はこの焼成曲線をＲを用いて分析してみます。
まず表をＣＳＶファイルで保存。
続いて、このファイルをＲで読み込みます。
wwwww2=read.csv(“Book777.csv”)
読み込んだ表をグラフにプロット。
plot(ondo~.,wwwww2,xlim=c(0,5000),ylim=c(0,900),ylab=””)
これで上記のグラフと同じものが現れます。
（グラフ描画だけなら、
xlim=c(0,5000),ylim=c(0,900),ylab=””は要らないのですが、
後でこのグラフを再利用するので。）
さて、ここからがＲの本領発揮です。
まずこのグラフに沿った曲線のモデルを見つけます。
前回は単回帰分析で直線のモデルを作りましたが、
今回は曲線です。
曲線などの非線形は単回帰分析では不向きなので、
他の方法を使います。
いくつか種類がありますが、今回は「指数曲線」と
「ゴンペルツ曲線」でモデルを作ってみましょう。
まず、指数曲線です。
y=a*b^x+c
上記式のa,b,cの値を見つけます。
これはＲのnls関数を使用します。
nls(ondo~SSasymp(neturyou,Asymp,R0,lrc),data=wwwww2)
Asymp,R0,lrc,それぞれの値が出ました。
続いてこれらの値をy=a*b^x+cのa,b,cに直します。
一つずつ計算するのは面倒なので、当方がオリジナルで作った
zenkinsisu関数を使います。
（zenkinsisu関数は文末に掲載しておきます。
ご自由に使ってください。）
引数にはそれぞれAsymp,R0,lrcの値を入れます。
zenkinsisu(Asymp,R0,lrc)に数値を入れます。
zenkinsisu(1399.449,31.928,-8.588)
これでenterキー。
“-1367.521 * 0.999813688986062 ^ x + 1399.449”
指数曲線の数式が返ってきます。
y=-1367.521 * 0.9998 ^ x + 1399.449ってことです。
こいつを先ほどのグラフに描画します。
par(new=T)
curve(-1367.521 * 0.9998
^ x + 1399.449,xlim=c(0,5000),ylim=c(0,900))
どうでしょう。結構あてはまってますね。
もういっちょ。ゴンペルツ曲線も作ってみましょう。
同じくnls関数を使います。
nls(ondo~SSgompertz(neturyou,Asym,b2,b3),data=wwwww2)
Asym,b2,b3の値が返ってきます。
指数曲線と同じく、ゴンペルツ曲線の数式に直します。
ゴンペルツ曲線式は
y=a*b^exp(-c*x)
で表されます。
何とも小難しそうな式ですね。
こいつのa,b,cの値を出します。
これも一つずつ出すのが面倒なので、
当方がオリジナルで作ったgonpe関数を使いましょう。
（gonpe関数は文末に掲載しておきます。
ご自由に使ってください。）
引数にはそれぞれAsym b2 b3の値を入れます。
gonpe(Asym,b2,b3)ですから、
gonpe(914.4192,2.3795,0.9994)
で、enterキー。
“914.4192 * 0.092596864369849
^ exp (- ( 0.000600180072032461 ) * x )”
ゴンペルツ曲線の数式が返ってきます。
y=914.419 * 0.0926 ^ exp (- ( 0.0006 ) * x )ってことです。
こいつも先ほどのグラフに描画しましょう。
par(new=T)
curve(914.419 * 0.0926 ^ exp (- ( 0.0006 ) * x ),xlim=c(0,5000),ylim=c(0,900),col=”red”)
末尾のcol=”red”は、赤線で描画せいってことです。
どうでしょう。
指数曲線とゴンペルツ曲線の２つが描画されましたが、
どちらがあてはめがいいでしょうか。
赤のゴンペルツ曲線の方が合致具合がいいように見えますが、
チェックしてみましょう。
指数曲線とゴンペルツ曲線を適当な変数に入れます。
dat1=nls(ondo~
SSasymp(neturyou,Asymp,R0,lrc),data=wwwww2)
dat2=nls(ondo~
SSgompertz(neturyou,Asym,b2,b3),data=wwwww2)
dat1が指数曲線。dat2がゴンペルツ曲線です。
で、赤池情報量基準です。
AIC(dat1)
[1] 89.59759
AIC(dat2)
[1] 86.29048
数値の小さい方があてはまりがいいので、
ゴンペルツ曲線の方を選択しましょう。
ちなみにそれぞれの曲線の詳細を知りたければ、
summary()です。
summary(dat2)
これで「素焼き焼成曲線モデル」の完成です。
（あくまでもIsodaデータです。窯や焼成引っ張り時間に
よって変わってきます。ご参考までに）
y=914.419 * 0.0926 ^ exp (- ( 0.0006 ) * x )
xにneturyouを入れれば、温度が何度になるかがわかる。
あるいは、ある温度ｙの時のneturyouがわかる。
便利なモデルですね。
※但し、xが０の時はこの曲線では正確な数値は出ません。
予測式として使う場合はxの値が増加していくのを
想定しているので、これといった不都合は無いのですが。
０からの立ち上がり近辺は、ロジスティック曲線、
ゴンペルツ曲線ともに、予測値として使いづらいようです。
その場合は違うモデルを選択しましょう。
ちなみにneturyou５０００の時の温度をＲで出してみましょう。
ｘに５０００を代入します。
y=914.419 * 0.0926 ^ exp (- ( 0.0006 ) * 5000 )
> y
[1] 812.2614
812度と出ました。
さて、今回は陶芸の素焼きに応用しましたが、
このゴンペルツ曲線は商品の販売量、
人間の寿命なんかにも利用できる曲線のようです。
（もともとヒトの老化現象を解析するために
作られた曲線らしいです。）
言われてみれば、どんな画期的な商品だって
このゴンペルツ曲線のように、最初は知名度が無く
売れ行きもちょぼちょぼ。
やがて宣伝広告なんかの効果が出て右上がりで売れ、
一通り人々に行き渡ると限界点を迎えます。
最近だと地デジ対応テレビなんかがいい例ですよね。
ビジネスの現場でこの手法を使うのなら、
ある商品群の売上限界点を予測。
それらが頂点に達する前に、次の新製品の投入、
あるいは販路の拡大、価格の再検討ってところでしょうか。
いろんな所で、統計解析ってのは利用できますね。
※付録
zenkinsisu関数、 gonpe関数を掲載します。
ご自由に改変して使ってください。
zenkinsisu=function(Asymp,R0,lrc){
a=R0-Asymp
b=exp(-exp(lrc))
c=Asymp
siki=paste(a,”*”,b,”^”,”x”,”+”,c)
return(siki)
}
gonpe=function(Asym,b2,b3){
a=Asym
b=exp(-(b2))
c=-log(b3)
siki=paste(a,”*”,b,”^”,”exp”,”(-“,”(“,c,”)”,”*”,”x”,”)”)
return(siki)
}
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
２０１４年１月７日
『Ｒで統計学は最強なのか』
昨年度は統計関係の本がベストセラーになるなど、
統計学がちょっと気になる年でした。
その流れに乗って、今回紹介するのはこれ。
Ｒという統計解析ソフトです。
（画面背景の色は筆者仕様で、デフォルトは真っ白です。）
このＲはフリーソフトなので、公式サイトから
いつでもダウンロード可能。
公式サイト http://www.r-project.org/
興味のある方はぜひパソコンにインストールして
使ってみてください。
さて今回はこのＲを使って、身の回りの
いろんな確率を求めてみましょう。
まずは、パチンコや競馬などのいわゆるギャンブルと
呼ばれるもので複数回勝負した時の勝率を求めてみよう。
使用する関数は筆者がある数式をもとに
オリジナルで作ったgame_vic関数。
この関数には引数として、４つのパラメータを与えます。
勝つ回数、試行回数、勝つ確率、負ける確率。以上の４つ。
例えば、勝率０．５８のあるゲームを１３０回試行して、
８０勝以上する確率を知りたいと。
その時は、game_vic（勝つ回数、試行回数、勝つ確率、負ける確率）で、
game_vic(80,130,0.58,0.42)と入力。
すると、０．２３…と、結果を返してくれる。
80勝以上する確率約２３％のようです。
ちなみにこのgame_vic関数は
あくまでも「お遊び」として試行しています。
それでは遊び半分でいろいろやってみましょう。
まずは、競馬１日１２レースとして、単勝一番人気に全部賭けた場合、
6回以上勝てる確率は？
（この場合の「勝つ」の定義は「投資金額を１円でも上回る」
としましょう。）
で、１レースに勝つ確率は０．３としましょう。
（単勝一番人気なので、一応こんなもんかと。根拠の無い数字です。）
これは結構きついですね。
次に庶民の娯楽、パチンコはどうかと。
パチンコの１回の勝てる確率はいくつにしましょうか。
学生時代はよくやったのですが、最近はご無沙汰です。
当時の記憶で、３～４回に１回は勝てるとしましょうか。
勝率は０．２８位にしておきます。
週１回、年５６回として、そのうち半分の３０回
勝てる確率はどうでしょう。
1.951946e-05
0.0000195…。
かなり厳しいですね。
それでは、大数の法則チェック。
10000回の試行で2800回の勝ちは。
ん～ん。50パーに近づきますね。
てことは…。
ちなみに宝くじもこのgame_vic関数でやってみたのですが、
あまりにも夢のない数値が出たので、掲載は止めときます。
そう、宝くじは夢を買うのです…。
さて、どうでしょう。統計ソフトＲ。
本来なら統計学を駆使して、「～検定」やら「～分析」なんかに
使うソフトなんですが、こんな感じで入り口は楽しんで
使おうではありませんか。
えっ、統計解析ソフトなんだから、分析して欲しいって。
ん～ん。じゃひとつだけやってみましょう。
データは何がいいでしょうか。
そうですね…。
「陶芸作品の完成度と焼成回数」なんてどうでしょう。
やはり創り屋ですから、芸術は外せません。
まずはデータをエクセルのシートに書き込みます。
Ｒにもデータを書き込むエディタがあるので、
それを使ってもいいのですが、データ管理は
エクセルで行っているユーザーも多いと思います。
そのためエクセルのデータを使う方法で読み込んでみましょう。
まず、焼成回数（syousei）と作品完成度（kanseido）の
データを記入してＣＳＶファイルで保存。
※この「作品完成度」という数値にしづらい変数は、
数量化理論Ⅰ類を使って数値に置き換えます。
簡単に言えば、いい作品は１、今ひとつの作品は０で置き換えます。
その後、作品数に乗じて数値を出しました。
次に、保存したＣＳＶデータをＲで読み込みます。
Ｒコマンドを使います。
data11=read.csv(“Book1.csv”)
まず、この２つの数値に相関があるかチェック。
cor(data11)
ありそうですね。１に近い数値が出ています。
続いて、グラフにプロット。
Ｒはグラフの描画も簡単にできます。
plot(data11)
ん～ん。右上がりの典型的な相関図です。
これは正の相関ありですね。
お次は最小２乗法による単回帰分析です。
zzz=lm(kanseido~syousei, data11)
（kanseido値をy = b + ax の yにもっていきたいので、
kanseido~syousei を入れました。
syousei~kanseidoだとsyousei値がyに入ります。）
続いて、詳細チェック
summary(zzz)
いろいろ書いてありますが、大事なポイントを数点。
インターセプト 8.0455。
係数 3.8455。
これで、次の単回帰式が導かれます。
y = 8.0455 + 3.8455x
この数式の意味はこんな感じです。
作品の完成度 = 8.0455 + 3.8455 * 焼成回数
この式が出てしまえば、予測値に使えます。
例えば、焼成回数が１４回だと、作品完成度はどのくらいなのか。
焼成回数に１４を入れて、計算すればいいのです。
( いわゆる外挿になるので、あくまでも参考予測値です )
Ｒで計算してみましょう。
> ggg= 8.0455 + 3.8455*14
> ggg
[1] 61.8825
１４回目の焼成で、約６１の完成度と出ました。
陶芸の場合、やはり「適度な作品数」と、「焼成回数」が多い人ほど、
「作品の完成度」も高くなるという結果が出ています。
(Isoda調べ)
このように、データに相関関係があれば、数式に直して
将来の予測に利用できる。
非常に便利な方法ですね。
ちなみに、決定係数、Ｐ値なんかも大事な数値ですが、
これらについては、また後日取り上げてみましょう。
さて最後に、先ほど導いた式をグラフに書き込みましょう。
abline(ｚｚｚ)
これでＯＫです。
プレゼンに使えそうなほど立派な資料が、
ものの数十分で完成してしまいました。
こういった「～分析」をビジネスの現場で生かすのなら…。
例えば、店舗面積と売上高、サービスメニューと利益率など、
ｘ軸に何をもってくるかを中心に考えればいいと思います。
こいつは、いろんな用途に使えそうですね。
さて、どうでしょう。
統計学を使いこなすには、まずこのＲ。
こいつを使いこなそうではありませんか。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
Ｒでのアドオン率から実質年率を算出する関数作りました。
ご自由に改変してお使いください。
引数には、アドオン率x、期間nを入れる。
nenritu=function(x,n){
#アドオン率の式をちょい変換
y=function(r){
(r/12*(1+r/12)^n*n/((1+r/12)^n-1))-(1+x)
}
#近似解を求める。応用が利くようにニュートン法で。
#c(,)内の引数変更可（初期値-1～10)
ans=uniroot(y, c(-1,10))
#結果を返す 実質年率
return(ans)
}
nenritu関数の使い方
引数にはアドオン率x、期間nを入れる。
例
アドオン率0.02728978
期間12ヶ月
xに0.02728978を入れる。
nに12を入れる。
nenritu(0.02728978,12) でenterキーを押す。
Ｒでの表示
$root
[1] 0.04999692
$f.root
[1] -1.689646e-06
$iter
[1] 8
$estim.prec
[1] 6.103516e-05
最初に表示された0.04999692が近似解です。
実質年率 約0.05です。
***********************************
adon3関数
Ｒでのローン支払いの際のアドオン率と
利子、毎月支払額を調べる関数を作りました。
単純に方程式に代入すればもっと簡単なコードになりますが、
応用が利くようにニュートン法なんぞも使っています。
ご自由に改変してください。
引数には実質年率r、支払回数n、借入額tを入れる。
adon3=function(r,n,t){
#アドオン率の式をちょい変換
y=function(x){
(r/12*(1+r/12)^n*n/((1+r/12)^n-1))-(1+x)
}
#近似解を求める。応用が利くようにニュートン法で。
#c(,)内の引数変更可（初期値0～10)
ans=uniroot(y, c(0,10))
#ans[1]を型変換
ans2=as.numeric(ans[1])
#利子を算出
risi=t*ans2
#毎月の支払額を算出
kekka=(t+risi)/n
#結果を返す アドオン率、利子、毎月支払額
return(list(ans2,risi,kekka))
}
adon3関数 使い方
引数には実質年率r、支払回数n、借入額tを入れる。
例
借入額4000000円
実質年率3.5％
１０年間の毎月払い
上記のアドオン率、利子、毎月支払額を知りたい。
adon3(r,n,t)
rに0.035
nに120
tに4000000
を入れる。
adon3(0.035,120,4000000)で、enterキー
Ｒでの表示
[[1]]
[1] 0.1866304
[[2]]
[1] 746521.6
[[3]]
[1] 39554.35
アドオン率 18.66％
利子 746522円
毎月支払額 39554円と出ました。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
追記 2014年1月10日
Ｒでの「反復試行の確率」関数作りました。
ソースを掲載しておきます。ご自由に改変して使ってください。
hanpuku_sikou=function(r,n,p,q)
{
ans=(gamma(n+1)/(gamma(r+1)*gamma((n-r)+1)))*p^r*q^(n-r)
return(ans)
}
例題
マークシートの４択問題で適当に○付けして、
全１０問中、６問正解する確率は？
hanpuku_sikou(正解する数,試行回数,当たる確率,当たらない確率）
hanpuku_sikou(6,10,0.25,0.75)
解 0.016222
1.6パーセントと出ました。
ご参考までに。
※ちなみに階乗計算しているので
試行回数180前後が限界のようです。
＊＊＊＊＊＊＊＊＊＊＊＊
追記 ２０１４年１月１７日
Ｒによる
ブラックショールズ式、オプションプレミアムの算出
ソースコード作りました。
bl=function(s,k,r,sg,t){
#日数を３６５日で割る
t2=t/365
#d1、d2をそれぞれ算出
d1=(log(s/k)+(r+(sg^2/2))*t2)/(sg*sqrt(t2))
d2=d1-sg*sqrt(t2)
#コール、プットを算出
call=s*pnorm(d1)-k*exp(-r*t2)*pnorm(d2)
put=call-s+exp(-r*t2)*k
#結果を返す
return(list(call,put))
}
使用法
下記の５つの引数を入れます。
s 現在株価
k 行使価格
r 安全資産利子率
sg ボラティリティ
t 期間（日数）
例
現在株価 15000円
行使価格 15000円
安全資産利子率 1％
ボラティリティ 20％
期間（日数） 30日
bl(15000,15000,0.01,0.2,30)を入力します。
すると、１番目はコール、２番目はプットが返ってきます。
[[1]]
[1] 349.1287
[[2]]
[1] 336.805
制作後記
※Ｒはpnorm関数があるから楽ですね。
これがＶＢ6だと結構な手間がかかります。
＊＊＊＊＊＊＊＊＊＊＊＊＊
Ｒによる
オプションデルタの算出（ブラックショールズ式利用）。
ご自由に改変して使ってください。
delta=function(s,k,r,sg,t){
t2=t/365
d1=(log(s/k)+(r+(sg^2/2))*t2)/(sg*sqrt(t2))
cdelta=pnorm(d1)
pdelta=pnorm(d1)-1
return(list(cdelta,pdelta))
}
使用法
下記の５つの引数を入れます。
s 現在株価
k 行使価格
r 安全資産利子率
sg ボラティリティ
t 期間（日数）
例
現在株価 15000円
行使価格 15000円
安全資産利子率 1％
ボラティリティ 20％
期間（日数） 30日
delta(15000,15000,0.01,0.2,30)を入力します。
１番目はコールデルタ、２番目はプットデルタが返ってきます。
[1]]
[1] 0.5171507
[[2]]
[1] -0.4828493
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊