【親ページ >>> Q&A (初級者コース)

初心者のための R および RjpWiki に関する質問コーナー

新規投稿はできません。

過去の記事のアーカイブ



MacOS10.5でR2.9.0、日本語の表示ができません。

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")
}

回帰をグループごとに実行し、同データフレームに出力したいのですが。。。

頭バクハツです、、、もぅ、、、 (2009-04-23 (木) 06:34:03)


非常に初歩的な質問かもしれないのですが、
私にとってはとても難しく、悩んでも答えが出ませんでした。
どなたか、ぜひ、ご教示ください。

以下のようなデータフレームXXがあります。
変数A,Bを用いて回帰回帰分析を行い、最終的に変数Xの予測値、PREDICT.Xを予想することを目的とします。
これだけであればよかったのですが、次のことをする必要に迫られています。(データは簡略化し、かつ適当な数字です。)

1958年とそれ以降のデータを分け、別個に回帰する。
あるいは、変数label内で、同じラベルを張られているものは、同じグループとして回帰する。
(1958年以前のラベルはすべて1ですので、上記の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つを使用しました”となってしまいます。

どのようにコードすれば、データフレームXXに予測値PREDICT.Xを計算することが出来るのでしょうか。
どうか、よろしくお願いします。

pls.ldaの使い方について

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からBまで)←←こう出力したいのです。

}


example(1,2) →→→→ グラフのタイトルを ”1から2まで”と出力したい。
example(3,4) →→→→ グラフのタイトルを ”3から4まで”と出力したい。

chisq.test()の出力について

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"となる場合があります。これらはどのように違うのでしょうか? また、後者はどのようなときに出力されるのでしょうか?

EPSファイルへのTeXコマンド埋め込み(2.9.0 for Winでは失敗)

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を見ても、この問題に対応する仕様変更は見つけられませんでした。対応方法をご存知の方はいらっしゃいませんでしょうか?

R2.9.0のインストール

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",  : 
  パッケージをインストール出来ませんでした 。

boxplotの計算方法について

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 [シート1$]" )
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

pairwise.t.test()について

ビクス (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という引数が
ないのでしょうか?
宜しくお願いします。

geoRを用いて自分のデータを使う場合

ど素人 (2009-04-06 (月) 10:35:31)

バリオグラムを描きたいのですが、自分のデータをエクセルからcsvファイルに変換してバリうグラムを描きたいのですがどうしたらよいのでしょうか?
libraryに入っているデータならなんとかなるのですが、新しいデータを用いて活用したいのです。

ACF(自己相関係数)の計算結果のファイルへの出力

サザエ (2009-04-04 (土) 02:21:50)

頑張ってRのこと、勉強しております。
質問です。ACF(自己相関係数)の計算結果をExcelで読めるようにファイルに移したいんです。そのやり方を教えてください。
よろしくお願いします。

function と引数について

迷い猫 (2009-04-03 (金) 21:47:08)

全くの初心者です。
やれることはやったのですが、解決策が分からないのでどうか、どうか、よろしくお願いします。

データフレームとして、Excelから読み込んだJapan, U.S, U.Kが準備されています。
必要に応じたデータフレームを使用したいので、SASで行うマクロのような関数を作りたいと思い、以下のようなプログラムを作りました。

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関数に、データフレームを直接引数として渡す方法は存在するのでしょうか。

どうか、よろしくお願いします。

system関数の不具合?

レッドクリフ (2009-04-03 (金) 08:49:56)

system関数の不具合?で困っています.
通常はsystem関数で実行した結果を返すのに、Rでメモリを大量に消費しているときに、コマンドの結果がエラーで返ってきます.

> sessionInfo()
R version 2.8.0 (2008-10-20)
x86_64-unknown-linux-gnu
64ビットのLinux機、メモリは8G積んでいます.

具体的な症状は以下の通りです.

  • 成功する例
    > 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
    • このように5Gのメモリを消費して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:
    • セグメンテーションフォルトで落ちます・・・
  • このような現象の回避方法はございませんでしょうか.よろしくお願いいたします.

Windows Vista で R2WinBUGS をバッチ実行する方法

初心者 (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)

サンプリング間隔が不等な時系列データを、内挿等によって等間隔のデータを
作成するにはどうすればよいのでしょうか。お教え下さい。

特異値分解(singular value decomposition)の実行とエラーについて

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


質問は,以下の3つです。それぞれの質問について,下で補足します。
1.「数千行×数千列」の行列をRで特異値分解する方法はないか。
2.「数万行×数万列」の行列をRで特異値分解する方法はないか。
3.svd()関数がどのようなスクリプトによって定義されているか見る方法はないか。

「メモリ」と「svd()」をキーワードに過去ログやWEBを探しました。
質問1は,OSやメモリの扱いによる解決を期待した質問です。
質問2は,速度を犠牲にしてもスクリプトがストップせずに特異値分解を実行できるアルゴリズムや方策を期待した質問です。
質問3は,自分で特異値分解のスクリプトをかけばもっと省メモリでやれるかもしれないという期待からお尋ねしています。

何卒よろしくお願いいたします。

set.seedの使い方

カーネルおじさん (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)の位置はどちらでもいいのではと思いますが,Rはプログラム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では,#1〜10000=-124.08, #10001〜20000=-121.98, #20001〜30000=-115.20, #30001〜40000=-113.31
プログラム2では,#1〜10000=-121.43, #10001〜20000=-120.52, #20001〜30000=-115.89, #30001〜40000=-116.48



プログラム1
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

tkplotのグラフサイズ・背景色の変更方法は?

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)

図3 福田総理の所信表明演説文における名詞の共起ネットワークマップ
の様なネットワーク図を出力しています。

質問は
・「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) などとしていています。
 もっとましな方法がありそうですが

宜しくご教授ください。

csvデータのインポート

初学者 (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で用意されているのか知りたいと思い、質問させて頂きました。

以上、ご教授の程よろしくお願いします。

chplotで観測点の色の指定

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"))

としてみましたが、凸包輪郭線の色のみ反映し、観測点の色が変更できません。

samplenanbuwks20090312.png

どのようにしたらいいでしょうか?

levelpotでcolorkeyを対数表示する方法

u-kun (2009-03-12 (木) 15:42:40)

latticeパッケージのlevelplotを使ってグラフを作成しようとしています。
格子状の測定点での測定値の分布を可視化したいのです。
しかし、z軸に対応するデータが10^(-10)から10^(-3)まで幅があるので、
colorkeyの色の変化(目盛り)を対数にしたいのですがうまくいきません。
例えばplotであればlog="x"などとすればx軸が対数になるように、
colorkeyを対数にするにはどうすれば良いのでしょうか?

boxplotで観測点を消す方法

nanbuwks (2009-03-12 (木) 11:49:19)

x <- c(0,1,1,2,2,3,5)
boxplot(x)

とすると、箱ひげ図とその範囲外の観測点が表示されますが、この観測点を表示させない方法ってあるのでしょうか?

行の逆向き累乗積総和の計算

カーネルおじさん (2009-03-11 (水) 10:32:21)
お世話になります.Rの基本的なことがまだまだわかっておりません
しかし,Rを楽しんでます.
今般,i行のそれまでの行の累乗積の総和をi行に返すプログラムに苦労しています.
つまり,行列Xから行列Zの変換で,新しい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が以下の通りだとして,

> 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行列が求まります.

z
    [,1] [,2] [,3]
[1,]    2    3    1
[2,]    6    6    5
[3,]   36    8   10
[4,]   80   36   45

Xを単純な累乗積Yに変換するプログラムは,以下の通りに書けましたが,
i行から逆向きに累乗積を計算するのが難しく,
どうもうまくいきません.プログラムの外側にもうひとつFor文(j in 1:4)を挟み,
初期値i=1(i<2の箇所)をi=jなりに変え,行列Yで得られた行をZに落していけば
ワークしそうなのですが,実力不足を実感しております.

> 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,]}
}

scatterplotのhelp

nanbuwks (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です。よろしくお願いします。

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ですが、あまり必要なさそうなので結局使いませんでした。言いだしっぺということで、ご参考までに。より簡潔で美しいコードをご教示いただければ、もちろん乗り換えます(笑)。みなさま、よろしくお願いします。

sample.png

部分代入不可の場合

るる (2009-03-05 (木) 19:42:16)

R初心者です。自分の解析結果を入れたファイルが容量が多く、119行あって上の1:47が見れません。

それで、
>ファイル名[1:47]
と打ち込んだら、部分代入不可というエラーが出ました。
部分代入不可ならどうやって、47行目までの入るの内容を読めばよいのでしょうか。

3次元標本値に対応するカーネル密度値を求めたいのですが・・・

カーネルおじさん (2009-03-04 (水) 08:18:39)

R初心者です.お世話になります.約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のグラフの軸名と目盛りを省略し、
さらに余白を調節して3つのグラフを連結しようと考えました。
つまり以下のように入力します。

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軸のプロットする位置を自由に指定することは
できないでしょうか?また、ここでは3つのグラフを書きましたが、
グラフの数に依存せずに調節したいと考えています。
どなたかご教授宜しくお願いします。

相関係数からクラスター分析

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 D

1 A 0 1 2 2
2 B 1 0 2 2
3 C 2 2 0 2
4 D 2 2 2 0

data1.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"

plclustのHeight値の範囲を変更

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.314

Heightの目盛の下限を20に変更するには

plclust(hc, hmin=20)

太陽質量 M_\odot を表示

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

数値変換

初学者 (2009-02-16 (月) 14:21:57)

はじめまして。お世話になります。

使用環境は、Windows XPで、RのVersionは、下記の通りです。
R version 2.8.1 (2008-12-22)
i386-pc-mingw32

アンケートに逆転項目があるため、回答の数値を逆するため、下記のようなプログラムを記述しましたが、問題がありエラーになりました。

[プログラム例]
回答者1〜10の方のQ1のデータを、1は7に、2は6に、・・・7は1に変換するため、下記のプログラムを実行しました。(下記のデータでは、6や7の回答はありませんが、実データには含まれている可能性があるので、実データ6や7に対しても、データ変換を試みています)

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}

以上、アドバイスよろしくお願いいたします。

lines()が使えない

capetapeta (2009-02-14 (土) 10:06:54)

Macintosh、Windows、Linux(Win,Linuxは仮想マシン)いずれにしてもRを何回ダウンロードしてもlinesという命令文を書いても正しくQuartz2にUPされないんです。以前にも色々試してnew plot? abline?
という投稿をしましたが、結局線がプロットできないという結論に至りました。これを解決するためには何か有効な手だてはあるのでしょうか?或は線がプロットできるようになるようなパッケージがあるのでしょうか?

while等で、1日廻したいが、1時間でハングする。

独学中 (2009-02-14 (土) 09:15:58)

いろいろ検索してみましたが、いい案がなく、どなたか助言をお願いします。
やりたいこと
 測定データのファイルが、あるディレクトリに入るのをwhileで待って、
 入ってくれば、統計処理させたい。
 測定は、約3分ごと。
 データファイルの大きさは、70KBで、36本ほど。
困っていること
 一時間ほど、走らせると、Rがハングしている。
 
使用環境 winXP, R ver2.8.1 memory.limitは、1500MB
以下の工夫をしているが、効果なし。
 変数は、使い終われば、rm(),gc(),gc().
 sinkで、出力先を変更
 Rconsoleが、メッセージで一杯かと思い、10分ほどで、ESCしてみるが、
 メッセージは、表示されない。(何事もなく、>を表示。)

コックス回帰でのAICの求め方について

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が一応出てくることは確認しています。
必要な情報を失念している場合もあるかと思いますので、その場合はご指摘頂ければ幸いです。よろしくおねがいします。

高次多項式によるデータの近似式の作成

hoso (2009-02-05 (木) 16:13:39)

データとして(X,Y)のペアが100個ぐらいあり、高次多項式Y=f(X) を求めてカーブフィッティングしたいと思っています。(データの性質にもよると思いますが、割と滑らかで高次多項式でのフィッティングが妥当と思わせるデータです。)これだけだったらエクセルでも可能です。ところがエクセルでやりますと、高次項の係数が2e-9などと桁落ちになっております。それらを用いて計算すると、Xの値を大きくして高次項が利き始めると大きなエラーとなってしまいます。係数が2.3462e-9などと精度が高い場合、エラーが小さくなるものと思います。エクセルの分析ツール+回帰分析を行ってみると直線近似でした。高次の式(6次ぐらい)でもエクセルの分析ツールで処理できるのでしょうか。

このような問題をRで実行するとどのようになるでしょうか。具体的なやり方(データの読み込みは分かりますが)を教えて頂けると助かります。あるいは関連した実例集を指示して頂くのも有難いです。

よろしくお願いします。

par$maiをPDFに反映させたい

<ふ> (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)をインストールした時、下記のエラーが出てしまいインストールができません。既存のQ&Aを見ると、Rの再インストールやミラーサイトの変更などが出てきたため試してみましたが、未だに解決されません。
また、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)

開かれた URL 

downloaded 2.3 Mb

パッケージ 'Rcmdr' は無事に開封され、MD5 サムもチェックされました 
以下にエラー normalizePath(path) : 
 path[1]: 指定されたファイルが見つかりません。


library(Rcmdr)

以下にエラー library(Rcmdr) :  'Rcmdr' という名前のパッケージはありません。


パッケージのインストールでつまづくと、作業が進まず非常に困っています。
どなかたアドバイスいただけないでしょうか。お願いします。

plot.newとabline?

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版R 2.8.1でPDFを作成するときの日本語フォントが

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 freedom
      Residual deviance: 175.12 on 5 degrees of freedom
      AIC: 48.402

      Number of Fisher Scoring iterations: 2

      となりますが,その後にどのようにして95%信頼限界,R2,F値を求めたらよいかわかりません.
      すみませんが,よろしくお願いします.

Rのコンソールの内容をファイルに保存する

初心者 (2009-01-25 (日) 20:43:44)

R で長いプログラムを書いております.
Rのコンソールの内容(命令+実行結果+エラーメッセージ+警告メッセージ)
をファイルに保存するにはどうすればよいのでしょうか?
例えば,sinkだと実行結果しか保存できません.教えていただければと思います.
ちなみに R 2.8.1, OS は Ubuntulinux を使っております.

plotでylim

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軸を上にずらす以外のやりかたで方法がないでしょうか。

メモリ不足

(2009-01-23 (金) 15:17:52)

SVMを使いたいのですが、「エラー:サイズ 276.9 Mb のベクトルを割り当てることができません」となってしまい、実行できません。
PCのメモリは2GBで、57600×101(ヘッダー含む)のデータを判別したいのです。

メモリを変える以外になにか解決方法ありませんでしょうか?

クラスタ分析後の樹形図の作成

初学者 (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="ウォード法")

以上、アドバイスの程よろしくお願いいたします。

最小全域木

初心です (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)をインストールした時、下記のの用になりましたが、初心者で内容や、回復方法が分かりません。既存のQ&Aを見ましたが、見落としを含めて、見当たりません。回復方法を教えていただけると助かります。よろしくお願いします。

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)

開かれた URL 

downloaded 148 Kb

パッケージ 'tree' は無事に開封され、MD5 サムもチェックされました 
以下にエラー normalizePath(path) : 
 path[1]: 指定されたファイルが見つかりません。

library(tree)

以下にエラー library(tree) :  'tree' という名前のパッケージはありません

2変量混合正規分布に従う乱数

Rin (2009-01-16 (金) 15:50:03)

はじめまして。
プログラムについての質問です。
φを平均-1,分散1の正規分布に従う確率変数,ψを平均2,分散1の正規分布に従う確率変数として

   f(x,y)=0.6φ(x,y)+0.4ψ(x,y)~

に従う2変量混合正規乱数を発生させたいのですが、どのようにすればよいでしょうか?
現在、自分が書いたものではないのですが、他の方が書いたプログラムとして、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))
       }

1行目の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です。
よろしくお願い致します。

ちなみに例題の解は1にはならないような...

arulesで、listデータからtransactionの処理ができなくなった。

NewComer (2009-01-14 (水) 14:18:19)

R2.7.0、XPで、アッソシエーション分析をしていました。1月になって
arulesを再インストール(全パケージをインストールしたため)したところ
listからtransactionの処理が出来なくなりました、本を見ながらあちこち試しましたがNG、友人にも友人のPCで試してもらいましたがNG、初心者で
分かりません。以下、データを添付しますので、よろしくお願いします。

#テスト用に作成したデータ~
install.packages("arules",dependencies=TRUE)
library(arules)
{data<-list(c("A","B","C","D"), c("E","A","G","B"))}
data.tran<-as(data,"transactions")
class(data.tran)
#結果
> {data<-list(c("A","B","C","D"),
+ c("E","A","G","B"))}
> 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%

VineLinuxでソースからのビルド

当方 (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

R2.8.1におけるRエディタの背景色の自動変更について

Saito (2009-01-11 (日) 19:11:22)

いつもお世話になっております。
R2.8.1になって、編集→GUIプリファレンスと進みコンソール画面の背景色を変更するとRエディタの方でも背景色が自動で変更されるようなのですが、output text/User inputで入出力の色を変更してもRエディタの方は入力する文字の色が黒のままで変更されず、背景色との組み合わせによっては非常に見づらいです。何度もアンインストールとインストールと再起動をして確かめたのですがやはり直りません。どこかでRエディタの色を変えることができるのか調べたのですがわかりませんでした。このようなことになるのは私だけでしょうか?どなたか解決法をご存知でしたら教えていただけると幸いです。
なお使用環境はWindowsXpです。

Contourグラフ

EDI (2009-01-11 (日) 09:02:29)

Contnurグラフを作成しようと考えています。

xy平面上にplot()を用いて散布図を作成しました。

平面上に散布された点(因子)の密度を数値化し、その数値をz値に使用と考えています。
まだ初心者なので、調べ方も不十分なのかとは思いますが、多くのContnurグラフを作成するfunctionたとえばcontour()などはこちらで(x,y,z)を定義する必要があるので使えていません。

そこで、(x,y)を定義するだけで、私の考えているようなContnurグラフを作成させられる関数があれば、教えていただければ幸いです。
あるいは、xy平面上に散布されたの点の密度分布をすべての格子点で求められるような関数があれば、そこからz値を定義できるのでそういった関数でもかまいません。何かアドバイスいただけると幸いです。よろしくお願いします。

回帰線の95%信頼区間

(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)

エクセルへの書き出し

みち (2009-01-09 (金) 19:59:59)

Rでの実行結果をエクセルに書き出す方法を教えてください。

DPpackageのインストール

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を使用しています。
どなたか解決方法を教えて下さい。

minhash法

okinawa (2009-01-08 (木) 17:40:48)

いつもお世話になっております。okinawaです。jaccard係数を用いたminhash法をRで利用する方法をどなたかご存じないでしょうか?数式の解説やsasのサンプルは見つけたのですが、読みこなせなくて・・・。
環境は:windowsXPSP3,R2.7.2

繰り返し計算

ぽち (2009-01-06 (火) 08:59:17)

すみません。一通り探したのですがわかりませんでした。
以下のような流れでspreadsheetについて計算をしたいのです。これ自体は動くのですが、問題はこれを1000回繰り返し、1000個の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を使って呼び出すというところまでやってみたのですが・・・。

更なるご助言をお願いします。

凡例の文字のフォントについて

初心者です (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ゴシックのような半角文字で表示したいのです。
初心の質問で申し訳ありませんが、よろしくお願いします。

print.rpartの結果をオブジェクトとして保存できませんか?

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 

・・・・
どなたか方法をご存知でしょうか?
よろしくお願い致します。

BRugsでの切断正規分布

お正月 (2009-01-05 (月) 12:47:11)

ベイズ推定でBRugsを使っていて、事前分布において、あるパラメータが切断正規分布(<0)に従うのですが、切断正規分布の関数がないため、モデルの記述ができません。年末年始ずっと調べたり考えたりしたのですが、方法が見つかりせん。
どなたか方法をご存知でしょうか?
よろしくお願い致します。

オプション価格計算パッケージについて

しょうりゅう (2009-01-03 (土) 14:33:22)

 データ分析においてRのすばらしさが実感しておりますが、デリバティブ(たとえばアメリカンオプションのプレミアムやリスクファクターの計算)が扱うようなパッケージはないでしょうか?よろしくお願いします。

関数scatterplotを見つけることができませんでした

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関連は使えなくなっているのでしょうか?よろしくお願いします。

3D棒グラフについて

totoro (2008-12-28 (日) 20:14:54)

3D棒グラフ(縦)の作り方をご教示いただけないでしょうか。
Excelの3D棒グラフと同じようなものを作成しようと考えております。
persp関数も試したのですが、点を結ぶ形になってしまい上手くいきません。
よろしくお願いします。

Web 上の raw データを画像として扱う方法

師走 (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"となるような関数定義は
どのようにすれば良いのでしょうか.

gdbを使ったデバッグ

syou6162 (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) 

バージョン等で出力が違ったりするのでしょうか? -- syou6162 2008-12-24 (水) 02:54:40

Windowsメタファイル

(く) (2008-12-23 (火) 12:27:32)

クラスタ分析の樹形図をWindowsメタファイルとしてPowerPointに取り込もうとしています。(Windows版2.8.0)
画面上で日本語表示をさせて、それをメタファイルとしてPowerPointに取り込むとき:
1)クリップボード経由だと文字化けします。
2)ファイル経由で取り込むと一見表示できるように見えますが、メタファイルとして編集しようとすると文字化けしてしまいます。

おそらく、WMFに含まれていると思われるビットマップが表示される限りでは、日本語が正常に表示され、編集しようとすると図形の情報が使われて文字化けするのではないかと思います。
どのようにすれば日本語で編集可能なメタファイルが作成できるのでしょうか?

因子分析:変数の数と因子の数について

flyer (2008-12-20 (土) 15:21:53)

因子分析での因子数について

いつもお世話になります。
「すべてがわかるアンケートデータの分析」(菅民郎 現代数学社)で因子分析を勉強しています。
例題で10人の生徒の数学、国語、英語の試験データから文系能力、理系能力と2つの因子で分析
する問題が載っております。
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

txtファイルを区切る方法

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ファイルがあった場合、:や/で数値を区切るにはどうしたらよいでしょうか?

SASのLSMeans(調整済み平均値;Least-Square means)をRで

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

というものに変更したいとき、どのようにすればよいでしょうか?
よろしくお願いいたします。

plmにおけるNAの処理の仕方

take69 (2008-12-18 (木) 02:43:16)

Rの初心者です。

現在パネルデータを用いて、分析をしているんですが、パッケージ『plm』において、パネルデータに含まれる「NA」と表記されるデータがない場合の数値はどのようにパッケージ『plm』の中で処理されているんでしょうか?

パッケージの中身の内容なので、ご存知の方がいないかもしれませんが、よろしくお願いします。

BRugsでのAIC,DICの計算について

yuki (2008-12-17 (水) 11:57:15)

WinBUGSではDICの計算を行うコマンドがあるようですが,BRugsでもそのようなコマンドはあるのでしょうか?
できればAICとDICを計算したいのですが・・・

Rをサーバープログラムとして動かす事は可能でしょうか?

Rを使い始めた者 (2008-12-16 (火) 16:55:57)

<質問1>
クライアント・サーバー方式ですが、Rをサーバープログラムとして、動かす事は可能でしょうか?
<質問2>
C++からRを呼び出す事は可能でしょうか?
<質問3>
Rの連続したコマンド命令をスクリプトとして作成、スクリプトを自動実行する事は可能でしょうか?

パッケージのインストールが出来ません。

初めの一歩 (2008-12-16 (火) 03:37:25)

Rを使い始めたばかりの初心者です。
Rversion 2.8.0をインストールし、さらに必要なパッケージ複数をインストールしたいのですが、

le.create(f.tg) :

  fife 'C:\PROGRA~1\R\R-28~1.0/doc/html/packages.html' を作れません, 理由は'Permission denied'です。


と表示されます。
何がいけないのでしょうか。
ネットや本で検索をするのですが、解決策が見出されません。

使用PCはパナソニック 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でコサイナー法を行うためのパッケージや関数をご存知ないでしょうか.
インターネットで検索しても,ヒットしなくて困っています.
目的としては,血圧の日周変動の位相や振幅を求められるような解析をしたいと思っています.
よろしくお願いいたします.

RでX-12ARIMAを実行するには

初心者 (2008-12-15 (月) 15:53:02)

初心者で勉強不足で恐縮ですが、Rを使ってX-12ARIMAの季節調整を行う方法を、御存じの方がおられましたら、お教え下さりますようお願いいたします。

多変量自己回帰モデルの寄与率

R初心者 (2008-12-13 (土) 15:57:49)

mulnos()関数によりノイズ寄与率を求めるのは、多変量自己回帰モデルからノイズ寄与率(パワー寄与率)を求めていると考えて良いのでしょうか?
mulnos()関数で求めた寄与率を横軸に周波数、縦軸に寄与率(%)のグラフで表すにはどのようにすれば良いのでしょうか?
初心者なので出来るだけ具体的に教えていただけたら幸いです。

Random Forest

SE (2008-12-13 (土) 03:39:22)

R を使って、集団学習のひとつである Random Forest 法に挑戦しました。
データはCSVファイル(cd_rf.csv)で以下のように1行目がヘッダで第1〜6列がデータ、第7列(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.2

cdtrain<-as.matrix(read.csv("myData/cd_rf.csv",header=T))

is.matrix(cdtrain)
[1] TRUE

cd.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"

サイト1に紹介されている例(Rに組み込まれているspamデータを解析)は問題なく実行できましたので、どうもデータ読み込みなど基本的な部分で失敗している気がします・・・。

参考にしたサイト>
1)http://www1.doshisha.ac.jp/~mjin/R/0603_32.pdf
2)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

書籍の通りに進めてもエラーメッセージ

初心者 (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 149

lm(娘~父+母)

以下にエラー 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

4半期ごとに罫線を引きたいのですが

子羊 (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)
}

という関数を書いてみたのですが、罫線が出力されません。

どこが間違っているのでしょうか。どなたか御存じの方が
いらっしゃいましたら、お教え頂ければ幸いです。

VineLinuxにR

当方 (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

のように表記されるものがあります。
これらの結果の読み方を教えてください。

RODBCのインストールのやり方

文系で初心者です (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を使用しています。
非常に簡単な質問かもしれませんが、ご回答いただければ助かります。

n相主成分分析について

大学生 (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)
	)
)

としてみました.どのような形式ならうまくいくのか,できるだけ具体的に教えてください.よろしくお願いします.

hist( )関数について

tt (2008-12-01 (月) 03:34:05)

hist()関数で、y軸のfrequencyを%にするにはどうしたらよいのでしょうか?

ksvmの出力について

なお (2008-12-01 (月) 01:27:36)

kernlabパッケージの中のksvmについて教えてください。
結果を別のシステムに実装するために、結果のパラメータなどを書き出すことは可能でしょうか?
predict()で指定する引数の部分を別のプログラム言語で使える形で出力したいのです。
勉強不足で申し訳ありませんが、よろしくお願いします。

c関数について

ふじたま (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
のように、引数を整数、小数を混在させた時は
最初の引数が優先されるのでしょうか?

Rの画像の画質について

さと (2008-11-29 (土) 10:33:52)

表題の件について、質問いたします。
先日Rで作成したグラフを付けた論文を投稿し、アクセプトされたのですが、編集から画質が悪いので作り直すように指示がありました。(メッセージは以下の通りです)

Attached figures are not usable due to pixilated text and lines.

とりあえず Cairo packageをつかって、dpi=400で作り直しましたが、まだ画質が悪いといわれました。

作り直さねばいけませんが、何かよい方法をご存知の方がいらっしゃいましたら、ご教授いただければ助かります。

なにとぞよろしくお願いいたします。

軽くしました

R掲示板へ移動しました)

2つのデータを結合する方法は?

まいくろ (2008-11-23 (日) 00:10:23)

マイクロアレイのデータ解析をしようとしているのですが、アレイのデータに別のファイルからの数値を結合したいと思っています。
この際、まず最初に2つのファイルに存在する「遺伝子名」が同じものを探し、その上でrowとrowを結合したいのですが、どのようにしたら良いでしょうか?
遺伝子が30万個ほどあるので、自動的にできる方法があればよいのですが。
よろしくお願いします。

例えば、データセット1

	name	aa	bb	cc
	a	NA	6	4
	b	6	4	6
	c	4	4	1
	d	4	3	4

とデータセット2

	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

としたいのですが、できますでしょうか?

列の成分が同じものを削除した行列を作成

レオ (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で判定して削除する
というようにしていますが、効率が悪い気がします。

もっと単純な方法があれば教えてください。
よろしくお願いいたします

for文を使って規則性のある名前の行列を生成したい

エル (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" は存在しません

詳しい方ご教授願います。

積分を最尤法に組み合わせて最適値を求めたい

山口 (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です。
どうぞよろしくお願いします。

データフレームの特定の行だけ消去したい

心理学徒 (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        V4

1 16.01704  7.389199  8.523638  2.508838
これを下記のコマンドで4行10列の行列に拡張(Re2)して
Re2 <- matrix( Re1 , nrow=10 , ncol=4, byrow=T)

        V1       V2       V3       V4

1 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 : 二項演算子の引数が数値ではありません

何が問題でエラーが出るのかが分かりません。
おそらく基本的なミスをしているのだと思いますが原因が分かる方ご教示ください。

画像ファイルの読み込みについて

さと (2008-11-19 (水) 06:21:51)

jpeg, pngなどの画像ファイルを、Rに読み込む方法についてお尋ねします。

jpeg(), png()などのコマンドによって、Rで作った図を様々なファイル形式で保存することが可能ですが、逆に既に保存されている画像ファイルをRの中に読み込んで表示することは可能でしょうか。

具体的には、graphics windowをいくつかに分割して、その一部にはRで作った図を書き込み、他の部分にはすでに作成して保存しているpngなどのファイルの画像を張り込むことを考えています。

ご教授いただければ幸いです。

windowsパスのRでの自動置換

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です。
よろしくお願いします。


添付ファイル: filesample.png 1590件 [詳細] fileoodraw01.png 1731件 [詳細] filesamplenanbuwks20090312.png 1582件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2015-03-01 (日) 01:15:59