proc fcmp outlib = sasuser.funcs.logbase; function logbase(x, base); return (log(x) / log(base)); endsub; run; options cmplib = sasuser.funcs; data test; a = logbase(48, 7); b = logbase(49, 7); proc print; run;
OBS a b 1 1.98940 2
確率変数 $Y_i$ が互いに独立に canonical form の指数型分布族、
\[
f_{Y_i}(y_i)=\exp\left\{ [y_i \theta_i - b(\theta_i)] / \phi^2 - c(y_i, \phi) \right\},
\] に従う事を仮定する ($Y_i \overset{\rm{i.i.d.}}{\sim} f_{Y_i}(y_i)$)。
また、$b(\theta_i)$ および $\phi^2$ は、仮定した確率分布によって "決まっている" 事に注意する必要がある (例えば、ポアソン分布を仮定した場合、$\phi^2=1$)。
確率変数 $Y_i$ の期待値を $E[Y_i]=\mu_i$ と書く。
この期待値に対して、既知のリンク関数、$g(\cdot)$、を用いて
\[
g(\mu_i)={\bf x}_i^{\mathrm t}\boldsymbol{\beta},
\] を仮定する。
${\bf x}_i$ は各標本のデザインベクトル、$\boldsymbol{\beta}$ はパラメータベクトルである。
対数尤度関数は、
\[
l=\sum_{i=1}^{n} [y_i \theta_i - b(\theta_i)] / \phi^2 - \sum_{i=1}^{n} c(y_i, \phi),
\] と書くことができる。
$\theta_i$ をパラメータと見なし、スコア関数の期待値が $0$ になるという性質、
\[
\mathrm{E}\left[ \left\{ Y_i - \frac{\partial b(\theta_i)}{\partial \theta_i} \right\} \bigg/ \phi^2 \right] = 0
\] から、
\[
\mathrm{E}[Y_i]=\mu_i=\frac{\partial b(\theta_i)}{\partial \theta_i},
\] が分かる。
Fisher 情報行列の性質、
\[
\mathrm{Var}\left[ \left\{ Y_i - \frac{\partial b(\theta_i)}{\partial \theta_i} \right\} \bigg/ \phi^2 \right] = -\mathrm{E}\left[ -\frac{1}{\phi^2} \frac{\partial^2 b(\theta_i)}{\partial \theta_i^2} \right],
\] から、
\[
\mathrm{Var}(Y_i)=
\phi^2 \frac{\partial^2 b(\theta_i)}{\partial \theta_i^2} \equiv \phi^2 v(\mu_i),
\] であり、$v(\mu_i)$ は分散関数とも呼ばれている。
$\boldsymbol{\beta}$ に対するスコア関数は、
\[
\begin{align*}
\frac{\partial l}{\partial \boldsymbol{\beta}} &=
\frac{1}{\phi^2} \sum_{i=1}^{n} \left\{ y_i \frac{\partial \theta_i}{\partial \boldsymbol{\beta}} - \frac{\partial b(\theta_i)}{\partial \theta_i}\frac{\partial \theta_i}{\partial \boldsymbol{\beta}} \right\} \\ &=
\frac{1}{\phi^2} \sum_{i=1}^{n} (y_i - \mu_i) \frac{\partial \theta_i}{\partial \mu_i} \frac{\partial \mu_i}{\partial \bf{\beta}} \\ &=
\frac{1}{\phi^2} \sum_{i=1}^{n} (y_i - \mu_i) \frac{1}{v(\mu_i)} \left(\frac{\partial g(\mu_i)}{\partial \mu_i}\right)^{-1}{\bf x}_i \\ &=
\frac{1}{\phi^2} \sum_{i=1}^{n} (y_i - \mu_i) w_i \delta_i {\bf x}_i
\end{align*}
\] である。ただし、
\[
\delta_i = \frac{\partial g(\mu_i)}{\partial \mu_i}
\] であり、
\[
\frac{\partial \theta_i}{\partial \mu_i}=\left( \frac{\partial \mu_i}{\partial \theta_i}\right)^{-1}=\frac{1}{v(\mu_i)}
\] \[
\frac{\partial \mu_i}{\partial \boldsymbol{\beta}}= \frac{\partial \mu_i}{\partial g(\mu_i)}\frac{\partial g(\mu_i)}{\partial \boldsymbol{\beta}}= \left(\frac{\partial g(\mu_i)}{\partial \mu_i}\right)^{-1}\frac{\partial {\bf x}_i^{\mathrm t}\boldsymbol{\beta}}{\partial \boldsymbol{\beta}}= \left(\frac{\partial g(\mu_i)}{\partial \mu_i}\right)^{-1}{\bf x}_i
\] を用いた。
行列を用いて書くと、
\[
\frac{\partial l}{\partial \boldsymbol{\beta}}= \frac{1}{\phi^2}\bf{X}^{\rm{t}}\bf{W}\boldsymbol{\Delta}(\bf{y} - \boldsymbol{\mu})
\] であり、
\[
\frac{\partial l}{\partial \boldsymbol{\beta}} \bigg|_{\boldsymbol{\beta}=\hat{\boldsymbol{\beta}}} = \bf{0}
\] を満たす解、$\hat{\boldsymbol{\beta}}$、が最尤推定量である。
解析的に解けない場合の方がたぶん多い。
対数尤度関数の二階偏導関数は、
\[
\frac{\partial^2 l}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^{\rm{t}}} = -\frac{1}{\phi^2} \bf{X}^{\rm{t}} \bf{W} \boldsymbol{\Delta} \frac{\partial}{\partial \boldsymbol{\beta}^{\rm{t}}} (\bf{y} - {\boldsymbol \mu}) + \rm{\frac{1}{\phi^2}} \bf{X}^{\rm{t}} \frac{\partial \bf{W} \boldsymbol{\Delta}}{\partial \boldsymbol{\beta}^{\rm{t}}} (\bf{y} - \boldsymbol{\mu})
\] である。
これについて期待値をとり、マイナスを付けると
\[
I(\boldsymbol{\beta}) = -\mathrm{E}\left[ \frac{\partial^2 l}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^{\rm{t}}} \right]= \frac{1}{\phi^2} \bf{X}^{\rm{t}}\bf{W}\boldsymbol{\Delta}\frac{\partial \boldsymbol{\mu}}{\partial \boldsymbol{\beta}^{\rm{t}}} + \bf{0}= \rm{\frac{1}{\phi^2}} \bf{X}^{\rm{t}}{\bf W}\boldsymbol{\Delta}{\boldsymbol \Delta}^{-1}{\bf X} = \rm{\frac{1}{\phi^2}} \bf{X}^{\rm{t}}\bf{W}\bf{X}
\] が得られる。
間違ってたら怖いな、特にベクトル二階偏微分のあたりを見直そう。
あとは、quasi-likelihood についても復習しようかな。
これは覚えておくとお得かもしれません。
set data; by variable; if last.variable then output;
の仲間ですね。
data sample; input x y z @@; cards; 1 2 3 4 5 6 7 8 9 ; data sample_final; retain xx xy xz yy yz zz 0; set sample end = final; xx = xx + x * x; xy = xy + x * y; xz = xz + x * z; yy = yy + y * y; yz = yz + y * z; zz = zz + z * z; if final then output; proc print; run;
先日 Reference class に roxygen 用のコメントを入れる方法 という記事を書いたのですが、実は全然ドキュメントを作れていなかったことに今更気づいたので、せっかくなので roxygen2 用の場合を作ってみました。
基本的には変わっていませんが、@section
で fields 等の部分を記述していたり、[http://www.mail-archive.com/roxygen-devel@lists.r-forge.r-project.org/msg00091.html:title=@aliases]
の問題 (Roxygen-devel) があったり、個人的には解決しておきたい問題がありました。
パッケージの Rd ファイルは以下の様に書くことにしました。
#' testRefclassPkg: A test package for reference class documents using 'roxygen2' #' #' A Test Package for Reference Class Documents Using 'roxygen2' #' #' This package is an example of how to document reference class methods. #' #' @name testRefclassPkg-package #' @aliases testRefclassPkg #' @rdname testRefclassPkg-package #' @docType package #' @keywords documentation #' #' @author #' Author: Triad sou. #' Maintainer: Triad sou. \email{triadsou@@gmail.com} #' @seealso \code{\link[methods:ReferenceClasses]{ReferenceClasses}} #' @exportPattern '^[^\\.]' NULL
ポイント
\title
フィールド、\description
フィールド、\details
フィールドの三つを @
無しで略記できます。@name testRefclassPkg
: \name
フィールドを指定します。パッケージの roxygen2 コメントを書く場合、最後に続けて関数定義やクラス定義などを記入できないため、NULL
や {}
などを記入します。この場合、関数名やクラス名が存在しないため、\name
フィールドを自動的に作成できずエラーが発生してしまいますが、自分で指定しておけば問題ありません。@aliases testRefclassPkg
: help("keyword")
や ?keyword
でヘルプを検索するための検索ワードを指定しています。いちいち -package
を付けるのが面倒なので、個人的にパッケージ名を検索ワードとして登録しています。@rdname testRefclassPkg-package
: Rd ファイル名を指定します、好きな名前を指定して問題ありません。@docType package
: \docType
フィールドにパッケージ用の Rd ファイルであることを指定します。@keywords documentation
: \keyword
フィールドを指定します。@author
: \author
フィールドを指定します。@seealso
: \seealso
フィールドを指定します。ヘルプ間でリンクを張って、参考情報をたどれるようにするために指定しておきます。@exportPattern '^[^\\.]'
: NAMESPACE ファイルを自動生成するために、すべてにマッチする directive を書いておきます。#' My test class #' #' My Test Class (description) #' #' \code{test_refclass} class is a test reference class that has three methods (details). #' #' @name test_refclass-class #' @aliases test_refclass #' @docType class #' #' @section Fields: #' \describe{ #' \item{\code{f1}: }{numeric number of a test field.} #' \item{\code{f2}: }{numeric number of a test field.} #' } #' @section Contains: #' NULL #' @section Methods: #' \describe{ #' \item{\code{m1(arg1 = 1)}: }{\code{m1} method set a value to \code{f1}.} #' \item{\code{m2(arg1)}: }{\code{m2} method return \code{f1 + arg1 * 2}.} #' \item{\code{m3(arg1)}: }{\code{m3} method return \code{f1 + arg1 * 3}.} #' } #' @keywords documentation #' @family test_refclass test_refclass <- setRefClass( Class = "test_refclass", fields = c("f1", "f2"), methods = list( m1 = function(arg1 = 1) { f1 <<- arg1 }, m1 = function(arg1 = 1, arg2) { f1 <<- arg1 + arg2 }, m2 = function(arg1) { if (class(f1) == "uninitializedField") f1 <<- 2 f1 + arg1 * 2 }, m3 = function(arg1) { if (class(f1) == "uninitializedField") f1 <<- 3 f1 + arg1 * 3 } ) )
ポイント
@name test_refclass-class
: \name
フィールドを指定します。@aliases の問題で、@aliases test_refclass-class
と書いてしまうと、Rd ファイルの \alias
フィールドがダブルクォートで囲まれてしまい、Rcmd check
から警告をいただくことになります。ですので、@name
と @aliases
を逆に書いて、この問題に対処しています。@aliases test_refclass
: \alias
フィールドを指定します。@docType class
: クラス用の Rd ファイルであることを指定します。@section Fields:
: Bioconductor にある S4 クラスのドキュメントを参考にすると、\section{Slots}
フィールドや \section{Methods}
フィールドが用いられていることが分かります。それに合わせて、Reference class では \section{Slots}
の代わりに \section{Fields}
を作成してみました。ちなみに、@section XXX:
で任意のセクションをヘルプファイルに入れておくことができます。@section Contains:
: 不要かもしれませんが、個人的には \section{Contains}
フィールドを作って、継承に関する情報も併記してあった方が良いと思います。@section Methods:
: \section{Contains}
フィールドを指定します。@keywords documentation
: \keyword
フィールドを指定します。@family test_refclass
: roxygen2 の特殊キーワードで、@family グループ名
の形式で記入しておくと、同じグループ名を指定した Rd ファイルに対して \seealso
フィールドに相互リンクが自動的に生成されます。#' The \code{m1} method of \code{test_refclass} class #' #' \code{m1} method #' #' \code{m1} method set a value to \code{f1}. #' #' @name m1,test_refclass-method #' @aliases m1 #' @rdname test_refclass-m1 #' @docType methods #' #' @param arg1 numeric number of a test argument #' @usage \S4method{m1}{test_refclass}(arg1 = 1) #' @examples #' classObj <- test_refclass$new() #' print(classObj$m1(10)) #' @family test_refclass NULL #' The \code{m2} method of \code{test_refclass} class #' #' \code{m2} method #' #' \code{m2} method return \code{f1 + arg1 * 2}. #' #' @name m2,test_refclass-method #' @aliases m2 #' @rdname test_refclass-m2 #' @docType methods #' #' @param arg1 numeric number of a test argument #' @usage \S4method{m2}{test_refclass}(arg1) #' @examples #' classObj <- test_refclass$new() #' print(classObj$m2(10)) #' @family test_refclass NULL #' The \code{m3} method of \code{test_refclass} class #' #' \code{m3} method #' #' \code{m3} method return \code{f1 + arg1 * 3}. #' #' @name m3,test_refclass-method #' @aliases m3 #' @rdname test_refclass-m3 #' @docType methods #' #' @param arg1 numeric number of a test argument #' @usage \S4method{m3}{test_refclass}(arg1) #' @examples #' classObj <- test_refclass$new() #' print(classObj$m3(10)) #' @family test_refclass NULL
ポイント
@name m1,test_refclass-method
: \name
フィールドを指定します。@aliases の問題については同様のため。@aliases m1
: \alias
フィールドを指定します。@rdname test_refclass-m1
: Rd ファイル名を指定します、好きな名前を指定して問題ありません。@docType methods
: メソッド用の Rd ファイルであることを指定します。@param
: \arguments
フィールドを指定します。@usage
: \usage
フィールドを指定します。\S4method{generic}{signature_list}(argument_list)
を流用しています。@examples
: \examples
フィールドを指定します。@example /パッケージホームからの相対パス
で書くと、ファイルを指定して読み込むことも可能です。@example /inst/examples/xxxx.R
を用いている例が実際にいくつかありました。@family test_refclass
: 相互リンクを指定しておきます。この例ではクラスとメソッド全部に相互リンクを張っています。上記のスクリプトを *.R ファイルとして保存し、DESCRIPTION ファイルなど他の必要なファイルと一緒に適切にディレクトリに配置し、
require(roxygen2) roxygenize(path, copy.package = FALSE)
の様に roxygenize() 関数を適用すると、Rd ファイルが生成されますので、これらをまとめてパッケージ化すると完成です。
以下に、パッケージ化してインストール後に生成されたドキュメントを示しておきます。
roxygen 版の記事では @slot
云々の話を書いていたのですが、結局書いても無意味だった事に気づいていなかっただけで、roxygen2 で書いた方が良いということが分かりました!
ちなみに roxygen では @section
も利用できません。
以下のようにすると簡単に横棒付きの箱ヒゲ図を作成できます。
require("ggplot2") df <- data.frame( x = rep(1:5, each = 50), y = rt(2500, 3) ) ggplot(df, aes(x = factor(x), y = y)) + stat_boxplot(geom = "errorbar", stat_params = list(width = 0.5), geom_params = list()) + geom_boxplot() + theme_gray(25, "serif") ggsave("myplot.png", height = 7, width = 7)
この方法には一つ難点があって、stat_boxplot(geom = "errorbar", width = 0.25)
のように書いても横棒の長さを変更できませんでした。
理由がちょっと分からないので、情報をもらえると助かります。
kohske さんより頂いた情報で解決できました。
stat_params = list()
と geom_params = list()
は知っておくとお得かもしれません。
stat_qq()
の dparams = list()
みたいに投げたいところにパラメータを与えられるイメージですね。
stat_params = list(width = 0.25)
自分でスクリプトを書いてしまえば、もっと細かく設定できますね。
参考:How to add horizontal lines to ggplot2 boxplot?
1.0.1の次のバージョンでは
ggplot(df, aes(x = factor(x), y = y)) + stat_boxplot(geom = "errorbar", width = 0.5) + geom_boxplot() + theme_gray(25, "serif")
と書けるようになるみたい。
Rcmd Rd2dvi
で PDF マニュアルを作成することがあるのですが、いつの間にか失敗してしまうようになっていて困っていました。
しかし、単に Rcmd Rd2dvi
の実行には MikTeX が必要 (推奨) だったようで、インストールすると問題無く PDF マニュアルを作れるようになりました。
以下のようなログが出てしまう場合は MikTeX を入れてみると良いかもしれません。
Hmm ... looks like a package sed.exe: not found sed.exe: not found kpsewhich.exe: not found cat: not found Warning: running command '"texi2dvi.exe" --pdf "Rd2.tex" ' had status 1 Saving output to 'test-manual.pdf' ... Done You may want to clean up by 'rm -rf TEMP/RtmphKXMi6/Rd2pdf678e7c97'
sed.exe: not found
みたいなエラーだったので、パス関係かと思って色々いいじっていたのにこんな簡単に解決できるとは・・・
R で Excel などに入力されたちょっとしたデータを読み込む時、データ整形に結構手間がかかってめんどくさいと思っていましたが、
How to read space delimited data into a data frame from your script/document file? (Stack Overflow) ではコピペだけで実行できる方法が紹介されていました。
dat <- read.table(textConnection("person1 12 15 person2 15 18 person3 20 14"), stringsAsFactors=FALSE) dat
> dat V1 V2 V3 1 person1 12 15 2 person2 15 18 3 person3 20 14
header
とかが使えた... (2011/09/16) 少し編集し、変数名入りのものでも使えるようにしてみました。
dat1 <- read.table(textConnection(" person1 12 15 person2 15 18 person3 20 14 "), stringsAsFactors = FALSE, header = FALSE) dat1 dat2 <- read.table(textConnection(" name x1 x2 person1 12 15 person2 15 18 person3 20 14 "), stringsAsFactors = FALSE, header = TRUE) dat2
> dat1 V1 V2 V3 1 person1 12 15 2 person2 15 18 3 person3 20 14 > dat2 name x1 x2 1 person1 12 15 2 person2 15 18 3 person3 20 14
調べ切れていないだけで R base にこの機能入ってたりするのかな?
クリップボードから直接読み込む方法です。
id:dichika さん、ありがとうございます。
read.delim("clipboard") #windonws read.delim(pipe("pbpaste")) #Mac
Rcmdr で使用されているスクリプト。
read.table( "clipboard", header = TRUE, sep = "", na.strings = "NA", dec=".", strip.white = TRUE )