ホーム  •  ニュース  •  フォーラム  •  アカウント情報  •  サイト内検索  •  新規登録
 ログイン
ユーザ名:

パスワード:


パスワード紛失

新規登録
 カウントダウンブロック
カウントダウンイベントはありません
 メニュー
メイン
   地球科学のための R
     R の Tips		投稿するにはまず登録を

スレッド表示 | 新しいものから 前のトピック | 次のトピック | 下へ
投稿者 スレッド
森
投稿日時: 2004-10-26 12:52
登録日: 2004-7-29
居住地: 地球
投稿: 303
R の Tips
古いホームページに掲載していた R に関するちょっとしたことを投稿していきます。似たような Tips がありましたら、ここにどうぞ。
森
投稿日時: 2004-10-26 12:54
登録日: 2004-7-29
居住地: 地球
投稿: 303
多次元配列の作り方

  x <- 1:16
  dim(x) <- c(2,2,2,2)

とすると、

  x[,,1:2,1]

等として要素にアクセスできます。~

NCEP/NCARのデータは、三次元配列に変換するととっても便利です。
森
投稿日時: 2004-10-26 13:00
登録日: 2004-7-29
居住地: 地球
投稿: 303
プログラミングに関して-その1-関数の切り分け
** はじめに

プログラミングの質は、関数の作り方にかかっているといって良いと思います。
プログラミングの考え方に、「構造化プログラミング」「オブジェクト指向」といったものがありますが、どれも、主に関数の作り方に関する考え方です。それだけ、関数の作り方が大切である、ということです。

以下に、その例を書きます。例えば、常微分方程式を解くことを考えましょう。一般に、

  dx/dt = f(x)

という関数を積分する状況を考えるわけですが、f(x) は、もちろん、いろいろな場合がありえます。このような場合、どのように対処するのが良いでしょうか。次の二つを比較して下さい。

+ 毎度、プログラムの該当個所をさがして変更する。
+ f(x) を使っている場所がはっきりするようにプログラミングしておき、そこだけ修正する。

どちらが良いかは明らかです。使い捨てのプログラムなら、前者でもいいのですが、プログラムを別の場所に再利用したいと思ったら、後者であるべきですし、そのために、書き捨てのプログラムの場合でも普段から心がけるべきです。

同様に、積分のアルゴリズムの部分も、独立した関数にしておくことが望ましいのは、明らかです。他の問題にも使えるようになります。

さらに、こうした観点からプログラミングをすると、全体に機能別に関数が分離され、再利用性だけでなく、プログラミング自身が読みやすくなります。これも大変重要なことです。なぜなら、プログラミングにかかる時間の90%は、デバッグ(プログラムの問題点を見つけ出して修正・改善すること)だからです。

** ここまでのまとめ

まとめ
-- 上手に関数をつくろう。次のような点に気を配ろう。
+++ 変更する可能性がある場所はどこか
+++ 再利用性・可搬性・汎用性
+++ デバッグのしやすさ

** 実践例
桑子君が書いてくれている例では、例えば、

 Method.Euler <- function( dxdt, t, dt, x ){
     x <- x + dxdt(x,t) * dt
     t <- t + dt
 
     return( t=t, x=x )
 }

 DxDt.Simple <- function( x, t ){
     return( x )
 }

などとしておけばよいのです。すると、

 Method.Euler( DxDt.Simple, 0, 0.01, 1 )

なんてできる訳です。これを何回も繰り返すループの中に入れれば、あとは図にするにしても何にしても全体の見通しが大変良くなる訳です。アルゴリズムや、f(x) を変えるにしても、この行だけ変えればいいのですから。

是非、試みて下さい。そして、どこまで見通しがよいプログラミングができるか、挑戦してみて下さい。そして、それが、最初に作ったものとどれほど違うのか、よく見比べて下さい。

最後に、その成果をここに書いて下さい。
森
投稿日時: 2004-10-26 13:05
登録日: 2004-7-29
居住地: 地球
投稿: 303
プログラミングに関して-その1-関数の切り分け(例)
* プログラミングの例

たとえば、何も工夫をしない場合に、オイラー法で次のように積分していたとします。

 plot(0, type='n', xlim=c(0,10), ylim=c(0,2000))
 data <- list( x=1, t=0 )
 dt   <- 0.1
 for( i in 1:100 ){
   data$x <- data$x + data$x * dt
   data$t <- data$t + dt
   points( data$t, data$x )
 }

いろいろな問題点があります。
+ data$x が二回現れるがそれぞれの意味がわかりにくい。
+ 同じことであるが、どのような積分方法をとっているのか、と、dx/dt = f(x) の f(x)の部分がどうなっているのか、が分離できていないのでプログラム全体が読みずらい。
+ 積分方法を変える、f(x)を変える、というときにどのように工夫すればいいのかわかりずらい。
+ 再利用できない。
+ などなど。

そこで、次のようにします。


 #- オイラー法 ----------------------------------------------#
 
 Method.Euler <- function( dxdt, data, dt ){
   data$x <- data$x + dxdt( data$x, data$t ) * dt
   data$t <- data$t + dt
   return( data )
 }



 #- 1次元の関数の例 -----------------------------------------#
 
 DxDt.Simple <- function( data, time ){
   return( data )
 }



 #- 実行例 --------------------------------------------------#
 plot(0,type='n', xlim=c(0,10), ylim=c(0,2000))
 data <- list( x=1, t=0 )
 for( i in 1:1000 ){
   data <- Method.Euler( DxDt.Simple, data, 0.01 )
   points( data$t, data$x )
 }

このようにすることで、オイラー法の積分部分も独立して利用できます。~
プロットする部分は、積分法、f(x)、の他に、初期値によっても変更することが必要です。このようにたくさんの部分で頻繁に変更する部分は、関数に取り込む必要性は低いです。関数にすることで、かえってどの変数がどんな役割を担っているのかわかりにくくなる可能性があるからです。

さらに、今回のプログラミングでは、オイラー法の部分は多次元でもそのまま使えてしまう、という利点があります。2次元の問題では、Method.Euler() 部分はそのままで、次のようにできます。


 #- 2次元の例 -(単振動)--------------------------------------#
 
 DxDt.2d.Simple <- function( data, time ){
   return( c( matrix( c(0, -1, +1, 0 ), 2 ) %*% matrix(data,2) ) )
 }



 #- 実行例 --------------------------------------------------#
 plot(0,type='n', xlim=c(0,10), ylim=c(-2,2))
 data <- list( x=c(1,0), t=0 )
 for( i in 1:1000 ){
   data <- Method.Euler( DxDt.2d.Simple, data, 0.01 )
   points( data$t, data$x[1] )
 }


試してみて下さい。
森
投稿日時: 2004-10-26 13:08
登録日: 2004-7-29
居住地: 地球
投稿: 303
いっぱい図を出すときに止めながら見たい！

図を出すときに、ダーッと出力したいときがありますが、パカパカ画面が変わったのでは大事なところもきちんと見ることができません。そこで、図を止めながら見る方法を。

例えば、こんな感じです。data[] の各行の値を順番にプロットしたい場合：

 par(ask=interactive())
 
 for( i in 1:dim(data)[2] ){
   print(i)
   plot( data[i,] )
 }
 
 par(ask=F)

par() がポイント。~
はじめの par() で、リターンキーを押さないと絵を描かないようにし、二つ目の par() で、それを解除している。

2003-12-04 (木) 11:47:56
森
投稿日時: 2004-10-26 13:09
登録日: 2004-7-29
居住地: 地球
投稿: 303
postscript() を使うとき、大きさや向きをいちいち指定するのが面倒
そんなあなたに ps.options()。

 ps.options(horizontal=F,width=7,height=7)

としておけば、大きさは 7インチ四方、TeX にも取り込みやすい、縦向きの
指定ができます。

あとは、postscript("hogehoge.ps") と、ファイル名を指定するだけでOK

2003-12-04 (木) 18:35:42
森
投稿日時: 2005-1-13 12:57
登録日: 2004-7-29
居住地: 地球
投稿: 303
軸に描く値を変更したい！
グラフを描くときに、軸の値をどのように変更するか、について、整理しておきます。

1.標準的な軸をかかないようにする。frame.plot=F にすると、外側の枠を描かなくなる。
 plot( .....,  axes=F,  frame.plot=T )

2. 軸を追加する。
axis(1,c(3,4,5,6,8,10))
axis(2, 1/c(0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2 ),  c(0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2 ))

axis(1,... )で横軸を描く。2番目の引数は、数字を入れる場所の値を指定。
axis(2,...)で縦軸を描く。2番目の引数は、やはり、数字を入れる場所の本来の値を指定する。3番目の引数では、その場所に入れる値を指定する。どちらも、ベクトルで指定。
森
投稿日時: 2005-10-18 20:17
登録日: 2004-7-29
居住地: 地球
投稿: 303
東西南北の縮尺が(ほぼ)同じであるような地図を描く
関東地方近辺であれば、緯度と経度の１度あたりの長さの比は一定であると思って良い。そこで、この仮定に基づいて、東西と南北の縮尺が同じであるような地図を描いてみる。
###########################################################
#
#       どこでもミュージアムエコ用の地図
#
###########################################################
#
# 指定した点(base)から、北と東に max [km] ずつの領域を描く。
# ticlen km ごとに、軸にマークを入れる。
#
# 地図の描画は、
#       source("~/lib/R/MAP/map.R")
#       map()
# で行う。緯度経度の座標系で地図が描けるのであれば、
# 他のものでもよい。
#
###########################################################
source("~/lib/R/MAP/map.R")
#
# 原点と、幅の計算
#
base     <- c(138.5, 34.5)
max      <- 220
ticlen   <- 10                               # 10km ごとに tic
deglen   <- max / 6370 / pi * 180
latlen   <- deglen/cos((base[2]+deglen/2)/180*pi)

#
#   正方形   軸を描かない
#
par(pty="s", xaxs='i', yaxs='i', xaxt='n', yaxt='n')
plot(type='n',
     xlim=c(base[1],base[1]+deglen/cos((base[2]+deglen/2)/180*pi)),
     ylim=c(base[2],base[2]+deglen),
     xlab='', ylab='',
     x=0)
map()
#
#   軸を描き加える
#
par(xaxt='s', yaxt='s')
axis(1, at=seq(base[1], base[1]+latlen, latlen / max * ticlen ), pos=base[2])
axis(2, at=seq(base[2], base[2]+deglen, deglen / max * ticlen ), pos=base[1])
AM
投稿日時: 2005-10-18 20:19
管理人
登録日: 2004-7-27
居住地:
投稿: 16
不偏分散と標本分散
分散と一口に言うが、統計学的には、標本分散と不偏分散がある。通常、無限の数の母集団の中から標本を抽出し、分散を推定するという作業が行われるが、不偏分散を求めることで母集団の分散を推定できる。R では、通常、分散を求めるときには、不偏分散が使われる。

ところが、幾つかの地球科学的なアプリケーションの場合には、標本分散が想定されていることがあるので注意が必要である。例えば、分散の値を１に規格化して何かを行うようなプログラムでは、標準装備の関数 scale() を利用するのではなく、次のような関数を使うべきだ。

spscale <- function(x, center = TRUE, scale = TRUE){
  return ( scale(x,center,scale) / sqrt ( (length(x)-1) / length(x) ) )
}
スレッド表示 | 新しいものから 前のトピック | 次のトピック | トップ

投稿するにはまず登録を
 

WWW を検索 meteorology.jp を検索

Powered by XOOPS 2.0 © 2001-2006 The XOOPS Project, Maitained by A. Mori
FI Theme :: XOOPS 2 Theme by ImageSquare :: Costomized by matchan and A.Mori