統計分析ソフトRのミニソース集

かゆいところに手が届くようなRのミニソースを作成

【ver.1.00】パッケージ化しました(ggplotでF分布と棄却域を描く関数)

先月,ggplotでF分布と棄却域を描く関数について,記事にした。
mini-tool.hatenablog.com

この関数について,オプション引数を増やしてソースパッケージ化した。
ソースのダウンロードはコチラからどうぞ。
右クリックして「名前をつけてリンク先を保存」を選択すると良い。
http://iku-sawa.info/R/F.dist_100.txt

使い方(基本)

#読み込み
#ggplot2とgridパッケージを予めインストールしておく必要有り(パッケージの読み込みは自動)
source("F.dist_100.txt")

#必須の引数は,df1(第1引数)とdf2(第2引数)である。
#自由度(3,120)のF分布とその棄却域
F.dist(3,120)

f:id:mini-tool:20160124174847p:plain

オプション引数

  • alpha(第3引数,デフォルト=0.05)
  • fmax(第4引数,デフォルト=20)
    • 描画するF軸の最大値
  • arrow(第5引数,デフォルト=T)
  • arrowPoint(第6引数,デフォルト=F)
    • 塗り潰し範囲の下限値(矢印を指すポイントの変更)

オプションの使い方(例1)

#自由度(3,120)のF分布で1%水準の棄却域を塗りつぶす。F値の最大値は5とする。
F.dist(3,120,alpha=0.01,fmax=5)

f:id:mini-tool:20160124180046p:plain

オプションの使い方(例2)

#自由度(3,120)のF分布で,F>=1.44の範囲を塗りつぶす。F値の最大値は5とする。
F.dist(3,120,fmax=5,arrowPoint=1.44)

f:id:mini-tool:20160124180244p:plain

オプションの使い方(例3)

#自由度(3,120)のF分布で,塗りつぶしは行わない。F値の最大値は5とする。
F.dist(3,120,fmax=5,arrow=F)

f:id:mini-tool:20160124180557p:plain

ggplotでF分布と棄却域を描く自作関数

RでF分布と棄却域を描画したいとき,片所強先生が公開されている関数がたいへん便利である。

引数に"norm", "t", "chisq", "F"のいずれかを指定するだけで簡単に分布曲線が描かれます。描こうとする分布によっては、x軸の値が適当でないこともあるので、その場合は引数myself = TRUEとすることによって分布を描くための関数自体が返されます。*1

分布のプロット(正規分布、t分布、F分布、χ2分布)

片所先生の関数を使えば,例えば,df(2,100)のF分布と有意水準を5%とした場合の棄却域は,次のように書けば良い。

#※予め関数の定義を読ませておく
d <- my.dist("F", myself=TRUE)
#自由度を聞かれるので2と100を入力する
plot(d, 0, 10) #第2,3引数は,x軸の描画範囲
my.paint(d, qf(1-0.05, 2, 100), 10, col="grey") #第2,3引数は,塗りつぶしに対応するx軸の範囲

f:id:mini-tool:20151220020637p:plain

これで十分であるが,場合によっては,もう少し丁寧に描画したいことがある。
今回は,ggplotを使った関数を書いた。
以下に関数の定義コードを示す。

F.dist <- function(df1,df2,alpha=0.05,fmax=20){
my.F <- data.frame(X=x<-seq(0,fmax,len=10000), Y=df(x,df1,df2)) #F分布の確率密度関数データフレーム
p <- ggplot()
p <- p + geom_path(data=my.F, aes(x=X, y=Y), size=1.0)
p <- p + labs(x = "F value", y = "Probability density")

#以下,棄却域
upper <- qf(1-alpha, df1, df2) #臨界値
ymin <- par("usr")[3]
ymax <- par("usr")[4]
pos1 <- ((ymin+ymax)/100)*20*0.70 #矢印↓のy座標上限
pos2 <- ((ymin+ymax)/100)*20 #テキストのy座標

p <- p + geom_ribbon(data=subset(my.F, X>=upper), aes(x=X, ymin=0, ymax=Y), fill="gray")
p <- p + annotate("text",size=4.5,label=format(upper,digits=2,nsmall=2,scientific=F),x=upper,y=pos2, fill="white")
p <- p + annotate("segment",size=1, x=upper, xend=upper, y=pos1, yend=0.0,arrow=arrow())
p <- p + scale_x_continuous(expand = c(0, 0))
p <- p + scale_y_continuous(expand = c(0, 0))
p <- p + annotate("text", x = Inf, y = Inf, label = paste("F(",df1,",",df2,")\nα=",alpha,sep=""),
            hjust=1.1, vjust=1.1, col="black", cex=6,
            fontface = "bold", alpha = 0.8)

#以下,書式の調整
p <- p + theme(axis.title.x = element_text(vjust = 0), #vjustでラベル位置調整
		axis.title.y = element_text(vjust = 2, angle = 90), #vjustでラベル位置調整
		panel.background = element_rect(fill = "white", colour = NA),
		plot.background = element_rect(colour = "white"),
		axis.line = element_line(colour = "#000000", size = 1, linetype = "solid", lineend = "round"),
		axis.text = element_text(size = rel(1)),
		axis.text.x = element_text(vjust = 1,colour="#000000"),
		axis.text.y = element_text(hjust = 1,colour="#000000"),
		axis.ticks = element_line(size = 1,colour="#000000"),
		axis.ticks.length = unit(0.2, "cm"),
		axis.ticks.margin = unit(0.2, "cm"),
		panel.grid = element_blank(),
		text = element_text(face = "plain", colour = "black", size = 12, hjust = 0.5, vjust = 0.5, angle = 0, lineheight = 0.9))

plot(p)
return(paste("F分布(",df1,",",df2,")における,α=",alpha,"の棄却域はF>=",format(upper,digits=2,nsmall=2,scientific=F),"です。",sep=""))
}

この関数を使って,df(2,100)のF分布と有意水準を5%とした場合の棄却域を描くには,次のようにすれば良い。

#※予め関数の定義を読ませておく
library(ggplot2)
library(grid)
F.dist(2,100,fmax=10) #fmaxはx軸の描画上限値(省略時20)

f:id:mini-tool:20151220021505p:plain

有意水準を変えたい場合は,引数にalphaを追加すれば良い。

F.dist(2, 100, fmax=10, alpha=0.1)

f:id:mini-tool:20151220022335p:plain

一行で描画できるため,次のようなことを理解する際に良い。

  • 自由度を変えた場合に,F分布の形状がどう変化するか?
  • 有意水準を変えた場合に,棄却域はどう変化するか?

次回は,F分布に形状が近い,χ2乗分布の関数を作成する予定である。

*1:2015年12月19日閲覧

はてなブログから各種ブログランキングへのping送信を完全リアルタイム自動化(PHP+mySQL+cron)

はてなブログは,インターフェイスがシンプルでわかりやすく,たいへん使い勝手が良い。

一方で,pingの送信機能がないため,新しい記事を追加したことを各種ブログランキングに通知できないという問題がある。

この問題に対処するため,一般的には以下のような方法が考えられる。

  • PINGOO!のようなサービスを利用して,自動化する。

しかし,毎回手動で送信するのは手間であるし,PINGOO!は,無料アカウントの場合,RSSの更新確認が不定期(有料アカウントでさえ,60分に1回が上限),という問題がある。

この問題については,id:MoneyReportさんの記事(2015, Feb.)が詳しい。

はてなブログと名乗っていながら、ブログの基本機能であったはずのPing通知が使えないのは時に不便です。

PINGOO!を利用して、はてなブログで「にほんブログ村」と「人気ブログランキング」に自動でPING通知する方法! - マネー報道 MoneyReport

このような問題を解決するため,今回は,「当該はてなブログRSSを1分おきに更新確認し,更新されていれば自動でping送信する」という,PINGOO!が提供するサービスに近いミニプログラムを自作した。

具体的なPHPコードのロジック概要は次の通りである。
このPHPをcronで1分毎に実行している。

<?php
//(1) RSSの取得
$rss = simplexml_load_file('http://●●●●●.hatenablog.com/rss');
foreach($rss->channel->item as $item){
	$title = $item->title;
	$pubDateTime = date("Y-n-j H:i:s", strtotime($item->pubDate));
	$url = $item->link;
	$content = mb_strimwidth (strip_tags($item->description), 0 , 100000, "", "utf-8");
	$content = str_replace("'", "?", $content); 
	$content = str_replace("\"", "?", $content);  //Quotoが入るとSQLでエラーになる

	//(2) DBとの同期的処理
	//もし新着記事なら,記事情報をレコードに追加
	//記事urlが,ユニークキーになるはず
	//コード略
}

//(3) ping送信処理
//もし新着記事があったら,以下を実行
// 送信先の設定
$ping_servers = array();
$ping_servers[] = "http://ping.blogmura.com/xmlrpc/●●●●●"; //ブログ村
$ping_servers[] = "http://blog.with2.net/ping.php/●●●●●/●●●●●"; //人気ブログ
$ping_servers = array_unique($ping_servers);

// 自分のサイト名
// (日本語名が使用できない送信先サーバがあります)
$site_name = "●●●●●";
// 自分のサイトURL(RSSファイルのリンクがあるページ)
$site_url = "http://●●●●●.hatenablog.com/";

// 送信するweblogUpdates.pingのXMLデータ
$post_data = '<?xml version="1.0" encoding="UTF-8"?>
<methodCall>
<methodName>weblogUpdates.ping</methodName>
<params>
<param><value>'.$site_name.'</value></param>
<param><value>'.$site_url.'</value></param>
</params>
</methodCall>
';

$headers = array(
  'Content-Type: application/xml',
  'Content-Length: '.strlen($post_data)
);

$context = stream_context_create(
  array(
    'http'=>array(
      'method'=>'POST',
      'header'=>implode( "\r\n", $headers ),
      'content'=>$post_data
    )
  )
);

foreach($ping_servers as $ping_server){
  $http_response_header = null;
  $response = @file_get_contents($ping_server,false,$context);
  //echo $response;
}

//(4) メール送信処理
//お好みで足すと便利
//もしping送信したら,指定したアドレスにメール通知
//例)ブログ(●●●●●)に,新しい記事(●●●●●)が投稿されました。
//ブログランキングへの新着記事としての掲載を,要求しましたので,お知らせいたします。
?>


尚,ping送信については,Pentan.infoさんによる,以下の記事を参考にした。
サイトの更新情報をPINGサーバに送信する方法 - [サンプルコード/PHP] ぺんたん info


cronで上記PHPを1分毎に実行するようにした結果,新しい記事を公開してから,遅くとも10分後には,ブログランキングの新着記事リストに紹介されていることが確認できた。


ーーー

ご案内

需要があれば,今後,Webサービス化を検討しております。
サーバの負荷がどのくらいかかるか等を検証してみたいので,もし,自身のはてなブログping送信も自動化してほしいという人がいらっしゃいましたら,下記情報を添えて,コメント欄にてお知らせください*1。先着10名様くらいで,無料でサービスを提供させて頂きます(※)。

  • ご自身のはてなブログのURL※1
  • ご希望の新着記事更新確認頻度(1分毎/10分毎のどちらかをお選びください※2)
  • ping送信先URL※3(ブログ村のマイページからコピペしてください)
  • ping送信時のメール通知希望(希望する/希望しないのどちらかをお選びください)
  • 連絡先兼/ping送信通知先メールアドレス※3(@gmail.comからメールを受信できるよう設定してください)

※サービスを終了する予定は,当分の間ありませんが,突然終了する可能性もありますので,予めご容赦願います。
※1 少なくとも,これまでに記事を10件以上投稿されているブログについて,ご希望ください。
※2 1分ごとに更新確認する場合,誤って記事を投稿した場合でも,すぐにping送信されてしまう恐れがあります。送信後の取消はできませんので,ご注意願います。
※3 当ブログのコメントは,承認後公開される設定です。これらの内容を含むコメントは,承認されず,内々に処理いたします。

一両日中に設定後,連絡先メールアドレス宛に,ご連絡いたします。

*1:追記:TwitterのDMでも構いません