![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
【親ページ >>> Ｑ＆Ａ (初級者コース)】
初心者のための R および RjpWiki に関する質問コーナー
新規投稿はできません。
過去の記事のアーカイブ
H.I. (2009-04-23 (木) 18:55:52)
初心者です。
MacOS10.5、X11でR2.9.0を動かしています。
グラフの中に日本語の文字を入れることができません。
この日本語化掲示板の2009年2月12日の回答を参考にやってみましたがだめでした。
お教えいただけると助かります。sessionInfo()
R version 2.9.0 (2009-04-17)
i386-apple-darwin8.11.1
locale:
ja_JP.UTF-8/ja_JP.UTF-8/C/C/ja_JP.UTF-8/ja_JP.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
if (.Platform$pkgType == "mac.binary"){
options(device="quartz")
}
get("familyset_hook", pos="MacJapanEnv")
をやると，
function() { if(names(dev.cur())=="quartz") par(family="sans")}
が表示されるんですよね。 -- 2009-04-25 (土) 13:51:40頭バクハツです、、、もぅ、、、 (2009-04-23 (木) 06:34:03)
非常に初歩的な質問かもしれないのですが、
私にとってはとても難しく、悩んでも答えが出ませんでした。
どなたか、ぜひ、ご教示ください。
以下のようなデータフレームＸＸがあります。
変数Ａ，Ｂを用いて回帰回帰分析を行い、最終的に変数Ｘの予測値、PREDICT.Xを予想することを目的とします。
これだけであればよかったのですが、次のことをする必要に迫られています。（データは簡略化し、かつ適当な数字です。）
1958年とそれ以降のデータを分け、別個に回帰する。
あるいは、変数label内で、同じラベルを張られているものは、同じグループとして回帰する。
（1958年以前のラベルはすべて１ですので、上記の2文は全く同じグループ分けを指します。）data.frame XX year label X A B PREDICT.X 1951 1 NA 1 191 ? 1952 1 NA 2 9 ? 1953 1 NA 3 29 ? 1954 1 NA 4 349 ? 1955 1 4 5 2 ? 1956 1 3 5 9 ? 1957 1 NA 6 3 ? 1958 1 NA 6 495 ? 1959 2 2 7 692 ? 1960 2 3 8 94 ? 1961 2 3 9 95 ? 1962 2 NA 7 93 ? 1963 2 NA 6 7 ? 1964 2 NA 5 3 ? 1965 2 4 8 3 ? 1966 2 3 4 2 ?この為に、以下のようにコードしました。
if(XX$year<1958){ model.X<-lm(PREDICT.X~A+B) pred.X <-predict(model.X,XX)}が、エラー”条件が2以上になったので、最初の1つを使用しました”となってしまいます。
どのようにコードすれば、データフレームＸＸに予測値PREDICT.Xを計算することが出来るのでしょうか。
どうか、よろしくお願いします。
kuri (2009-04-22 (水) 17:36:13)
状況；パッケージplsgenomicsのpls.ldaで、クラス分類されたデータに対し て、モデルを構築しようとしております。
問題点；cross validationでモデルを構築したいのですが、
pls.lda(Xtrain, Ytrain, Xtest=NULL, ncomp＝1:20, nruncv=20, alpha=9/10, priors=NULL)
で、実行したところ同じデータセットを用いても、出力結果が二回すれば二回
とも異なることになってしまいます。
考えられる原因；cross validation を行うときにtest setとtraing setの選び方がランダムに選ばれている可能性があるのではないかと思います。
質問内容；結果を何回行っても同じようにするためにはどのようにcrossvalidationを行うべきでしょうか？教えて頂きたいと思います。
どうぞよろしくお願いいたします。
Wanko (2009-04-22 (水) 08:16:10)
R version 2.8.1 (2008-12-22)
i386-pc-mingw32
locale:
LC_COLLATE=Japanese_Japan.932;LC_CTYPE=Japanese_Japan.932;LC_MONETARY=Japanese_Japan.932;LC_NUMERIC=C;LC_TIME=Japanese_Japan.932
attached base packages:
[1] stats graphics grDevices utils datasets methods base
今日は、散布図のプロット座標変換の質問をさせて下さい。
たとえば
plot(x,y) として描いたプロットのX,Y軸のデータ範囲が
X: 最大値 8 最小値 -6
Y: 最大値 3 最小値 -2
となっているとします。
これを X, Y軸とも
最大値 -1 最小値 +1 の範囲に成るよう座標点のX,Y値を変更したいのです。
今は、数行のコードで変換していますが、Rにはそのような関数が用意されて
いないでしょうか？
よろしくお願いします。
ちっぽけな自分 (2009-04-21 (火) 23:01:16)
以下のように、functionの引数（数値型）として、例えばA,Bを使用しました。
function内でグラフを作成し、そのタイトルに数値型の引数A,Bを利用して
"AからBまで"と出力したいとします。
どのようにコードすればよいのでしょうか？？？？
長く時間をかけるのですが、未だに解決できません。
どなたか、解決策を教えて頂けませんか？。
よろしくお願いします。
example <- function(A,B){
kaiki<-lm(...)
plot(kaiki...)
title(AからＢまで)←←こう出力したいのです。
}
example(1,2) →→→→ グラフのタイトルを ”1から２まで”と出力したい。
example(3,4) →→→→ グラフのタイトルを ”3から4まで”と出力したい。
MH (2009-04-19 (日) 22:21:28)
初めて書かせて頂きます。書き方や書く場所が適切でなかったら、申し訳ありません。chisq.test()について、3つほど質問があります。非常に初歩的（かつ舌足らず）な質問で恐縮です。
1) 2x2よりも大きい分割表（例えば、2x3）を分析する際には、引数correctをTにしてもFにしても、同じ結果がかえってきます。また、Tにした場合も、出力の1行目に"with Yates' continuity correction"という文言が出ないため、イェーツの補正がかかっていないということでしょうか？
2) Rでのイェーツの補正は、X-squared=0の場合（実測値と期待値が全く同じ場合）にも補正をし、その結果としてx-squaredが0よりも大きい値となるようですが、これは統計的に言って妥当なのでしょうか？
3) chisq.test()を実行すると、出力の1行目が"Pearson's Chi-squared test"となる場合と"Chi-squared test for given probabilities"となる場合があります。これらはどのように違うのでしょうか？ また、後者はどのようなときに出力されるのでしょうか？
chisq.test(c(2,3,1,2,3,4)) chisq.test(matrix(c(2,3,1,2,3,4), 2, 3))
Hiro (2009-04-19 (日) 10:31:45)
WinXP pro SP2、R 2.9.0 for Winです。
RでEPSファイルを作成する時にTeXコマンドを埋め込み、LaTeXでPSfrag.styを経由して文書を作成しています。2.8.1までは、次のようなスクリプトを作り、LaTeXファイル側でpsfragscanonコマンドを使うことで、EPS作成時に埋め込んだTeXコマンドは問題なくインライン展開されていました。
ps.options(family="ComputerModern") plot(sin, main="\\tex[][]{\\LARGE Curve of $\\sin(x)$ and $\\cos(x)$}", xlab="\\tex[][]{\\large $x$}", ylab="\\tex[][]{\\large $\\sin(x)$ and $\\cos(x)$}") curve(cos, add=T) dev.copy2eps(file="sin.eps")
EPSファイルの中のTeXコマンド部分は次のようになっていたので、TeXも理解できていたようです。
265.27 467.43 (\\tex[][]{\\LARGE Curve of $\\sin\(x\)$ and $\\cos\(x\)$}) .5 0 0 t
しかし、同じ事を2.9.0で実行すると、出来上がるEPSファイルの中で、TeXコマンドに対応する文字列が、以下のように無残にぶった切られ、TeXコンパイルも当然の如く失敗します。
73.82 455.43 (\\tex[][]{\\LAR) 0 ta -0.434 (GE Curv) tb -0.434 (e of $\\sin\(x\)$ and $\\cos\(x\)$}) tb gr
このEPSファイルを、次のLaTeXファイルにかけると、
\documentclass{jarticle}
\usepackage{psfrag,graphicx}
\begin{document}
\begin{center}
\psfragscanon
\includegraphics[width=10cm]{sin.eps}
\end{center}
\end{document}
! Undefined control sequence.
<argument> \LAR
) 0 ta -0.434 (GE Curv) tb -0.434 (e of $\sin (x)$ and $\cos...
l.6 \includegraphics[width=10cm]{sin.eps}
と叱られます。
2.9.0のChangelogやNEWSを見ても、この問題に対応する仕様変更は見つけられませんでした。対応方法をご存知の方はいらっしゃいませんでしょうか？
yishii (2009-04-19 (日) 00:53:45)
VistsでR2.8.1を使っていましたが、今度R2.9.0をインストールしました。ところが、packagesのインストールができなくなりました。
環境は以下のとおり。sessionInfo()
R version 2.9.0 (2009-04-17)
i386-pc-mingw32
locale:
LC_COLLATE=Japanese_Japan.932;LC_CTYPE=Japanese_Japan.932;LC_MONETARY=Japanese_Japan.932;LC_NUMERIC=C;LC_TIME=Japanese_Japan.932
ローカルにダウンロードしたpackageAAAは2.9用はもちろん2.8用でも以下の警告が出るので困っています。
どうすればいいかどうかよろしくご教授ください。
install.packages(choose.files("", filters = Filters[c("zip", 中で警告がありました：'lib = "C:\Users\本人名\Local Settings\Application Data\R-core\R/R/win-library/2.9"' は書き込み可能ではありません 以下にエラー install.packages(choose.files("", filters = Filters[c("zip", : パッケージをインストール出来ませんでした 。
nanbuwks (2009-04-15 (水) 08:12:01)
boxplotについて、教えてください。
http://www.okada.jp.org/RWiki/?%A5%B0%A5%E9%A5%D5%A5%A3%A5%C3%A5%AF%A5%B9%BB%B2%B9%CD%BC%C2%CE%E3%BD%B8%A1%A7%C8%A2%B7%BF%BF%DE
に
「上下のひげは、それぞれ上(下)側四分位数の位置から、極値までの間に引かれます。極値とは、上(下)側四分位数から四分位範囲の 1.5 倍以内にあるデータのうちの最大(小)値です。極値よりも大きい、または小さい値は外れ値としてひげの先にプロットされます。」
と説明があります。
例えば、100人分の100点満点のテストのデータがあるとすると、データ値順に25人目の値が90点、50人目の値が60点、 75人目の値が50点とすると、極値とは105点〜45点の中にある最大最小値ということになるのでしょうか？
どうしてもわかりません。。 (2009-04-10 (金) 07:49:40)
R2.8.1を使っています。
OSはWindowsXPです。
例えば、以下のコードでエクセルを読み込み出力してみました。
しかし、途中までは読み込めているのですが、
途中からデータが存在するにも関わらずNAが表示されてしまいます。
何が考えられる原因なのでしょうか。
教えて頂けると助かります。
よろしくお願い致します。
library(RODBC)
original <- odbcConnectExcel("えくせる.xls")
sqlTables(original)
X<-sqlQuery(original,"select*from [シート１$]" )
odbcClose(X)
元データ（例）
A B C D E F G
1975 1 1 1744 85.24 164 8
1975 2 1 1618 79.08 187 9
1975 3 1 1651 80.69 170 8
1975 4 1 1572 76.83 136 6
1975 5 1 1603 78.35 102 4出力結果
A B C D E F G
1975 1 1 1744 85.24 NA NA
1975 2 1 1618 79.08 NA NA
1975 3 1 1651 80.69 NA NA
1975 5 1 1603 78.35 NA NA
library(RODBC)
original <- odbcConnectExcel("c:/ABCDEFG.xls")
sqlTables(original)
X<-sqlQuery(original,"select*from [ABCDEFG$]" )
Y<-X #Xは参照渡しなのでＹに値渡しする
odbcClose(original)
Yビクス (2009-04-08 (水) 20:21:24)
こんにちは。
pairwise.t.test()の引数pool.sdについて伺いたいことがあります。
これはpool.sd=Tとすると各群の標準偏差が等しいことを意味し、
つまり等分散を仮定することと同義だと思ったのですが、
等分散を仮定したt.test()と補正無し、pool.sd=Tのpairwise.t.test()の
p-valueが一致しません。たとえば以下のような3群xを検定します。
x <- c(c(1:20),c(21:40),c(41:60)) group <- c(rep(1,20),rep(2,20),rep(3,20)) t.test(x[1:40]~group[1:40], var=T) Two Sample t-test data: x[1:40] by group[1:40] t = -10.6904, df = 38, p-value = 5.168e-13 pairwise.t.test(x,group, pool.sd=T, p.adj="none") Pairwise comparisons using t tests with pooled SD data: x and group 1 2 2 3.1e-15 - 3 < 2e-16 3.1e-15
このようにt.test()では1群と2群はp-value=5.168e-13、pairwise.t.test()
では3.1e-15となります。
pool.sd=Tは等分散が等しいことを意味していないのでしょうか？
ヘルプを見ても意味が理解できません。そもそも、なぜvar=Tという引数が
ないのでしょうか？
宜しくお願いします。
> set.seed(141421356)
> x <- c(rnorm(20), rnorm(20, mean=0.5), rnorm(20, mean=1))
> group <- c(rep(1,20),rep(2,20),rep(3,20))
> t.test(x[1:40]~group[1:40], var.equal=TRUE)
Two Sample t-test
data: x[1:40] by group[1:40]
t = -2.4653, df = 38, p-value = 0.01832
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.3690472 -0.1344442
sample estimates:
mean in group 1 mean in group 2
0.06037503 0.81212077
> pairwise.t.test(x,group, pool.sd=FALSE, p.adj="none", var.equal=TRUE)
Pairwise comparisons using t tests with non-pooled SD
data: x and group
1 2
2 0.018 -
3 0.046 0.839
P value adjustment method: noneど素人 (2009-04-06 (月) 10:35:31)
バリオグラムを描きたいのですが、自分のデータをエクセルからcsvファイルに変換してバリうグラムを描きたいのですがどうしたらよいのでしょうか？
libraryに入っているデータならなんとかなるのですが、新しいデータを用いて活用したいのです。
サザエ (2009-04-04 (土) 02:21:50)
頑張ってRのこと、勉強しております。
質問です。ACF(自己相関係数)の計算結果をExcelで読めるようにファイルに移したいんです。そのやり方を教えてください。
よろしくお願いします。
迷い猫 (2009-04-03 (金) 21:47:08)
全くの初心者です。
やれることはやったのですが、解決策が分からないのでどうか、どうか、よろしくお願いします。
データフレームとして、Excelから読み込んだJapan, U.S, U.Kが準備されています。
必要に応じたデータフレームを使用したいので、ＳＡＳで行うマクロのような関数を作りたいと思い、以下のようなプログラムを作りました。
byage <- function(country,){
country <- substitute(country)
kakunin<-country
kakunin
if(country==japan){data=Japan} #引数にjapanを渡した時は、データフレームJapanを使いたいと思ったので。
ifelse(country==us){data=U.S}
else(country==uk){data=U.K}
data
}
byage(japan)
これを実行すると、次のエラーが返ってきます。
「以下にエラー byage(japan) : オブジェクト "japan" は存在しません」
何が悪くて、どう修正すればよいのか、さっぱり分かりません。。。
また、function関数に、データフレームを直接引数として渡す方法は存在するのでしょうか。
どうか、よろしくお願いします。
> func <- function(z) # 渡されたデータフレーム z の先頭6行を表示する関数
+ {
+ print(head(z))
+ }
> x <- data.frame(a=1:10, b=letters[1:10])
> y <- data.frame(x=rnorm(20), y=runif(20))
> func(x) # データフレーム x を渡す
a b
1 1 a
2 2 b
3 3 c
4 4 d
5 5 e
6 6 f
> func(y) # データフレーム y を渡す
x y
1 2.15205491 0.0295724
2 -1.32527541 0.3563355
3 0.17068233 0.2448835
4 -0.84898390 0.3565478
5 -0.01846693 0.5706420
6 1.78453651 0.4487693レッドクリフ (2009-04-03 (金) 08:49:56)
system関数の不具合？で困っています．
通常はsystem関数で実行した結果を返すのに、Ｒでメモリを大量に消費しているときに、コマンドの結果がエラーで返ってきます．> sessionInfo() R version 2.8.0 (2008-10-20) x86_64-unknown-linux-gnu ６４ビットのLinux機、メモリは８Ｇ積んでいます．具体的な症状は以下の通りです．
- 成功する例
> gc() used (Mb) gc trigger (Mb) max used (Mb) Ncells 111691 6 350000 18.7 350000 18.7 Vcells 120708 1 786432 6.0 456124 3.5 > x <- system("ls") Makefile blaslapack.o fblaswr.h work.0.txt workL2.txt workS3.txt > x [1] 0- 失敗する例
> dim(model_mat) [1] 383907 1886 > gc() used (Mb) gc trigger (Mb) max used (Mb) Ncells 503372 26.9 597831 32.0 506595 27.1 Vcells 737629983 5627.7 904433606 6900.3 748064327 5707.3 > x <- system("ls") > x [1] -1
- このように５Ｇのメモリを消費してsystem関数を使用した場合にうまく動きません．
- オブジェクトを消去して再度実行した場合
> rm(list=c("model_mat")) > gc() used (Mb) gc trigger (Mb) max used (Mb) Ncells 117672 6.3 597831 32.0 506595 27.1 Vcells 12803200 97.7 723546884 5520.3 748064327 5707.3 > x <- system("ls") Makefile blaslapack.o fblaswr.h work.0.txt workL2.txt workS3.txt > x [1] 0
- このようにうまく動きます．
- メモリを大量に消費してsystem関数にintern引数を与えた場合
> dim(model_mat) [1] 383907 1886 > gc() used (Mb) gc trigger (Mb) max used (Mb) Ncells 503477 26.9 667722 35.7 506705 27.1 Vcells 737630000 5627.7 820415969 6259.3 748064344 5707.3 > x <- system("ls", intern=T) *** caught segfault *** address (nil), cause 'memory not mapped' Traceback: 1: system("ls", intern = T) Possible actions: 1: abort (with core dump, if enabled) 2: normal R exit 3: exit R without saving workspace 4: exit R saving workspace Selection:
- セグメンテーションフォルトで落ちます・・・
- このような現象の回避方法はございませんでしょうか．よろしくお願いいたします．
> sort(sapply(ls(), function(x)object.size(get(x))))
a b c
48 104 8000200
> gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 181395 9.7 407500 21.8 350000 18.7
Vcells 1132760 8.7 1912818 14.6 1133141 8.7
> save(c,file="c.dat")
> rm(c)
> gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 181391 9.7 407500 21.8 350000 18.7
Vcells 132759 1.1 786432 6.0 1133141 8.7
> system("ls")
> load("c.dat")#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include <R.h>
#include <Rdefines.h>
// exec vfork
SEXP system_vfork(SEXP exCmd)
{
pid_t pid;
SEXP ret;
int iRet;
char *cmd_for_execv[] = {"/bin/bash", "-c", NULL, NULL};
char *cmdStr;
if( !isString(exCmd) ){
Rprintf("not CHAR.");
return R_NilValue;
}
cmdStr = (char *)CHAR(VECTOR_ELT(exCmd, 0));
printf("cmd str : %s\n", cmdStr);
cmd_for_execv[2] = cmdStr;
if( (pid=vfork()) < 0 ){
/* fork 失敗 */
printf("fail to fork\n");
return R_NilValue;
}else if( pid == 0 ){
/* 新しく起動した子プロセスに処理させたい内容 */
//sleep (3);
//iRet = execlp("sleep", "sleep", "5", (char *)NULL);
iRet = execv(cmd_for_execv[0], cmd_for_execv);
//sleep (30);
printf("ret value : %d\n", ret); /* (actually never reached)*/
exit (EXIT_SUCCESS); /* 子プロセス終了*/
}
int status = 0;
pid_t wait_pid;
Rprintf ("parents, child is %d\n", pid);
wait_pid = wait(&status);
if (wait_pid == -1){
err (EXIT_FAILURE, "wait error");
//exit(-1);
return R_NilValue;
}
//Rprintf("input char > ");
//scanf("%c", &c);
return R_NilValue;
}$ gcc -std=gnu99 -I/usr/local/lib64/R/include -I/usr/local/include \
-fpic -g -O2 -c fork.c -o fork.o
$ gcc -std=gnu99 -shared -L/usr/local/lib64 -o fork.so fork.o \
-L/usr/local/lib64/R/lib -lRdyn.load("fork.so")
is.loaded("system_vfork")
res <- .Call("system_vfork", "ls")
res <- .Call("system_vfork", "ls -la")
res <- .Call("system_vfork", "sleep 3; ls -la")初心者 (2009-04-03 (金) 07:25:25)
Windows Vista で R2WinBUGS を使用しております．Rのバージョンは2.7.2です．
R2WinBUGS を動かした場合，WinBUGSから画像ファイルやテキストファイルを
出力するため，Rを「管理者権限として実行」して起動する必要が
ありますが，バッチ実行する場合に「管理者権限として実行」する
方法が分かりませんでした…．
'C:\Program Files\R\R-2.7.2\bin\R.exe' --no-restore --no-save < 'C\temp\test.R'
のような命令をコマンドプロンプトから起動しているのですが，
何かオプション等はありますか？どうかよろしくお願い致します．
kf (2009-03-22 (日) 23:00:27)
サンプリング間隔が不等な時系列データを、内挿等によって等間隔のデータを
作成するにはどうすればよいのでしょうか。お教え下さい。
kushi (2009-03-22 (日) 00:38:04)
はじめて質問させて頂きます。よろしくお願い致します。
実行環境は以下でした。
Windows Vista, 32bit, 4GB RAM
R version 2.8.1 (2008-12-22)
なお，Rのショートカットは以下のように設定しています。
"C:\Program Files\R\R-2.8.1\bin\Rgui.exe" --max-mem-size=2047M
私は文系の博士課程大学院生です。それほど詳しいPCの知識はありません。
研究の必要上，数千行×数千列〜数万行×数万列の大きい行列に対して，特異値分解を行う必要が出てきました。
Rで実行しようとしたところ，以下のスクリプトを実行すると，以下のようなエラーがでます。> y = NULL > y = matrix(0.0, 5000, 4000) > svd_y = NULL > svd_y = svd(y) エラー： サイズ 366.4 Mb のベクトルを割り当てることができません > gc() used (Mb) gc trigger (Mb) max used (Mb) Ncells 114305 3.1 350000 9.4 350000 9.4 Vcells 20076602 153.2 90632978 691.5 172087953 1313.0
質問は，以下の３つです。それぞれの質問について，下で補足します。
１．「数千行×数千列」の行列をRで特異値分解する方法はないか。
２．「数万行×数万列」の行列をRで特異値分解する方法はないか。
３．svd()関数がどのようなスクリプトによって定義されているか見る方法はないか。
「メモリ」と「svd()」をキーワードに過去ログやWEBを探しました。
質問１は，OSやメモリの扱いによる解決を期待した質問です。
質問２は，速度を犠牲にしてもスクリプトがストップせずに特異値分解を実行できるアルゴリズムや方策を期待した質問です。
質問３は，自分で特異値分解のスクリプトをかけばもっと省メモリでやれるかもしれないという期待からお尋ねしています。
何卒よろしくお願いいたします。
> y = NULL
> y = matrix(0.0, 5000, 4000)
> svd_y = NULL
> svd_y = svd(y)
> gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 240738 6.5 467875 12.5 350000 9.4
Vcells 56183965 428.7 159348538 1215.8 148191804 1130.7すみません，Macのファイルのようで，どういったファイルなのか確認できませんでした…。
＞ 解釈
お答えになっているかわからないのですが…。
潜在意味インデキシング，もしくは潜在意味解析という手法を用いたいと考えています。つまり，コーパスに基づいた数千×数千〜数万×数万の語句‐文書の共起行列を，SVDで分解された３つの行列のうち，特異値を元の数から300次元程度に減らして，もう一度行列を復元することで，人間の知識モデルとしたいのです。
どうぞよろしくお願いいたします。-- kushi 2009-03-22 (日) 22:21:53
カーネルおじさん (2009-03-21 (土) 12:23:44)
いつもお世話になります．確率を与えてブートストラッピングを行い，時系列を発生させることを行ってます．
単純化して，1本のサンプルパスが3期(i=3)で，10,000本(j=10000)のパスを発生させるとします．
再試ができるように，set.seedで特定のseedを与えますが，今回はj回ごとにseedを変動させることを考えました．
ところが，はじめはプログラム1のようにset.seed(a+j-1)の位置をfor(i)の次に書きました．
結果をみるとバイアスが大きく，x行列は各列同じu値しか出現していないので変だと思いました．
プログラム2のようにfor(j)とfor(j)の間にset.seed(a+j-1)の位置を変更したら，期待通りになりました．
初歩的な質問ですが，
1)set.seed(a+j-1)の位置はどちらでもいいのではと思いますが，Ｒはプログラム1をどう読んでいるのでしょうか？
2)たとえばset.seed(100）と入力すると，以降のsampleは，他のseedに変えるか，Rを再起動しない限り，
ずっと100のseedで乱数を発生し続けるという理解でよろしいですか？
(そう観察したので，j回ごとseedを変えた方がランダム性が高まると考えました）
3)seedの種類はこの事例のように40,000個もあるのでしょうか？
2),3)はhelp(set.seed)や乱数TIP大全を見てもわからなかったので，質問させていただきます．よろしくお願いします．
下記の確率設定ですと，すべてのパスの合計値の期待値は-120になりますが，
それぞれのaに数値(1,10001,20001,30001)を代入したところ合計値は以下の通りとなって相違があります．
プログラム１では，#1〜10000=-124.08, #10001〜20000=-121.98, #20001〜30000=-115.20, #30001〜40000=-113.31
プログラム２では，#1〜10000=-121.43, #10001〜20000=-120.52, #20001〜30000=-115.89, #30001〜40000=-116.48
プログラム１
x<-matrix(0,nr=3,nc=10000)
p<-c(0.3,0.5,0.2)
u<-c(-0.02,0,0.01)
a<-1
for (j in 1:10000){
for (i in 1:3){
set.seed(a+j-1)
x[i,j]<-sample(u,1,prob=p,replace=T)
}
}
y<-sum(x)
y
プログラム2
x<-matrix(0,nr=3,nc=10000)
p<-c(0.3,0.5,0.2)
u<-c(-0.02,0,0.01)
a<-1
for (j in 1:10000){
set.seed(a+j-1)
for (i in 1:3){
x[i,j]<-sample(u,1,prob=p,replace=T)
}
}
y<-sum(x)
y
Wanko (2009-03-20 (金) 14:35:48)
始めて投稿させていただきます。
R version 2.8.1 (2008-12-22)
i386-pc-mingw32
locale:
LC_COLLATE=Japanese_Japan.932;LC_CTYPE=Japanese_Japan.932;LC_MONETARY=Japanese_Japan.932;LC_NUMERIC=C;LC_TIME=Japanese_Japan.932
attached base packages:
[1] stats graphics grDevices utils datasets methods base
Rは使い始めて間もないです。
WinXP で VBAでの連携処理しています。
こちらの金先生の資料を参考に
http://www1.doshisha.ac.jp/~mjin/R/200808_61.pdf
下記のコマンドを実行して
kekkaDF<-NgramDF(targetText, type=1, N=3, pos="名詞")
sortlist<-order(kekkaDF[,3],decreasing = TRUE)
fwn<-kekkaDF[sortlist,]
wng<-graph.data.frame(fwn[1:100,])
tkplot(wng, vertex.label=V(wng)$name, layout=layout.fruchterman.reingold, vertex.size=1)
図３ 福田総理の所信表明演説文における名詞の共起ネットワークマップ
の様なネットワーク図を出力しています。
質問は
・「Graph plot 1」などのウインドウサイズは事前或いは事後に変更するのは
可能でしょうか？ → ソースでは450×450で固定のように思えますが
ユーザー関数に作り直す？事無く、変更できませんか？
・ウインドウ背景色が灰色で定年前の年寄りには見にくいです。
→ 白色とか薄いベージュなどに変更したいのですが？
・VBAからはCall oShell.Run(Environ("COMSPEC") & " /c " & sRpath & " --no-save < " _ & txtR_workfolder & "\RelatedTerm.R", vbNormalFocus, False)sRpath：C:\Program Files\R\R-2.8.1\bin\R.exe
と言った形で行っていますが、実行後にDOSプロンプトが閉じてしまい、グラフも
同時に消え去ります。
現在、コマンド最後に Sys.sleep(1800) などとしていています。
もっとましな方法がありそうですが
宜しくご教授ください。
初学者 (2009-03-16 (月) 17:31:58)
12万件のデータが、csv形式であります。
そのファイルをインポートしようと下記のプログラムを実行しました。
RDT1<-read.table("ローデータ.csv", header=T, sep=",", dec=".", nrows=65536, skip=1, na.strings="")
すると、以下のようなエラーメッセージが出ます。以下にエラー scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings, : '15' 行目には，74 個の要素がありません
変数の数は74ですが、分岐している質問があるため、null値が非常に多く、それが悪さをしているのかと推測し、エクセルでnull値を-99に置き換え、インポートしてみると読み込まれました。
しかしながら、データファイル取得ごとに、これらの処理をしなければならないので、前処理に時間を要してしまいます。
質問としては、上記のようなnull値が多いファイルを、そのまま読み込む良い方法、または違ったプログラムオプションが、Rで用意されているのか知りたいと思い、質問させて頂きました。
以上、ご教授の程よろしくお願いします。
nanbuwks (2009-03-12 (木) 19:14:12)
アドオンパッケージ chplot の chplot ライブラリを使ってみたのですが、色の指定が思うようにできません。library(chplot) data(iris) chplot(Petal.Length ~ Petal.Width | Species, data = iris, legend = list(cex = 0.6),plot.points=TRUE,pch = 18, cex = 0.5, col=c("black","orange","yellow"))としてみましたが、凸包輪郭線の色のみ反映し、観測点の色が変更できません。
どのようにしたらいいでしょうか?
u-kun (2009-03-12 (木) 15:42:40)
latticeパッケージのlevelplotを使ってグラフを作成しようとしています。
格子状の測定点での測定値の分布を可視化したいのです。
しかし、z軸に対応するデータが10^(-10)から10^(-3)まで幅があるので、
colorkeyの色の変化(目盛り)を対数にしたいのですがうまくいきません。
例えばplotであればlog="x"などとすればx軸が対数になるように、
colorkeyを対数にするにはどうすれば良いのでしょうか？
> range.z <- pretty(range(z)) > levelplot(z, colorkey = list(labels = list(at = range.z, labels = log (range.z) ) )
> log.z <- pretty(range(log10(z))) > levelplot(log10(z) ~ x * y,col.regions=topo.colors(16), colorkey = list(labels = list(at = log.z, labels = log.z)))おかげさまで、ほぼ私のイメージ通りのグラフが描けました。
> (label.z <- paste(10, log.z, sep = "^")) > levelplot(log10(z) ~ x * y, colorkey = list(labels = list(at = log.z, labels = substitute(label.z, list(label.z = label.z))) ) )leveplot()で実際に描画を行っているのは panel.levelplot()関数（そしてdraw.colorkey()関数）なので，このヘルプを参考にしてください．
nanbuwks (2009-03-12 (木) 11:49:19)
x <- c(0,1,1,2,2,3,5) boxplot(x)とすると、箱ひげ図とその範囲外の観測点が表示されますが、この観測点を表示させない方法ってあるのでしょうか?
カーネルおじさん (2009-03-11 (水) 10:32:21)
お世話になります．Ｒの基本的なことがまだまだわかっておりません
しかし，Ｒを楽しんでます．
今般，i行のそれまでの行の累乗積の総和をi行に返すプログラムに苦労しています．
つまり，行列Xから行列Ｚの変換で，新しいi行は，
z[i,]=x[i,]*x[i-1,]+x[i,]*x[i-1,]x[i-2,]+...+x[i,]*x[i-1,]x[i-2,]...x[1,]
を計算する必要に迫られました（Excelでやってますが，2GMパソコンが悲鳴をあげています．それにできればRで自己完結したいです）
具体例では行列Ｘが以下の通りだとして，
> x
[,1] [,2] [,3]
[1,] 2 3 1
[2,] 3 2 1
[3,] 4 1 5
[4,] 2 4 3
#y...z
[,1] [,2] [,3]
[1,] 2 3 1
[2,] 6 6 1
[3,] 12+24 2+6 5+5
[4,] 8+24+48 4+8+24 15+15+15
という計算をへて次のＺ行列が求まります．
z
[,1] [,2] [,3]
[1,] 2 3 1
[2,] 6 6 5
[3,] 36 8 10
[4,] 80 36 45
Ｘを単純な累乗積Yに変換するプログラムは,以下の通りに書けましたが，
i行から逆向きに累乗積を計算するのが難しく，
どうもうまくいきません．プログラムの外側にもうひとつFor文(j in 1:4)を挟み，
初期値i=1(i<2の箇所）をi=jなりに変え，行列Ｙで得られた行をＺに落していけば
ワークしそうなのですが，実力不足を実感しております．
> y
[,1] [,2] [,3]
[1,] 2 3 1
[2,] 6 6 1
[3,] 24 6 5
[4,] 48 24 15
#行列xの各行の累乗積を行列yへ返す
x<-matrix(c(2,3,4,2,3,2,1,4,1,1,5,3),nc=3)
y<-matrix(nr=4,nc=3)
for (i in 1:4)
{ if (i<2) {y[i,]<-x[1,]}
else
{y[i,]<-y[i-1,]*x[i,]}
}
> (x <- matrix(c(2, 3, 4, 2, 3, 2, 1, 4, 1, 1, 5, 3), ncol=3))
[,1] [,2] [,3]
[1,] 2 3 1
[2,] 3 2 1
[3,] 4 1 5
[4,] 2 4 3
> y <- matrix(0, nrow=nrow(x), ncol=ncol(x))
> for(j in 2:nrow(x)) {
+ n <- matrix(0, nrow=nrow(x), ncol=ncol(x))
+ n[j-1,] <- x[j-1,]
+ for(i in j:nrow(x)) n[i,] <- n[i-1,] * x[i,]
+ if(j != 2) n[j-1,] <- 0
+ y <- y + n
+ }
> y
[,1] [,2] [,3]
[1,] 2 3 1
[2,] 6 6 1
[3,] 36 8 10
[4,] 80 36 45> x <- matrix(c(2, 3, 4, 2, 3, 2, 1, 4, 1, 1, 5, 3), ncol=3)
> y <- matrix(0, nrow=nrow(x), ncol=ncol(x))
> for(i in 2:nrow(x)) y[i,] <- (y[i-1,] + x[i-1,]) * x[i,]
> y[1,] <- x[1,]
> y
[,1] [,2] [,3]
[1,] 2 3 1
[2,] 6 6 1
[3,] 36 8 10
[4,] 80 36 45> a <- apply(x, 2, function(d) sapply(1:4, function(i) sum(cumprod(rev(d[1:i]))[2:i])))
> a[1,] <- x[1,]
> a
[,1] [,2] [,3]
[1,] 2 3 1
[2,] 6 6 1
[3,] 36 8 10
[4,] 80 36 45> bx <- function(x,i) { ifelse(i==1,x[1],ifelse(i==2,x[1]*x[2],(x[i-1]+bx(x,i-1))*x[i])) }
> z <- apply(x,2,bx,1)
> for (i in 2:nrow(x)) z <- rbind(z,apply(x,2,bx,i))
> z> bx <- function(x,i) { ifelse(i==1,x[1],x[i]*(bx(x,i-1)+1)) }
> z <- apply(x,2,bx,1)
> for (i in 2:nrow(x)) z <- rbind(z,apply(x,2,bx,i))
> znanbuwks (2009-03-08 (日) 15:40:16)
scatterplotのオプションを調べたいと思い、> help(scatterplot)しましたが、
No documentation for 'scatterplot' in specified packages and libraries:you could try 'help.search("scatterplot")'
となります。helpやマニュアルを調べるにはどこを参照すればいいでしょうか？
なお、sessionInfo()はR version 2.6.2 (2008-02-08) i486-pc-linux-gnu locale: LC_CTYPE=ja_JP.UTF-8;LC_NUMERIC=C;LC_TIME=ja_JP.UTF-8;LC_COLLATE=ja_JP.UTF-8; LC_MONETARY=ja_JP.UTF-8;LC_MESSAGES=ja_JP.UTF-8;LC_PAPER=ja_JP.UTF-8;LC_NAME=C; LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=ja_JP.UTF-8;LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods baseです。ちなみに、2変数の箱ヒゲ図を描きたいと考えています。scatterplot では軸に箱ヒゲ図
の表示ができますが、これを散布図を描かずに箱ヒゲ図のみxy座表内に描くことができないか
と考えました。
よろしくお願いします。
nami (2009-03-05 (木) 22:40:10)
一定時時間ごとに計測したデータの日平均、最高、最低を得たいのですがうまくいきません。例えばdatAというオブジェクトの中身がDate Time A B C 2009/1/1 06:00:00 1 15 200 2009/1/1 12:00:00 2 20 300 2009/1/1 18:00:00 3 25 400 2009/1/2 00:00:00 4 30 500 2009/1/2 06:00:00 5 35 600 2009/1/2 12:00:00 6 40 700 2009/1/2 18:00:00 7 45 800 2009/1/4 06:00:00 12 70 1300 2009/1/4 12:00:00 13 75 1400・・・
のようになっているとします(たまに欠測があったりします)。これにcombine.time <- paste(datA$Date, datA$Time) date.time <- as.POSIXct(combine.time) datA["Date2"] <- date.timeのようにして日付型データを組み込むまではできたのですが、このデータの日付ごとの平均値、最高値、最低値のデータセットのつくり方がわかりません。最終的には時系列のグラフを作りたいと考えてております。どなたかお教えいただけ無いでしょうか。
環境はWindowsXP sp3、R2.8.1です。よろしくお願いします。
threenum <- function(x){
summary(x)[c(1,4,6)]
}
tapply(datA$C,
list(format(datA$Date2, "%Y-%m-%d")),
threenum)datB <- tapply(datA$C, list(format(datA$Date2, "%Y-%m-%d")),threenum)として、
datB$ObsNumber <- 1:3
matplot(datB$ObsNumber,datB[1])
以下にエラー matplot(datB$ObsNumber, datB[1]) :
'x' と 'y' は同じ数の行を持たなければなりません
と出ます。なお、
datB $`2009-01-01` Min. Mean Max. 200 300 400 $`2009-01-02` Min. Mean Max. 500 650 800 $`2009-01-04` Min. Mean Max. 1300 1350 1400 $ObsNumber [1] 1 2 3となっています。ご指導ください。
datB <- tapply(datA$C, list(datA$Date),threenum)
datB <- as.data.frame(do.call("rbind",datB))
datB["Date"] <- as.Date(row.names(datB))
row.names(datB) <- c(1:nrow(datB))
names(datB) <- c("Min","Mean","Max","Date")
attach(datB)
days <- c(Date[nrow(datB)]-Date[1])
plot(Date,Mean,type="b",xaxt="n",ylim=c(0,1500))
axis.Date(1, at=seq(Date[1],Date[days],"1 days"), format="%m月 %d日",las=2)
segments(Date,Min,Date,Max)
detach()
というように、かなーり強引な感じがします(データフレームに変換してしまったあたり)。それと、最初に組み込んだDate2ですが、あまり必要なさそうなので結局使いませんでした。言いだしっぺということで、ご参考までに。より簡潔で美しいコードをご教示いただければ、もちろん乗り換えます(笑)。みなさま、よろしくお願いします。
るる (2009-03-05 (木) 19:42:16)
R初心者です。自分の解析結果を入れたファイルが容量が多く、１１９行あって上の１：４７が見れません。
それで、
＞ファイル名[1:47]
と打ち込んだら、部分代入不可というエラーが出ました。
部分代入不可ならどうやって、４７行目までの入るの内容を読めばよいのでしょうか。
カーネルおじさん (2009-03-04 (水) 08:18:39)
Ｒ初心者です．お世話になります．約200個の3次元ベクトル標本から，kde3d(in library misc3d)を使って，カーネル密度dをcsvファイルで出力しました．ただ，出力結果を見ると，グリット点(50x50x50)に対応する密度の数値が表示されており，関連のhelpを読んだところでは，標本ベクトル値に対応する密度数値は，出力できないようです．
おそらく，グリット点にいちばん近いところの密度を読み取ることになると思いますが，大変そうです．
密度dのExcelスプレッドシート出力を見ると，x軸は縦に50行ですが，横はy,z軸混ぜて，50x50=2500列から構成されています．このままではExcelのvlookupが使えず，なんらかの工夫をしない限り，200回，目で探してくることになりそうです．
それに，腕力でやるよりできるだけ，Rの中で自己完結したいという思いもあり，投稿いたしました．
質問をまとめますと，
1.kde3dは標本ベクトル値に対応するカーネル密度値は返さないのでしょうか?
2.グリット点にいちばん近いところを探してくるRでの方法のヒントをいただければ幸いです．
(なお，将来の構想として，標本値のカーネル密度をウエイトにして，ブートストラッピングを行いたいので，
標本の密度が必要という事情があります．長文になってすみません)
R1 (2009-03-02 (月) 21:27:51)
お世話になります。A,B,Cの3群がそれぞれx,yという値をもっており、mfrow=c(1,3)で連続してstripchart()でY軸方向に散布図を書きます。
ここでY軸の単位が共通なので、BとCのグラフの軸名と目盛りを省略し、
さらに余白を調節して３つのグラフを連結しようと考えました。
つまり以下のように入力します。
F <- factor(c("X","Y")) x1 <- rep(1:6,5) y1 <- rep(7:12,5) x2 <- rep(9:14,5) y2 <- rep(15:20,5) x3 <- rep(21:26,5) y3 <- rep(27:32,5) par(mfrow=c(1,3), mar = c(5, 4, 4, 0)) COL <- c("red","blue") stripchart(rbind(x1,y1) ~ F, vertical = T,col=COL, xlab="A") par(mar = c(5, 0, 4, 0)) stripchart(rbind(x2,y2) ~ F, vertical = T, yaxt = "n", col=COL,xlab="B",ylab="") par(mar = c(5, 0, 4, 2)) stripchart(rbind(x3,y3) ~ F, vertical = T, yaxt = "n", col=COL,xlab="C",ylab="")しかし、中央のBのグラフのx軸が長くなるのと、各グラフのx軸の両端に
値がプロットされてしまいなってしまい不格好です。
x軸の長さをそろえ、かつx軸のプロットする位置を自由に指定することは
できないでしょうか？また、ここでは３つのグラフを書きましたが、
グラフの数に依存せずに調節したいと考えています。
どなたかご教授宜しくお願いします。
TM (2009-02-22 (日) 16:33:45)
超初心者です。
マトリックスの相関係数からクラスター分析をしたいと考えています。相関係数rをd=2(1-r)の距離に直したマトリックスをエクセルで作成しています。既に距離になっているはずなのですが、これをクラスター分析しようと、入力したところ、
data1<-read.table("ox008.csv",header=TRUE,sep=",")
data1
data1.hc<-hclust(data1)
plot(data1.hc)
data1<-read.table("ox008.csv",header=TRUE,sep=",")
data1
X A B C D1 A 0 1 2 2
2 B 1 0 2 2
3 C 2 2 0 2
4 D 2 2 2 0data1.hc<-hclust(data1)
以下にエラー if (n < 2) stop("must have n >= 2 objects to cluster") : 引数の長さが0ですplot(data1.hc)
以下にエラー plot(data1.hc) : オブジェクト "data1.hc" は存在しませんのようなエラーメッセージが出ます。
何が問題なのでしょうか？
ge (2009-02-20 (金) 10:13:35)
コンソール画面でタイプミスをすると、エラーのダイアログが出てきますが、かなり反応（数秒以上）が遅いです。
設定か何か工夫して反応を早くできないものでしょうか。
R version 2.8.1 (2008-12-12)
i386-pc-mingw32
です。
plott(0,1)
"Error: could not find function "plott"
sh (2009-02-18 (水) 12:20:14)
plclustのHeight値の範囲を変更することは可能でしょうか？ 例えば、以下の例でHeightの目盛の上限を200や300に変更することは可能でしょうか？hc <- hclust(dist(USArrests), "ave") plclust(hc)この場合、Heightの最大値は
max(hc$height) [1] 152.314Heightの目盛の下限を20に変更するには
plclust(hc, hmin=20)
> plot(as.dendrogram(hc), ylim=c(0,300)) 以下にエラー plot.default(0, xlim = xlim, ylim = ylim, type = "n", xlab = xlab, : 仮引数 ”ylim” が複数の実引数にマッチしました
junkoda (2009-02-17 (火) 08:26:06)
はじめまして。
太陽質量 (solar mass) の記号（latex で書くと M_\odot、M の添字に丸書いて点）をグラフのラベルに書きたいのですが、出来ますでしょうか。
一般化すると、Unicode や Harshey Vector Font に含まれるけれども plotmath にはない記号を添字として表示するよい方法はないものでしょうか。
添字は plotmath で書けます。> plot(1, xlab=expression(M[O]))太陽記号は Hershey Vector Font で書けます
> par(family="HersheySerif") > plot(1, xlab="M\\SO")しかし、Hershey と plotmath が並用できません。
太陽記号は unicode にはありますが（"\u2299"）、Adobe Symbol にはありません。
最終的には eps ファイルとして欲しいので、> plot(1, xlab=expression(M["\u2299"])) > dev.copy2eps(file="out.eps")などと、強引にやってみますと
Warning messages: 1: In dev.copy(file = "out.eps", device = function (file = ifelse(onefile, : font metrics unknown for Unicode character U+2299以下、10個の警告が出されて、太陽記号は ... に変換されてしまいました。
助言などありましたら助かります。最後に sessionInfo を書いておきます。
R version 2.8.1 (2008-12-22)
i686-pc-linux-gnu
locale:
LC_CTYPE=ja_JP.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=C;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C
#!/bin/bash
TMPTEX=$$
TEX=/usr/bin/latex
DVIPS=/usr/bin/dvips
cd /tmp
cat <<EOF > $TMPTEX.tex
\documentclass{article}
\usepackage[dvips]{graphicx}
\usepackage[scanall]{psfrag}
\begin{document}
\thispagestyle{empty}
\includegraphics[scale=.8]
EOF
echo {$1} >> $TMPTEX.tex
cat <<EOF >> $TMPTEX.tex
\end{document}
EOF
$TEX $TMPTEX
$DVIPS -E $TMPTEX
echo You got /tmp/$TMPTEX.ps
rm $TMPTEX.tex $TMPTEX.dvi $TMPTEX.log $TMPTEX.aux $TMPTEX.pfg
cd -
使い方は、\tex{}を埋め込んだRグラフィックのepsをこのスクリプトの引数に指定します。初学者 (2009-02-16 (月) 14:21:57)
はじめまして。お世話になります。
使用環境は、Windows XPで、RのVersionは、下記の通りです。
R version 2.8.1 (2008-12-22)
i386-pc-mingw32
アンケートに逆転項目があるため、回答の数値を逆するため、下記のようなプログラムを記述しましたが、問題がありエラーになりました。
［プログラム例］
回答者１〜10の方のＱ１のデータを、１は７に、２は６に、・・・７は１に変換するため、下記のプログラムを実行しました。（下記のデータでは、６や７の回答はありませんが、実データには含まれている可能性があるので、実データ６や７に対しても、データ変換を試みています）
sid<-c(1,2,3,4,5,6,7,8,9,10)
q01<-c(2,4,5,2,1,3,2,3,4,5)
RDT<-data.frame(SID=sid, Q01=q01)
RDT$Q01a<-NA
if(RDT$Q01 == 1) {RDT$Q01a <- 7}
if(RDT$Q01 == 2) {RDT$Q01a <- 6}
if(RDT$Q01 == 3) {RDT$Q01a <- 5}
if(RDT$Q01 == 4) {RDT$Q01a <- 4}
if(RDT$Q01 == 5) {RDT$Q01a <- 3}
if(RDT$Q01 == 6) {RDT$Q01a <- 2}
if(RDT$Q01 == 7) {RDT$Q01a <- 1}
以上、アドバイスよろしくお願いいたします。
> sid <- 1:10 > q01 <- c(2, 4, 5, 2, 1, 3, 2, 3, 4, 5) > RDT <- data.frame(SID=sid, Q01=q01) > > RDT$Q01a <- NA > RDT$Q01a[RDT$Q01 == 1] <- 7 > RDT$Q01a[RDT$Q01 == 2] <- 6 > RDT$Q01a[RDT$Q01 == 3] <- 5 > RDT$Q01a[RDT$Q01 == 4] <- 4 > RDT$Q01a[RDT$Q01 == 5] <- 3 > RDT$Q01a[RDT$Q01 == 6] <- 2 > RDT$Q01a[RDT$Q01 == 7] <- 1 > RDT SID Q01 Q01a 1 1 2 6 2 2 4 4 3 3 5 3 4 4 2 6 5 5 1 7 6 6 3 5 7 7 2 6 8 8 3 5 9 9 4 4 10 10 5 3
capetapeta (2009-02-14 (土) 10:06:54)
Macintosh、Windows、Linux（Win,Linuxは仮想マシン）いずれにしてもRを何回ダウンロードしてもlinesという命令文を書いても正しくQuartz2にUPされないんです。以前にも色々試してnew plot? abline?
という投稿をしましたが、結局線がプロットできないという結論に至りました。これを解決するためには何か有効な手だてはあるのでしょうか？或は線がプロットできるようになるようなパッケージがあるのでしょうか？
独学中 (2009-02-14 (土) 09:15:58)
いろいろ検索してみましたが、いい案がなく、どなたか助言をお願いします。
やりたいこと
測定データのファイルが、あるディレクトリに入るのをwhileで待って、
入ってくれば、統計処理させたい。
測定は、約３分ごと。
データファイルの大きさは、70KBで、３６本ほど。
困っていること
一時間ほど、走らせると、Rがハングしている。
使用環境 winXP, R ver2.8.1 memory.limitは、1500MB
以下の工夫をしているが、効果なし。
変数は、使い終われば、rm(),gc(),gc().
sinkで、出力先を変更
Rconsoleが、メッセージで一杯かと思い、１０分ほどで、ESCしてみるが、
メッセージは、表示されない。（何事もなく、＞を表示。）
> while(1) {
if (file.exists("test.R")) {
print("Oh!, I got it")
# do somethig
file.remove("test.R")
}
}R初心者7 (2009-02-08 (日) 22:42:22)
はじめまして。お世話になります。
使用環境は、WindowsXPで、Rのパッケージは下記のようになっています。
R version 2.6.2 (2008-02-08)
i386-pc-mingw32
locale:
LC_COLLATE=Japanese_Japan.932;LC_CTYPE=Japanese_Japan.932;LC_MONETARY=Japanese_Japan.932;LC_NUMERIC=C;LC_TIME=Japanese_Japan.932
attached base packages:
[1] stats graphics grDevices utils datasets methods base
私は、樹木の生存時間に効く要因を求めようとしています。
複数の要因を組み合わせたモデルを作成し、どのモデルが生存を良く説明しているのか、AICを使ってモデル選択をしようとしています。
例えば下記のようなデータでコックス回帰し、AICを求めようとすると、
testdata
time event density height light
5 1 3.05 48 2.7
5 0 6.36 44 6.3
3 0 5.38 107 6.9
8 0 2.57 39 5.8
2 1 0.38 150 2.7
9 1 1.99 194 3.1
#density, height, lightが要因
library(KMsurv)
要求されたパッケージ survival をロード中です 要求されたパッケージ splines をロード中ですrequire(survival)
data <- read.delim("testdata.txt",T)
model.01 <- coxph( Surv (time, event) ~ density + height , data=data)
AIC(model.01)
以下にエラー UseMethod("logLik") : "logLik" に適用可能なメソッドがありません
となり、AICが求められません。
どなたか、解決方法をご教示頂ければ幸いです。
ちなみにstepを使えば、AICが一応出てくることは確認しています。
必要な情報を失念している場合もあるかと思いますので、その場合はご指摘頂ければ幸いです。よろしくおねがいします。
> library(survival) > data <- data.frame( + time = c(5, 5, 3, 8, 2, 9), + event = c(1, 0, 0, 0, 1, 1), + density = c(3.05, 6.36, 5.38, 2.57, 0.38, 1.99), + height = c(48, 44, 107, 39, 150, 194), + light = c(2.7, 6.3, 6.9, 5.8, 2.7, 3.1) + ) > model.01 <- coxph( Surv (time, event) ~ density, data=data) > model.02 <- coxph( Surv (time, event) ~ density + height, data=data) > model.03 <- coxph( Surv (time, event) ~ density + height + light, data=data) Warning message: In fitter(X, Y, strats, offset, init, control, weights = weights, : Ran out of iterations and did not converge > extractAIC(model.01) [1] 1.000000 6.075746 > extractAIC(model.02) [1] 2.000000 7.280983 > extractAIC(model.03) [1] 3 6
hoso (2009-02-05 (木) 16:13:39)
データとして(X,Y)のペアが100個ぐらいあり、高次多項式Y=f(X) を求めてカーブフィッティングしたいと思っています。（データの性質にもよると思いますが、割と滑らかで高次多項式でのフィッティングが妥当と思わせるデータです。）これだけだったらエクセルでも可能です。ところがエクセルでやりますと、高次項の係数が2e-9などと桁落ちになっております。それらを用いて計算すると、Xの値を大きくして高次項が利き始めると大きなエラーとなってしまいます。係数が2.3462e-9などと精度が高い場合、エラーが小さくなるものと思います。エクセルの分析ツール＋回帰分析を行ってみると直線近似でした。高次の式（６次ぐらい）でもエクセルの分析ツールで処理できるのでしょうか。
このような問題をＲで実行するとどのようになるでしょうか。具体的なやり方（データの読み込みは分かりますが）を教えて頂けると助かります。あるいは関連した実例集を指示して頂くのも有難いです。
よろしくお願いします。
d <- data.frame(
x = c(1.6, 2.7, 4.2, 5.6, 6.1, 6.7, 7.9, 8.8, 9.3, 9.7),
y = c(6.9, 6.1, 5.6, 6.3, 7.2, 8.2, 7.9, 8.5, 8.4, 8.0))
plot(y~x, data=d, ylim=c(1,9))
a <- lm(y~x+I(x^2)+I(x^3)+I(x^4)+I(x^5)+I(x^6)+I(x^7)+I(x^8), data=d)
x2 <- seq(1, 10, length=1000)
d2 <- data.frame(x=x2)
y2 = predict(a, newdata=d2)
lines(y2~x2)＜ふ＞ (2009-02-02 (月) 13:52:56)
R2.6.1を使っています。
グラフを書かせるにあたって、maiで余白を設定し（c(1,1,1,1)）ちょうどいい具合になったので、PDFに出力しようとしたところ、maiの設定が反映されませんした。mfrow=c(2,2)のような設定もききません。
PDFに出力する場合の、このようなパラメータの設定はどこでやればいいでしょうか。
よろしくお願いいたします。
新参 (2009-02-02 (月) 12:01:40)
お世話になります。
R2.8.1、XPにパッケージ(Rcmdr）をインストールした時、下記のエラーが出てしまいインストールができません。既存のＱ＆Ａを見ると、Ｒの再インストールやミラーサイトの変更などが出てきたため試してみましたが、未だに解決されません。
また、zipファイルをローカルに落として、"ローカルにあるパッケージからインストール"を選んでも同様のエラーが出ます。
またRcmdrだけではなく、他のパッケージでも（例：glmmMLやfarawayなど）同様のエラーが出ますが、大部分のパッケージでは問題なくミラーサイトからインストールできます。
回復方法を教えていただけると助かります。よろしくお願いします。
utils:::menuInstallPkgs()
--- このセッションで使うために、CRANのミラーサイトを選んでください --- URL 'http://cran.md.tsukuba.ac.jp/bin/windows/contrib/2.8/Rcmdr_1.4-7.zip' を試していますContent type 'application/zip' length 2389722 bytes (2.3 Mb)
開かれた URLdownloaded 2.3 Mb
パッケージ 'Rcmdr' は無事に開封され、MD5 サムもチェックされました 以下にエラー normalizePath(path) : path[1]: 指定されたファイルが見つかりません。
library(Rcmdr)
以下にエラー library(Rcmdr) : 'Rcmdr' という名前のパッケージはありません。
パッケージのインストールでつまづくと、作業が進まず非常に困っています。
どなかたアドバイスいただけないでしょうか。お願いします。
capetapeta (2009-02-01 (日) 21:14:59)
x <- c(1,3,4,6,8,9,12)
y <- c(5,8,6,10,9,13,12)
par(mfrow = c(1,3))
plot(x,y,ylim = c(0,16))
abline(0,0.68)
以下にエラー int_abline(a = a, b = b, h = h, v = v, untf = untf, ...) : まだ plot.new が呼ばれていません
というエラーが発生しました。質問項目やヘルプを探しましたがエラーを改善する方法を見つける事ができませんでした。昨日Rを始めたばかりで右も左もわからない状態です。宜しくお願いいたします。
初心者 (2009-01-28 (水) 14:27:40)
こんにちわ。
「pgmm」を用いてダイナミックパネルデータ分析を行い、summaryで結果を確認すると、R-squareが表示されませんでした。
いろいろ自分で試行錯誤してみましたが、わからなかったので、ご存知の方がいらっしゃいましたら、やり方を教えてください。
Vista使い (2009-01-27 (火) 00:51:14)
Vista版RでPDFを作成するときの日本語フォントはどうすればよいのでしょう？
setHook(packageEvent("grDevices", "onLoad"),function(...) grDevices::ps.options(family="Japan1"))
ですが、"Japan1"でも、"Japan1"のところをいろいろ変えても文字化けしてしまいます。
どうかコメントおねがいします。
setHook(packageEvent("grDevices", "onLoad"),function(...)
grDevices::ps.options(family="Japan1"))
setHook(packageEvent("grDevices", "onLoad"),function(...)
grDevices::pdf.options(family="Japan1"))
とすることで、文字化けしないPDFが出来ました。以前だったらps.optionsだけ設定していればPDFの方にも反映されてたのですが、Vistaはpdf.optionsも設定しないといけないんですね。なかまさま、どうもありがとうございました。-- Vista使い 2009-01-28 (水) 09:16:53
kato (2009-01-26 (月) 11:11:44)
こんにちは．
step関数を用いて重回帰分析を行った後，得られたそれぞれの説明変数の非標準化係数の95％信頼区間と決定係数(R2)，F値を得たいのですが，それらを調べるには，どのような関数を使ったら良いのでしょうか．
検索で調べてみましたが，95％信頼区間はグラフ上に示す方法などは見つかりましたが，数値として出力する（上限と下限）ような記載は見つかりませんでした．
例えば以下のようなデータで，glm()後にstep()を行うと，data
身長 体重 年齢1 179 82 26
2 170 67 31
3 188 89 36
4 180 88 39
5 183 89 35
6 181 94 26
7 184 83 27
summary(result.glm)
Call:
glm(formula = 体重 ~ 身長, data = data)
Deviance Residuals:1 2 3 4 5 6 7
- 0.4665 -4.4154 -4.5175 4.3056 1.6220 9.0777 -5.6059
Coefficients:
Estimate Std. Error t value Pr(>|t|)(Intercept) -137.3270 78.1518 -1.757 0.1392
身長 1.2279 0.4323 2.840 0.0362 *
- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 35.02454)
Null deviance: 457.71 on 6 degrees of freedomResidual deviance: 175.12 on 5 degrees of freedom
AIC: 48.402
Number of Fisher Scoring iterations: 2
となりますが，その後にどのようにして95％信頼限界，R2，F値を求めたらよいかわかりません．
すみませんが，よろしくお願いします．
初心者 (2009-01-25 (日) 20:43:44)
R で長いプログラムを書いております．
Rのコンソールの内容（命令＋実行結果＋エラーメッセージ＋警告メッセージ）
をファイルに保存するにはどうすればよいのでしょうか？
例えば，sinkだと実行結果しか保存できません．教えていただければと思います．
ちなみに R 2.8.1, OS は Ubuntulinux を使っております．
nanbuwks (2009-01-24 (土) 11:53:03)
こんにちは。plot(1:10,type="l",ylim=c(0,10))とすると、y軸の目盛が0,2,4,6,8,10と振られますが、y軸の底(x軸と交わる点)は0になりません。
大体プロット領域は-0.5から10.5ぐらいになっているようです。これを0から10までにするには
どうしたらいいでしょうか。
なお、barplotと重ね合わせをする目的なので、x軸を上にずらす以外のやりかたで方法がないでしょうか。
plot(1:10,yaxs="i")plotを重ねるならpoints()とか
barplot(matrix(1:10,2)) points(1:5)検索すればいろいろでてくると思います。デモデータとかあると試せますが… -- akira 2009-01-25 (日) 10:20:05
(2009-01-23 (金) 15:17:52)
SVMを使いたいのですが、「エラー：サイズ 276.9 Mb のベクトルを割り当てることができません」となってしまい、実行できません。
PCのメモリは2GBで、57600×101(ヘッダー含む)のデータを判別したいのです。
メモリを変える以外になにか解決方法ありませんでしょうか？
‘--max-mem-size=N’ (Windows only) Specify a limit for the amount of memory to be used both for R objects and working areas. This is set by default to the smaller of 1.5Gb and the amount of physical RAM in the machine, and must be between 32Mb and the maximum allowed on that version of Windows.
"C:\Program Files\R\R-2.8.0\bin\Rgui.exe" --max-mem-size=1800Mとすることで、使用できるメモリが増えました．ただしXPでは１プロセスで使用できるメモリが２Ｇまでなので、ご注意を．
> x <- matrix(0.0,57600,101*20)
エラー： サイズ 887.7 Mb のベクトルを割り当てることができません
> gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 105955 2.9 350000 9.4 350000 9.4
Vcells 74584 0.6 786432 6.0 353875 2.7> x <- matrix(0.0,57600,101*20)
> gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 105950 2.9 350000 9.4 350000 9.4
Vcells 116426576 888.3 128697668 981.9 116426905 888.3初学者 (2009-01-21 (水) 18:25:08)
はじめまして。お世話になります。
使用環境は、Windows XPで、RのVersionは、下記の通りです。
R version 2.8.1 (2008-12-22)
i386-pc-mingw32
locale:
LC_COLLATE=Japanese_Japan.932;LC_CTYPE=Japanese_Japan.932;LC_MONETARY=Japanese_Japan.932;LC_NUMERIC=C;LC_TIME=Japanese_Japan.932
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] cluster_1.11.12 MASS_7.2-45
loaded via a namespace (and not attached):
[1] tools_2.8.1
質問ですが、クラスター分析後に樹形図を作成しましたが、（縦型）樹形図の最下部にある表示が数字で表示されています。その部分を、たとえば、各サンプルの名前など、データに存在するサンプル名に表示できないでしょうか？
たとえば、下記のプログラム例ですと、（縦型）樹形図の最下部に数字が表示されます。irisデータで、Speciesの表記、setosa、versicolor、virginicaを表示させることは可能でしょうか？
※実際には、「単一の名前」を用いたデータを利用しますが、説明のためirisデータにより投稿いたしました。
（プログラム例）
DST.d <- dist(iris)
DST.hc1<-hclust(DST.d, method="ward")
plot(DST.hc1, hang=-1, main="ウォード法")
以上、アドバイスの程よろしくお願いいたします。
DST.d <- dist(iris) DST.hc1<-hclust(DST.d, method="ward") plot(DST.hc1, hang=-1, main="ウォード法",labels=iris$Species)
初心です (2009-01-19 (月) 23:30:10)
最小全域木のようなグラフを描きたいです
matplot(x,type="p",pch=20)で点を書きsegments(7,70,8,78)で一本づつ線を描けばかけるのですが、数値を変更すると線は書き直しになり、細かく線を書いていくと時間がかかります。何かアドバイスを頂ければ幸いです。よろしくお願いします。
初級者 (2009-01-18 (日) 19:17:48)
お世話になります。
R2.8.1、XPにパッケージ(tree）をインストールした時、下記のの用になりましたが、初心者で内容や、回復方法が分かりません。既存のＱ＆Ａを見ましたが、見落としを含めて、見当たりません。回復方法を教えていただけると助かります。よろしくお願いします。install.packages("tree")
--- このセッションで使うために、CRANのミラーサイトを選んでください --- URL 'http://cran.md.tsukuba.ac.jp/bin/windows/contrib/2.8/tree_1.0-26.zip' を試していますContent type 'application/zip' length 152558 bytes (148 Kb)
開かれた URLdownloaded 148 Kb
パッケージ 'tree' は無事に開封され、MD5 サムもチェックされました 以下にエラー normalizePath(path) : path[1]: 指定されたファイルが見つかりません。library(tree)
以下にエラー library(tree) : 'tree' という名前のパッケージはありません
Rin (2009-01-16 (金) 15:50:03)
はじめまして。
プログラムについての質問です。
φを平均-1，分散1の正規分布に従う確率変数，ψを平均2，分散1の正規分布に従う確率変数としてf(x,y)=0.6φ(x,y)+0.4ψ(x,y)~に従う2変量混合正規乱数を発生させたいのですが、どのようにすればよいでしょうか？
現在、自分が書いたものではないのですが、他の方が書いたプログラムとして、２変量正規乱数を発生させるための以下のようなプログラムがあります。# nc:発生させる乱数のデータ数~ #r:相関係数~ gendat <- function(nc, r) ~ {~ # 仮のデータ行列を作る。この時点では変数間の相関は近似的に0である。 z <- matrix(rnorm(2*nc), ncol=2) # 主成分分析を行い，主成分得点を求める。この時点で変数間の相関は完全に0である。 res <- eigen(r2 <- cor(z)) coeff <- solve(r2) %*% (sqrt(matrix(res$values, 2, 2, byrow=T))*res$vectors) z <- t((t(z)-apply(z, 2, mean))/sqrt(apply(z, 2, var)*(nc-1)/nc)) %*% coeff # コレスキー分解の結果をもとに，指定された相関係数行列 を持つように主成分得点を変換する。 z %*% chol(matrix(c(1, r, r, 1), ncol=2)) }~このプログラムを、2変量混合正規乱数を発生させるためのプログラムにしようと思い、改変したものが以下のプログラムです。
gendat2 <- function(nc, r)~ { z <- matrix((rnorm(2*nc)-1)+3*(runif(2*nc) > 0.6), ncol=2)#~ res <- eigen(r2 <- cor(z)) coeff <- solve(r2) %*% (sqrt(matrix(res$values, 2, 2, byrow=TRUE))*res$vectors) z <- t((t(z)-apply(z, 2, mean))/sqrt(apply(z, 2, var)*(nc-1)/nc)) %*% coeff z %*% chol(matrix(c(1, r, r, 1), ncol=2)) }
１行目のzの部分を変えてみました。
z <- matrix((rnorm(2*nc)-1)+3*(runif(2*nc) > 0.6)
の部分は、rjpwikiの『ヒストグラムと密度の推定』の、『f(x) に従う乱数を生成する関数』の項目を参考にさせて頂きました。
しかし結果を見てみると、どうもf(x,y)=0.6φ(x,y)+0.4ψ(x,y)に従うような乱数が発生していないようです。
コレスキー分解の部分を理解していないので、きっとその辺りに原因があると思うのですが、このソースで間違えている部分が分かる方教えてください、お願いします。
Saito (2009-01-14 (水) 23:49:00)
いつもお世話になっております。
ある点の座標とその座標からある直線までの距離がわかっているとき、ある点から直線に垂直に降ろした線分が直線と交わる点の座標を求めるにはどうすればよいでしょうか？以下にサンプルを示します。x <- 0;y <- sqrt(2)#座標の作成 plot(x, y, xlim=c(0,3), ylim=c(0,3), col="red")#座標の表示 abline(a=0, b=1)#直線式の定義(ここではy=x) d <- 1#座標から直線までの距離(既知)このようなとき、つまり(x,y)=(0,sqrt(2))という座標からy=xまでの垂直距離(今回は距離d=1)がわかっているときの交点の座標を求めたいのです。ちなみに答えは(sqrt(2)/2,sqrt(2)/2)になるはずです。もちろん手計算でやればでますが、どのようにプログラミングしてよいのかがわかりません。
最初にoptimを試して垂直距離が最小になるような座標を推定したところ、大体うまく推定できるようなのですがたまに初期値が悪いとうまく推定できないようです。solveで連立方程式が解けるらしいので勉強したのですがどうもよくわからず使いこなせません。
どなたか交点の座標が分かる方がおりましたらご教授願えませんでしょうか。
なお、環境はWindowsXp、R-2.8.1です。
よろしくお願い致します。
x <- 0 ; y <- sqrt(2)/2 #座標の作成
plot(x, y, xlim=c(0,3), ylim=c(0,3), col="red") #座標の表示
intercept <- 0 ; slope <- 1 # 直線の定義 (y = x)
abline(a = intercept, b = slope) # 直線の表示
slope2 <- -1/slope # 直交する直線の傾き
intercept2 <- -slope2 * x + y # 直交する直線の切片
abline(b=slope2, a=intercept2) # 直交する直線の表示もしておくか
cross.p <- solve( # 連立方程式の解(交点座標) = これが求めたい解
matrix(
c(-slope, -slope2, 1, 1), 2, 2),
matrix(
c(intercept, intercept2)))[,1]
(d <- sqrt(sum((cross.p - c(x, y))^2))) # ついでに距離の計算ちなみに例題の解は1にはならないような．．．
NewComer (2009-01-14 (水) 14:18:19)
R2.7.0、ＸＰで、アッソシエーション分析をしていました。１月になって
arulesを再インストール（全パケージをインストールしたため）したところ
listからtransactionの処理が出来なくなりました、本を見ながらあちこち試しましたがＮＧ、友人にも友人のＰＣで試してもらいましたがＮＧ、初心者で
分かりません。以下、データを添付しますので、よろしくお願いします。＃テスト用に作成したデータ~ install.packages("arules",dependencies=TRUE) library(arules) {data<-list(c("Ａ","Ｂ","Ｃ","Ｄ"), c("Ｅ","Ａ","Ｇ","Ｂ"))} data.tran<-as(data,"transactions") class(data.tran) ＃結果 > {data<-list(c("Ａ","Ｂ","Ｃ","Ｄ"), + c("Ｅ","Ａ","Ｇ","Ｂ"))} > data.tran<-as(data,"transactions") 以下にエラー validObject(.Object) : 不正なクラス "ngCMatrix" オブジェクト： row indices are not sorted within columns ←エラーの意味がわからない。＃昨年末arulesは改版されているので、これと関係あるのだろうか？
Package: arules
Version: 0.6-7
Date: 2008-13-12
Title: Mining Association Rules and Frequent Itemsets
#library(arules)を実行すると、下記の様なメッセージが出るようになったと思うが
関係あるのだろうか？要求されたパッケージ Matrix をロード中です 要求されたパッケージ lattice をロード中です 次のパッケージを付け加えます: 'Matrix' The following object(s) are masked from package:stats : xtabs The following object(s) are masked from package:base : colMeans, colSums, rcond, rowMeans, rowSums 次のパッケージを付け加えます: 'arules' The following object(s) are masked from package:base : %in%
当方 (2009-01-11 (日) 19:28:21)
Linux の中には実行ファイルが用意されていないものがあります．その場合は Unix 版 R をインストールするのと同じ方法でソースからビルドするかと思います． 当方，VineLinux 4.2 を使っており，以下の手順でソースからのビルドを試みております．
$ apt-get update $ apt-get install gcc-g77 $ apt-get install readline $ apt-get install gcc $ apt-get install gcc-c++
$ ./configure $ make $ make install
configure: error: --with-readline=yes (default) and headers/libs are not available
./configure --with-readline=no --with-x=no
とすると，何とかインストールできました．ありがとうございます！ちなみに，「--with-readline=no」だけだと以下のエラーが出るので「--with-x=no」も付けました．-- 当方 2009-01-11 (日) 22:03:09
configure: error: --with-x=yes (default) and X11 headers/libs are not available
$ apt-get install XOrg-devel $ apt-get install readline-devel
Saito (2009-01-11 (日) 19:11:22)
いつもお世話になっております。
R2.8.1になって、編集→GUIプリファレンスと進みコンソール画面の背景色を変更するとRエディタの方でも背景色が自動で変更されるようなのですが、output text/User inputで入出力の色を変更してもRエディタの方は入力する文字の色が黒のままで変更されず、背景色との組み合わせによっては非常に見づらいです。何度もアンインストールとインストールと再起動をして確かめたのですがやはり直りません。どこかでRエディタの色を変えることができるのか調べたのですがわかりませんでした。このようなことになるのは私だけでしょうか？どなたか解決法をご存知でしたら教えていただけると幸いです。
なお使用環境はWindowsXpです。
EDI (2009-01-11 (日) 09:02:29)
Contnurグラフを作成しようと考えています。
xy平面上にplot()を用いて散布図を作成しました。
平面上に散布された点(因子)の密度を数値化し、その数値をz値に使用と考えています。
まだ初心者なので、調べ方も不十分なのかとは思いますが、多くのContnurグラフを作成するfunctionたとえばcontour()などはこちらで(x,y,z)を定義する必要があるので使えていません。
そこで、(x,y)を定義するだけで、私の考えているようなContnurグラフを作成させられる関数があれば、教えていただければ幸いです。
あるいは、xy平面上に散布されたの点の密度分布をすべての格子点で求められるような関数があれば、そこからz値を定義できるのでそういった関数でもかまいません。何かアドバイスいただけると幸いです。よろしくお願いします。
ま (2009-01-10 (土) 20:09:47)
glmで誤差分布をポアソン分布とし，一般化線形回帰モデルを使用しています。その結果として出てくるモデルの”回帰線”の95%信頼区間を求めたいのですが，どなたか良い方法を教えていただけれませんでしょうか？
参考として，疑似データでの解析したスクリプトを添付いたします。
#疑似データ作成
x1<- c(1:10)
y1<- c(0,0,5,15,10,30,40,200,290,1000)
plot(x1,y1,ylim=c(0,1500))
#作図（回帰直線）
m <- glm(y1~x1,family=poisson(link=log))
summary(m)
x.val <- seq(min(x1),max(x1),0.1)
lines(x.val,predict.glm(m,newdata=data.frame("x1"=x.val),type="response"),lwd=2)
data(cars) library(simpleboot) car.lm <- lm(dist~speed, data=cars)#線形モデルで回帰 car.boot <- lm.boot(car.lm, R =2000)#ブートストラップ plot(car.boot, xlab="speed", ylab="dist")ポアソンだとこれでは上手くいかないようですが、参考までに。-- Saito 2009-01-12 (月) 16:34:09
みち (2009-01-09 (金) 19:59:59)
Rでの実行結果をエクセルに書き出す方法を教えてください。
tonton (2009-01-08 (木) 22:38:15)
CRANのhttp://cran.r-project.org/web/packages/DPpackage/index.htmlでDPpackage_1.0-6.tar.gzをダウンロードしました。
これをRでインストールしようとしたのですが失敗します。
以下のようになります。
install.packages("C:/DPpackage_1.0-6.tar.gz",contriburl = NULL)
以下にエラー gzfile(file, "r") : コネクションを開くことができません 追加情報: Warning messages: 1: In zip.unpack(pkg, tmpDir) :~ zip ファイルから抽出中にエラー 1 が起こりました 2: In gzfile(file, "r") :~ 圧縮されたファイル 'DPpackage/DESCRIPTION' を開くことができません, 理由は 'No such file or directory' です
WindowsXPでR2.8.0を使用しています。
どなたか解決方法を教えて下さい。
okinawa (2009-01-08 (木) 17:40:48)
いつもお世話になっております。okinawaです。jaccard係数を用いたminhash法をＲで利用する方法をどなたかご存じないでしょうか？数式の解説やsasのサンプルは見つけたのですが、読みこなせなくて・・・。
環境は：windowsXPSP3,R2.7.2
ぽち (2009-01-06 (火) 08:59:17)
すみません。一通り探したのですがわかりませんでした。
以下のような流れでspreadsheetについて計算をしたいのです。これ自体は動くのですが、問題はこれを１０００回繰り返し、１０００個のAB値を求めたいのです。si <-data.frame(sample(data[2:11],10,replace=T))~ a <-data.frame(rowMeans(si[,1:5]))~ b <-data.frame(rowMeans(si[,6:10]))~ AB <-a/b~どなたかお教えいただけないでしょうか？
sapplyでできました。
sapply( 1:10 , function(i)
{si <-data.frame(sample(data[2:11],10,replace=T))
a <-data.frame(rowMeans(si[,1:5]))
b <-data.frame(rowMeans(si[,6:10]))
AB <-a/b} )
forの場合、どうやら変なことをしてしまっているようで動かないのですが、どなたか間違いを指摘していただけないでしょうか？
set.seed(101)
m <- 10
replicate10 <- numeric(m)
for (i in 1:m) replicate500[i] <-function(i)
{si <-data.frame(sample(data[2:11],10,replace=T))
a <-data.frame(rowMeans(si[,1:5]))
b <-data.frame(rowMeans(si[,6:10]))
AB <-a/b}
また、関数を定義しようとしてみたのですが、これもどうやらおかしなことをしてしまっているようです。仮に、
ABC <-function(i)
{si <-data.frame(sample(data[2:11],10,replace=T))
a <-data.frame(rowMeans(si[,1:5]))
b <-data.frame(rowMeans(si[,6:10]))
AB <-a/b}
とした際に、ABCとした際に、sourceを使って呼び出すというところまでやってみたのですが・・・。
更なるご助言をお願いします。
ABC <- function(i)
{si <-data.frame(sample(i[2:11],10,replace=T))
a <-data.frame(rowMeans(si[,1:5]))
b <-data.frame(rowMeans(si[,6:10]))
AB <-a/b
}
ABC(i=data)
でどうでしょうか？ -- akira 2009-01-07 (水) 12:48:09初心者です (2009-01-06 (火) 00:45:24)
legend(0,0.6,c("R","TPR","ATP"),col=1,pch=c(4,0,16),lty=1:3)で凡例を表示させました。グラフに表示させると（R,TPR,ATP)が半角で描かれているようにみえません。全角のようにみえます。たぶんMSPゴシックで表示されていると思うのですが、MSゴシックのような半角文字で表示したいのです。
初心の質問で申し訳ありませんが、よろしくお願いします。
rpart (2009-01-06 (火) 00:13:24)
rpartオブジェクトのprintの結果をオブジェクトとして保存して加工したいです。
下記のようなテキストのツリー図が標準出力に表示されますがオブジェクトとして保存する方法が見つかりません。1) root 150 102.1683000 5.843333 2) Petal.Length< 4.25 73 13.1391800 5.179452 4) Petal.Length< 3.4・・・・
どなたか方法をご存知でしょうか？
よろしくお願い致します。
お正月 (2009-01-05 (月) 12:47:11)
ベイズ推定でBRugsを使っていて、事前分布において、あるパラメータが切断正規分布（＜０）に従うのですが、切断正規分布の関数がないため、モデルの記述ができません。年末年始ずっと調べたり考えたりしたのですが、方法が見つかりせん。
どなたか方法をご存知でしょうか？
よろしくお願い致します。
しょうりゅう (2009-01-03 (土) 14:33:22)
データ分析においてRのすばらしさが実感しておりますが、デリバティブ（たとえばアメリカンオプションのプレミアムやリスクファクターの計算）が扱うようなパッケージはないでしょうか？よろしくお願いします。
nanbuwks (2009-01-03 (土) 11:48:21)
こんにちは。Rを試そうとして、いくつかのサイトを参考にインストール・操作をしているのですが、scatterplot.matrix( ~ID+月曜日+火曜日+水曜日, reg.line=lm, smooth=TRUE, span=0.5, diagonal = 'density', data=X )
エラー： 関数 "scatterplot.matrix" を見つけることができませんでした
scatterplot( USD~seq, reg.line=lm, smooth=TRUE, labels=FALSE, boxplots='xy', span=0.5, data=X )
エラー： 関数 "scatterplot" を見つけることができませんでした
となります。
MS-WindowsXP Professional上でバージョンはR2.8.1+Rcmdr1.4-6、およびパッケージとしてrglをインストールし、
"R User Configuration" http://androids.happy.nu/doc/r-tips#USER-CONFIG
を適用しています。
このエラーのためか、Rコマンダーのグラフメニューに「散布図行列」の項目が表示されません。
これは何かの不具合あるいはインストールミスなのでしょうか？それとも現バージョンではscatterplot関連は使えなくなっているのでしょうか？よろしくお願いします。
totoro (2008-12-28 (日) 20:14:54)
3D棒グラフ(縦)の作り方をご教示いただけないでしょうか。
Excelの3D棒グラフと同じようなものを作成しようと考えております。
persp関数も試したのですが、点を結ぶ形になってしまい上手くいきません。
よろしくお願いします。
> library(latticeExtra)
> cloud(VADeaths, panel.3d.cloud = panel.3dbars,
col.facet = "grey", xbase = 0.4, ybase = 0.4,
screen = list(z = 40, x = -30))
師走 (2008-12-26 (金) 11:49:50)
Web 上より以下のように画像を取得して、画像として表示などの処理をする方法ご教示願えないでしょうか？
> library(RCurve) > x<-getBinaryURL("http://www.okada.jp.org/RWiki/Rlogo.jpg") > class(x) [1] "raw" > summary(x) Length Class Mode 3173 raw raw
上の"raw"形式を画像にする方法が分りません。rgdal も試してみましたが、ディスク上に保存しているものにしか対応できないようです。
hidee (2008-12-25 (木) 19:08:22)
> a <- 5 > (function(x){x})(a) [1] 5 >
ですよね．
ここで戻り値が"a"となるような関数定義は
どのようにすれば良いのでしょうか．
> a <- 10
> (function(x){substitute(x)})(a)
a
> class((function(x){substitute(x)})(a))
[1] "name"
> (function(x){deparse(substitute(x))})(a)
[1] "a"
> class((function(x){deparse(substitute(x))})(a))
[1] "character"
>
であることがわかりました．-- hidee 2008-12-25 (木) 23:50:24syou6162 (2008-12-23 (火) 13:46:15)
gdbを使ってRオブジェクトがどのように表わされるか調べようとしています。
http://www.is.titech.ac.jp/~mase/mase/R-exts.jp/R-exts.jp.html#SEC50
を参考にしているのですが、「中断点を do_get に置き、R のプロンプトに get("DF") とタイプすると、 DF のメモリ中のアドレスを得ることが出来る」の部分がよく分かりません。Value returned is $1 = (SEXPREC *) 0x40583e1c
と書いてあるのですが、これをどうやって得ればよいのかが分かりませんでした。
途中までの操作は以下のようになっています。
/Users/yasuhisa% R -d gdb
GNU gdb 6.3.50-20050815 (Apple version gdb-768) (Tue Oct 2 04:07:49 UTC 2007)
Copyright 2004 Free Software Foundation, Inc.
GDB is free software, covered by the GNU General Public License, and you are
welcome to change it and/or distribute copies of it under certain conditions.
Type "show copying" to see the conditions.
There is absolutely no warranty for GDB. Type "show warranty" for details.
This GDB was configured as "i386-apple-darwin"...Reading symbols for shared libraries ..... done
(gdb) r
Starting program: /Library/Frameworks/R.framework/Versions/2.8/Resources/bin/exec/R
Reading symbols for shared libraries ++++........... done
Reading symbols for shared libraries . done
Reading symbols for shared libraries . done
R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Rはフリーソフトウェアであり、「完全に無保証」です。
一定の条件に従えば、自由にこれを再配布することができます。
配布条件の詳細に関しては、'license()'あるいは'licence()'と入力してください。
Rは多くの貢献者による共同プロジェクトです。
詳しくは'contributors()'と入力してください。
また、RやRのパッケージを出版物で引用する際の形式については
'citation()'と入力してください。
'demo()'と入力すればデモをみることができます。
'help()'とすればオンラインヘルプが出ます。
'help.start()'でHTMLブラウザによるヘルプがみられます。
'q()'と入力すればRを終了します。
[以前にセーブされたワークスペースを復帰します]
Reading symbols for shared libraries ............................... done
Reading symbols for shared libraries . done
> DF <- data.frame(a = 1:3, b = 4:6)
>
Program received signal SIGINT, Interrupt.
0x9037c5e2 in select$DARWIN_EXTSN ()
(gdb) b do_get
Breakpoint 1 at 0x3bbbdf: file envir.c, line 1600.
(gdb) c
Continuing.
Program received signal SIGINT, Interrupt.
0x9037c5e2 in select$DARWIN_EXTSN ()
(gdb) c
Continuing.
^D
Save workspace image? [y/n/c]: c
> get("DF")
Breakpoint 1, do_get (call=0xde0, op=0x1016c28, args=0x18529f0, rho=0x1852a98) at envir.c:1600
1600 SEXP rval, genv, t1 = R_NilValue;
(gdb)
Current directory is ~/
GNU gdb 6.3.50-20050815 (Apple version gdb-768) (Tue Oct 2 04:07:49 UTC 2007)
Copyright 2004 Free Software Foundation, Inc.
GDB is free software, covered by the GNU General Public License, and you are
welcome to change it and/or distribute copies of it under certain conditions.
Type "show copying" to see the conditions.
There is absolutely no warranty for GDB. Type "show warranty" for details.
This GDB was configured as "i386-apple-darwin".
/Users/yasuhisa/63250: No such file or directory.
Attaching to process 63250.
Reading symbols for shared libraries . done
Reading symbols for shared libraries ....................... done
0x9037c5e2 in select$DARWIN_EXTSN ()
(gdb) b do_get
Breakpoint 1 at 0x3bbbdf: file envir.c, line 1600.
(gdb) finish
Run till exit from #0 0x9037c5e2 in select$DARWIN_EXTSN ()
0x004b864a in R_SelectEx (n=1, readfds=0x5cea40, writefds=0x0, exceptfds=0x0,
timeout=0x0, intr=0x4b8bc0 <handleInterrupt>) at sys-std.c:152
152 val = select(n, readfds, writefds, exceptfds, timeout);
(gdb) c
Continuing.
Breakpoint 1, do_get (call=0xde0, op=0x1012628, args=0x184a268, rho=0x184a310) at envir.c:1600
1600 SEXP rval, genv, t1 = R_NilValue;
(gdb) finish
Run till exit from #0 do_get (call=0xde0, op=0x1012628, args=0x184a268, rho=0x184a310) at envir.c:1600
do_internal (call=0x11c9d18, op=0x1010238, args=0x184a268, env=0x184a310) at names.c:1143
1143 if (flag < 2) R_Visible = flag != 1;
Value returned is $1 = (struct SEXPREC *) 0x192a4a8
(gdb) p *$1
$2 = {
sxpinfo = {
type = 19,
obj = 1,
named = 2,
gp = 0,
mark = 0,
debug = 0,
trace = 0,
spare = 0,
gcgen = 0,
gccls = 1
},
attrib = 0x14f208c,
gengc_next_node = 0x192a488,
gengc_prev_node = 0x192a4c8,
u = {
primsxp = {
offset = 2
},
symsxp = {
pname = 0x2,
value = 0x0,
internal = 0x1a91c38
},
listsxp = {
carval = 0x2,
cdrval = 0x0,
tagval = 0x1a91c38
},
envsxp = {
frame = 0x2,
enclos = 0x0,
hashtab = 0x1a91c38
},
closxp = {
formals = 0x2,
body = 0x0,
env = 0x1a91c38
},
promsxp = {
value = 0x2,
expr = 0x0,
env = 0x1a91c38
}
}
}
(gdb) p R_PV($1)
$3 = void
(gdb) p R_PV($1->attrib)
$4 = void
(gdb)
アドレスは分かったのですが、R_PVで見てみようとするとvoidが返ってきてしまいます。また、*$1の内容も若干違うようです(vecsxpがない辺り)。sessionInfoは以下の通りです。
> sessionInfo() sessionInfo() R version 2.8.0 (2008-10-20) i386-apple-darwin9.5.0 locale: ja_JP.UTF-8/ja_JP.UTF-8/C/C/ja_JP.UTF-8/ja_JP.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods base
バージョン等で出力が違ったりするのでしょうか? -- syou6162 2008-12-24 (水) 02:54:40
(く) (2008-12-23 (火) 12:27:32)
クラスタ分析の樹形図をWindowsメタファイルとしてPowerPointに取り込もうとしています。（Windows版2.8.0）
画面上で日本語表示をさせて、それをメタファイルとしてPowerPointに取り込むとき：
１）クリップボード経由だと文字化けします。
２）ファイル経由で取り込むと一見表示できるように見えますが、メタファイルとして編集しようとすると文字化けしてしまいます。
おそらく、WMFに含まれていると思われるビットマップが表示される限りでは、日本語が正常に表示され、編集しようとすると図形の情報が使われて文字化けするのではないかと思います。
どのようにすれば日本語で編集可能なメタファイルが作成できるのでしょうか？
$ pstoedit -f svm tmp.eps tmp.svmとsvm形式（OpenOfficeのメタファイル形式）に変換し、次に、OpenOffice Drawでファイルを読み込んで、コンテキストメニューの「切り離し」を実行すれば、自由に編集できます。その後、PowerPointに持って行けばいかが？オブジェクトを切り離すのがコツです。 -- 谷村 2008-12-24 (水) 16:52:46
TT Arial : plain TT Arial : bold TT Arial : italic TT Arial : bold&italic次の4行に変更しました
TT MS Gothic : plain TT MS Gothic : bold TT MS Gothic : italic TT MS Gothic : bold&italicすると、メタファイルが文字化けせず編集できるようになりました。
flyer (2008-12-20 (土) 15:21:53)
因子分析での因子数について
いつもお世話になります。
「すべてがわかるアンケートデータの分析」（菅民郎 現代数学社）で因子分析を勉強しています。
例題で10人の生徒の数学、国語、英語の試験データから文系能力、理系能力と２つの因子で分析
する問題が載っております。
Rにてこの問題を行おうとしたのですが、
"以下にエラーfactanal(seiseki, factors = 2) : 2 factors is too many for 3 variables"とエラーメッセージが表示されます。
factorsを"1"で指定すると、factanalを実行できますが、求めたい数字にはなりません。
3つの変数（ここでは3教科）では2つの因子を指定することは出来ないのでしょうか。
よろしくお願いします。> seiseki <- matrix(c(2,3,3,1,4,5,2,2,3,3,2,1,5,4,5,4,4,4,8,5,4,6,3,4,7,6,7,4,5,6),10,3,byrow=TRUE) > colnames(seiseki) <- c("数学","英語","国語") > rownames(seiseki) <- c("001","002","003","004","005","006","007","008","009","010") > seiseki 数学 英語 国語 001 2 3 3 002 1 4 5 003 2 2 3 004 3 2 1 005 5 4 5 006 4 4 4 007 8 5 4 008 6 3 4 009 7 6 7 010 4 5 6 > (seiseki.fac <- factanal(seiseki, factors=2)) 以下にエラーfactanal(seiseki, factors = 2) : 2 factors is too many for 3 variables > (seiseki.fac <- factanal(seiseki, factors=1)) Call: factanal(x = seiseki, factors = 1) Uniquenesses: 数学 英語 国語 0.639 0.005 0.239 Loadings: Factor1 数学 0.601 英語 0.998 国語 0.872 Factor1 SS loadings 2.117 Proportion Var 0.706 The degrees of freedom for the model is 0 and the fit was 0.1292
- 無理です。1因子のときの結果で、自由度が0とでていますね（最終行）。つまり、2因子では自由度が足りないということです。もう10年も前の書籍のようなので、主因子法前提なのかもしれません。Rでは最尤法という比較的新しい計算方法を採用しているのですが、主因子法に比べて、計算を走らせるための縛りがきびしいという特徴があります。 -- 林 2008-12-21 (日) 00:29:35
- 林さん、どうもありがとうございます。最尤法を易しく説明している参考書などあるでしょうか。 -- flyer 2008-12-21 (日) 17:33:14
- 統計学を勉強したいのでなければ、分析手法として因子分析を解説しているものの方がよいと思います。次の書籍は文系向けに書かれていてわかりやすいと思います。http://www.amazon.co.jp/多変量データ解析法―心理・教育・社会系のための入門-足立-浩平/dp/4779500575/ref=sr_1_2?ie=UTF8&s=books&qid=1229856307&sr=8-2 -- 林 2008-12-21 (日) 19:50:49
- 林さん、参考にさせていただきます。ありがとうございました。 -- flyer 2008-12-23 (火) 00:24:19
tt (2008-12-20 (土) 12:15:11)
00:18:37.25 13/08/2008 137.199799 35.094155 423.5 39.3 0.60 7.4
00:18:37.50 13/08/2008 137.199799 35.094154 423.5 39.3 0.52 12.0
00:18:37.75 13/08/2008 137.199800 35.094156 423.5 39.3 0.61 15.5
00:18:38.00 13/08/2008 137.199800 35.094156 423.5 39.3 0.30 15.0
00:18:38.25 13/08/2008 137.199799 35.094156 423.5 39.3 0.39 10.0
00:18:38.50 13/08/2008 137.199800 35.094157 423.5 39.3 0.55 5.0
00:18:38.75 13/08/2008 137.199800 35.094157 423.5 39.3 0.65 356.5
00:18:39.00 13/08/2008 137.199800 35.094157 423.5 39.3 0.68 348.0
00:18:39.25 13/08/2008 137.199800 35.094157 423.5 39.3 0.72 340.5
というようなtxtファイルがあった場合、:や/で数値を区切るにはどうしたらよいでしょうか？
$ sed -e 's/[:/]/ /g' tmp.txt 00 18 37.25 13 08 2008 137.199799 35.094155 423.5 39.3 0.60 7.4 00 18 37.50 13 08 2008 137.199799 35.094154 423.5 39.3 0.52 12.0 00 18 37.75 13 08 2008 137.199800 35.094156 423.5 39.3 0.61 15.5 00 18 38.00 13 08 2008 137.199800 35.094156 423.5 39.3 0.30 15.0 00 18 38.25 13 08 2008 137.199799 35.094156 423.5 39.3 0.39 10.0 00 18 38.50 13 08 2008 137.199800 35.094157 423.5 39.3 0.55 5.0 00 18 38.75 13 08 2008 137.199800 35.094157 423.5 39.3 0.65 356.5 00 18 39.00 13 08 2008 137.199800 35.094157 423.5 39.3 0.68 348.0 00 18 39.25 13 08 2008 137.199800 35.094157 423.5 39.3 0.72 340.5
SASとRの両刀使い (2008-12-19 (金) 15:36:28)
本質問をRjpWikiで質問するのは大変恐縮なのですが、
SASのLSMeans（調整済み平均値；Least-Square means）を
Rで計算する関数はあるのでしょうか？> RSiteSearch("lsmeans")で調べてみても、上手い関数が無さそうで、Google検索をしてみても
良さそうな記事がありません。「multcomp」パッケージで計算出来るかも
という記事があったのですが、「multcomp」パッケージの中の関数には、
それらしきものがありませんでした。
よろしければお教えください。
> set.seed(7777)
> y <- c(rnorm(10),1+rnorm(10))
> DATA=data.frame(Y =y,
+ GROUP=c(rep(1,10),rep(2,10)),
+ COV1 =rbinom(20,1,0.5),
+ COV2 =jitter(y,amo=3))
> DATA$GROUP <- as.factor(DATA$GROUP)
> DATA$COV1 <- as.factor(DATA$COV1)
> head(DATA)
Y GROUP COV1 COV2
1 -1.8779550 1 1 -2.6682028
2 0.3264572 1 0 1.0997025
3 -0.3018570 1 1 -1.7436962
4 1.1816399 1 0 0.3136283
5 -0.1637332 1 1 -0.7099280
6 1.2964903 1 1 3.1079893
データ「DATA」に対して，とりあえず関数 lm() で回帰分析（SASではGLMプロシジャまたはmixedプロシジャ）を行います．
> options(contrasts = c("contr.treatment","contr.treatment"))
> result <- lm(Y ~ GROUP + COV1 + COV2, data=DATA)
> summary(result)
Call:
lm(formula = Y ~ GROUP + COV1 + COV2, data = DATA)
Residuals:
Min 1Q Median 3Q Max
-1.694067 -0.494184 0.002832 0.253098 2.225632
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4867 0.3604 1.351 0.19562
GROUP2 0.5943 0.5444 1.092 0.29110
COV11 -0.3924 0.4754 -0.825 0.42131
COV2 0.3774 0.1147 3.290 0.00461 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
この結果より，回帰直線の推定式は（GROUP(2)：GROUP==2のときに1でそれ以外は0，COV1(1)：COV1==1のときに1でそれ以外は0）
y=Intercept+GROUP(2)*0.5943+COV1(1)*(-0.3924)+0.3774*COV2
と推定されます．SASでは，この回帰直線の推定式からLSMeans（調整済み平均値；Least-Square means）を計算しています．
> mean(DATA$COV2) # COV2の平均（GROUPを無視して全ての平均） [1] 0.8944106 > 0.4867+0.0000-0.3924/2+0.3774*0.8944106 # GROUP==1 の y の LSMeans [1] 0.6280506 > 0.4867+0.5943-0.3924/2+0.3774*0.8944106 # GROUP==2 の y の LSMeans [1] 1.222351
上記の結果はSASの値と一致しています．中澤先生の方法で，正にSASと同じ方法が得られます！
> result2 <- dummy.coef(result)
> result2$GROUP + result2$"(Intercept)" + mean(result2$COV1) + result2$COV2*mean(DATA$COV2)
1 2
0.6280501 1.2223804
ちなみに，ザブさんから教わったJohn Fox氏のeffectsパッケージですと，SASと若干結果が異なるようです．
> library(effects)
> effect("GROUP", result)
GROUP effect
GROUP
1 2
0.6084311 1.2027614
ただ，交互作用モデルになると，ちょっと計算が手間になるので，適当なパッケージがあればなぁと思った次第です．
> options(contrasts = c("contr.treatment","contr.treatment"))
> result <- lm(Y ~ GROUP + COV1 + COV2 + GROUP*COV1 + GROUP*COV2 , data=DATA)
> summary(result)
Call:
lm(formula = Y ~ GROUP + COV1 + COV2 + GROUP * COV1 + GROUP *
COV2, data = DATA)
Residuals:
Min 1Q Median 3Q Max
-1.47336 -0.59849 -0.04832 0.24913 2.44732
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.59860 0.42080 1.423 0.1768
GROUP2 0.26255 0.92717 0.283 0.7812
COV11 -0.65934 0.66949 -0.985 0.3414
COV2 0.39928 0.16734 2.386 0.0317 *
GROUP2:COV11 0.58255 1.02201 0.570 0.5777
GROUP2:COV2 -0.02241 0.24532 -0.091 0.9285
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> mean(DATA$COV2)
[1] 0.8944106
> 0.5986006+0.000000-0.6593449/2+0.3992763*0.8944106+0.0000000/2+0.00000000
[1] 0.6260451
> 0.5986006+0.262553-0.6593449/2+0.3992763*0.8944106+0.5825505/2- 0.02241451*0.8944106
[1] 1.159826
>
> result2 <- dummy.coef(result)
> result2$ GROUP + result2$"(Intercept)" + mean(result2$COV1) + result2$COV2*mean (DATA$COV2)+
+ c(mean(result2$"GROUP:COV1"[c(1,3)]), mean(result2$"GROUP:COV1"[c(2,4)])) +
+ result2$"GROUP:COV2"*mean(DATA$COV2)
1 2
0.626045 1.159826
う〜ん，何時間も調べてみたのですが，適当なパッケージは無さそうですねぇ・・・．どうもありがとうございました！-- SASとRの両刀使い 2008-12-21 (日) 10:00:25
カホ (2008-12-18 (木) 21:09:27)
いつも参考にさせていただいております。
さっそくですが、質問させていただきます。
今、ベクトルが3つあり、それぞれ
S1=[,1][,2][,3][1,] 5 20 3
S2=[,1][,2][,3][1,] 4 20 8
S3=[,1][,2][,3][1,] 5 20 16
であるとき、
S1=[,1][,2][,3][1,] 5 20 3
S2=[,1][,2][1,] 4 8
S3=[,1][1,] 16
というものに変更したいとき、どのようにすればよいでしょうか？
よろしくお願いいたします。
take69 (2008-12-18 (木) 02:43:16)
Rの初心者です。
現在パネルデータを用いて、分析をしているんですが、パッケージ『plm』において、パネルデータに含まれる「NA」と表記されるデータがない場合の数値はどのようにパッケージ『plm』の中で処理されているんでしょうか？
パッケージの中身の内容なので、ご存知の方がいないかもしれませんが、よろしくお願いします。
yuki (2008-12-17 (水) 11:57:15)
WinBUGSではDICの計算を行うコマンドがあるようですが，BRugsでもそのようなコマンドはあるのでしょうか？
できればAICとDICを計算したいのですが・・・
Rを使い始めた者 (2008-12-16 (火) 16:55:57)
＜質問1＞
クライアント・サーバー方式ですが、Rをサーバープログラムとして、動かす事は可能でしょうか？
＜質問2＞
C++からRを呼び出す事は可能でしょうか？
＜質問3＞
Rの連続したコマンド命令をスクリプトとして作成、スクリプトを自動実行する事は可能でしょうか？
[name@zzzz ~]$ ps aux | grep Rserve name 5002 0.0 0.0 29132 3332 ? Ss May07 20:00 /usr/local/lib/R/bin/Rserve name 15142 0.0 0.0 4740 752 pts/31 S+ 18:34 0:00 grep Rserve
初めの一歩 (2008-12-16 (火) 03:37:25)
Ｒを使い始めたばかりの初心者です。
Ｒversion 2.8.0をインストールし、さらに必要なパッケージ複数をインストールしたいのですが、
le.create(f.tg) :fife 'C:\PROGRA~1\R\R-28~1.0/doc/html/packages.html' を作れません, 理由は'Permission denied'です。
と表示されます。
何がいけないのでしょうか。
ネットや本で検索をするのですが、解決策が見出されません。
使用ＰＣはパナソニック Let'note CF-Y7で、OSはWindows Vista Businessです。sessionInfo()
R version 2.8.0 (2008-10-20)
i386-pc-mingw32
locale:
LC_COLLATE=Japanese_Japan.932;LC_CTYPE=Japanese_Japan.932;LC_MONETARY=Japanese_Japan.932;LC_NUMERIC=C;LC_TIME=Japanese_Japan.932
attached base packages:
[1] stats graphics grDevices utils datasets methods base
どなたか、解決方法をご教示頂ければ幸いです。
必要な情報を失念している場合もあるかと思いますので、その場合はご指摘頂ければ幸いです。
よろしくお願いします。
学生 (2008-12-15 (月) 17:12:16)
Rでコサイナー法を行うためのパッケージや関数をご存知ないでしょうか．
インターネットで検索しても，ヒットしなくて困っています．
目的としては，血圧の日周変動の位相や振幅を求められるような解析をしたいと思っています．
よろしくお願いいたします．
初心者 (2008-12-15 (月) 15:53:02)
初心者で勉強不足で恐縮ですが、Rを使ってX-12ARIMAの季節調整を行う方法を、御存じの方がおられましたら、お教え下さりますようお願いいたします。
R初心者 (2008-12-13 (土) 15:57:49)
mulnos()関数によりノイズ寄与率を求めるのは、多変量自己回帰モデルからノイズ寄与率(パワー寄与率)を求めていると考えて良いのでしょうか？
mulnos()関数で求めた寄与率を横軸に周波数、縦軸に寄与率(％)のグラフで表すにはどのようにすれば良いのでしょうか？
初心者なので出来るだけ具体的に教えていただけたら幸いです。
SE (2008-12-13 (土) 03:39:22)
R を使って、集団学習のひとつである Random Forest 法に挑戦しました。
データはCSVファイル（cd_rf.csv）で以下のように１行目がヘッダで第１〜６列がデータ、第７列（status）が教師シグナル(-1 OR 1)になっています。
MAP Linda,S23 IgG,Red Star IgG,Bakers IgA,Red Star IgA,CA IgA,E faecalis,status
0.255048271,0.061732594,0.025085013,0.103566194,0.012575701,0.051569017,0.449583816,-1
0.269702681,0.100519062,0.106106155,0.332437714,0.022432332,0.078892899,0.528922137,-1
0.363834506,0.106774655,0.177443201,0.074255007,0.010397355,0.069219808,0.245495733,-1
：
このデータを使って以下のように解析を始めたのですが、classification でなく regression のモードになってしまいます。もし解決法（classification モードにする方法）をご存知でしたらご教示ください。
library(randomForest)
randomForest 4.5-25
Type rfNews() to see new features/changes/bug fixes.
Warning message:
package 'randomForest' was built under R version 2.6.2cdtrain<-as.matrix(read.csv("myData/cd_rf.csv",header=T))
is.matrix(cdtrain)
[1] TRUEcd.rf<-randomForest(status~., data=cdtrain,na.action="na.omit")
Warning message:
In randomForest.default(m, y, ...) :
The response has five or fewer unique values. Are you sure you want to do regression?cd.rf$type
[1] "regression"
サイト１に紹介されている例（Rに組み込まれているspamデータを解析）は問題なく実行できましたので、どうもデータ読み込みなど基本的な部分で失敗している気がします・・・。
参考にしたサイト>
１）http://www1.doshisha.ac.jp/~mjin/R/0603_32.pdf
２）http://cran.r-project.org/web/packages/randomForest/randomForest.pdf
OS: Windows XP version 5.1 (Build 2600.xpsp_sp2_gdr.080814-1233, Service pack 2)
R version 2.6.1 (2007-11-26)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
cdtrain <- read.csv("myData/cd_rf.csv",header=T)
cdtrain$"status" <- as.factor(cdtrain$"status")
cd.rf <- randomForest(status~., data=cdtrain,na.action="na.omit")
でどうかと。 -- akira 2008-12-15 (月) 13:04:28初心者 (2008-12-12 (金) 13:39:12)
最近Rを使い始めたばかりで初歩的な質問で申し訳ありません．
まずマニュアル本を使って，Rを触ってみているのですが，本の通りにデータを入力し，全く同じ方法で進めているはずなのに，エラーが出てしまいその理由がわかりません．
「Rによるやさしい統計学」のp.291からの重回帰分析でエラーが出るのですが，それ以前のT検定や分散分析，多重比較などは，本の通りに進めたらうまくいきましたので，この解析でいきなりなのです．
本の通りに進めると，
重回帰データ
娘 父 母1 155 159 167
2 158 168 156
3 163 173 153
4 151 153 153
5 157 169 158
6 158 169 152
7 155 163 149
8 155 163 158
9 155 163 151
10 154 161 155
11 156 164 151
12 155 160 161
13 153 162 153
14 157 167 152
15 157 172 153
16 156 166 148
17 157 163 150
18 157 161 145
19 160 161 149lm(娘~父+母)
以下にエラー eval(expr, envir, enclos) : オブジェクト "娘" は存在しません
ということになります．
Rの使用環境は以下の通りです．
何かおかしいところがあれば，教えていただけましたら幸いです．
sessionInfo()
R version 2.8.0 Patched (2008-12-07 r47101)
i386-pc-mingw32
locale:
LC_COLLATE=Japanese_Japan.932;LC_CTYPE=Japanese_Japan.932;LC_MONETARY=Japanese_Japan.932;LC_NUMERIC=C;LC_TIME=Japanese_Japan.932
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] sem_0.9-14
loaded via a namespace (and not attached):
[1] tools_2.8.0
lm(娘 ~ 父 + 母，data = x)
なり、
lm(娘$x ~ 父$x + 母$x)
とする必要があります。 -- 2008-12-12 (金) 14:21:54
子羊 (2008-12-12 (金) 11:06:37)
時系列データのプロットで、見やすくするために4半期ごとに縦の罫線を引こうと思って、
quarterly.lines <- function(x){
x.ts<- as.ts(x,start=c(1900,1),frequency=4)
plot(x.ts,ylab=" ",type="l")
par(new=T)
abline(v=1:23*4,lty=3)
}
という関数を書いてみたのですが、罫線が出力されません。
どこが間違っているのでしょうか。どなたか御存じの方が
いらっしゃいましたら、お教え頂ければ幸いです。
当方 (2008-12-09 (火) 01:23:21)
VineLinux4.2を使っています．VineにRを入れようとしているのですが，apt-getでインストールできません．
etc/apt/source.listの中身を見ても「VinePlus extra *1」に関するコメントがなく，どう編集していいか分かりません．どうか御教授いただければとおもいます．
peko (2008-12-07 (日) 11:00:23)
Steel-Dwass法の結果の読み方を教えてください。t p 1:2 3.866170 0.00006392122のような表記の場合と、
t p 1:2 4.3593873 7.692932e-05のように表記されるものがあります。
これらの結果の読み方を教えてください。
文系で初心者です (2008-12-06 (土) 11:15:20)
ビスタを使っています。
＞ install.packages(choose.files("", filters = Filters[c("zip", 中で＞警告がありました：
＞ 'lib = "C:/PROGRA~1/R/R-28~1.0/library"' は書き込み可能ではありませ＞ん
＞ 以下にエラー install.packages(choose.files("", filters = Filters[c＞("zip", :
＞ パッケージをインストール出来ませんでした
というのが出てきて使えません。
パソコンのセキュリティの問題かもと言われたのですが、よくわかりません。どうすればいいでしょうか？
初心者です。 (2008-12-04 (木) 12:04:11)
最近、初めてRを使い始めた初心者です。
Rを立ち上げた時に必ず出てくるエラーメッセージの意味が分からないで困っています。
メッセージは『関数 "etHook" を見つけることができませんでした』です。
現在、使用しているPCはMacPro(IntelMac)でOSはleopardを使用しています。
非常に簡単な質問かもしれませんが、ご回答いただければ助かります。
大学生 (2008-12-01 (月) 22:23:21)
R cranにてPTAkパッケージをダウンロードして，その中のPCAn(Principal component analysis on n modes)を利用しようと思ったのですが，以下にエラー if (length(X$met[[d]]) == dim(X$data)[d]^2) { : 引数の長さが0ですとエラーが出てしまいます．
関数の引数（データ形式）はlist型で行うようだったので，ex <- list( met = list( met1 = matrix(c(sample(1:25)),5,5), met2 = matrix(c(sample(1:25)),5,5) ), data = list( data1 = matrix(c(sample(1:25)),5,5), data2 = matrix(c(sample(1:25)),5,5) ) )としてみました．どのような形式ならうまくいくのか，できるだけ具体的に教えてください．よろしくお願いします．
line 167
PCnam <- outer(PCnam, 1:dim[q], FUN = "paste", sep = "")
たぶんこれの間違い
PCnam <- outer(PCnam, 1:dim[t], FUN = "paste", sep = "")
そこを直してみて，第１引数を以下のように定義して，PCAn を実行するととにかく何か結果らしきものは出る。正しいかどうかわからないけど。
ex <- array(matrix(24), dim=c(2,2,2,3)) PCAn(ex, dim=c(2,2,2,3))どうなんでしょうね？
tt (2008-12-01 (月) 03:34:05)
hist()関数で、y軸のfrequencyを％にするにはどうしたらよいのでしょうか？
a <- hist(x, ylab="%", yaxt="n") scale <- sum(a$counts)/100 labels <- pretty(0:max(a$counts)/scale) axis(2, at=labels*scale, labels=round(labels, 1))
なお (2008-12-01 (月) 01:27:36)
kernlabパッケージの中のksvmについて教えてください。
結果を別のシステムに実装するために、結果のパラメータなどを書き出すことは可能でしょうか？
predict()で指定する引数の部分を別のプログラム言語で使える形で出力したいのです。
勉強不足で申し訳ありませんが、よろしくお願いします。
ふじたま (2008-11-30 (日) 06:55:34)
本当に初歩的な質問で申し訳ないのですが、x<-c(1.1,1)
x
[1] 1.1 1.0
is.double(x)
[1] TRUE
is.integer(x)
[1] FALSE
のように、引数を整数、小数を混在させた時は
最初の引数が優先されるのでしょうか？
> x <- c(1.1, 1) > x [1] 1.1 1.0 > is.double(x) [1] TRUE > is.integer(x) [1] FALSE最初の要素を整数，二番目の要素を実数にして，同じことをやってみる。
> y <- c(1, 1.1) > y [1] 1.0 1.1 > is.double(y) [1] TRUE > is.integer(y) [1] FALSEどっちも同じ結果だ。
> x <- c("1.2", 1.1, 1)
> class(x)
[1] "character"
> x
[1] "1.2" "1.1" "1"> typeof(1.1) [1] "double" > typeof(1) [1] "double" > typeof(1L) [1] "integer"そもそも, 最初に整数として扱ってないので.
さと (2008-11-29 (土) 10:33:52)
表題の件について、質問いたします。
先日Rで作成したグラフを付けた論文を投稿し、アクセプトされたのですが、編集から画質が悪いので作り直すように指示がありました。（メッセージは以下の通りです）
Attached figures are not usable due to pixilated text and lines.
とりあえず Cairo packageをつかって、dpi=400で作り直しましたが、まだ画質が悪いといわれました。
作り直さねばいけませんが、何かよい方法をご存知の方がいらっしゃいましたら、ご教授いただければ助かります。
なにとぞよろしくお願いいたします。
（R掲示板へ移動しました）
まいくろ (2008-11-23 (日) 00:10:23)
マイクロアレイのデータ解析をしようとしているのですが、アレイのデータに別のファイルからの数値を結合したいと思っています。
この際、まず最初に２つのファイルに存在する「遺伝子名」が同じものを探し、その上でrowとrowを結合したいのですが、どのようにしたら良いでしょうか？
遺伝子が３０万個ほどあるので、自動的にできる方法があればよいのですが。
よろしくお願いします。
> (d1 <- data.frame(a=c(3,2,4), b=c(4,5,1), c=c(7,3,9))) a b c 1 3 4 7 2 2 5 3 3 4 1 9 > (d2 <- data.frame(e=c(4,3,1,6), f=c(3,5,8,7), g=c(4,2,5,9), b=c(3,4,5,7))) e f g b 1 4 3 4 3 2 3 5 2 4 3 1 8 5 5 4 6 7 9 7 > merge(d1, d2) b a c e f g 1 4 3 7 3 5 2 2 5 2 3 1 8 5 > merge(d1, d2, all=TRUE) b a c e f g 1 1 4 9 NA NA NA 2 3 NA NA 4 3 4 3 4 3 7 3 5 2 4 5 2 3 1 8 5 5 7 NA NA 6 7 9
例えば、データセット１
name aa bb cc a NA 6 4 b 6 4 6 c 4 4 1 d 4 3 4
とデータセット２
Name ee ff gg b NA 6 4 c 6 4 6 a 4 4 4
があった場合、求めたいものとして、
name aa bb cc ee ff gg a NA 6 4 4 4 4 b 6 4 6 NA 6 4 c 4 4 1 6 4 6 d 4 3 4 NA NA NA
としたいのですが、できますでしょうか？
こんな風になるんだけど？それでもまだあなたのやりたいこととは違うというなら， 例題が悪いんだね〜 > a name aa bb cc 1 a NA 6 4 2 b 6 4 6 3 c 4 4 1 4 d 4 3 4 > b name ee ff gg 1 b NA 6 4 2 c 6 4 6 3 a 4 4 4 > merge(a, b, all=TRUE) name aa bb cc ee ff gg 1 a NA 6 4 4 4 4 2 b 6 4 6 NA 6 4 3 c 4 4 1 6 4 6 4 d 4 3 4 NA NA NA
レオ (2008-11-22 (土) 14:59:49)
いつも拝見させていただいております。
さて、
「ある行列において列の成分が完全に重複するもの(順番は問わない)を1つにまとめる」
のに良い方法はないでしょうか？
例えば、以下のようなデータでは、
A= [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 2 2 5 3 2 2
[2,] 5 3 2 2 4 1
[3,] 3 4 3 5 3 5
とあった場合に
B= [,1] [,2] [,3]
[1,] 2 2 2
[2,] 5 3 1
[3,] 3 4 5
という行列が返ってくるようなことを望みます。
今のところ、毎iに、A[,i]の成分をsortし、
すべてのiで、A[,i]と同じものはallで判定して削除する
というようにしていますが、効率が悪い気がします。
もっと単純な方法があれば教えてください。
よろしくお願いいたします
> a <- matrix(c(2,5,3,2,3,4,5,2,3,3,2,5,2,4,3,2,1,5), 3)
> t(unique(t(apply(a, 2, sort))))
[,1] [,2] [,3]
[1,] 2 2 1
[2,] 3 3 2
[3,] 5 4 5エル (2008-11-22 (土) 14:07:23)
R歴3ヵ月の初心者です。
Cからの乗換で、Rにあったプログラムがなかなか書けていません。
行列を自動で複数個作りたいのですが、どのようにすればいいでしょうか。
イメージとしては以下のような感じです。
for(i in 1:3){ A[i] <- matrix(nrow=i,ncol=5) }
とすると
A1 = [,1] [,2] [,3] [,4] [,5] [1,]
A2 = [,1] [,2] [,3] [,4] [,5] [1,] [2,]
A3 = [,1] [,2] [,3] [,4] [,5] [1,] [2,] [3,]
このようなマトリックスを作るようなプログラムを作成したいです。
上記のプログラムだと以下のようなメッセージが出てきます。
以下にエラー A[i] <- matrix(nrow = i, ncol = 5) :
オブジェクト "A" は存在しません
詳しい方ご教授願います。
A <- NULL
for(i in 1:3){
A[[i]] <- matrix(NA, nrow=i, ncol=5)
}
~
やりたいこととはちょっと違うのかな？ --トマソン 2008-11-22 (土) 14:15:19山口 (2008-11-21 (金) 10:49:28)
いつも拝見させていただいております。
ある関数式がありまして、その関数をパラメータが未知のまま(-1.96,∞)まで積分したものを、最尤推定によってパラメータを推定するにはどうすればよいでしょうか。integrate()、optim()、function()など組み合わせてみたのですが、私の理解が足りないために全く動いてくれません。以下がそのプログラムです。
set.seed(1);th1<-rnorm(10000) myfunk<-function(x,sigma,u){ integrate((x/sqrt(2*pi*(sigma^2)))*exp(-((x-u)^2)/(2*(sigma^2)), -1.96, Inf)/ integrate((1/sqrt(2*pi*(sigma^2)))*exp(-((x-u)^2)/(2*(sigma^2)), -1.96, Inf) #sigmaとuがパラメータ } optim(c(1,1),myfunk,method="L-BFGS-B",th1)どなたか分かる方がいらっしゃいましたらご教示願います。
なお、環境はWindowsXp,R-2.8.0です。
どうぞよろしくお願いします。
integrate(f, lower, upper, subdivisions=100,
rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol,
stop.on.error = TRUE, keep.xy = FALSE, aux = NULL, ...)
ですので、【f】のところは関数オブジェクトを与えてあげなければなりません．なので
func1 <- function(x, sigma, u){
(x/sqrt(2*pi*(sigma^2)))*exp(-((x-u)^2)/(2*(sigma^2))
}
func2 <- function(x, sigma, u){
(1/sqrt(2*pi*(sigma^2)))*exp(-((x-u)^2)/2*(sigma^2))
}
こんな感じにしましょう．また【optim】で推定したいパラメータはベクトルとしてひとつの変数かつ第一引数として与えなければならないので
myfunk<-function(par, x){
sigma <- par[1]
u <- par[2]
integrate(f=func1, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value/
integrate(f=func2, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value
#sigmaとuがパラメータ
}
【par】がoptimで動かす変数になります．でも、まあ、いまのままだと【th1】がどういう意味を持ってるのかさっぱりなので、これでも動かないとは思いますが・・・．これを参考にがんばってみてください・・・myfunk<-function(par){
sigma <- par[1]
u <- par[2]
integrate(f=func1, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value/
integrate(f=func2, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value
#sigmaとuがパラメータ
}
なるのかな・・・．これでc(1,1)を入れてみると
> myfunk(c(1,1)) [1] 1.005001となりますね．これでoptimを動かそうと思ったら
optim(c(1,1),myfunk,method="L-BFGS-B")こうかな．でも、これだと推定したい値が、どこまでも動いてしまって「myfunk」でエラーが出るときがあると思うので、「optim」の
lower, upper引数を使って、その範囲を指定してあげてください．でもその前に、func1とfunc2が山口さんの思ってる式になっているか確認してみてくださいね・・・．
integrate(f=func1, lower=-1.96, upper=Inf, u=u)$value/ integrate(f=func2, lower=-1.96, upper=Inf, u=u)$valueはある定数値を取ることを事前に仮定していまして、ある定数値を与えたときのsigmaとuの値を最尤推定したいのです。私は単純に
myfunk<-function(par,spl){
sigma <- par[1]
u <- par[2]
spl<-integrate(f=func1, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value/
integrate(f=func2, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value
#sigmaとuがパラメータ
}
optim(c(1,1),myfunk,method="L-BFGS-B",lower=1, upper=5,spl=4)#仮にある定数値を4とする。
としてやったのですが、何とか動いてくれるもののlowerの値に左右されてうまくいっていないように思います。この部分は私が自分でできると思い質問しなかったのですが、逆に二度手間になってしまい申し訳ありません。これはlowerの値に左右される性質のものなのでしょうか。ちなみにfuncは正規分布の式にxをかけたものが分子でそのままのものが分母になると考えているので、これで合っていると思います。
どうぞよろしくお願い致します。 -- 山口 2008-11-28 (金) 17:05:52> myfunk<-function(par, spl){
+ sigma <- par[1]
+ u <- par[2]
+ (spl-integrate(f=func1, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value/
+ integrate(f=func2, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value) 2
+ #sigmaとuがパラメータ
+ }
> optim(c(1,1),myfunk,method="L-BFGS-B",lower=0.01, upper=25, spl=2, control=list(maxit=10000, trace=10))#仮にある定数値を4とする。
final value 0.000000
converged
$par
[1] 0.950756 1.999649
$value
[1] 4.374987e-18
$counts
function gradient
6 6
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"integrate(f=func1, lower=-Inf, upper=Inf, u=u)$value/ integrate(f=func2, lower=-Inf, upper=Inf, u=u)$valueとしたときに、それは正規分布の平均値を表していることを確認します。
func11 <- function(x){
(x/sqrt(2*pi*2^2))*exp(-((x-3)^2)/(2*(2^2)))#平均3標準誤差2とします。
}
func22 <- function(x){
(1/sqrt(2*pi*2^2))*exp(-((x-3)^2)/(2*(2^2)))
}
integrate(func11,-Inf,Inf)$value/
integrate(func22,-Inf,Inf)$value
[1] 3#ちゃんと3という値が返ってきます。
ですので、
integrate(func11,tau,Inf)$value/ integrate(func22,tau,Inf)$value#今まではtau=-1.96としていました。はtauに影響を受けて平均からずれた、平均値が返ってくることになります。今tau=-1.96とすると
integrate(func11,-1.96,Inf)$value/ integrate(func22,-1.96,Inf)$value [1] 3.03709#若干3からずれました。この値が定数値splとなります。私がやりたいのはtauが与えられているときに得られたsplを使ったuとsigmaの最尤推定ですので(ずれた後の平均値が分かっているときにずれる前のパラメータを推定したい)、
func1 <- function(x, sigma, u){
(x/sqrt(2*pi*(sigma^2)))*exp(-((x-u)^2)/(2*(sigma^2)))
}
func2 <- function(x, sigma, u){
(1/sqrt(2*pi*(sigma^2)))*exp(-((x-u)^2)/(2*(sigma^2)))
}
myfunk<-function(par, spl){
sigma <- par[1]
u <- par[2]
(spl-integrate(f=func1, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value/
integrate(f=func2, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value)^2
#sigmaとuがパラメータ
}
optim(c(1,1),myfunk,method="L-BFGS-B",lower=0.01, upper=25, spl=3.03709,
control=list(maxit=10000, trace=10))
final value 0.000000
converged
$par
[1] 1.100345 3.037075
$value
[1] 1.397589e-20
$counts
function gradient
6 6
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
と、平均値は大体上手く推定できているような気がしなくもないですが、sigmaの精度はよくないようです。正解例を示したかったのですが、sigmaが全くうまくいきませんでした。申し訳ありません。 -- 山口 2008-11-28 (金) 19:26:05func1 <- function(x, sigma, u){
(x/sqrt(2*pi*sigma))*exp(-(x-u)^2/2*(sigma)^2)
}#Before
から
func1 <- function(x, sigma, u){
(x/sqrt(2*pi*(sigma^2)))*exp(-((x-u)^2)/(2*(sigma^2)))
}#After
のように、最初のルートの中のsigmaを二乗したのと、それに伴い括弧を付けたし、後半のexpの中も括弧で計算の順序を整えました。func2については
func2 <- function(x, sigma, u){
(1/sqrt(2*pi*sigma))*exp(-(x-u)^2/2*(sigma)^2)
}#Before
から
func2 <- function(x, sigma, u){
(1/sqrt(2*pi*(sigma^2)))*exp(-((x-u)^2)/(2*(sigma^2)))
}#After
とfunc1と同様、sigmaを二乗してから括弧を付けたし、expの中も括弧で整えました。要は正規分布を表す式となっていなかったので、それを正しく直しました。 -- 山口 2008-11-28 (金) 20:49:55> optim(c(1,1),myfunk,method="L-BFGS-B",lower=0.01, upper=25, spl=3.03709,
control=list(maxit=10000, trace=10))
final value 0.000000
converged
$par
[1] 1.100345 3.037075
$value
[1] 1.397589e-20
$counts
function gradient
6 6
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"integrate(func11,-Inf,Inf)$value/ integrate(func22,-Inf,Inf)$value [1] 3で平均値が返ってくるので問題ないかと思っておりましたが…これが原因なのでしょうか？？ -- 山口 2008-11-28 (金) 21:38:14
func11 <- function(x){
(x/sqrt(2*pi*2^2))*exp(-((x-3)^2)/(2*(2^2)))#平均3標準誤差2とします。
}
func22 <- function(x){
(1/sqrt(2*pi*2^2))*exp(-((x-3)^2)/(2*(2^2)))
}
と定義しておりまして、このときに平均値を3、標準誤差を2と定めています。すみません、計算式の横に書いてあったのが不適切でした。ですから、3という正解が返ってくる以上、負の値をfunc1がとっても問題ないかと思いました。また、このように定義していますのでoptim()で最後に求まる値はsigma=2とu=3であってほしいのです。 -- 山口 2008-11-28 (金) 21:44:00func1 <- function(x, u){
(x/sqrt(2*pi*(2^2)))*exp(-((x-u)^2)/(2*(2^2)))
}
func2 <- function(x, u){
(1/sqrt(2*pi*(2^2)))*exp(-((x-u)^2)/(2*(2^2)))
}
myfunk<-function(par, spl){
u <- par[1]
(spl-integrate(f=func1, lower=-1.96, upper=Inf, u=u)$value/
integrate(f=func2, lower=-1.96, upper=Inf, u=u)$value)^2
#uのみがパラメータ
}
optim(1,myfunk,method="L-BFGS-B",lower=0.01, upper=25, spl=3.03709,
control=list(maxit=10000, trace=10))
final value 0.000000
converged
$par
[1] 3.000000
$value
[1] 6.01349e-16
$counts
function gradient
7 7
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
となり、きれいなu=3という推定値が導けるようです。もう少しデザインがどうにかならないか考えてみます。ありがとうございました。 -- 山口 2008-11-28 (金) 23:08:43山口さんが書いたこと： splを使ったuとsigmaの最尤推定ですので(ずれた後の平均値が分かっているときにずれる前のパラメータを推定したい)これを読んで感じたのは「tau=-1.96」でtauをずらした場合平均値も分散もずれるんじゃないかなと言うことです．それなのに「spl=3.03709」とずれたあとの平均だけを指定して、元の平均と分散を求めようとしている？このあたりにヒントがありそうな気がしますが、いかがでしょう？？？（まったくの思い付きです）
func33 <- function(x){
((x-3)^2)*(1/sqrt(2*pi*2^2))*exp(-((x-3)^2)/(2*(2^2)))
}#分散を確率密度を使って表現、平均3標準誤差2
integrate(func33,-Inf,Inf)$value/
integrate(func22,-Inf,Inf)$value
[1] 4
#確かに2^2の値が返ってきます(標準誤差を2と定めているので)
sqrt(integrate(func33,-1.96,Inf)$value/
integrate(func22,-1.96,Inf)$value)
[1] 1.953467
#これがtauによって2からずれた標準誤差の値
ということを確認して、以下が作成したプログラムです。
func1 <- function(x, sigma, u){
(x/sqrt(2*pi*(sigma^2)))*exp(-((x-u)^2)/(2*(sigma^2)))
}
func2 <- function(x, sigma, u){
(1/sqrt(2*pi*(sigma^2)))*exp(-((x-u)^2)/(2*(sigma^2)))
}
func3 <- function(x, sigma, u){
((x-u)^2)*(1/sqrt(2*pi*(sigma^2)))*exp(-((x-u)^2)/(2*(sigma^2)))
}
myfunk<-function(par, spl, spl2){
sigma <- par[1]
u <- par[2]
(spl-integrate(f=func1, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value/
integrate(f=func2, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value)^2+
(spl2-sqrt(integrate(f=func3, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value/
integrate(f=func2, lower=-1.96, upper=Inf, sigma=sigma, u=u)$value))^2
#sigmaとuがパラメータ,splにずれた平均値を,spl2にずれた標準誤差を設定します
}
optim(c(1,1),myfunk,method="L-BFGS-B",lower=0.01, upper=25, spl=3.03709,
spl2=1.953467, control=list(maxit=10000, trace=10))
final value 0.000000
converged
$par
[1] 2.000000 3.000000#きれいに標準誤差と平均値が推定できました。
$value
[1] 1.235012e-14
$counts
function gradient
7 7
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
おかげさまで当初の目的を達成することができました。皆様本当にありがとうございました。そしてななしさん、適切なアドバイスに心から感謝申し上げます。 -- 山口 2008-11-30 (日) 16:29:08心理学徒 (2008-11-20 (木) 15:52:44)
こんにちは。
下のような感じのデータフレームで多変量解析を試みたのですが，多変量外れ値を示すオブザベーションがあったので，その行だけ消去して新しいデータフレームを作りたいのですが，どうやればよいでしょうか？
どうぞよろしくお願いします。sex age id x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 1 0 19 8055 3 3 2 3 2 2 2 1 1 1 2 1 19 8045 5 5 4 3 4 5 2 4 5 4 3 0 19 8061 4 4 1 1 1 1 2 1 5 4 4 0 18 8084 5 4 2 4 5 4 4 2 2 5 5 1 19 8051 3 2 3 4 4 2 2 1 1 2 6 0 18 8047 4 4 1 4 2 1 1 1 1 2
北京原人 (2008-11-19 (水) 10:39:21)
Rの初心者です。
下記のデータを「Re1」で取り込んで、V1 V2 V3 V41 16.01704 7.389199 8.523638 2.508838
これを下記のコマンドで４行10列の行列に拡張（Re2）して
Re2 <- matrix( Re1 , nrow=10 , ncol=4, byrow=T)
V1 V2 V3 V41 16.01704 7.389199 8.523638 2.508838
2 16.01704 7.389199 8.523638 2.508838
3 16.01704 7.389199 8.523638 2.508838
4 16.01704 7.389199 8.523638 2.508838
5 16.01704 7.389199 8.523638 2.508838
6 16.01704 7.389199 8.523638 2.508838
7 16.01704 7.389199 8.523638 2.508838
8 16.01704 7.389199 8.523638 2.508838
9 16.01704 7.389199 8.523638 2.508838
10 16.01704 7.389199 8.523638 2.508838
行列計算しようとするとエラーが出てしまいます。
Re2+Re2
以下にエラーRe2 + Re2 : 二項演算子の引数が数値ではありません
何が問題でエラーが出るのかが分かりません。
おそらく基本的なミスをしているのだと思いますが原因が分かる方ご教示ください。
> Re2 <- matrix( data.matrix(Re1) , nrow=10 , ncol=4, byrow=T)
> Re2+Re2
[,1] [,2] [,3] [,4]
[1,] 32.03408 14.77840 17.04728 5.017676
[2,] 32.03408 14.77840 17.04728 5.017676
途中省略
[9,] 32.03408 14.77840 17.04728 5.017676
[10,] 32.03408 14.77840 17.04728 5.017676さと (2008-11-19 (水) 06:21:51)
jpeg, pngなどの画像ファイルを、Rに読み込む方法についてお尋ねします。
jpeg(), png()などのコマンドによって、Rで作った図を様々なファイル形式で保存することが可能ですが、逆に既に保存されている画像ファイルをRの中に読み込んで表示することは可能でしょうか。
具体的には、graphics windowをいくつかに分割して、その一部にはRで作った図を書き込み、他の部分にはすでに作成して保存しているpngなどのファイルの画像を張り込むことを考えています。
ご教授いただければ幸いです。
windows (2008-11-18 (火) 10:16:24)
いつもお世話になっております。
windowsパスをRで読み込む際に"\"を"/"に置き換えなくてはいけないことが
わかりましたが、手作業ではなく自動で置き換えることを検討しています。
下記のようなコマンドを作成しましたが、エラーが出て止まってしまいます。
解決策をご教授いただけないでしょうか。
よろしくお願いいたします。
a <- gsub("\", "/", "C:\Documents and Settings\user\デスクトップ\hoge")
Saito (2008-11-17 (月) 20:45:38)
いつもお世話になっております。
地形の変化などによって調査区内において予測地点間が均質でない場合、どうやれば不均一なグリッドデータ(緯度経度の座標)が作成できるでしょうか。例えば以下のような場合です。library(geoR) a<-pred_grid(c(10,20), c(10,20),by=1)#予測グリッドの作成 plot(a,xlab="long",ylab="lat") b<-c(10:20,20:10) c<-sample(b,11) c2<-c(c,numeric(11)) d<-data.frame(b,c2) points(d,col=2,type="l") polygon(d,col=2)#赤色で塗り潰されて"いない"ところが予測したい場所。この場所だけを含む緯度経度データを作成したい。似たような質問で"gstatパッケージでのクリギングについて"という記事の中で、expand.grid()を用いたグリッドデータの作成方法が紹介されています。ですが、この方法だと長方形のグリッドデータが用意されてしまいます(上記のpred_gridと一緒)。予測したい地点の緯度経度のみを算出するにはどうすればよいでしょうか。
ちなみに環境はWindowsXP、R-2.8.0です。
よろしくお願いします。
> library(geoR)
> a<-pred_grid(c(10,20), c(10,20),by=1)#予測グリッドの作成
> plot(a,xlab="long",ylab="lat")
> b<-c(10:20,20:10)
> set.seed(1); c<-sample(b,11) # 再現性のため乱数種を指定
> c2<-c(c,numeric(11))
> d<-data.frame(b,c2)
> points(d,col=2,type="l")
> polygon(d,col=2)
>
> J <- NULL
> d1 <- d[,1]; d2 <- d[,2]
> for (i in seq(nrow(d))) {
for (j in seq(nrow(a))) {
if (d1[i] < a[j,1] & d2[i] == a[j,2])
J <- c(J,j)
}}
> a <- a[unique(J),]
> dim(a)
[1] 51 2cc<-rep(c(10,11,13,14,15,16,17,18,19,20),c(5,3,6,7,4,5,1,2,10,8))#x軸の作成
cc2<-c(16,17,18,19,20,18,19,20,15,16,17,18,19,20,
14,15,16,17,18,19,20,17,18,19,20,16,17,18,19,20,
20,19,20,11,12,13,14,15,16,17,18,19,20,13,14,15,16,17,18,19,20)#y軸の作成
a2<-data.frame(cc,cc2)
a2
と一致するはずなのですが、私の理解が間違っているでしょうか。
また、「不規則領域に限るのは予測をしてからにした方が面倒でないと思います」と言って下さったので、私がやりたかったことを言いますと、クリギングによって予測された地点の値の合計です。具体例としましては、
library(geoR)
a<-pred_grid(c(10,20), c(10,20),by=1)
plot(a,xlab="long",ylab="lat")
b<-c(10:20,20:10)
set.seed(1); c<-sample(b,11)
c2<-c(c,numeric(11))
d<-data.frame(b,c2)
points(d,col=2,type="l")
polygon(d,col=2)
cc<-rep(c(10,11,13,14,15,16,17,18,19,20),c(4,2,4,4,4,5,1,2,5,4))
cc2<-c(16,17,18,20,18,20,15,17,19,20,
14,15,16,20,17,18,19,20,16,17,18,19,20,
20,19,20,11,17,18,19,20,13,15,17,20)
cc3<-rep(c(10,11,13,14,15,16,17,18,19,20),c(5,3,6,7,4,5,1,2,10,8))
cc4<-c(16,17,18,19,20,18,19,20,15,16,17,18,19,20,
14,15,16,17,18,19,20,17,18,19,20,16,17,18,19,20,
20,19,20,11,12,13,14,15,16,17,18,19,20,13,14,15,16,17,18,19,20)
par(mfrow=c(1,3))
plot(a,main="予測グリッド")
plot(cc,cc2,main="実際の観測点")
plot(cc3,cc4,main="予測値を合計したい地点")
set.seed(2);cc5<-abs(rnorm(35))
cc6<-data.frame(long=cc,lat=cc2,value=cc5)
library(gstat)
library(lattice)
vario<-variogram(object=value~long+lat,
locations=~long+lat,cutoff=5,width=1,data=cc6)
plot(vario$dist,vario$gamma)
model2<-fit.variogram(vario,model=vgm(0.03,"Sph",2.5,0.42))
cc7<-gstat(id="val",formula=value~long+lat,locations=~long+lat
,data=cc6,model=model2)
a<-data.frame(long=a$Var1,lat=a$Var2)
cc8<-predict.gstat(cc7,a)#予測グリッドにおける内挿
levelplot(val.pred~long+lat,cc8)
sum(cc8$val.pred)#合計したくないところまで値が合計されてしまう:plot(cc3,cc4)で 表される地点の合計値のみが欲しい
このような感じになります。つまり合計したい値は任意の緯度より上の値だけなのです。
library(gstat) data(meuse.grid) plot(meuse.grid$x,meuse.grid$y)で表されるように、クリギングをする前には最初から予測グリッドを任意の地点のみに指定しておいたほうがあとで合計するときなどに便利だとと思い質問させていただいたのです。私の場合は(イメージとして)山の裾野の境界の緯度データがありまして、平地だけでクリギングを行いたいと考えています。境界データはものすごく長いので、手作業で予測グリッドを作ると、どこかで呪われそうなのです。そういうわけで、なんとか境界の位置データがわかっている場合に、そこを末端とする予測グリッドを作成したかったのです。しかし、予測したあとに不規則領域を区切ったほうが面倒でないということは、何か上手いやり方があるのでしょうか？もし何かご存知でしたら、どうかご教授願えないでしょうか？ -- Saito 2008-11-18 (火) 23:41:47
set.seed(1);b<-c(10:20,20:10);c<-sample(b,11)
S2<-NULL
R2<-NULL
c3<-20-c
for(i in 1:length(c)){
ifelse(-c[i]==-20,S<-0,
S<-seq(-20,-c[i]-1))
ifelse(c3[i]==0,R<-0,
R<-rep(c(9+i),c3[i]))
R2<-c(R2,R)
S2<-c(S2,S)
}
a<-data.frame(R2,S2)
a<-subset(a,a["R2"]!=0)
a$S2<--a$S2
c2<-c(c,numeric(11))
d<-data.frame(b,c2)
plot(a)
points(d,col=2,type="l")
polygon(d,col=2)
どうもありがとうございました。 -- Saito 2008-11-19 (水) 17:37:04
samplenanbuwks20090312.png
oodraw01.png
sample.png
