自由度 $\gamma$ の $t$ 分布の確率密度関数は、 \[ f(T=t\mid \gamma)= \frac{1}{\sqrt{\gamma}B\left(\frac{1}{2}, \frac{\gamma}{2}\right)}\left(1+\frac{t^2}{\gamma}\right)^{-\frac{1+\gamma}{2}} \] である。 ただし、$B(a, b)=\Gamma(a)\Gamma(b)/\…

R 2.10.x にも大分慣れましたが、ヘルプファイルの形式が text (R 画面内に出てくるやつ) だったので、 html に切り替えました。 よく忘れてしまうので、またメモ。 options(help_type="html") これを、".Rprofile" にでも入れておけばいいみたいです。 ".Rp…

Subversion Featuring Ely* - Alright

なに procedure だっけ、と思って調べたのでメモ。 data d1; input x1 x2 @@; diff = x1 - x2; cards; 10 12 20 21 30 21 15 18 10 22 ; proc univariate data = d1; var diff; output out = t1 t=t probt=probt msign = msign probm = probm signrank = sig…

$n$ 個の独立な確率変数 $\mathbf{Y}$ が、$p$ 個のパラメータ $\boldsymbol{\theta}$ を持つ指数型分布に従うとすると、$\mathbf{Y}$ の確率密度関数は、 \[ \begin{align*} f(\mathbf{Y}=\mathbf{y} \mid \boldsymbol{\theta}) &= \left\{\prod_{i=1}^n f(…

ちょっとだけ勉強したのでメモしておこう。 一般化線型モデル まず、確率変数ベクトル $\mathbf{Y}$ の期待値ベクトルを \[ \boldsymbol\theta=\mathrm{E}(\mathbf Y) \] と書いておきます。 一般化線型モデルと呼ばれている統計モデルでは、 \[ g(\boldsymb…

SAS Excel一括変換するマクロ に書いたものを改良したものです。 最近は R ばかり使っていたので、作ったのを忘れていました (2007/07/02に作った物を多少弄ってあります)。 Macro Name: SAS to Excel (xls / csv) batch conversion macro Author: Triad Sou…

ggplot2 は非常に良いパッケージですね、R をグラフィックスで推すときに、説得力のある実例になるんじゃないかと思いました。 Kaplan-Meier plot の実装例がない様なので、自作してみました。 かなり整理されていなくて申し訳ないのですが、 library(ggplot…

久々に W32TeX をアップデートしてみたら Undefined control sequence... ! Undefined control sequence. \Gm@rmargin ->\Geom@rmargin l.52 ...ext margin left=1em,text margin right=1em} ? ! Emergency stop. \Gm@rmargin ->\Geom@rmargin l.52 ...ext m…

服部 哲弥. 統計と確率の基礎. 学術図書出版社, 2006 2ch で服部先生のルベーグ積分の問題と回答集がいいって紹介されていたので、「統計と確率の基礎」も分かりやすいかなと思って買ってみました。 理学部向けっぽい感じがするので、数式を読むのが苦じゃな…

重み付き最小二乗法において、残差分散 \[ V_{e}= \frac{\sum_{i=1}^nw_i(Y_i-\hat{Y}_i)^2}{\phi_e},~ \phi_e=n-p \] が $\s^2$ の不偏推定量である事、すなわち \[ \mathrm{E}\left(\sum_{i=1}^n w_i e_i^2\right)=(n-p)\s^2 \] が成り立つ事を示したい。 …

アドオン FireGestures (マウスジェスチャ) NoScript (Javaスクリプト, Javaなどをブロック) Adblock Plus (広告抹消) Drag Drop Upload (そのまんま) Conquery (文字列選択して右クリ→カスタムサーチ、xml を直接編集すると任意のサーチエンジンを登録でき…

RjpWiki (REvolutionR は連邦の新型か) を読んでちょっと触ってみました。 環境は、CPU: Core 2 Duo E8400 (3.0 GHz), Memory: 3.3 GB, Windows XP 32bit SP3。 R-benchmark-25.R を5回ずつ実行してみた。 R-org: R-2.10.0 デフォルト ATLAS: ATLAS (C2D, CR…

パッケージ情報が定義してあるPACKAGESファイルと、必要なパッケージファイルを同じディレクトリに配置しておきます。 上のファイル群は、同じタイミングでダウンロードしておかないと、バージョンがずれて失敗するので最悪手で修正します。 Rってパッケージ…

あまり好きなグラフではないのですが、axis 関数を使って各時点のサンプルサイズを付加する事ができるようなので、やってみました。 Rには gplots::plotmeans と Rcmdr::plotMeans という関数がありますが、一部のグラフオプションが使えなかったりするので…

新年初投稿、ちょっと調べ物をしたのでメモ。確率変数 $X$ の期待値 $\theta=\mathrm{E}(X)$ について、モンテカルロ法による近似を考える。 何らかの方法で、$X$ と同じ確率分布に従う $n$ 個の独立な乱数列 $X_i$ ($i=1, 2, \ldots, i, \ldots, n$) が得ら…

Lemieux C. Monte Carlo and Quasi-Monte Carlo Sampling (Springer Series in Statistics). Springer, 2009. とりあえずこれを読みました。 Chapter 4 の Variance reduction techniques がまとまっていて読みやすかったです。 quasi-Monte Carlo methods …

やろうと思っていること 研究 線型モデル (大分すすんだ) 特殊なパラメータ変換をしたときの conditional exact test (積多項分布モデル) リサンプリング・ブートストラップを使った方法 ベイズ流の解析方法、MCMC EMアルゴリズム ネットワークアルゴリズム …

平岡 和幸, 堀 玄. プログラミングのための線形代数. オーム社 2004. 自分が読んだ線型代数の本でピカイチ。 これを読んで、行列の変換としての意味が完璧に理解できた。 学部1年では東京大学出版会の線型代数入門 (緑のやつ, 解析学は青) が使われていた。 …

パッケージをもりもり作っていたら、 * checking Rd cross-references ...が表示されたまま止まってしまう。 自分の環境だと何故か Sys.setenv( ) 関数を .Rprofile に書いておいても働いてくれない、これも不思議。 しょうがないので check.pl の240行目ぐ…

ご要望があったので、機能追加版を作ってみました。 プロットに用いたデータを、グラフを描画したシートに移動するように変更 箱の幅を指定できるように変更 箱の幅指定に関しては、Excel 2007の致命的なバグで、書式設定画面から幅を指定できないため、VBA…

Rcmdr package自体のインストールは、car packageを追加した上で、「ローカルから...」を選べばよいので簡単でした。 しかし、Rcmdrの初回起動時に、依存関係があるパッケージをディレクトリを指定して一括でインストールする画面がでます。 ここで、 must i…

せっかくなので触ってみました。 data logistic; call streaminit(88745623); a = 0; b = 1; do i = 1 to 200; do x = 0 to 1; logit = a +b * x; p = exp(logit) / (1 + exp(logit)); y = rand("Bernoulli", p); output; end; end; keep x y; run; ods list…

意外に色々あるんですね。 ods listing close; ods pdf file = "c:\out.pdf" style = Journal; ods graphics on; proc xxxx; run; ods graphics off; ods pdf close; ods listing; その他使いそうな形式。 ods csvall file = "c:\out.csv"; ods rtf file = "…

09/11/29版「W32TeX」のdvipdfmxの挙動について (TeX Q&A) 私もほぼ同じ症状が出てしまいました。 数式の一部と、文章中の斜体 (要は斜体の一部?) が歯抜けになって表示できない状態でした。 斜体のアルファベットとギリシャ文字は歯抜けして、subscriptの大…

関数を作ってみました。 HS.tapply <- function(x, index, type) { # create index indicate matrix index <- as.numeric(factor(index)) mindex <- NULL for(i in 1:length(unique(index))) { mindex <- cbind(mindex, index-(i-1)) } if (type == 5 || typ…

うーん、こりゃ参ったねｗ n <- 20 rep <- 100000 d <- matrix(rnorm(n*rep), nrow=n, ncol=rep ) # apply Rprof() varA <- apply(d, 2, var) Rprof(NULL) print(a <- summaryRprof()$sampling.time) # matrix operation (1) Rprof() varB <- NULL for(i in …

前回書いた結果をもとにして、ちゃんと近似できているのか試してみることにしました。 今回は Welch (1938) の (2) 式の統計量 \[ v=\frac{\bar{X}_1-\bar{X}_2}{\sqrt{\frac{V_1}{n_1}+\frac{V_2}{n_2}}} \] を導出する過程で出てくる分布の近似を考えてみ…

[1] の説明が分からなかったので、考えたり調べたりしていたら [2] を見つけることが出来た。 ある尤度関数 $p(y \mid \theta)$ に対して、無情報事前分布を構成する事を考える。 どんなパラメータ $h(\theta)$ に対して一様な事前分布を考えるかによって、$…

Ntzoufras I. Bayesian Modeling Using WinBUGS (Wiley Series in Computational Statistics). New York: John Wiley & Sons, Inc. 2009 ぱらぱらっと読ませてもらいました。 うーん、欲しいけど高い。