「統計学関連なんでもあり」の過去ログ--- 048

No.22600　NaNを1に置き換えたい　　【赤岳】　2018/08/21(Tue) 17:30

いつもお世話になっております。
次の文をRで指定しましたが，logの値がマイナスになり，エラーとなります。
logの中の乱数がマイナスやゼロの場合，1を返すようにしたいのですが，うまく文をかけません。
ご教示いただけないでしょうか。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊

> N <-100
> P<- rnorm(N, 5, 1)
> dP <- c(1:N)*0
> EP <- c(1:N)*0
> count<- 0
> for (i in 1:N) { 
+ dP[i] <- P[i] - log(rnorm (1, 1, 0.5))
+ EP[i] <- exp (dP[i])/exp(P[i])
+ if (EP[i] >0.5) 
+ count <- count+1
+ }
 log(rnorm(1, 1, 0.5)) で警告がありました:   計算結果が NaN になりました 
 if (EP[i] > 0.5) count <- count + 1 でエラー: 
   TRUE/FALSE が必要なところが欠損値です 
> count/N
[1] 0.41

No.22601　Re: NaNを1に置き換えたい　　【青木繁伸】　2018/08/21(Tue) 22:23

本筋ではないところをコメント
dP <- c(1:N)*0 は dp <- numeric(N) のほうが適切でしょう。
rnorm(1, 1, 0.5) が 0 以下のとき，log(rnorm(1, 1, 0.5)) を 1 にするのですか？

for ループを使っている場合は，特に何の工夫も必要ないです
N <-100
P <- rnorm(N, 5, 1)
dP <- c(1:N)*0
EP <- c(1:N)*0
count<- 0
for (i in 1:N) {
	x <- rnorm (1, 1, 0.5) # ここから
	if (x <= 0) {
		x <- 1
	} else {
		x <- log(x)
	}
	dP[i] <- P[i] - x # ここまで追加（変更）
	EP[i] <- exp (dP[i])/exp(P[i])
	if (EP[i] >0.5)
		count <- count+1
}
count/N
for ループをなくしてみましょう。ベクトル演算します。
N <- 100
p <- rnorm(N, 5, 1)
x <- rnorm(N, 1, 0.5) # まとめて N 個生成
x <- ifelse(x <= 0, exp(1), x) # x <= 0 のとき log(x) を 1 にするなら，x = exp(1) にしておく
dP <- P - log(x) # x が正なら log(x)，正でなければ log(exp(1)) = 1 になる
EP <- exp(dP) / exp(P) # EP の計算
mean(EP > 0.5) # count = sum(EP > 0.5); count/N は，mean(EP > 0.5)
注: ifelse(x <= 0, 1, log(x)) としても，一度全ての x に log(x) を適用して，x <= 0 のときに 1，そうでないときに log(x) の値を取るように動作するので，エラーが出る。
そこで，上のようにへんな（？）ことをする。

No.22602　Re: NaNを1に置き換えたい　　【赤岳】　2018/08/21(Tue) 22:53

青木先生
ご教示ありがとうございました。> mean(EP > 0.5) # count = sum(EP > 0.5); count/N は，mean(EP > 0.5)この mean がこんな風に使われるのは初めて知りました。
また，確かに ifelse(x <= 0, 1, log(x)) のようにしましたが，エラーメッセージはクリアーできませんでした。
解決しました。ありがとうございました。
（p <- rnorm(N, 5, 1) は，P <- rnorm(N, 5, 1) と修正しました。）