統計解析：CV%の合成（？）について

質問

ある分析法AのCV%がx、別の分析法BのCV%がy、と判っている時、それらの分析法で得られたデータを合成した値のCV%は、どう計算すれば良いでしょうか？

例1）
分析法Aで得られたデータをa、分析法Bで得られたデータをbとする時、「a/b」のCV%は？

例2）
分析法Aで得られたデータをa、分析法Bで得られたデータをbとする時、「a*b」のCV%は？

...よろしければ、こういう事が判りやすく書かれている実践的な内容の書籍についても、アドバイスを頂ければ助かります。

回答No.5

grothendieck · ベストアンサー

申し訳ありませんん。No3,4の回答は誤りのようです。
　　b = (1+ε)mb
より
　1/b ≒ (1-ε)/mb、　1/b^2 ≒ (1-2ε+ε^2)/mb^2
としましたが、1/b^2をε^2のオーダーまで計算するには1/bもε^2のオーダーまでとる必要がありました。そこで
　1/b ≒ (1-ε+ε^2)/mb、1/b^2 ≒ (1-2ε+3ε^2)/mb^2
に修正します。
　σ^2[a/b] = E[(a/b - E[a/b])^2]
　 = E[a^2]・E[1/b^2] - (ma・E[1/b])^2
に
　E[a^2] = (1 + CVa^2)ma^2
　E[1/b^2] = (1 + 3 CVb^2)/mb^2
　E[1/b] = (1 + CVb^2)/mb
を代入すると
σ^2[a/b] = (CVa^2 +CVb^2 + 3CVa^2・CVb^2)ma^2/mb^2

よって
CVa/b =√(CVa^2 +CVb^2 + 3CVa^2・CVb^2) /(1+CVb^2)

すると下のようにシミュレーションとの一致が良くなりました。（CVb<CVaという条件は必要ありません）

ma　sda　mb　sdb　乱数　　上の式
500　9　1000　3　0.0183　0.0182
500　15　1000　5　0.0305　0.0304
500　20　1000　5　0.0403　0.0403
1000　3　200　 9　0.0453　0.0450

間違えやすい点を多々含んでおり、（E[1/b]=1/E[b]とかついやってしまいそう）これを考えるのは良い演習問題になると思います。上の計算に使ったRの関数を次に示します。

cvawb <- function(ma,sda,mb,sdb){
a <- rnorm(500000,mean=ma,sd=sda)
b <- rnorm(500000,mean=mb,sd=sdb)
ab <- a/b
cab <- sd(ab)/mean(ab)
ca <- sda/ma
cb <- sdb/mb
cv <- sqrt(ca^2 +cb^2 +3*ca^2*cb^2)/(cb^2 +1)

return(list(cab,cv))
}

シミュレーションでやるのがどんな場合でも間違いなく確実にできるとも言えます。

お礼

再度の詳細な御説明、誠に有り難う御座います。大変助かっております。

素人目には、「3CVa^2」の「3」が、そういう数字が出るのだなと面白く感じられました。
また「シミュレーションでやるのが...」は、なるほどなと思わされました。今は、どこにでも「パーソナル」コンピュータがあり、こういう計算やシミュレーションも個人で行えてしまう、大変良い時代なのですね。

締め切ってしまうと、重ねての御質問がしにくくなり、惜しい気がするのですが...　ともあれ、誠に有り難う御座いました。また、どうぞ宜しく御願い致します。

回答No.4

grothendieck

CVa/b = √(CVa^2 + 1)/(CVb^2 + 1) - 1)

は

　CVa/b = √( (CVa^2 + 1)/(CVb^2 + 1) - 1)

に訂正させて頂きます。

お礼

再び詳細に御回答頂き、大変助かります。有り難う御座いました。
素人ながら薄々思っていたのですが、「a/bは一般には複雑になると思います」との事、やはりそうなのか...と感じました。簡単な計算では求められないのですね。
この中の「かつCVb<CVaのときの近似式」という点についてですが、この近似式はCVb>CVaだとどこを変えなくてはいけないのでしょうか。そもそも、お示し頂いている計算の過程が良く判っていないのでお恥ずかしいのですが、この式（最終の公式や、計算の過程）の中でCVb<CVaの場合とCVb>CVaの場合では何が違うのか、知っておかなくてはと思いました。もしも御時間が頂けましたら、恐縮ですがお教え下さい。

回答No.3

grothendieck

「CVab = √(CVa^2・CVb^2 + Cva^2 + CVb^2)」
は変動係数を
　Cv=(s/Xa)　sは標準偏差、Xaは平均
としたときのabの変動係数です。CV%だと
「CVab%=√(CVa%^2・CVb%^2/10000 + Cva%^2 + CVb%^2)」
となります。
a/bは一般には複雑になると思います。簡単には下の回答に書いた様に乱数を使ってシミュレートして評価することができると思います。CVbが１に比べて小さく、
かつCVb<CVaのときの近似式を次ぎに導きます。E[b]=mbとして確率変数bがとる値を

　b = (1+ε)mb
とすると
　1/b ≒ (1-ε)/mb
この近似の下でE[1/b]=1/mbであり、
　　E[a/b] = ma/mb
　σ^2[a/b] = E[(a/b - ma/mb)^2]
　 = E[a^2]/E[b^2] - (ma/mb)^2

(正確には E[1/b]=1/mb でない)などより

　CVa/b = √(CVa^2 + 1)/(CVb^2 + 1) - 1)

前と同様に乱数でシミュレートした結果とこの式の比較例は

ma　sda　mb　sdb　乱数　　上の式
500　9　1000　3　0.0182　0.0177
500　15　1000　5　0.0305　0.0296
500　20　1000　5　0.0403　0.0397

お礼

大変有り難う御座いました。
（次項に続きます）

回答No.2

grothendieck

無理なことなんかない。確率変数aとbが独立でE[a]=ma, E[b]=mb とします。このとき

　E[ab] = ma・mb
　σ^2[ab] = E[(a・b - ma・mb)^2]
　 = E[a^2]E[b^2] - (ma・mb)^2

これから、aの変動係数をCVa、bの変動係数をCVb、abの変動係数をCVabとすると

　CVab = √(CVa^2・CVb^2 + Cva^2 + CVb^2)

aとbが正規分布に従うとして、平均と標準偏差を与えた時、乱数で計算したabの変動係数と上の式で計算したCVabを比較する関数を統計ソフトR
http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html
で作りました。

cvab <- function(ma,sda,mb,sdb){
a <- rnorm(500000,mean=ma,sd=sda)
b <- rnorm(500000,mean=mb,sd=sdb)
ab <- a*b
cab <- sd(ab)/mean(ab)
ca <- sda/ma
cb <- sdb/mb
cv <- sqrt(ca^2*cb^2 + ca^2 + cb^2)
return(list(cab,cv))
}

aの平均ma、標準偏差sda、bの平均mb、標準偏差sdbを与えたときの計算例を示します。

ma　sda　mb　sdb　　乱数　　　上の式
10　 1　 20　 2　0.1414726　0.1417745
10　 1　 20　 4　0.2249844　0.2244994
10　 3　 20　 2　0.3179832　0.3176476
10　 3　 20　 4　0.3651448　0.3655133

なお統計に関する質問はここでするより
http://aoki2.si.gunma-u.ac.jp/lecture/mb.html
でした方が良いと思います（私は質問したことはありません）

お礼

御回答、大変有り難う御座いました。
済みません、素人なものでまだ理解が不足しているのですが、ここで
「CVab = √(CVa^2・CVb^2 + Cva^2 + CVb^2)」
としてお示し頂きましたのは「a*b」のCV%の計算式であるという理解でよろしいでしょうか。...まだ、これを応用出来るまでの力が無いもので、「a/b」のCV%ならどういう計算式になるか、併せて御教示頂けましたら大変助かります。

統計ソフトRや、「統計学関連のメッセージボード」の件、気になってはいながらそれぞれに敷居が高かったのですが、今後活用して行きたいと思います。
こういう分野は、困った時に助けを求められる相手が周囲にいませんので、grothendieckさんにこうして教えて頂ける様な事があると、ウェブというのは有り難いものだなと思います。

回答No.1

endlessriver

Cv=(s/Xa)*100
sは標準偏差、Xaは平均ですから無理と思います。

お礼

ありがとうございました。
当方は素人ですので、こういう計算・確認は荷が重いものです。