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

新規投稿はできません



丸め方向の指定

abc? (2004-11-30 (火) 19:25:46)

windows版のRを使用しているのですが、丸めモードの変換を行い、上向けに丸めたときの結果と下向きに求めた時の結果の比較を行いたいのですが、
どのようにすればよろしいのでしょうか?教えてください。
よろしくお願いします。

  • かなり意味不明。 -- 2004-11-30 (火) 20:22:47
  • 「丸めモードの変換」ってのがよく分かりませんが、丸めた結果は一緒じゃないんですか?例を挙げてもらわないと分からないです。 -- 2004-11-30 (火) 20:25:21
  • 浮動小数点数を_controlfp(_RC_DOWN,_MCW_RC); -- 2004-11-30 (火) 20:47:12
  • もしかして、round, floor, ceil 関数の比較という意味ですか。もしそうなら JIS,ISO式四捨五入 なども参考になるかも知れません。 -- 2004-12-01 (水) 16:57:26
  • 何を意図しているかはわかりませんね。誠意のないところを見ると,即刻削除がそうとうかと。 -- 2004-12-01 (水) 23:22:50
  • 確かに何を意図している質問かわかりませんが,あまりにも初心者に対する回答の態度が高圧的ではありませんか?初心者にとっては,正確に質問すると言うことは,とても難しいことです.また,質問することもとても勇気がいることです.初心者用掲示板なのですから,1度目の質問では,適度な易しさも必要かと -- 2004-12-02 (木) 01:24:59
  • R(LAPACKも)はIEEE演算の丸めを基に開発されていますので,CPUの丸めモードを変更(CPUのモード変更ですよね)すると言うのはR的(LAPACK的にも)意味の無い結果をもたらします.よって,丸めモードの変更を行って何か調べたいのなら,言語的にもCで調べるのが適切ではないかと考えます.初心者に文句をたれる方はここでの書き込みは卒業して頂いて上級編へどうぞ.初心者(色んな種類の初心者がいるでしょうけど)を罠に^H^H導くのがこの掲示板になれば良いと思います。:-) -- なかま 2004-12-02 (木) 10:16:52

エクセルへ値コピー

ABA? (2004-11-29 (月) 23:47:56)

Rで数値データを一列に出力させたとき、値の横に一行目には[1,]、
2行目には[2,]・・・と表示されるのですが、出力されたデータの値のみを
エクセルにコピー&ペーストするには、どうすればよいのでしょうか。
色々と試行錯誤したのですが、分からないので教えてください。
宜しくお願いします。

  • Windows環境だとすれば,Rで数値データがベクトルXに入っているとして,writeClipboard(as.character(X))としてExcelに行ってペーストするというのはどうですか? -- 中澤? 2004-11-30 (火) 00:18:41
  • Macな人は,apply(X, 1, function(i) cat(paste(i, "?t"), "?n")) などとやってみるとよいかも。画面には思っていた通りには出力されないかもしれないが,数値部分をコピーしてExcel上にペーストするとあら不思議。ちゃんとセルに値が入ってくれます。このままだと NULL という行が最後に表示されるのでその部分はコピーしないか,全体を invisible( ) で囲ってやるか, junk に付値するか。気に入ったら関数化するとか。
    d.list <- function(d)
    {
        invisible(apply(d, 1, function(i) cat(paste(i, c(rep("?t", length(i)-1), "?n")))))
    }
    質問を読み返してみたら,数値データはベクトルに入っていることを想定していて,Excelには縦に一列にペーストしたいのだろうか。だとしたら,y <- rnorm(10) などになっているなら,cat(paste(y, "?n")) -- 青木繁伸 2004-11-30 (火) 12:54:30
  • ありがとうございます。中澤先生と青木先生の両方の手法で解決できました。ご丁寧な説明、大変感謝です。ちなみに私はWINDOWSでした。今度からはきちんと書きますので宜しくお願いします。 -- ABA? 2004-12-01 (水) 19:32:52
  • ご丁寧な対応,ありがとうございます。わかっている人には,中澤先生のやりかたで十分なことは自明です。 -- 青木繁伸 2004-12-01 (水) 23:25:34

標準化回帰係数

B? (2004-11-26 (金) 13:58:42)

重回帰をした際に標準化回帰係数を知りたいのですが、summaryでは表示されません。手計算で出すことができると思うのですが、出力する関数があれば教えてください

  • ここにあります -- 中澤? 2004-11-26 (金) 17:59:29
  • 発想の転換で,scale を使って標準化したデータセットで同じ分析を行う。計算時間はもったいないけど,複雑なモデルだと間違いは少ないかも。
    d がデータフレームだとして,
    ds <- data.frame(scale(d)) # 標準化する
    lm(y ~ x1+x2+x5, d)
    lm(y ~ x1+x2+x5, ds) # 第二引数だけが違う
    # 厳密に言えば,後者のモデルの最後に -1 を付けておくのがいいのかもしれない
    前者は回帰係数,後者は標準化回帰係数が得られる(^_^) -- 青木繁伸 2004-11-26 (金) 19:20:51
  • ありがとうございます。両方ともやってみます -- B? 2004-11-27 (土) 13:01:18

rimage インストール時のエラー

ささかま? (2004-11-26 (金) 00:10:45)

R2.0.1 mac版を使っています。OS10.3.6 Xcode 1.1を入れてあります。 rimageをインストールするため、Finkを使ってfftwとlibjpegを入れました。その後R上からrimageをインストールしようとしましたが、

checking jpeglib.h usability... no
checking jpeglib.h presence... no
checking for jpeglib.h... no
configure: error: Sorry, can't find jpeglib header

と出てしまい、インストールできませんでした。fftwをcheckするときにも同様のメッセージが出ましたが、以前教えて頂いたように(「macでソースファイルからパッケージをインストールする方法」という題です)/usr/local/include/にシンボリックリンクを作ることでエラーは回避できました。しかし、jpeglib.hについては、シンボリックリンクを作ってもエラーが出てしまいます。別の場所にリンクを作らなければならないのでしょうか?どなたかご存知の方ご教授願います。 この掲示板の使い方がよくわからず上手く書き込めませんでしたが、どなたか直して頂いたようで、ありがとうございました。

  • BTW(by the way)ですが,Finkのページによると,Xcode の最新は 1.5 で,しかしそれについている gcc にはパッチが入ったとのこと。 -- 青木繁伸 2004-11-27 (土) 15:56:36
  • レスありがとうございます。Xcodeを1.5にしてgccのパッチも当てました。しかしFinkの設定か何かが悪いらしく上手く行きませんでした。そこでlibjpegとfftwをFinkを使わずにインストールしたところ、無事にrimageをインストールできました。ちなみにfftw3.0.1ではだめで、fftw2.1.5を入れたら上手くいきました(3は2以前と互換性がないらしいです)。どうもお騒がせしました。 -- ささかま? 2004-11-27 (土) 23:26:19
  • こんにちは。ただいま、ささかまさんと同じく、Mac版R2.1にrimageをインストールしようとしているのですが、うまくインストールできません.
    finkを使ってfftwとjpeglibはインストールしました。
    Rでパッケージをソースからコンパイルするとき、/usr/includeを見にいっているようなのですが、これをfinkでインストールした/sw/includeに変える方法をご存知の方はいらっしゃらないでしょうか? -- setoyama? 2005-12-28 (水) 16:38:30

統計的有意の記号

y.a.? (2004-11-20 (土) 22:37:54)

例えば,回帰分析の結果に対して summary.lm()を signif.stars=T
オプション付きで実行するとt値の各有意水準に応じて星をつけてくれます.
この星の水準を変更するにはどうすればよいでしょうか?
具体的には,10%, 5%, 1% 水準でそれぞれ,*, **, ***
をつけたいのです(経済学ではこれがデフォルトです).

summary.lm をみると,係数の行列は ans$coefficients ですが,
これに星をつけている箇所はよく分かりませんでした.

  • 私も不思議(経済学の標準も不思議ですが、経済データでは 1% 有意でもびっくりものということでしょうか)だったので r-help で尋ねてみたら、早速反応がありました。以下に紹介します。 あんまり簡単ではなさそうですね。-- 間瀬茂 2004-11-21 (日) 07:18:02
    For example, calling summary(lmObject) dispatches on method summary.lm()
    hwich creates an object of class "summary.lm". The latter is printed by
    method print.summary.lm() which calls printCoefmat().
    
    The stars are hard-coded there, and I don't think anybody is going to change
    that. I suggest to turn of the printing of siginificant codes by specifying
    print(summary(.....), signif.stars = FALSE) or by setting the corresponding
    option().
    
    Uwe Ligges
    
    It would be possible to re-define 'printCoefmat' privately so as to change
    the lines
    
          cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
          symbols = c("***", "**", "*", ".", " ")) 
    
    towards the end of its code into whatever you prefer, e.g.
    
          cutpoints = c(0, 0.01, 0.05, 0.1, 1),
          symbols = c("***", "**", "*", " "))
    
    or
    
          cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
          symbols = c("****", "***", "**", "*", " "))
    
    (both of which are compatible with your description of what is needed).
    
    The most straightforward way of redefining it is to copy the code for
    'printCoefmat' into a file, e.g.
    
      sink("printCoefmat.R")
      printCoefmat
      sink()
    
    and then edit that file.
    NOTE that the code written to the file does not include
    the name of the function, i.e. it starts
    
      function (x, digits = max(3, getOption("digits") - 2),....
    
    so the first modification has to be 
    
      printCoefmat<-function(x, digits = .... )
    
    Then, when you want your private version, simply do
    
     source("printCoefmat.R")
    
    and it will overlay the original version. (Experts will have to advise
    whether this clashes with any "namespace" issues.
    On my reading of the code, it doesn't seem to; but I'm no expert!)
    
    If your friend wants to use this new definition all the time,
    then one way to arrange this is to put the revised function
    definition (as in the edited file) into his .Rprofile,
    or put the command
      source("printCoefmat")
    into that file.
    
    Best wishes,
    Ted.
  • 早速のお返事ありがとうございます. -- y.a.? 2004-11-21 (日) 11:34:31 経済学(や他の社会科学?)では,5%水準で有意であればだいたいヨシとすると思います.以前,工学では基準が0.1%と聞いて驚愕したことがあります.
    ところで,上記の上書きする方法を試したのですが,どうもうまくいきません… 以下の例では,printCoefmat が上書きされているにも関わらず,結果はオリジナルのprintCoefmatの通りに出力されます.
    > tail(printCoefmat, 15)
    74                 Signif <- symnum(pv, corr = FALSE, na = FALSE,         
    75                   cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),         
    76                   symbols = c("***", "**", "*", ".", " ")) 
    (以下略)            
     
    > data(swiss)
    > summary(lm(swiss$Fertility ~ swiss$Education + swiss$Examination))
    
    (前略)
    Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
    (Intercept)        85.2533     3.0855  27.630   <2e-16 ***
    swiss$Education    -0.5395     0.1924  -2.803   0.0075 ** 
    swiss$Examination  -0.5572     0.2319  -2.402   0.0206 *  
    ---
    Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
    (略) 
    
    > # 上記の方法でprintCoefmatを上書きする
    > source(file="printCoefmat.R")
    > tail(printCoefmat, 15)  # 上書きできているか確認
    74                 Signif <- symnum(pv, corr = FALSE, na = FALSE,          
    75                   cutpoints = c(0, 0.01, 0.05, 0.1), symbols = c("***", 
    76                     "**", "*"))                                         
    (略)
    
    > summary(lm(swiss$Fertility ~ swiss$Education + swiss$Examination))
    (前略)
    Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
    (Intercept)        85.2533     3.0855  27.630   <2e-16 ***
    swiss$Education    -0.5395     0.1924  -2.803   0.0075 ** 
    swiss$Examination  -0.5572     0.2319  -2.402   0.0206 *  
    ---
    Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
    (略)
  • 私の理解を越えていますので単なる推測ですが第二回答者が書いている namespace issue というのが曲者ですね。実際 printCoefmat 関数のコードを表示すると最後に <environment: namespace:stats> という行が現れます。この意味がわからないので全くの憶測ですが、print.summary.lm 関数が参照する printCoefmat 関数はある特殊な仕方で記録され(そしてそこから参照される)ているのでは? ですから表面的にそれを変えても何も変わらない? -- 2004-11-21 (日) 18:51:59
  • 今回のパズルは前にも首をかしげた R のメカニズムに関係するようですね。まず lm 関数の結果はクラス "lm" とされ、総称的関数 summary はそれにたいしてメソッド関数 summary.lm を起動(派遣, dispatch)する。次に summary.lm 関数は結果をクラス "summary.lm" とし、summary.lm 関数の返り値が(付値されないで)表示(その際、暗黙のうちに print 関数が用いられる)されるときに、総称的関数 print のクラス "summary.lm" オブジェクト用のメソッド関数 print.summary.lm を暗黙のうちに起動(dispatch)される。その際幾つかの追加の加工がされる、...等。奥が深い。後は namespace の謎ですが? とりあえず経済畑の方には手作業が必要なようです。-- 間瀬茂 2004-11-21 (日) 21:06:45
  • ありがとうございます.どうもprint.summary.lmが呼び出しているprintCoefmatが別モノではないか(namespace関係で?)と推測しているのですが,print.summary.lmの中身が見られないのでよく分かりません.とりあえず,今回は諦めます.なお,分野によって何%を有意水準にするかは,それぞれ慣習が違うと思いますので, そこをhard-codedしないでオプションすると親切だと思いました.間瀬先生,わざわざR-helpへ代理投稿して頂いてありがとうございました.-- y.a.? 2004-11-21 (日) 23:18:25
  • r-help への回答記事の続きです。-- 間瀬茂 2004-11-22 (月) 08:04:36 第二回答者のやりかたでは駄目との U. Ligges 氏(第一回答者)コメント:
    Ted, it "clashes"! Functions in the namespace are looked up at first.  
    そしてこうすれば良いという G. Grothendieck 氏(r-help の常連回答者) からの有望な提案(やはり namespace とやらを操作する必要があるそうです、私はうまくいくかどうか未確認)。こうしたいわば R の聖地にある関数を変更した場合、何か思わぬ副作用があるかも知れません。
    True, but one can still get the effect by using assignInNamespace. For
    example, run these two lines (the body(...) <- line is just for illustration
    here.  You want to ultimately replace that line with your redefined
    printCoefmat:  printCoefmat <- function... as discussed by Ted.)
    
    body(printCoefmat) <- parse(text = "cat('Greetings from printCoefmat!!!')")
    ## 次の命令が namespace を操作する魔法の呪文らしい?! 
    assignInNamespace("printCoefmat", printCoefmat, "stats")
    
    Now running summary.lm as shown below displays the desired Greetings line:
    
    R> example(lm)...snip...
    R> summary(lm.D90)
    Call:lm(formula = weight ~ group - 1)
    Residuals:    Min      1Q  Median      3Q     Max -1.0710 -0.4938  0.0685  0.2462  1.3690 
    Coefficients:Greetings from printCoefmat!!!
    Residual standard error: 0.6964 on 18 degrees of freedom
    Multiple R-Squared: 0.9818,     Adjusted R-squared: 0.9798 F-statistic: 485.1 on 2 a
  • 例えば print.summary.lm 等の関数コードがそのままでは表示されない関数のコードを見るには getS3method("print", "summary.lm") とするようです。 -- 2004-11-22 (月) 14:53:23
  • 他にも, stats:::print.summary.lm とすると内部の名前空間を呼び出せます. -- なかま 2004-11-22 (月) 16:30:21
  • G. Gro氏の提案でうまくいきますね。勉強になりました。 -- 間瀬茂 2004-11-22 (月) 17:22:54
    T. H氏のやりかたで plotCoefmat 関数(の一部)を次のように変えてみました:
    cutpoints = c(0, 0.01, 0.05, 0.1, 1),
    symbols = c("TottemoYuui!", "KanariYuui!", "MaamaaYuui!", "YuuijyaNai!"))
    次に
    > assignInNamespace("printCoefmat", printCoefmat, "stats")  
    > example(lm)
      --- snip ---
    > summary(lm.D9)
    Call:
    lm(formula = weight ~ group)
    Residuals:
        Min      1Q  Median      3Q     Max 
    -1.0710 -0.4938  0.0685  0.2462  1.3690 
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)             
    (Intercept)   5.0320     0.2202  22.850 9.55e-15 TottemoYuui!
    groupTrt     -0.3710     0.3114  -1.191    0.249 YuuijyaNai! 
    ---
    Signif. codes:  0 `TottemoYuui!' 0.01 `KanariYuui!' 0.05 `MaamaaYuui!' 0.1  `YuuijyaNai!' 1 
    Residual standard error: 0.6964 on 18 degrees of freedom
    Multiple R-Squared: 0.07308,    Adjusted R-squared: 0.02158 
    F-statistic: 1.419 on 1 and 18 DF,  p-value: 0.249 
    注意:こうして printCoefmat 関数を変更すると rm(printCoefmat) としても、今度は逆に元には戻らないようです。仕方がないからまた同じ手順で既定動作に戻すはめになりました。あまり面白がってあれこれ変更しない方が良さそうです。
  • getAnywhere(print.summary.lm) 等としてもコードが得られるようです。複数ある時は getAnywhere(foo)[2] 等とすると2番めに見つかったもののコードが得られます。。 -- 2004-11-22 (月) 18:54:21
  • namespace メカニズムはアドオンパッケージ等で、他のパッケージ等の同名の関数とオブジェクト名の衝突が起きないようにするために使われているようですね。 -- 2004-11-22 (月) 21:53:08
  • Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
    という行は,分野の違いにかかわらず,共通認識だと思います。
    記号の付け方の慣習が違うというのは私はみたことがありませんでしたので,びっくり仰天でした。
    そもそも,記号で有意性を示すのは前近代的であり,P値を明示すれば何の問題もないことです。 -- 青木繁伸 2004-11-22 (月) 22:33:26
  • namespace,奥が深いですね… どうも,ありがとうございました.なお,経済学における統計的有意の前近代的な記法については,なんでも掲示板?に少し書いておきました. -- y.a.? 2004-11-23 (火) 01:18:26

ラベルの縦書き

学生? (2004-11-19 (金) 19:30:47)

日本語版R 2.00(Win版)を使用しています.

グラフを描くに,x軸やy軸のラベルを縦書きにしたいです.
以前,掲示板に少し書かれていますが,具体的に,どこをどう設定すれば,実現するのかが分かりません.
どなたか,ご存知の方,ご助言頂ければ幸いです.

  • こんなのじゃ駄目ですよねぇ・・・ -- 2004-11-20 (土) 12:03:44
    > plot(1:10, type="n", ylab="")
    > text(3,3, "A?nB?nC?nD?nE")
  • 基本的にはこうした日本 specific な工夫は労多くして益少ないのであまりお勧めではない(R は所詮英語文化圏作品ですから)ですが、どうしても大々的にやりたければ画像を pictex に落して、それを手作業でいじることでしょうか。私はとてもやる気は起きませんが、コツさえ飲み込めばいかようにもいじれるはず。--> 勉強半分にテストをしてみました(意外に簡単かも)。 y 軸マークはお遊びです。-- 2004-11-20 (土) 12:28:48 まず R で
    x <- rnorm(10)
    pictex("test.tex")
    plot(x)
    dev.off()
    これで作られたファイル test.tex の最初に次を加える
    ?documentclass[a4paper]{jarticle}
    ?usepackage{pictex}
    ?begin{document}
    当然最後に次を加える
    ?end{document}
    次に test.tex からx,y軸ラベルに相当するところを見つけ編集、つまり
    ?put {Index}  [rB] <0.00pt,0.00pt> at 169.59 5.17
    ?put {x}  [lB] <0.00pt,0.00pt> at 9.72 154.73 
    をたとえば
    ?put {インデックス}  [lB] <0.00pt,0.00pt> at 169.59 5.17
    ?put {?shortstack[r]{ス??[-7mm] ク??[-7mm] ッ??[-7mm] エ}}  [rB]
           <0.00pt,0.00pt> at 2 154.73 
    に変更する。at 以下はテキストをおく位置の x,y 座標ですから、これをいじくれば見栄えの良い位置になります。?shortstack 内のテキストはさかさまにします。??[-7mm] は垂直間隔を調整するためです。さらに微妙な調整をするには pictex の命令を勉強する必要がありそうです。画像ファイルが欲しければ dvips test.dvi で ps ファイルにし、さらに Imagemagick 等の画像ソフトで最終的に欲しいフォーマットに変換する、でしょうか。特に複雑な画像はひょっとすると latex コンパイル時にトラブルが起きるかも知れません。参考画像は latex の一頁分そのままです(でっかく、見苦しくて済みません)。x 軸マークはお遊びです、どうやったかはもうおわかりだと思います。
    test.png
  • 丁寧なご解説ありがとうございました.MS Wordしか使ったことがなかったので,TeXの勉強から始めます. -- 学生? 2004-11-21 (日) 20:34:48
  • TeXの勉強から始めるくらいなら,PowerpointとかOpenOffice?.orgのDrawとかに貼り付けて編集する方が早いんではないかと思いますがどうでしょう? -- 2004-11-21 (日) 23:22:56
  • おかげさまで,emf形式で保存して,MS Power pointで,「挿入」→「図」→「ファイルから」と押して,保存したemfを挿入した後,画像をダブルクリックして「MS office 描画オブジェクト」に変換することで,ラベルの縦横を変更できました.少し疑問なのですが,MS wordで同じ処理をすると,なぜか文字化けしてしまいますね. -- 学生? 2004-11-22 (月) 01:47:02
  • 文字化けのWORDのバージョンにもよると思いますが,ツール+オプション+編集?の図の編集は何になっていますか?,選択肢が他にあるならそれに変更してみてください.私も組合せの関係もあるのではっきりとは断言出来ませんが,わーどのバグではないかと. -- なかま 2004-11-22 (月) 13:21:31
  • Word2000を使用しています.ツール+オプション+編集と日本語入力の「図の編集ツール」を 「Microsoft Word」から「Microsoft Draw 98」に変更することで文字化けが消えました.ありがとうございました. -- 学生? 2004-11-22 (月) 13:31:22
  • 「Drawとかに貼り付け」とか「emf形式で保存」とは MSW だけで通用する技でしょうか? Linux の OpenOffice? でも使える技でしょうか? -- 2004-11-22 (月) 20:33:05
  • emf自体はOpenOffice? for Linuxでも編集可能です. -- なかま 2004-11-22 (月) 23:46:34

Decisiontreeはどうすればできるのでしょうか?

MI? (2004-11-19 (金) 17:51:39)

RでDecisiontreeができると人づてに聞きました。
調べるとtreeのような感じですが、よくわかりません。
どなたか、ご存知ないでしょか?

  • 使ったことはありませんが,rpartというパッケージでできるらしいです。Maindonald and Braun "Data analysis and graphics using R"のChapter 10に詳しい記載があります。 -- 中澤? 2004-11-19 (金) 18:25:25
  • Jin's Pageにも解説があります. -- 2004-11-19 (金) 19:20:17
  • Jin's Pageっていいですね。 -- 2004-11-19 (金) 22:47:38
  • 読みやすいですよね. -- 2004-11-19 (金) 23:03:17
  • 忘れてましたが,英語に抵抗がなければ,http://www.liacc.up.pt/~ltorgo/DataMiningWithR/PDF/DataMiningWithR.pdf も参考になると思います -- 2004-11-20 (土) 00:41:20
  • ありがとうございます。Jin'sPageは丁寧にまとめてあり初級者には非常に参考になります。 -- 2004-11-20 (土) 10:36:38

Scatterplot3Dについて

みち? (2004-11-15 (月) 16:12:57)

色々調べてみても分からなかったので、こちらに投稿させて頂きました。

1.Scatterlot3D で Type=h として、プロットした点から底面に垂直線を引かせた後、点を滑らかな線で繋ぐにはどうしたら良いのでしょうか。

2.Scatterlot3Dで軸ラベルを回転させたり移動させたりするにはどうすればよいのでしょうか。

宜しくお願いします。

  • type=h でプロットした後、par(new=TRUE) として、今度は type=l で重ね書きすれば一応それらしくなりますね。見栄え良くするには線の太さや色を変える必要があるかも知れません。いずれにしても type = l は点の間を直線で結だけです。回転・移動は結構不自由だと思います。angle パラメータを変えて再実行すれば、疑似的に回転できますが、お気にめすかどうか。本格的な回転機能を持つグラフィックス関数には rgl (グラフィックス参考実例集参照)がありますが、これはイメージ図を対象にしているので scatterplot3d のような三次元棒グラフを画くには相当苦労が必要な気がします。 -- QDU? 2004-11-15 (月) 17:47:33
  • ありがとうございます。やってみましたら、それらしくなりました。暫くこれでやってみます。 -- みち? 2004-11-15 (月) 19:49:45

複数条件下でのベクトルの取り出し

P? (2004-11-13 (土) 18:14:19)

初歩的な質問と思いますが,あれこれと調べてどうしても分からなかったので,質問させていただきます.

例えば,dataという配列もしくはデータフレームがあったとして,A列が1,B列が2を満たすC列を取り出したい場合,
data[,3][data$A==1][data$B==2]
では上手く行かないのですが,どのようにすればよいのでしょうか.
ご教示いただければ幸いです.

  • たとえば行列なら -- QDU? 2004-11-13 (土) 18:28:11
    > x = matrix(c(2,1,2,2,2, 3,4,3,4,3, 5,6,7,8,9), 5,3)
    > x     
         [,1] [,2] [,3]
    [1,]    2    3    5
    [2,]    1    4    6
    [3,]    2    3    7
    [4,]    2    4    8
    [5,]    2    3    9
    > x[x[,1]==2 & x[,2]==3, ]     
         [,1] [,2] [,3]
    [1,]    2    3    5
    [2,]    2    3    7
    [3,]    2    3    9
    > x[x[,1]==2 & x[,2]==3, 3]
    [1] 5 7 9
    データフレームなら
    > x = data.frame(C1=c(2,1,2,2,2), C2=c(3,4,3,4,3), C3=c(5,6,7,8,9))
    > x
      C1 C2 C3
    1  2  3  5
    2  1  4  6
    3  2  3  7
    4  2  4  8
    5  2  3  9
    > x[x$C1==2 & x$C2==3,] # 結果はデータフレームになる
      C1 C2 C3
    1  2  3  5
    3  2  3  7
    5  2  3  9
    > x[x$C1==2 & x$C2==3,]$C3
    [1] 5 7 9
  • 大変参考になりました. -- P? 2004-11-13 (土) 20:33:49

macでソースファイルからパッケージをインストールする方法

ささかま? (2004-11-12 (金) 20:33:39)

R2.0.0 mac版を使用しています。OSは 10.3.6です。xcode1.0を入れてあります。
lpSolveというパッケージを導入したいのですが、CRANにはmac用のバイナリはありません。Package Installer を使ってソースからインストールしようとしましたが、
In file included from commonlib.c:17:
commonlib.h:6:20: malloc.h: No such file or directory
make: *** [commonlib.o] Error 1
ERROR: compilation failed for package 'lpSolve'
というエラーが出てしまい、出来ませんでした。
ソースを手動でダウンロードし、ターミナルから
R CMD INSTALL lpsolve
と実行しても同様のエラーが出てしまいました。ルートになってもだめでした。
Rを使い始めて数週間ですし、UNIXは素人です。色々と検索しては見たのですが、解決できませんでした。どなたか解決法をご存知の方、教えていただけないでしょうか?

  • /usr/include/malloc/malloc.h を /usr/local/include にコピーあるいはシンボリックリンクしておけば、コンパイルは通ります。ターミナルからだと、ln -s /usr/include/malloc/malloc.h /usr/local/include/ でしょうか。 -- 2004-11-12 (金) 23:19:05
  • ありがとうございます。教えて頂いたようにしましたら無事インストールできました。不可視ファイルなのでFinderからは見えなかったのですね。他にもmac用バイナリの無いパッケージ(rimageやtsiriesなど)があるようですので、UNIXの勉強をして入れてみようと思います。また質問するかも知れませんが、よろしくお願いします。 -- ささかま? 2004-11-13 (土) 17:58:55

RExcelのインストール

NK? (2004-11-10 (水) 14:44:56)

R(D)COM Server V1.2を今まで使っていて、V1.2をアンインストールして残っていた関係したファイルも全て削除してからV1.35をインストールしましたが,インストールの途中で”C:WINNT?System32?msvcirt.dll Restart Replace failed Move FileEx? Failed code 5"というエラーが出てしまい正常にインストールできません。どのように対応すれば良いでしょうか?宜しくお願いします。

  • ロード済のDLLを置き換えられないと言っているので該当ファイルをリネーム後に再度インストールすれば良いかと. -- なかま 2004-11-13 (土) 23:35:14
  • お忙しいところ有り難うございました。お蔭様でインストールできました。 -- NK? 2004-11-15 (月) 10:49:42

出力画像のサイズを変更したい

fujim? (2004-11-04 (木) 22:15:36)

初めまして。近頃Rをはじめようと思い、R-tips,データサイエンス入門
などを読みつついろいろと試しているものです。
環境:Windows2000,R2.0.0jp 実行方法:test.Rファイルを作成し.batファイルより R.exe CMD BATCH --slave --vanilla test.Rで実行
分からない点:例えばグラフィックス参考実例集barplotの項のbarplot03.R等において、出力画像.pngのサイズをどのように変更すればよいのでしょうか?
調べた点:parパラメータを以下のように変更してみた。
1) par(pin=c(5,10))
2) par(fin=c(5,10))
いずれも.outファイルにエラーが出て画像が出力できなくなります。
エラーは以下の通りです。
Error in plot.new() : Figure region too large
Execution halted
Warning message:
calling par(new=) with no plot
宜しくお願いします。

  • 見つけました。済みません。png("outfilename.png",width=600)とかするんですね。盲点でした失礼しました。 -- fujim? 2004-11-04 (木) 22:30:10
  • 便乗質問ですが、pngの場合はwidthやheightの単位がpxになっていますが、postscriptの場合はinchです。pngの出力指定単位がpxと知らずにいつもどおりinchのつもりで指定して半日費やしたことがあります。どこかにデフォルトの設定単位の一覧表はありませんか。 -- 2004-11-06 (土) 13:11:43
  • こんな感じでしょうか.
関数単位default width & heightdefault pointsize
bitmap()inch6‐(適当)
bmp()pixel48012
dev2bitmap()inch6
jpg()pixel48012
pdf()inch6
pictex()inchw=5, h=4なし
png()pixel48012
postscript()inch
quartz()inch512
windows()inch712
win.graph()inch712
win.metafile()inch712
win.print()inch712
x11()inch712
X11()inch712
xfig()inch

本に載せるいいネタを戴きました.ありがとうございます! -- 舟尾? 2004-11-06 (土) 13:53:27

スクリプト実行

taku? (2004-11-04 (木) 14:19:56)

perlやrubyの様に、一行目に次のようなRのパスを書いて、
一連のコマンドをスクリプトとして実行したいのですが、(ファイル名hello)

#!/usr/bin/R

ARGUMENT './hello' __ignored__
というメッセージが出て、普通にRのプロンプトが出てしまいます。

Rプログラムをスクリプトとして実行するにはどうしたらいいのでしょうか?

  • スクリプトとして実行なんかできるんですか?
    まずもって,R --help で,オプション確認 -- 2004-11-04 (木) 14:27:08
  • この話は r-help 等で繰返し議論されています。一度 batch 等のキーワードで検索してみて下さい。 -- 2004-11-04 (木) 17:13:12
  • いや,だから。バッチと「Rプログラムをスクリプトとして実行する」というのは違うでしょうということです。
    taku さんは,どういうスクリプトファイルを用意したのですか?hello の内容が
    #!/usr/bin/R
    x <- rnorm(2000)
    print(mean(x)
    などとして,それを chmod 755 hello なんかしても,./hello では実行できませんよね。
    例えば prog というファイルにやりたいことを書いておいて
    $ cat prog
    x <- rnorm(2000)
    print(mean(x))
    そこから入力して R を実行する
    $ R --vanilla --quiet < prog
    > x <- rnorm(2000)
    > print(mean(x))
    [1] -0.007029891
    こういうことをやりたいのでしょう?
    quiet を slave にすれば,入力ファイルからの入力行は出力されない。とにかく,そのようなことなら,R --help すべし。 -- 2004-11-04 (木) 17:29:02
  • PerlスクリプトからRを起動するなどがヒントになるでしょうか。 -- 2004-11-04 (木) 21:52:56
  • いろいろありがとうございます。 どうやら、#!/usr/bin/Rというのはできないみたいですね。 perlやrubyを使う身としては、実行のたびにいちいち'R --vanilla...'とタイプしたくはないので、 このあたりは何とかして欲しいものですが… 結局、今のところはshを使うのがよさそうです。
    #!/usr/bin/sh
    R --vanilla --slave <<EOF
    x <- rnorm(50)
    y <- rnorm(x)
    plot(x,y)
    EOF
  • 面倒ならば,
    alias rb        'R --vanila --slave ?!*'~
    しておいて,
    rb < ファイル名
    で動くようにしておいたらいかが。
    sh のなかでいちいち
    #!/usr/bin/sh
    R --vanilla --slave <<EOF
        なんたら
    EOF
    なんて書く方が面倒だ。 -- 2004-11-05 (金) 11:12:45
  • ああ、こういう人もいるんじゃないかと思ってはいましたが。根本的な解決ではない、という点では同じでしょ。#!/usr/bin/Rとできない時点で面倒は避けられないんですよ。 -- taku? 2004-11-05 (金) 23:51:53
  • どうぞ。おすきなように。Cプログラムを記述してコンパイルして実行するときも同じような手順を踏みたいですか? -- 2004-11-06 (土) 19:04:24 内容が
    R --vanilla --slave < $1
    な実行ファイルを/usr/bin/Rsとでもしておけば,
    #!/usr/bin/Rs
    で動く. -- 2004-12-01 (水) 11:15:44
  • 今は Rscript があります。#!/usr/bin/Rscript -- 2009-02-19 (木) 04:28:12

関数の中身を見る関数

akosuke? (2004-10-29 (金) 10:45:01)

自分で作った関数で名前はわかるけど中身を忘れてしまった、、、というような場合に関数の中身を見るための関数はありますでしょうか?
ここでも調べてみたのですが調べ方が悪いのか見つかりませんでした。
くだらない質問で申し訳ありませんがどうぞよろしくお願いします。

  • コンソールにその関数の名前だけを入力してみる.... これは,自作のものに限らない.... ls() ではなく,ls って入力してみそ.... -- 2004-10-29 (金) 11:00:27
  • おわかりになりましたか? -- 2004-11-04 (木) 14:36:30

AICの計算

AIC? (2004-10-26 (火) 21:50:37)

The R Bookの第5章を読んでいる者です。p.109の上から2行目で「ARIMA分析で行ったようにプログラムを組む・・・」と書いてあるのですが、このARIMAモデルではq=0としてpのみ動かすという形にしていますよね。ところがGARCHでパラメタp,qを2つ動かしてAICを最小にするp,qを求めるにはどのようにプログラムを組んだら良いのかが分からないので、教えていただけると幸いです。
また、AICのパラメータ数はこの場合npram<-2と定義していますが、GARCHの
場合、式(5.14)にあるようにω、α、βの3と定義するべきではないのでしょうか?
それから(長くてすみません)GARCHの派生モデルを"fSeries"というパッケジで分析できると聞き、garchと同じ要領で書けばよい(例えばGJRなら
gjr(data,order(p,q,γ))??)と思ってやってみたのですができません。
ご存知でしたら教えていただけないでしょうか。

  • 本来ならば「The R bookご意見」でいただくべきご質問ですが、矢野に余裕がないのでこちらで返答させていただきます。 -- 矢野? 2004-10-27 (水) 09:56:47
  • まず1番目のご質問は以下のようにしてください。
    for (p in 1:3){
      for ( q in 1:3){
        tmp <- garch(dnikkei, order=c(p,q))
        cat("p=", p, "q=",q, ", AIC=", (tmp$n.likeli+p+q),"?n")
      }
    }
    warningなどが出て少し見づらいですが、とりあえず。-- 矢野? 2004-10-27 (水) 09:57:52
  • 2番目のご質問に対しては、AICはモデルの当てはまりの良さに対する相対的な値なので、ωの分(つまり+1)は必要ありません。たとえば、上のプログラムでは(tmp$n.likeli+p+q)としてありますが、これを(tmp$n.likeli+p+q+1)としなくても一貫していればAIC同士を比較するに当たっては特に違いはないということです。それとよく使われるAICの定義は2*(tmp$n.likeli+p+q)なのですが、これも次数選択(モデル選択)という観点からはそうしなくても特に支障はありません。 -- 矢野? 2004-10-27 (水) 10:12:26
  • 3番目のご質問のfSeriesに関しては詳しくないのですが、aparchFitでoptのgammaを設定すればいけるのではないかと思います。helpをよく読んでチャレンジしてみてください。 -- 矢野? 2004-10-27 (水) 10:15:07
  • ありがとうございます。上のプログラムで実行した結果、1番目と2番目の質問は解決しました。ご丁寧な説明に感謝です。AICに関しては確かに1はなくてもいいですよね。 -- AIC? 2004-10-27 (水) 22:25:11
  • 3番目の質問に関しては、fSeriesをインストールしたあと、どのようなコマンドで実行するのかが分かりません。help("fSeries")をみてもGarchModelling?()とありますが、これを実行してもObject "garchmodelling" not foundとエラーが出るので。。 -- AIC? 2004-10-27 (水) 22:33:29
  • fSeriesを用いたGJR推定は以下のようになります。 The R book第5章のデータを使った例(GJR(1,1))を示します。
    library(fSeries)
    tmp <- aparchFit(x=dnikkei,
           order = list(alpha.lags = 1, beta.lags = 1, delta = 2),
           opt = list(gamma = TRUE, delta = FALSE,
           disparm = FALSE), doprint = FALSE)
    tmp
    ポイントは上でgamma=TRUEとなっている部分です。この 部分がGJR定義式でのγを示します。gamma=FALSEとすると ただのGARCHになります。GJRで日経225を推定すると非対称性があることが 分かります。aparchFitについてより詳しくは
    help(aparchFit)
    としてhelpの内容を参照してください。-- 矢野? 2004-10-28 (木) 18:51:27
  • 懇切丁寧な説明、ありがとうございます。頑張ってみます。 -- AIC? 2004-10-28 (木) 22:53:18

0を埋める

(2004-10-25 (月) 18:52:25)

とってもFAQだと思うのですが…
1 -> 001,42 -> 042
のように有効桁数(この例だと3ケタ)を決めて,
それに満たない場合は0で埋めるにはどうすればいいでしょうか?

  • sprintfでこのサイトを検索してみて下さい。 -- 谷村 2004-10-25 (月) 19:02:46
  • ありがとうございました.formatC という関数がありました:
    x <- c(1:20)
    formatC(x,width=2,flag="0")
  • 関数の紹介ありがとうございます。
    > x <- c(1:20)
    > formatC(x,width=2,flag="*")
    flag にいろいろな文字を指定するとおもしろい結果が出てきます。
    本来は,無効な文字ははじかなくてはならないのでしょうが。
    ベクトル化されているという点では sprintf より使いでがあるのかもしれません。 -- 2004-10-25 (月) 21:34:15

変数名に別の変数の内容を含めたい

質問者? (2004-10-23 (土) 15:43:00)

大変初歩的な質問と思いますが、うまくいかないのでよろしくお願いします。

month <- c("Jan","Feb","Mar","Apr","May","Jun",
           "Jul","Aug","Sep","Oct","Nov","Dec")
for (i in 1:12) paste("Hoge.",month[i],".total",sep="") <- 0

として、例えばHoge.Jan.totalに0が格納されるようにしたいのですが、paste()の結果がクオートされているせいか、Error: Target of assignment expands to non-language objectとエラーになります。ご指導をお願いします。

  • RプログラミングTips大全 に定石構文が紹介されています。-- 2004-10-23 (土) 17:06:03
    > month <- c("Jan","Feb","Mar","Apr","May","Jun",
                 "Jul","Aug","Sep","Oct","Nov","Dec") 
    > for (i in 1:12) eval(parse(text=paste("Hoge.",month[i],".total <- 0",sep="") ))
    > Hoge.Apr.total
    [1] 0
    ついでですが次のようにしても良いわけです。
    > for (i in 1:12) eval(parse(text=paste("Hoge",month[i],"total <- 0",sep="."))) 
    また R は月名を表す組み込みの文字列ベクトルを最初から持っています。
    >month.abb
    [1] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep" "Oct" "Nov" "Dec"
    > month.name
    [1] "January"   "February"  "March"     "April"     "May"       "June"
    [7] "July"      "August"    "September" "October"   "November"  "December"
  • ありがとうございました。evalも使って試行錯誤もやってのですが、やり方が悪かったみたいです。また月名のベクトルは参考になりました。 -- 質問者? 2004-10-23 (土) 19:56:52

データフレームの分割処理について

R-Fan? (2004-10-18 (月) 14:38:21)

このサイトでRを知りました。いつも有益な情報をありがとうございます。
表題についてですが、大きなデータをキー項目(下記例ではdate)ごとに分割して処理するプログラムを作っています。subset()によるデータ抽出は大変便利ですがデータが大きくなると処理にかかる時間も増えてしまいます。この処理を高速化したいと考えておりますが、データフレームをキー項目で分割するような組み込み関数はないでしょうか。それともデータが大きくなる場合は、前処理段階でデータファイルを分割しておくのが良いでしょうか。全体のデータサイズは20〜30万件で、取り出すデータは1000〜2000件です。

(案1)
for (key in sort(unique(maindata$date))){
 subdata <- subset(maindata, date == key)
 #subdataの処理
}

自己レスになりますが、subset()を使用するかわりにキー項目に対応する行の索引を作成することで処理時間を4割程度削減することができました。

(案2)
idx <- tapply(1:nrow(maindata), maindata$date,
              function(subidx){return (subidx)}) 
for (key in sort(unique(maindata$date))){
  subdata <- maindata[idx[[toString(key)]],]
  #subdataの処理
}

もっと根本的に改善する方法がありましたらご教授いただければ思います。 どうぞよろしくお願いいたします。

  • ケースにもよるでしょうけど,unique(maindata$date)をlevels(maindata$date)にmaindata[idxtoString(key)?,]をmaindata[levels(maindata$date)==key,]とかではいかがでしょう. -- なかま 2004-10-18 (月) 20:45:47
  • ご回答いただきましてありがとうございます。後者のご提案についてですが、下記のようなコードでは因子1を持つ行を取り出すということはできないようです。想定された使い方と違っているでしょうか。しばらくヘルプと格闘してみますが、チューニングについてはあまり無理をせず使いながら勉強したほうがいいのかもしれませんね。 -- R-Fan? 2004-10-19 (火) 14:16:26
    > d <- data.frame(x = factor(c(1,1,1,2,2,3)), y = 10:15)
    > d[levels(d$x) == 1,]
      x  y
    1 1 10
    4 2 13
  • 次のようになるのでは?
    > d[levels(d$x)[1]==d[,1],]
      x  y
    1 1 10
    2 1 11
    3 1 12
    > 
    > d[levels(d$x)[2]==d[,1],]
      x  y
    4 2 13
    5 2 14
    > d[levels(d$x)[3]==d[,1],]
      x  y
    6 3 15
    データファイルが分析のたびに付け加わっていくのでなければ,事前にファイル分割しておいた方がよいような気はします。
    そうそう,思い出した。by という関数,速いか遅いか知りませんが,記述は簡潔。 -- 青木繁伸 2004-10-19 (火) 14:29:52
  • levels()関数の例は大変よく分かりました。ご紹介いただいたby()関数はまさに今回の処理にうってつけです。最終的に以下の形になりましたが、簡潔なだけでなくパフォーマンスも良好で(案2)より若干速いようです。これ以上は望めないような気がします。いろいろとありがとうございました。 -- R-Fan? 2004-10-20 (水) 02:03:11
    (案3)
    by(maindata, maindata$date, 
      function(subdata){
      #subdataの処理
      }
    )

関数の読み込み

恐縮です? (2004-10-17 (日) 00:24:05)

Rを使い始めて間もない者です。
Rで色々な関数をinstall.packages()でインストールするのですが、ワークスペースを閉じるたびにまたインストールし直しになります。
一度インストールした関数をそのまま使い続けるにはどう保存すればよいのでしょうか?宜しくお願いします。

  • library(パッケージ名) という使い方ができないということですか。もう少し,状況を説明しないと,答えをもらえないと思いますよ。 -- 青木繁伸 2004-10-17 (日) 08:25:30
  • すみません。library(パッケージ名)を実行して使えるようにはなるのですが、Rを終了すると、たとえ終了時にワークスペースの保存をしたとしても、つぎにRを立ち上げた時にもう一度install.packages()とlibrary()を実行しなければなりません(私の場合)。そうしなくても良い方法をお聞きしたいということです。背景説明が至らず、申し訳ありませんでした。宜しくお願いいたします。 -- 恐縮です? 2004-10-17 (日) 20:50:47
  • library は毎回,必要なときにやらないといけないでしょう。それがいやなら,自前の関数のように扱って,その環境を保存しておけばいいわけではないでしょうか。 -- 青木繁伸 2004-10-17 (日) 21:17:15
  • OSを書かれていませんので、Linuxを前提に説明するとinstall.package()は一般ユーザではなく、管理者権限(root)で行う必要があります。Windowsの場合はadministratorでしょうか。 -- 2004-10-17 (日) 21:18:33
  • library()を毎回実行するのが面倒なら、青木先生の助言のように関数を作成するか、作業ファイル(例えばhoge.r)のテンプレートにlibrary()を書いておいてエディタで作業する際に自動的に読み込むようにするとか。 -- 2004-10-17 (日) 21:22:05
  • 普通,ライブラリのインストール時には管理者権限でやるのでは? -- 青木繁伸 2004-10-17 (日) 21:22:07
  • .First()関数に初期化処理を定義しておく方法はだめでしょうか? -- R-Fan? 2004-10-18 (月) 14:57:22
  • 参考URL http://www.okada.jp.org/RWebRef/RManualJP/Startup.jp.html -- R-Fan? 2004-10-18 (月) 15:02:50
  • R 2.0 からは、たしかパッケージは実際に使う際になって初めて読み込まれる仕組み (Lazy Loading) に代わったので、初期化ファイルで使いそうなものはすべて library 関数で読み込むことにしておいても、支障はなさそうです(実際に確かめたわけではありません)。 -- 2004-10-20 (水) 10:01:39

積集合(INTERSECT)と和集合(UNION)の中間的なものはありますでしょうか?

撃墜王? (2004-10-14 (木) 15:39:55)

R言語には、INTERSECTとUNIONといったコマンドが存在します。
ただ、入力されるリストが非常に多く存在した場合に、「それら全てのINTERSECT」を取ると、最終的なアウトプットがほぼ0になってしまいます。
また、「それら全てのUNION」を取ると、当然ながらアウトプットが膨大なものになってしまいます。。。どちらを採択するにしても大変悩ましく感じています。
そこで、質問させて頂きたいのは、INTERSECTとUNIONの中間的のようなもの、もしくは考え方などがあったら是非アドバイスいただければと思います。

なお、自分の中のイメージとしては

入力されるリスト 100個
1リストあたりに含まれる変数名(Aomori,Miyagi,Hokkaido...etc) 100個
全てのINTERSECTを取った場合 最終的に残る変数名 3個(少なすぎる。。。)
全てのUNIONを取った場合、 最終的に残る変数名 1024個(多すぎる。。。)

そこで、単純なUNIONではなくて、1つのリストにしか含まれないような変数名は除くようにしたいと思っています。
具体的な例をさらにあげますと

list1 Kanagawa,Tokyo,Saitama
list2 Chiba,Kanagawa,Tokyo
list3 Tokyo,Irabaki,Gunma

Intersect = Tokyo,Union=Kanagawa,Tokyo,Saitama,Chiba,Ibaraki,Gunma
となってしまいますが、上のようなリストの場合Kanagawaはlist1,2に共通して含まれ、唯一list3には含まれないような形なので、これは採用するといったようなUnionの考え方とIntersectの考え方をいっしょくたにしたような方法と言うのは実現可能でしょうか。色々調べたのですが、よい情報がみつからなかったので、諸先輩方アドバイスの程どうぞよろしくお願いいたします。

  • ある要素が,いくつのリストに含まれているかを勘定し,リストの何パーセント以上に含まれているものという指示で,条件を満たす要素を返せばいいのかな?
    > top.n <- function(lst, n)
    + {
    + result <- table(unlist(lst))
    + result[result >= n]
    + }
    > lst <- list(no1=c("aa", "bb", "cc"),  
                   no2 =c("bb", "dd", "ff", "hh"),
                   no3=c("cc", "ff", "xx", "yy"), 
                   no4=c("aa", "cc"))
    > top.n(lst, 3)
    cc 
     3 
    > top.n(lst, 2)
    aa bb cc ff 
     2  2  3  2 
    探すのがめんどいから作っちゃった。関数にするまでもないが。仕様が詳細に決まれば,お化粧はいかようにも。 -- 青木繁伸 2004-10-14 (木) 16:57:25
  • 青木先生有難うございます。先生の手に掛かるとあっというまですね。本当に頭の下がる思いです。さっそく試してみたいと思います。その他、先生の独自スクリプト以外にも何か良いアイデア、考え方などありましたら示唆いただけると幸いです。m(_ _)m -- 撃墜王? 2004-10-15 (金) 13:47:39

CSVファイルのインポートで errorが。

竹内? (2004-10-10 (日) 19:42:29)

こんにちは。Rを使い始めてまだ5日目の初心者です。使用環境は Mac OS X(10.3)で Rの versionは 1.91(英語版)です。 R-Bookの p59にある read.csv("ファイル名 ")で CSVファイルを読み込ませようとしても次のようなメッセージが出てしまいます。

Error in file(file, "r") : unable to open connection
In addition: Warning message:
cannot open file `test.csv'

R-Bookや googleで一応は調べた(つもり)なのですが、 errorが何を意味しているかわかりません。どなたかご教授をお願いします。

  • ファイルが開けないといっているのです。どこにそのファイルがあるかRがわからない。というか,Rが見ているところにそのようなファイルはないということですね。
    Tools --> Change Working Directory で,そのファイルがあるディレクトリを指定しましょう。デスクトップかどこかにいつも使う作業用のディレクトリを作っておいて,Preferences --> Misc でいつもそのディレクトリを使うように設定しておいてもいいですね。-- 青木繁伸 2004-10-10 (日) 19:46:02
  • さっそくのご回答ありがとうございました。教えていただいたようにやりましたところ、前述のような errorは出なくなりました。が、別の error(ファイルのヘッダーがおかしい)がでました。この件についてはもう少し自分で試行錯誤しようと思っています。もし、自分で解決できなかった場合にはまた相談させていただくかもしれません。 -- 竹内? 2004-10-12 (火) 07:53:49

frequencyについて

R初心者? (2004-10-07 (木) 23:49:37)

Rで株価の時系列分析をしようとしている者ですが、frequencyについて、年次のデータであれば1を、月次のデータであれば12とすればよいと分かるのですが、株価のように日次データでしかも日曜日の値がない場合、frequency=365とするとずれが生じて横軸の値が正確でなくなるのですが、どのように解決すればよいのでしょうか?

  • 土日が無いなら freq = 5 とする他無いのでは。祭日は NA でしょうね。-- 2004-10-08 (金) 00:06:01
    > ts(1:10, frequency = 5)
    Time Series:
    Start = c(1, 1)  # 始点は第一週の第1日と解釈
    End = c(2, 5)    # 終点は第二週の第5日と解釈
    Frequency = 5  
    [1]  1  2  3  4  5  6  7  8  9 10
  • ご質問は株価データが日次データ(土日祝日データなし)が何年分も並んだ場合の扱いに関するご質問かなぁと推測します。実はあまり根本的な解決方法はなくて、frequencyを250から255あたりにしておくとある程度はごまかせます(というか矢野はいつもそうやってごまかしています) -- 矢野? 2004-10-08 (金) 00:47:10
  • 御返事有難うございます。矢野先生のおっしゃる通り株価が何年分も並んだデータに関する場合です。frequencyを250あたりにすると確かにある程度は正確になりますが、その理由はなんなのでしょうk  -- R初心者? 2004-10-13 (水) 22:15:28
  • えーっと、まず「先生」というのはやめてください。先生ではないですし偉くもありませんので。 -- 矢野? 2004-10-14 (木) 09:50:32
  • それでfrequencyを250にすると良い理由ですが、それは簡単で日本で株式市場が開かれている日数が250日弱(正確には246日から248日くらい。年によって少し違います)だからです。要は365日から土日祝日と年末年始を除いた営業日が250日弱なのです。 -- 矢野? 2004-10-14 (木) 09:54:16
  • どうもありがとうございます。納得しました。ただ何年か前までは年間営業日が280日くらいだったみたいで1980年くらいから2004年までのデータだとfrequency=270くらいがちょうど良いですね。勉強してみます。 -- R初心者? 2004-10-15 (金) 23:08:50

GARCHモデルのパラメータ推定

ピエロ? (2004-10-07 (木) 22:18:58)

RでGARCHモデル分析する際、パラメータ推定は何の方法を用いているのでしょうか?
擬似最尤法とかいうやつでしょうか?

あと、Rでの時系列分析の細かな手法を扱っている本あるいはサイト等ございましたら、教えて頂けると幸いです。
宜しくお願い致します。

  • ソースを御覧になるより仕方がないのでは。こうした質問にすぐ答えられるほど日本の R ユーザーは多くありませんよ(と私はおもう)。頑張って本家の r-help で聞かれてはいかがでしょうか。 -- 2004-10-08 (金) 00:12:05
  • help(garch)とすると"maximum-likelihood estimates of the conditionally normal model"とでますね。実際の計算方法は参考文献にあるとおもいます。 -- 矢野? 2004-10-08 (金) 00:28:39
  • それと時系列解析の大まかな手法でよければ、矢野が書かせていただいたthe R bookの第5章が参考になるかと思いますが、細かい手法は以下の慶応大学の小暮先生のサイトがとてもよいと思います -- 矢野? 2004-10-08 (金) 00:31:56
    時系列解析法:http://web.sfc.keio.ac.jp/~kogure/courses/2003/ts/ts.html
  • それと上のコメントの方がおっしゃるとおり、Rの日本でのユーザはまだあまり多くない(特に時系列解析)ように感じることが多いので、みんなで一緒に頑張りましょう -- 矢野? 2004-10-08 (金) 00:43:17
  • それとRjpの以下のページは非常に便利ですばらしいです(作成された方に心より感謝いたします)。 -- 矢野? 2004-10-08 (金) 00:58:23
    Rの基本パッケージ中の時系列オブジェクト一覧
  • ↑を作成されているのは間瀬先生ですね。ご尽力ありがとうございます! -- 2004-10-08 (金) 12:44:56
  • 間瀬先生でしたか。改めてありがとうございます。 -- 矢野? 2004-10-08 (金) 15:43:43
  • 実は私の先生は時系列の専門家(藤井光昭先生)なのですが、私は時系列には全く興味が持てなかった不肖の弟子です。あまり関係のないコメントですが、要するに「編集内容に責任持ち(て)ません」。 -- 間瀬茂 2004-10-09 (土) 14:35:01
  • 矢野先生の書かれた本の第5章、読ませて頂きました。5.7.2節のGARCHモデルのプログラムで、plot(predict(ng))とありますが、これは何を予想しているのでしょうか?縦軸のseries1,2の意味も分からないのですが。あと、GARCH以外の手法(EGARCHとかGJRとか)の分析手法はあるのでしょうか? -- ピエロ? 2004-10-13 (水) 22:33:50
  • こちらも同じで、「先生」というのは勘弁してください -- 矢野? 2004-10-14 (木) 09:56:40
  • plot(predict(ng))というのはパラメータを決定したgarchモデルを使って将来の条件付標準偏差を予測しています。縦軸のseries 1, 2は予測の+値と-値がそれぞれ入っています(なぜそういう仕様になっているかはちょっと僕にも分からないです。helpにそういう仕様だと明確に書いてあるので、何か理由はあるのではないかと思いますが)。また、garchの派生モデルをRで見たことはないです。ご存知の方がおられたら教えてください。 -- 矢野? 2004-10-14 (木) 10:00:38
  • 御返事有難うございます。矢野さん。上記のplot(predict(ng))を今一度出力したのですが、将来の条件付標準偏差というのは、何日先の標準偏差のことなのでしょうか?縦軸は2003年1月から12月までになっているので、元データと変わらない気がしますが・・。お教え頂けると幸いです。 -- ピエロ? 2004-10-15 (金) 21:54:51
  • ちょっとは,自分で調べましょう。 -- 2004-10-15 (金) 23:27:49
  • 私の勉強不足で初歩的な質問になってしまい、申し訳ありませんでした。 -- ピエロ? 2004-10-17 (日) 23:26:16
  • 縦軸が2003年1月から12月までになってしまうのは昔からあるバグだと思います。そのうちに直るだろうと思って放っているのですが、まだ直っていないようです。本来ならばr-helpか開発者に連絡すべきなのですが、余裕がなくそのままになっています。 -- 矢野? 2004-10-27 (水) 10:23:56

pairsのプロット範囲

ぷーすか? (2004-10-07 (木) 20:17:03)

散布図行列で、プロット点をドットに変え、y=xの直線を書きもうとして、

pairs(x,panel=function(x,y){points(x,y,pch=".");abline(0,1,col=2)})

でうまくできました。
ここでさらに各図のプロット範囲を同じにしようとxlim、ylimの指定を追加しようと思いましたが、

pairs(x,xlim=c(-3,3), ylim=c(-3,3), panel=function(x,y){points(x,y,pch=".");abline(0,1,col=2)})

ではエラーが出てしまいます。
プロット範囲以外のカスタマイズはあきらめ、

pairs(x,xlim=c(-3,3), ylim=c(-3,3))

とすれば、warningは出ますが、一応プロット範囲は制御できました。

プロット内容とプロット範囲のカスタマイズを両立させる方法はないのでしょうか?

  • エラーメッセージはパネル関数は追加パラメータを受け付けないといっているわけです。追加パラメータを指定する定石 ... 引数を加えるとうまくいきます。なお、この問題の解決策は CRAN のキーワード検索で pairs と panel をキーワードに検索したらすぐに見つかりましたよ。-- QDU? 2004-10-07 (木) 23:53:06
    data(USJudgeRatings)
    pairs(USJudgeRatings[1:3], 
          panel=function(x,y, ...) {points(x, y, pch="."); abline(0,1,col=2)}, 
          xlim=c(5, 10), ylim=c(5, 10) )  
    pairspanel.png
  • ありがとうございます。panelに引数の追加を指定しつつも、xlim、ylim自体はpanelの外側に書くとは!  -- ぷーすか? 2004-10-12 (火) 10:04:16

eclipseの使い方

(2004-10-04 (月) 19:18:27)

eclipse (eclipse-SDK-3.0-linux-gtk.zip)をダウンロードし、3.0 用 Plug-inを入れました。eclipseを起動させたのですが、はて、使い方が分かりません。Window->Open Perspective->Other->StatET perspectiveを選択してみましたが、新規プロジェクトでRを選ぶこともできず、最初から躓いています。Eclipseを用いたRの使用について分かりやすく解説したサイトはないでしょうか?

  • 自己フォローです。File->New->R Script fileでエディタが現れました。実行はRun->External Tools->External Tools->Run Rのようです。簡単なものを書いて実行してみましたが、実行できたのか、実行結果はどこに現れるのか、訳が分からない状態です。 -- 質問者? 2004-10-04 (月) 19:28:53

エクセルファイルのデータをRに取り込む方法は?

初心者? (2004-09-27 (月) 11:51:31)

Rの初心者ですが、どうやってエクセルファイルのデータをRに取り込むんでしょうか?よろしくお願いいしマース〜〜

  • このページの下の方に,「CSVファイルのインポート」と「Rcmdrでグラフの変数選択欄」があります.これがお役に立つと思います. -- HOJO? 2004-09-27 (月) 12:33:24
  • 出来たら初心者の為のデータインポートとか初心者が躓きやすいノウハウを新しい頁にまとめていただけると良いのですが。:-) -- なかま 2004-09-27 (月) 13:09:51
  • R の公式マニュアルの 「R Data Import/Export」に他システムから(へ)のデータのやりとりがまとめられています。マニュアルの日本語訳がリンク集から得られます。 -- 2004-09-27 (月) 14:28:50
  • 最近の なんでも掲示板? の青木さんの関数も使えるかも知れませんね。 -- 2004-09-27 (月) 14:36:53
  • みなさんありがとうございました。何とかなりそうです。でも、ファイルが大きすぎかな?エクセルの2万5000ぐらいのレコード取り込み大丈夫でしょうか -- 初心者? 2004-09-27 (月) 15:19:34
  • なぜか繰り返される質問と言うのは現状のドキュメントでは不足しているか,機能が不十分か,質問者が何も考えてないかのどちらかではないかと。最後尾以外であれば改善したいのですが...レコード数だけでは判断出来ませんが,やってみるのが一番手っ取り早いでしょう。 -- なかま 2004-09-27 (月) 17:37:36
  • 「よろしくお願いいしマース〜〜」などと甘まったれた文句を平気で書く御仁は、どうせまず良く調べることなぞしないから、何をしても無理かも(独り言)。 -- 2004-09-27 (月) 19:31:48
  • 「2万5000くらいのレコード」ということは,2万5000個の数値ということではないので,最大2万5000×256個の数値ということになるので(本当か?レコードというものの概念がちゃんと理解されていればだが)そんなものを Excel で扱っているのも大変ご愁傷様です。
    質問者の常だが,「できるのでしょうか」なんて聞くより,やってみる方がはやい。やってみて,「できなかったのですが(できたけど時間がかかりすぎたのですが)なんとかなりませんか」というような質問ならまだ答えがいがあるような。
    初心者には優しくしましょうね(^_^)(発言タイトルはなおしておいたのです) -- 青木繁伸 2004-09-27 (月) 21:30:11
  • ごめんなさい!いろいろ試してみたんだけど、だめだったので、皆さんに聞こうかなと思ったんです。 -- 初心者? 2004-09-28 (火) 13:29:34
  • いつもこのようなメッセジーが出て、困っています。Warning messages: 1: [RODBC] ERROR: Could not SQLDriverConnect? ;2: ODBC connection failed in: odbcDriverConnect?(con) -- 初心者? 2004-09-28 (火) 13:38:35
  • 初心者の一人ですが,なかまさんのアドバイスに従って,「エクセルデータのcsvファイルによるインポート」というページを追加しました.初めてですので,フォローのほどお願いいたします. -- HOJO? 2004-10-04 (月) 14:32:14
  • なるほど,MsysのシェルにPATH追加してRとか打ってるようではわからん使いかたですね。(^_^; -- なかま 2004-10-04 (月) 15:02:48


ヒストグラムX軸目盛について。

tofu? (2004-09-21 (火) 11:11:26)

度数分布をヒストグラムにしたときののX軸の目盛の数字、目盛り幅を変えるのはどのようにすればいいのでしょうか?区間の境界値をいれたいのですが…教えてください。

  • axis 関数で,任意の場所に任意の文字列を書くことができます。 -- 青木繁伸 2004-09-21 (火) 12:11:21
  • せっかく グラフィックス参考実例集 があるのですから、それをまず一通り見て下さい。初等的な疑問への回答は大抵見つかるはずです。 -- 2004-09-21 (火) 14:30:26
  • すみません…いろいろ勉強不足なもので。アドバイスありがとうございました。 -- tofu? 2004-09-21 (火) 14:59:56
  • 実例集のページもずいぶん大きくなったのですね。見るのも大変は大変かも。右も左もわからない初心者ならなおさらかも。
    簡単に一例を示すと,以下のようになるか(縦軸は自分で考えてね)。 -- 青木繁伸 2004-09-21 (火) 15:47:49
    x <- rnorm(2000)
    hist(x, axes=F)
    axis(1, at=c(-1.96,0,1.96), labels=c("-1.96 SD", "Zero", "1.96 SD"))
    hist3.png
  • 丁寧な解説ありがとうございました。Rとともにこのサイトの見方も慣れていきたいと思います。 -- tofu? 2004-09-22 (水) 09:56:14

plotでドットだけを打ちたい

KabuTaro? (2004-09-20 (月) 04:02:13)

plot(c(1,2,3),c(4,5,6),type="p")
とすると(1,4),(2,5),(3,6)に丸印でプロットされますが
これを単にドットだけにするにはどうしたらよいでしょうか?
pchを弄くっても形が変わるだけでうまくいきません。
よろしくお願いいたします。

  • ドットだけの意味が分かりませんが,http://cran.r-project.org/doc/contrib/rdebuts_en.pdf のp.35にpchの細かい図が載ってます.plot(c(1,2,3),c(4,5,6),type="p",pch=16)やplot(c(1,2,3),c(4,5,6),type="p",pch=".")で如何でしょう? -- 2004-09-20 (月) 09:01:11
  • pch="."ですね。できました。ありがとうございました。これはピリオドですが、pch="・"だと何故かドットは表示されないようです。Rには座標平面に文字を打つ以外に通常の意味でのドットを打つことはできないんでしょうか? -- KabuTaro? 2004-09-21 (火) 01:07:38
  • plot(1,pch="A",cex=34) と plot(1,pch=".",cex=34) を比べてもらえると御理解いただけますでしょうか。 -- なかま 2004-09-21 (火) 08:14:28
  • 勝手注: 中間さんのコメントの通り pch="." は真正のドット(指定位置にちょうど置かれる1ピクセルの点)で、文字拡大率に影響されません。これはドットをつなげて曲線を描く等のため他の文字記号指定とはことなり、特別な仕様になっているようです。-- 2004-09-21 (火) 09:04:02
    pchtest.png
  • なるほど、pch="."には特別な意味があるんですね。てっきりピリオドだからピリオドの文字が表示されているかと思ってました。ありがとうございました。 -- KabuTaro? 2004-09-28 (火) 02:37:33

package(gregmisc)のheatmap.2の引数設定について

ちょろべぇ? (2004-09-14 (火) 18:57:31)

gregmiscパッケージのheatmap.2を用いてヒートマップを作図しております。
heatmap.2は様々なパラメータ設定が可能であり、自分としては、distfun,hclustfunを下記のように設定してHeatmapを書きたいと思っております。
(希望)
?distfun→距離関数だと思われるので、ユークリッド距離を計算したい
?hclustfun→average linkageで結合したい
?できればdivisiveではなくagglomerativeでクラスタしたい

ただ、

heatmap.2(data,distfun=dist(method="euclidian"),
          hclustfun=hclust(dist(method="euclidian"),method="average"),
          dendrogram="both")

のように設定しても、hclustfun部分でエラーが帰ってきてしまいます。(T_T)
かなり愚かな質問だとは思いますが、アドバイスありましたらどうぞよろしくお願いいたします。

  • こうした質問にはエラーメッセージをそえるのが礼儀作法です。まず "euclidian" => "euclidean"。 関数の引数に関数を指定する場合、引数関数のオプション引数はダイレクトには指定できない?いずれにしても dist 関数の既定手法は "euclidean" ですからわざわざ指定する必要はありません。hclustfun 引数は例えば次のようにすればどうなりますか? -- 2004-09-14 (火) 19:27:51
Hclust <- function(d) hclust(d, method = "average", members=NULL)
heatmap.2(data, hclustfun=Hclust, dendrogram="both")
# もしくは一度に済ませるには
heatmap.2(data, hclustfun=function(d) hclust(d, method = "average", members=NULL), 
          dendrogram="both")
  • こんな風なエラーメッセージが出てしまいますね --ちょろべぇ? 2004-09-14 (火) 22:43:16
    Error in as.matrix(x) : Argument "x" is missing, with no default 
  • euclidian?euclidean? どちらが正しいのでしょうか?
    ちょっと違いますがこんな感じでどうでしょう-- ○田五郎?2004-09-15 (水) 00:31:59
    my.dist <- function(x) dist(x, method="binary")
    my.hclust <- function(d) hclust(d, method="ward")
    hm <- heatmap(blablabla, distfun=my.dist, hclustfun=my.hclust)
  • 確認しないでコメントしました(最初のコメントのバグを修正しました)が、実際に確かめてみました。次の example(heatmap.2) 中の例を使ったコードは「一応」動きますね。dist 引数は dist() 関数を既定で使いますから省略してあります。 教訓は引数として関数を与える(そしてその関数の引数を特に指定する)場合はそれなりの指定方法があるということでしょう。-- 2004-09-15 (水) 07:16:09
    library(gregmisc)
    data(mtcars)
    x  <- as.matrix(mtcars)
    rc <- rainbow(nrow(x), start=0, end=.3)
    cc <- rainbow(ncol(x), start=0, end=.3)
    my.hclust <- function(d) hclust(d, method = "average")
    hv <- heatmap.2(x, col = cm.colors(256), scale="column",
                    hclustfun=my.hclust,
                    RowSideColors = rc, ColSideColors = cc, margin=c(5,10),
                    xlab = "specification variables", ylab= "Car Models",
                    main = "heatmap(<Mtcars data>, ..., scale = ?"column?")",
                    tracecol="green", density="density", dendrogram="both")
    # 形式 hclustfun=function(d) hclust(d, method = "average") でも可

    #ref(): The style ref(filename,pagename) is ambiguous and become obsolete. Please try ref(pagename/filename)

  • 諸先輩方、親切なアドバイスありがとうございました。ご教授痛み入ります。今後もRがんばっていこうとおもいます。m(_ _)m -- ちょろべぇ? 2004-09-15 (水) 09:35:50

エラーバー

Kei? (2004-09-11 (土) 16:54:52)

実験のデータをエラーバーで表現したいですが、R言語でどのふうすればいいのかよく分かりません。一応、自分で作ってみました。内容は以下のようです:

mybar <- function(X) {
  n <- nrow(X)
  df <- n-1
  X.mean <- apply(X,2,mean)
  X.var <- apply(X,2,var)
  d.error<-qt(0.975,df)*sqrt(X.var/n)
  d.error.x <- (0.2 + 1) * 1:5 - 1/2   # エラーバーを表示する位置(x座標)
  barplot(X.mean, names.arg=colnames(X), 
              ylim=c(0,ceiling(max(X.mean+(d.error)))),
              col=c("blue", "green", "red", "yellow"))
  arrows(d.error.x, X.mean, d.error.x, X.mean+ d.error, angle=90)
}
mybar(分析するデータ[,列:列])

軸の調整方についてちょっと質問したいですが、教えてください。

  • 自分で作れるというのは素晴らしいことですが、こうした問題はまず CRAN の検索機能でキーワード検索するというのが、定石です。 実際 Hmisc というパッケージにはそのものずばり errbar という関数があります。こうした関数を試した上で、なお不満があればそれらをお手本に改良する方が人生長く生きられます。もっとも青木先生のように、探すより自分で作った方が早いという猛者もいらっしゃいますが。 -- QDU? 2004-09-11 (土) 17:25:36
  • 縦軸をどのように調整したいのでしょうか?
    描画される範囲を,必要・十分に確保するというとことでしょうか。
    僭越ながら,私の書いた関数などが,もしかしたら参考になるかもしれません。
    要点は,描画する前にある軸(今回は縦軸)において,描画対象の範囲をまず探索しておく。座標軸だけを適切に描くために,以下のように plot 関数を呼ぶ(実際には,描画範囲を決め,軸だけを描く)。その後おもむろに,lines とか points で実際の描画を行う。
    なお,必要なら,描画範囲を決めるだけで,軸も実際には描かず,後で axis 関数でちゃんとした軸を描く,などなど。
    plot(c(x軸の最小値, 同最大値), c(y軸の最小値, 同最大値),
              type="n", # 描画しない
              xaxt="n", # x 軸を描かない(必要なら)
              yaxt="n", # y 軸を描かない(必要なら)
             その他のパラメータ)
      :
    points(...)
    lines(...)
    axis(...)
    その他諸々の描画関数
    健康のために使いすぎに注意。 -- 青木繁伸 2004-09-11 (土) 17:35:13

棒グラフの重ね合わせ

TETSU? (2004-09-10 (金) 11:51:50)

2つの棒グラフをbarplot関数を用いて描画したのですが、これらのグラフを比較するために、積み上げではなく重ね合わせで描画したいと思っています。色々と調べてみたのですがうまい方法が見つかりませんでした。どなたか描画法をご存知でしょうか?

  • barplot の beside という引数について,よく調べてください。
    それでもよくわからなかったら,グラフィックス参考実例集:棒グラフのページを見ましょう。
    このサイトで何かを調べようと思ったら,単語検索すればいいのですよぉ。 -- 青木繁伸 2004-09-10 (金) 14:06:32
  • ありがとうございます!ただ、besideをTrueにするとグラフは並置されてしまい重ならないのです。やはり重ねることは難しいのでしょうか?背景を透明にして、ふたつのグラフをパワーポイント上で強引に重ねるという荒技?も考えています。 -- TETSU? 2004-09-11 (土) 00:43:14
  • 「色々調べてみた」-> まず ?barplot でたくさんある引数を吟味してみましたか。普通思いつきそうなオプションは用意してあります。それでもなければ、「普通でない」ということですから、あまりお勧めではないということかも。barplot の space 引数をチェックして下さい。以下は example(barplot) の最後の例を少し変えてみたものです。あまり素晴らしいとは思えませんが。 -- QDU? 2004-09-11 (土) 08:35:45
    barplot(VADeaths, border = "dark blue", space=c(-0.5))
    bp.png
  • 青木先生、QDUさん、ありがとうございます。全ては私の日本語力不足のせいです。イメージとしてはアグレスティ著「カテゴリカルデータ解析入門」(渡辺、菅波ほか訳)の7ページにあるような図が書きたいのですが、、、皆様から頂いた知恵をもとに頑張ってみたいと思います。それにしても素晴らしいWEBページですね! -- TETSU? 2004-09-11 (土) 10:45:59
  • polygon 関数を使って,自分でヒストグラムを描く関数を書けば,何でもできます。 -- 青木繁伸 2004-09-11 (土) 11:25:51

層別の相関値

HOJO? (2004-09-10 (金) 11:39:45)

 欠損値を含む二つの変数について,相関値をcor(x,y,use=pairwise complete observations)で計算しました.そうしたところ,全体では0.919だったのに,上位30件では,0.910,次の70件では,0.395, 残りの100件のデータでは,0.045となってしまいました.このようなことはどうして起こるのでしょうか.層別の方法に何か問題があるのでしょうか.

  • 上位30件というのはどのように決めたのでしょうか。
    それはともかく,データの散布図を描けば問題解決の手がかりがつかめるでしょう。
    下の図の右上30個ほどの点は相関係数0.910,次の70個ほどは0.395,残りは0.045です。全体は0.877程度にしかなっていませんけど。
    plot.png
    R はグラフを描くのがとっても簡単なので,数値解析する前に図を描いてみる習慣をつけてはいかがでしょう。 -- 青木繁伸 2004-09-10 (金) 13:05:20
  • 青木先生いつも有り難うございます。散布図の特徴は,(1)原点付近に100件程度のデータが散らばっている。(2)座標(10000,600)付近に他と離れたデータが1件ある。(3)X座標で1000~4000,Y座標で50~300あたりに30件程度のデータが散らばっている。 一般に全体の相関係数は,層別した相関係数より大きくはならないと思うのですが,どうなのでしょうか。 -- HOJO? 2004-09-10 (金) 15:18:14
  • >一般に全体の相関係数は,層別した相関係数より大きくはならないと思うのですが
    そんなことありませんよ。例えば,原点付近と(10,10)を中心とした各変数の標準偏差1で,相関0の二次元正規分布データ10個ずつくらあるとき,20個のデータ点全体の相関を計算することを想像すればいいでしょう。
    > gendat2 # 任意の相関係数を持つデータを生成する関数
    function(nc, r)	
    {
    	z <- matrix(rnorm(2*nc), 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))
    }
    > d1 <- gendat2(10,0)
    > d2 <- gendat2(10,0)+10
    > cor(rbind(d1,d2))
              [,1]      [,2]
    [1,] 1.0000000 0.9615385 # プールすると相関は0.962
    [2,] 0.9615385 1.0000000
    > cor(d1)
                  [,1]          [,2]
    [1,]  1.000000e+00 -1.010639e-16 # プールする前は無相関
    [2,] -1.010639e-16  1.000000e+00
    > cor(d2)
                [,1]        [,2]
    [1,] 1.00000e+00 3.27047e-16 # プールする前は無相関
    [2,] 3.27047e-16 1.00000e+00
    極端な場合,層別では負の相関,全体では正の相関ということだってあるんですよぉ!
    >(10000,600)付近に他と離れたデータが1件
    というのは,すさまじいような気がしますが。その1つのデータを除いて同じように相関係数を計算してみませんでしたか?
    ともかく,HOJOさんが今扱っているデータではそのように計算されているわけではあるが,それが変だと思うなら,その変な解析結果の出てくる原因を突き止めないといけないですね。聞いた限りにおいては,「そのデータではそうなって当たり前かも」という意見しかないですが。
    散布図をここにペーストしてくれると,後学のためにもなるかとも思いますが。このスレッド「広い意味で R に関係すること」に入りますかねぇ? -- 青木繁伸 2004-09-10 (金) 15:32:11
  • >そんなことありませんよ。例えば,原点付近と(10,10)を -- HOJO? 2004-09-11 (土) 06:43:08
  • >そんなことありませんよ。例えば,原点付近と(10,10)を...これは考えつきませんでした.このデータは,単に企業の売上などを順に上位から200社程度ならべただけに過ぎないのですが,最上位の会社が巨大企業なので,他と飛び離れた値となってしまいます.  散布図を作成してみました.初めて貼り付けるのでうまくいくとよいのですが. #ref(D:?My Documents?40joho?scatterplot1.bmp) -- HOJO? 2004-09-11 (土) 06:50:38
  • #ref(D:?My Documents?40joho?scatterplot1.bmp) -- 2004-09-11 (土) 06:55:08
  • この wiki に画像を張り付けるには、jpg か png 形式で、ページ末にあるボタンの5つ目を押して、画像のアップロードをしなければ駄目です。もしくは自分のホームページに置いた画像の URL へのリンクを張ります。 -- 2004-09-11 (土) 09:07:11
  • すこし極端な例を。全体では相関 0, 層別すると相関が 1, -1 になります。理由は添付画像を見れば自明ですね。-- QDU? 2004-09-11 (土) 09:16:57
    cor.png
  • #ref(scplot1.jpeg) -- HOJO? 2004-09-11 (土) 20:48:50
    scplot1.jpeg
  • 「編修」を使ってようやくできました。 アドバイスありがとうございました。 -- HOJO? 2004-09-11 (土) 21:06:09
  • この図から見るにやはり「外れ値」とみなすべき例外的データがあるようですね。それらを取り除かなければ、相関係数を計算してもほとんど無意味でしょう。もしくはこうしたスケールが極端に異なるデータが混在する時は(x and/or y 値の)対数を取るというのが一つの定石です。こうすることにより、原点付近に密集したデータをより散らばせることもできます。ただしそうしたからといって相関係数が高くなるかどうかはやってみないとわかりません。まず散布図を書いて直線状かどうか見極めるのが安全です。 -- 2004-09-12 (日) 11:14:11
  • 青木先生,QDUさん他の皆様方アドバイス有り難うございました.そうすると以下のようなことがいえるのでしょうか.(1)層別相関係数と全体相関係数の数値について,特段にいえることはない.(2)外れ値を含むデータについて,相関係数を単純に求めると,(その外れ値を含まないデータ全体の)傾向を見誤る可能性がある.今後ともよろしくお願いいたします. -- HOJO? 2004-09-13 (月) 09:21:32

単純主効果

(2004-09-10 (金) 02:07:41)

2元配置の分散分析(完全無作為)をする場合に
交互作用がある場合,
1.「単純主効果の検定」,
2.「単純主効果の多重比較」
をしたいのですが,Rではどのよなプログラムを書けば良いのでしょうか.

参考までに, SASでの単純主効果についての情報と,そのデータを以下にリンクします.
例:http://koko15.hus.osaka-u.ac.jp/~kano/lecture/faq/interaction.doc
データ:http://mat.isc.chubu.ac.jp/R/tech.html#CRFpq

ご存知の方,教えて頂ければ幸いです

  • 単純主効果という言葉の意味が一寸わからないのですが、R の aov 関数の出力では得られないものなのでしょうか。R の多重比較用の関数については Rの基本パッケージ中の古典的検定関数一覧 を見るとそれらしい関数がありますが、これでは役にたちませんか。なお、私は分散分析も多重比較もほとんど使ったことのないものですから、的外れのコメントかも知れません。 -- 2004-09-11 (土) 17:34:39
  • aov関数では得られないようです.単純主効果の説明は,http://www.aichi-gakuin.ac.jp/~chino/anova/chapter1/sec1-4-3.html が詳しいですがRで出力できる関数があるかどうかがわかりません. -- 2004-09-11 (土) 20:30:13
  • CRAN の検索では simple main effect + anova では一件も引っかからないところを見ると、かなりマイナーな手法なのでしょうか。幾つかの解説を見たところ、要するにある要因の水準を固定して分散分析を行なう意味だと感じられますが?それならその水準を固定した副データフレームを取り出し分散分析すればいいのかな? -- 2004-09-12 (日) 11:23:06

Rcmdrでグラフの変数選択欄

HOJO? (2004-08-31 (火) 16:06:53)

たびたびですみません.RcmdrでIndexPlot?などを用いてグラフを作画しようとしても,変数選択欄に必要な変数名称が全ては現れず,作画することができません.またこの内容は作画すべきグラフの種類によっても異なりますがこれはどうしてでしょうか.なおPCは,W2Kを利用しています.

  • 変数選択欄に出るのは数値として認識出来た列のみになるかと。従って必要な変数が現れないのはデータを誤認識(因子としてデータフレームに取り込まれた)してしまったかなんかだと思います。読み込んだデータが意図した状態か確認してみて下さい。 -- なかま 2004-08-31 (火) 20:23:12
  • なかまさん,何度もお手を煩わして申し訳ありません.因子かどうか区別する関数は知らないのですが,Rコマンダーの「要約」でみると,変数欄に表示されない変数については,すべて出現頻度らしきものが表示されていますので,なかまさんのご指摘のとおりだと思います.ただこれを単なる数値に戻すのはどうすればよいのでしょう. -- HOJO? 2004-09-01 (水) 10:36:06
  • 変数毎にas.vectorを行うと,文字データになってしまいます。次いでas.numericを行うと,数値データにはなりますが,最初の30個がNAになってしまいます。 -- HOJO? 2004-09-01 (水) 11:31:31
  • こうした質問はもとのデータ(長過ぎるなら問題になりそうな箇所を含む)を一緒に紹介しないと、回答者に憶測をさせることになり、時間のむだでもあり、失礼でもあります。 -- 2004-09-01 (水) 11:35:21
  • どうしても意図した通によみこみたいならscanを使ってください。read.tableへの理解が深まると思います。 -- なかま 2004-09-01 (水) 12:37:46
  • 数値そのものに,あまり意味があるとは思えないので掲載しなかったのですが。数字は企業の売上(単位:億円)を上から順に並べただけです。最初の16個は以下のようになります.[1] 8,467.05 7,570.33 3,784.72 3,473.94 3,443.77 3,418.74 3,094.17 2,661.70 -- HOJO? 2004-09-01 (水) 12:53:30
  • scanを使用して解決できました。有り難うございました。 -- HOJO? 2004-09-01 (水) 14:05:01
  • もとのデータがただの数値で、単純に因子化されているならば、as.numeric() でもとの数値に戻るはずです。NA が返されるということは、そうはなっていないという意味でしょうか。因子のままの形で出力して、なにかおかしなものが混ざっていないか点検されたらいかがでしょうか。とくに水準がどう表されているか。 -- 2004-09-01 (水) 15:21:41
  • 詳細に記載すると,以下の処理をしてからscanを施し解決しています。 -- HOJO? 2004-09-02 (木) 09:38:17
  • 詳細に記載すると,csvファイルを作成する際に事前に以下の処理をしています。(1)3桁毎のカンマ表示を除く,(2)%表示ではなく,小数点表示とする。(3)空白欄にはすべてNAという文字を代入する。(2)は,%表示の変数がすべて因子と解釈されたからです.(1)については,NAに変換された最初の30個はすべて,カンマを含む数値(つまり1,000以上)でしたから.(3)は,そうしないと正しく全データを読みとることができなかったからです.このデータに対して,scanをかけて解決しました.ちなみに,このcsvファイルをread.tableで読み込むと全てnumあるいはintで読みとってくれます.結局,エクセルでcsvファイルを作成する際には,上記の(1)(2)(3)等を施してからでないと,read.tableが正確に働かないように思えます.皆様のおかげで,解決でき,大変有り難うございました. -- HOJO? 2004-09-02 (木) 09:51:20
  • (1)は,カンマ区切りの数値は引用符で囲まれるから(そうしないと,数値間のカンマと桁区切りのカンマが区別できない。もっとも,タブ区切りだろうが,空白区切りだろうが同じ扱いになってしまうが)。
    (2)は文字として % が実際にファイルに書き出されているため,R は文字列と認識するだろう。
    (3)空セルは文字通り,なんにも出力されないので,ケースあたりの数値の個数が不足するため。
    全て,Excel に罪はないとも言える(制御のためのオプションを用意してくれればいいが)。
    いちいち手作業で(1),(2),(3)の前処理を行うのも面倒なので(間違いも生じやすいので)前処理スクリプトか,CSV読みとりの専用プログラムを作るのがよさそう。 -- 青木繁伸 2004-09-02 (木) 10:11:11
  • 解説有り難うございました.相変わらず,やみなべ(R)をつついていますが,皆様のアドバイスでなんとか前に進むことができます. -- HOJO? 2004-09-02 (木) 16:06:50

CSVファイルのインポート

HOJO? (2004-08-31 (火) 11:26:52)

エクセルで作成した200行,20列程度の数値データのみのCSVファイルをRcmdrで読み込んだ時に次のエラーメッセージが出ます。RはV.1.9.0です.
1. Error: unable to open connection
2. data is not a data frame and cannot be attached
これはどのような点に気をつければよいのでしょうか。

  • 多分ですが,日本語化してないRで日本語ファイル名のファイルを読もうとしたのでは?tcltk上でUTF8で認識されてそのままRにUTF8ファイル名が渡って読めないと言う話ではないかと。推測のなになのであれですが、日本語ファイル名を使わないと幸せかと。 -- なかま 2004-08-31 (火) 12:47:09
  • なかまさん,いつも有り難うございます.ご指摘のとおりでした.フォルダ名が日本語でしたので,修正したらうまく行きました. -- HOJO? 2004-08-31 (火) 13:41:51

非心パラメータ

(2004-08-24 (火) 19:12:06)

正確な用語がわからないのですが
SASではTNONCT関数で算出できる,
「t分布の非心パラメータ(引数は,x, df, prob)を求める関数」は,
Rにないでしょうか?
http://support.sas.com/91doc/getDoc/lrdict.hlp/a000245960.htm

例えば,
x <- 3.1; df<- 58; prob <- .975;
?(x,df,prob)

1.04844
と求めたいのです.

  • そうした関数はないですが、要するに pt(3.1, 58, ncp) = 0.975 を ncp に付いて解けば良いわけですから、次のような関数で求められます。uniroot() 関数は一変数関数の零点を求める関数です。interval 変数はその中で解を求める範囲で、問題に応じて変える必要があるかも知れませんが、help(dt) によれば制限 ncp <= 37.62 と書いてあります。-- QDU? 2004-08-24 (火) 20:37:44
    > TNONCT <- function(x, df, prob, interval = c(-30, 30)) {
                         temp <- function(ncp) pt(x, df, ncp) - prob
                         uniroot(temp, interval)[[1]] }
    > TNONCT(3.1, 58, 0.975)
    [1] 1.048444
  • 同じようなものになりましたが,参考までに。 -- 青木繁伸 2004-08-24 (火) 22:53:56
    > TNONCT <- function(x, df, prob, interval=c(0, 37.62))
    + {
    +     uniroot(function(ncp) pt(x, df, ncp)-prob, interval)$root
    + }
    > TNONCT(3.1, 58, 0.975)
    [1] 1.048444
  • 非心分布など普段意識して使ったことがないので、いまひとつ良く分からないのですが、ncp が負になるケースは考える必要は無いのでしょうか?少なくとも pt(), qt() 関数は負の ncp パラメータを受け付けるようですね。 -- QDU? 2004-08-25 (水) 00:59:31
  • ありがとうございます.おかげさまで,目的の計算ができました. -- 2004-08-25 (水) 12:10:52
  • 私も,非心度を求めるようなことはやったことがありませんでした。非心t分布は必要になることもありますが。
    > TNONCT(3.1, 58, 0.975)
    [1] 1.048444
    > TNONCT(-3.1, 58, 0.025)
    [1] -1.048444
    という関係があるようです。interval=c(-37.62, 37.62) としておいた方が無難ですね。
    質問者が示した URL での解説から,SAS でも同じような関数を定義して使っていそうですね。
    google で site:r-project.org ncp pt で,もろ回答がありましたね(^_^;)
    R Help 2001 これでは,c(-1e7,1e7) としていますが,QDU さんの言うように,ncp <= 37.62 ですから,あんまり大きい値を指定するとまずいかな?と思って,,,やってみたら別に問題ないみたいでした。 -- 青木繁伸 2004-08-25 (水) 12:02:16

行列のグラフィカル表示image()の色に任意の値を割り当てるには

川村健介 (2004-08-18 (水) 13:14:01)

お世話になります。
行列をグラフィカルに表示するimage()関数の色調を任意に割り当て,さらにその色を凡例(棒)で表示するにはどうすればよいのでしょうか。
例えば配列の値が0〜0.88の間にあった場合,0.0〜0.5に白を割り当て,0.5〜0.9の間を0.1間隔でGreyの色調を与えたいのですが。

  • 色見本付きのイメージ図は example(filled.contour) 参照。.色の指定は最近のこのページで話題になりましたね。col=gray(c(0,0,0,0,0,0.6,0.7,0.8,0.9)) でどうでしょう。もしかすると col=gray(c(1,1,1,1,1,0.4,0.3,0.2,0.1)) かな?-- 2004-08-18 (水) 13:53:07
  • ありがとうございます。col=grey(c(1,1,1,1,1,0.4,0.3,0.2,0.1))で,目的の図を得ることができました。カラーテーブルを勉強しなくてはなりませんね。ところで,「S-PLUSによる統計解析」伊藤・大津・戸瀬・中東 訳(シュプリンガー・フェアラーク東京)のp66,表3.1の中にあるimageの説明で「levelplotの方を推奨」とあります。これはどういった理由なのでしょうか? -- 川村健介 2004-08-19 (木) 07:19:11
  • 数値と濃淡が正しく対応しているか慎重に確認して下さい。なお、当方は貧乏なので S-PLUS のことは知りません。 -- 2004-08-19 (木) 08:18:23

plot()で軸の数値をなくしたい

(2004-08-18 (水) 00:22:26)

10ほど並べた散布図で、最初の1つだけ軸の数値を表示し、残りの9個は軸の数値を省略したいとのですが、ラベルはxlab="",ylab=""でOKですが、数値の方はどうすればよいか分かりません。よろしくお願いします。

par(mfrow=c(2,5))
plot(sin)
plot(sin, xaxt = "n", yaxt = "n", xlab="", ylab="")
plot(sin, xaxt = "n", yaxt = "n", xlab="", ylab="")
# あと 7 回繰り返し

上のような感じで如何でしょうか?ご参考になれば幸いです.-- 舟尾? 2004-08-18 (水) 00:40:00

  • ありがとうございます。目的の図を得ることができました。 -- 2004-08-18 (水) 11:21:11

Rで扱える統計手法

石川誠? (2004-08-16 (月) 01:26:14)

お世話になります。
Rで扱える統計手法を知りたいのですが、何処かに日本語で一覧のようなものはないでしょうか?

  • まさにこのサイトでしょう。stats(R 統計)パッケージ中のオブジェクト一覧?が,最初に見るべきものでしょう。
    ほかにどのようなページがあるかは,ページの一覧を見ればよいし,単語検索によって必要なページを検索することもできますよ。 -- 青木繁伸 2004-08-16 (月) 10:10:57
  • R 本体のオブジェクト数は千以上、アドオンパッケージは少なくとも400あまり、といえば、そんなものは無いことはことはお分かりになるでしょう。具体的な手法をあげれば、回答がえられるかも知れません。library(help=stats) で調べて下さい。英語になれるのが一番です。 -- 2004-08-16 (月) 14:15:00
  • R に含まれる統計手法全部を知りたいというわけでもないでしょうし,全部知らないと R で統計学を使えないというわけでもないのですから,まずは,存在する資料を見るのがよろしいでしょう。また,舟尾さんの労作(pdfファイル)も参考になると思いますよ。 -- 青木繁伸 2004-08-16 (月) 15:33:43
  • Rで扱える統計手法でしたら,まず見なければいけないのは Rによる統計処理 でしょう.一般的な統計手法はほぼ全て網羅されていますし,詳細な解説が付いておりますので,非っ常〜〜に役に立ちますよ. 2004-08-16 (月) 20:33:43
  • え〜〜,本人から申し上げますが,あそこはまんまり最初っから見ない方がいいです(というか,最後まで見なくて済むなら,その方が幸せでしょう) -- 青木繁伸 2004-08-16 (月) 21:26:40
  • またまた〜、ご謙遜されてますねぇ〜。「R」+「統計」=「Rによる統計処理」という位、充実しているサイトですよ>石川さん。 -- 2004-08-17 (火) 00:27:35

行列の値上位5つの要素の位置を知りたい

川村健介 (2004-08-14 (土) 09:05:51)

初めて投稿させていただきます。
630*630の配列の値の中から,上位5位までの要素の位置(x,y)をしるいい方法はありませんか?
which.maxだと最も高い値の位置,例えば39726と帰ってきます。
できれば,上位から[326,52],[555,156]...といった感じで位置を知りたいのですが。

  • 欠損値 NA がなければもう少し簡単になるのですが。(PS. もし値にタイがあるとランクは整数になりませんから、うまく行きませんね。) -- QDU? 2004-08-14 (土) 10:01:10
    > y        [,1]      [,2]      [,3]
    [1,]         NA 0.7122784 0.7761505
    [2,] 0.47935820 0.4700606 0.4547836
    [3,] 0.02593524 0.8849352 0.3914027
    > yy = matrix(length(y)+1-rank(y, na.last=F), dim(y), byrow=FALSE)
    > y  [,1] [,2] [,3]
    [1,]    9    3    2
    [2,]    4    5    6
    [3,]    8    1    7
    > lapply(1:5, function(i) which(yy==i, arr.ind=T))
    [[1]]     
             row col
      [1,]   3   2
    [[2]] 
           row col
      [1,]   1   3
    [[3]]
           row col
      [1,]   1   2
    [[4]]     
           row col
      [1,]   2   1
    [[5]]     
          row col
      [1,]   2   2
  • すばやいレスありがとうございます。うまくいきました。ちなみに欠損値NAがなければもう少し簡単になるということですが,rank(y, na.last=T)ということですか?もっと簡単にする別の方法がありましたら,参考までに教えていただけないでしょうか。 -- 川村健介 2004-08-14 (土) 11:16:54
  • ここはやはりベクトルの n 番目に大きな値を求める関数を定義すべきでしょう。タイや欠損値がある場合も当然考えておかないと思わぬ落し穴になるかも。 (PS. 蒙古草原での実習なんて素敵ですね。ところでこの後期に岐阜大で「工学のためのデータサイエンス入門」の著者の一人神保先生による講義が予定されています。この本をテキストにつかうそうです。もっとも学部レベルでしょうが。)-- QDU? 2004-08-14 (土) 11:43:26
    # ベクトル(行列、配列も可)の n 番目に大きな値を求める関数
    # タイ や NA があっても使える(はず?) 
    > nthmax          
    function(x, n) {
                if ( !(n %in% 1:length(x)) )
                  stop( "n is not in 1:length(x)" )
                if (n == 1) return(max(x, na.rm=TRUE))
                  else
                    for ( i in 1:(n-1) ) {
                      x <- x[x != max(x, na.rm=TRUE)]
                      if ( length(x) == 0 )
                        stop("no 'n'-th max element in x")
                    }
                return(max(x, na.rm=TRUE))
              }
    
     > x = matrix(runif(9), 3,3)
     > for (i in 1:length(x)) print( which(x==nthmax(x,i), arr.ind=T) )
         row col
    [1,]   1   1      
         row col
    [1,]   3   1
         row col
    (途中略)
         row col
    [1,]   1   2
    
    > x[1,2] <- x[1,1] # タイを入れる
    > for (i in 1:length(x)) print( which(x==nthmax(x,i), arr.ind=T) )
         row col
    [1,]   1   1       # <- 最大値が二つある
    [2,]   1   2
         row col
    (途中略)
         row col
    [1,]   1   3
    Error in nthmax(x, i) : no 'n'-th max element in x
    
    > x[1,2] <- NA # 欠損値を入れる
    > for (i in 1:length(x)) print( which(x==nthmax(x,i), arr.ind=T) )
         row col
    [1,]   1   1
         row col
    (途中略)
         row col
    [1,]   1   3
    numeric(0)         # 最後がエレガントでない(欠損値のせい)
    Warning message:
    no finite arguments to max; returning -Inf
  • 初等的なやりかたは見つけた最大値を -Inf に置き換えて繰り返すことでしょうか。たとえば -- 2004-08-14 (土) 12:16:31
    > x = matrix(runif(9), 3,3)
    > which(x==max(x), arr.ind=T)
         row col
    [1,]   2   3
    > x[which(x==max(x), arr.ind=T)] <- -Inf
    > which(x==max(x), arr.ind=T)
         row col
    [1,]   1   3
    > x[which(x==max(x), arr.ind=T)] <- -Inf
    > which(x==max(x), arr.ind=T)
         row col
    [1,]   3   2
    -----
    > x[ print(which(x==max(x), arr.ind=T)) ] <- -Inf # 一行で済ませるにはこうする
  • わかりやすい解説ありがとうございます。岐阜大で講義があるのですね。それはいい情報をいただきました。ありがとうございます。 -- 川村健介 2004-08-14 (土) 15:39:09
  • そんなに実行時間の短いやり方が必要なんですか。スローライフを楽しみましょう(^_^;)
    作りつけの速い関数を求めるのもいいが,それを問い合わせる時間がもったいない。どんなにのろいアルゴリズムでも,数時間かかるようなことはない。それを,何万回も繰り返すのなら別だが。だけど,何万回繰り返すとしても,それを人間時間で処理することを考えると,そんなに時間の節約にもなっていないのではないかと思うわけです。
    作りつけの関数でできないことも,いろいろ付け加えられますしね。
    こういう考え方を許容できない人は無視すればいいだけです(^_^;)
    > x <- matrix(runif(630*630),630) 
    > 
    > max.at <- function (x, top=10)
    + {
    + 	result <- matrix(0, top, 3)
    + 	for (loop in 1:top) {
    + 		mx <- x[1,1]
    + 		for(i in 1:630) {
    + 			for  (j in 1:630) {
    + 				if (x[i,j] > mx) {
    + 					mx <- x[i,j]
    + 					mxi <- i
    + 					mxj <- j
    + 				}
    + 			}
    + 		}
    + 		result[loop,] <- c(mxi, mxj, mx)
    + 		x[mxi, mxj] <- -Inf
    + 	}
    + 	result
    + }
    > 
    > system.time(max.at(x))
    [ 1] 78.30  4.31 89.11  0.00  0.00
    >  max.at(x)
          [,1] [,2]      [,3]
     [1,]  466   95 0.9999976
     [2,]  239   31 0.9999970
     [3,]  205  620 0.9999952
     [4,]  595  183 0.9999935
     [5,]  500  609 0.9999923
     [6,]   37  556 0.9999901
     [7,]  292    3 0.9999889
     [8,]  515  278 0.9999818
     [9,]  539  237 0.9999802
    [10,]  483  331 0.9999786
    速くしたいなら,C ででもプログラムすればいいだけです。プログラムするにはアルゴリズムが必要ですが,これはアルゴリズムというような大げさなものでもないし。-- 青木繁伸 2004-08-14 (土) 19:55:23
  • こんなのはいかが?
    max.at2 <- function (x, top=10)
    {
    	n <- nrow(x)*ncol(x)
    	col <- col(x)
    	row <- row(x)
    	o <- order(x)
    	cbind(row[o], col[o], x[o])[n:(n-top+1),]
    } 
    > system.time(max.at2(x))
    [1] 2.52 0.18 2.80 0.00 0.00 # すんばらしく,速いわ!
    > max.at3(x)
          [,1] [,2]     [,3]
     [1,]    2    3       NA # NA はこの位置に入る
     [2,]   18  178 99.99992
     [3,]  549  252 99.99988
     [4,]  208  115 99.99891
     [5,]  182  140 99.99866
     [6,]  308  597 99.99823
     [7,]  207  156 99.99820
     [8,]  616  112 99.99811
     [9,]  499  249 99.99801
    [10,]  355   79 99.99794
    本当は,以下のようにしたいのだが,decleasing=TRUE が効かない??
    max.at3 <- function (x, top=10)
    {
    	n <- nrow(x)*ncol(x)
    	col <- col(x)
    	row <- row(x)
    	o <- order(x, decreasing = TRUE)
    	cbind(row[o], col[o], x[o])[1:top,]
    }
    いずれにせよ,
    max.at4 <- function (x)
     {
     	o <- order(x)
     	cbind(row(x)[o], col(x)[o], x[o])
     } 
    とすれば,小さい順の結果は全部得られるので,必要なところだけ(いくらでも)使えばよい。NA も考える必要はない?(考えて処理すればよい) -- 青木繁伸 2004-08-14 (土) 21:21:49
  • 青木先生,ありがとうございました。大変,勉強になりました。なるほど。上位5つを取り出すため,小さい順に並び替えて,必要なところだけを抜き出すように考えれば,NAを考える必要もなく上位からも下位からも利用できますね。 -- 川村健介 2004-08-16 (月) 15:19:37

statsライブラリ内のheatmap関数について

○田五郎? (2004-08-10 (火) 01:45:09)

こんばんわ。○田です。
この度はstasライブラリ内のheatmap関数について相談させていただきたく。
statsライブラリをロードして、heatmapがかけることはわかったのですが、
?heatmap内のデモスクリプトを書くなどしてみましたが、
heatmap関数の引数であるcol引数の種類が良くわかりません。。。
自分としては正の値が赤系、マイナスの値を青系で表示したいのですが
デモのcol = cm.colors(256)はレインボー系、、、
col = topo.colors(16)は緑と薄いこげ茶系と、自分の望んでいる色が出せないんですよね。。ググってもあまり良い情報も無く、もしみなさんどこか素敵なソース(情報源)等お知りでしたら、ご教授ください。
どうぞよろしくお願いしますm(_ _)m。

  • こんなのありますよー http://gdds.pharm.kyoto-u.ac.jp/lecture/n006/dry_report.ppt -- ちょろちゃん? 2004-08-10 (火) 01:47:47
  • アドオンパッケージ RColorBrewer? には様々な配色のパレットを選べる brewer.pal() という関数がありますので試してみたら。既定のカラーマップ
    rainbow(n, s = 1, v = 1, start = 0, end = max(1,n - 1)/n, gamma = 1)    
    heat.colors(n)  
    terrain.colors(n)     
    topo.colors(n)     
    cm.colors(n)
    と同じように使えます。一般的には個人で配色を決めるよりはこうした慎重に選択された配色を使う方が賢いです。-- 2004-08-10 (火) 08:05:42
       
    ## help(brewer.pal) より  
    mypalette<-brewer.pal(7,"Greens")    
    image(1:7,1,as.matrix(1:7),col=mypalette,xlab="Greens (sequential)",
          ylab="",xaxt="n",yaxt="n",bty="n")
  • グラフィックス参考実例集:カラーパレット を見て下さい。topo.colors() による色調はは海は青、陸は土色となるはずですが。 -- 2004-08-10 (火) 20:04:58
  • 書きたいことは,ちゃんと書きましょう。主張したいことがよくわからないのですが。 -- 青木繁伸 2004-08-10 (火) 21:36:26
  • 元記事投稿者が topo.colors では青い色がでないといっているように読んだので、青も出ますよと注意しました。そうではなく、海色部分はOKだが、草・土色部分が気に入らないというのであれば、次のような仕方もあるかも知れません。カラーマップは単に色名をある方針で並べただけの文字ベクトルですから、一部分をツギハギすることもできるわけです。個人的にはあまり好ましいとは思いませんが、しかし趣味の問題でしょう。-- 2004-08-10 (火) 23:14:57
    data(volcano)
    x <- 10*(1:nrow(volcano))
    y <- 10*(1:ncol(volcano))
    image(x, y, volcano, col = c(topo.colors(12)[1:4],
               heat.colors(8)[c(1,3,5,8)]), axes = FALSE)
    topoheat.png

mvtnorm内で帰ってくる値が小数点4桁以下がその都度異なる?仕様?バグ?

ちょろちゃん? (2004-08-08 (日) 15:34:16)

お世話になっております。
Rのmultcompライブラリのsimtest、method=dunnettを
利用しております。
同じデータセットを入力した際に、simtestのダネットの結果中、
T値は常に一定の値が帰ってくるのですが、P値の小数点4桁よりも小さな値は
計算するたびに異なっております。
どこで値が異なってしまったか追跡したところ、mulcompから呼び出されている
mvtnorm内でP値を計算するときにどうやら異なった値が帰ってきているようでした。
これは仕様なのでしょうか、それともmvtnormのバグなのでしょうか。
もしご存知の方がいらっしゃいましたらアドバイスいただければ幸いです。
どうぞよろしくお願いいたします。

  • どこかで乱数でも使っていない限り、仕様にせよバグにせよ、同じ計算の結果が異なるというのは妙ですね。こうした質問では、結果が再現できるような実際のデータか、人口的なデータ(を作る R 手続き)でも一緒に紹介しないと、一般的には答えようがないですね。 -- 2004-08-08 (日) 17:37:18
  • mvtnorm って,multcomp から呼ばれるんですか。よくわかりません。
    R の関数には,uniroot だっけを呼び出して解を求めるなどという関数もあり,そのような場合には,毎回似て非なる解が得られて驚くと言うこともあるようです。
    あなたが利用した関数の定義を参照すれば,毎回異なった解が得られるからくりもわかるかも。
    P値などというものは,0.05や0.01と比べて大きいのか小さいのかがわかればよいのだという,極端な意見もありますので(誰がそう言っているのか?),小数点以下4桁目が違うとか違わないとかはどうでもいいじゃないかという意見もあるかもしれない(誰がそう言うのか?)
    いろいろな統計量に関するP値も,近似値なので(統計ソフト間で比較すればよくわかる),小数点以下何桁目まで正しいかというのは,よくわかっていない(保証されていない?)ことなのではあります(どうでもいいことでもあります)。
    いずれにしろ,近似値だと言うことで。そんなに目くじらたてる必要もないでしょう。 -- 青木繁伸 2004-08-08 (日) 19:58:18
  • 自分がテストに使ったデータは青木先生のページに記載してあったデータです。http://aoki2.si.gunma-u.ac.jp/R/dunnett.html。simtestのdunnettや青木先生のページにあった関数を用いても、やはり最終的なP値が変わってしまったんですよね。。。青木先生の関数の途中結果を出力していたところ、1-pmvt(lower=-x, upper=x, delta=rep(0, a-1), df=df, corr=corr, abseps=0.0001)のところで毎回計算が異なっていることがわかりました。ただ青木先生のおっしゃるように、統計量に関するP値ももともと近似値ですので、たしかに小数点4桁以下の値に目くじらをたてることはないかもしれませんね。(^_^;) -- ちょろちゃん? 2004-08-08 (日) 20:10:55
  • 灯台もと暗しでしたか。mvt 関数が疑惑の張本人ですね。主要部分はFORTRANで記述された関数のようですが,解を得るために内部で乱数を呼んでいるような雰囲気です。
    ? pmvt を読みますと(めんどくさいので斜め読み。従って精度は不詳)
    Randomized quasi-Monte Carlo methods are used for
    the  computations.
    とか
    Value:
    error: estimated absolute error
    とか読めますが。
    また,pmvt の引数には,abseps = 0.001, releps = 0 というのがあり,許容誤差はユーザが制御できるようです。いずれにせよ,まずは ? pmvt のようです。 -- 青木繁伸 2004-08-08 (日) 20:41:18
  • なるほど。青木先生ありがとうございますm(_ _)m。感謝感激雨霰。私も?pvmtを見てみたいと思います。これからもRがんばります(^^) -- ちょろちゃん? 2004-08-08 (日) 21:27:21
  • 青木さんのおっしゃる通り、乱数が絡んでいるようですね。試しに実行前に乱数シードを毎回設定(例えば set.seed(222))とすると同じ結果が得られるようですね。こんなこともあるのですね、勉強になりました。 -- 2004-08-08 (日) 22:49:08

繰り返しのt検定の結果を表示させるには

はまちゃん? (2004-08-03 (火) 18:12:58)

お世話になっております。

a<-matrix(rnorm(96),c(16,6))
for (i in 2:6) t.test(a[,1],a[,i])

のように、繰り返しでt検定を行う場合に、結果が表示されません。
なにか方法があったら、教えてお願いできますか?
ほかの検定でも同じことが起こるのでしょうか?宜しくお願いします。

  • オブジェクトを書き出すには print を使うのが基本であるということ。
    a<-matrix(rnorm(96),c(16,6))
    for (i in 2:6) print(t.test(a[,1],a[,i]))
    これで,うまくいくでしょう。
    統計学的な観点から言うと,このようなt検定の繰り返しが正しいかどうかは,あなたのデータ・目的にもよるのだが。具体的に言うと,多重比較を行うべきではないかと言うことですが。 -- 青木繁伸 2004-08-03 (火) 18:37:15
  • -- MKR? 2004-08-03 (火) 20:28:05
    x = sapply(2:6, function(i) {a<-matrix(rnorm(96),c(16,6)); 
                                  t.test(a[,1],a[,i])} )
  • sapply は実行結果の正確な出力が得られないのでお勧めではないです。
    結果は
               [,1]                      [,2]                      [,3]                     
    statistic   -0.3054952                0.3751029

添付ファイル: filemultiplot.png 1134件 [詳細] filehist2.png 1298件 [詳細] filepairspanel.png 2315件 [詳細] filepchtest.png 2370件 [詳細] filescplot1.jpeg 2446件 [詳細] fileplot.png 2755件 [詳細] fileoutput.png 1299件 [詳細] filehist3.png 2406件 [詳細] filecor.png 2282件 [詳細] filetest.png 2249件 [詳細] filebp.png 2388件 [詳細] filescatterplot3d.png 1394件 [詳細] filetopoheat.png 2368件 [詳細] filehm2.png 1565件 [詳細]

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