以下の行列を考えます.
行列X(M行N列,成分の値は分散1のガウスノイズ)
行列Xの分散共分散行列S(M行M列),
行列Sの逆行列Y(M行M列)
N=M+1のときはN≠M+1のときに比べて行列Yの各成分の大きさがかなり大きくなる,という結果が得られました.
(計算にはMATLABの関数covとpinvを用いました.)
これは数学的に正しいことなのでしょうか.
また,正しいとしたらどうしてこのようなことが起こるのでしょうか.
行列Yの計算結果の例を以下に示します.(行列Yを500回算出し,それらの平均値を示しています.)
http://wisteria.orz.ne.jp/download/pinvcovX.jpg
縦・横方向の軸は行列Yのインデックスを,高さ方向の軸は行列Yの成分の値を示しています.
上述した現象の原因についてご教授いただけると幸いです.
No.2ベストアンサー
- 回答日時:
Rで試してみたら同じような結果になりました。
その結果が添付画像です。
それぞれのz軸の上限は「z:上限」で表しています。
繰り返し回数を増やしたら各成分が小さくなると思っていたのですが、逆に大きくなりました。
Nが大きくなるにつれ、単位行列に近づいていくようです。
これ自体は分散共分散行列が単位行列に近づくことから当然ですか。
試したコード
func1 <- function(x, M)
{
y <- matrix(x, ncol = M)
z <- cov(y)
c(as.vector(z), as.vector(solve(z)))
}
func2 <- function(M = 20, N = 21, R = 500)
{
x <- matrix(rnorm(M*N*R), nrow = R)
y <- apply(x, 1, function(t) func1(t, M))
z <- rowMeans(y)
list(cov = matrix(z[1:(M*M)], ncol = M), covinv = matrix(z[(M*M+1):(2*M*M)], ncol = M))
}
x <- func2(20, 21, 10000)
y <- func2(20, 21, 500)
z <- func2(20, 22, 500)
w <- func2(20, 100, 500)
par(mfrow = c(2, 2))
persp(x$covinv, zlim = c(-0.5, 10^10), xlab = "x", main = "M = 20, N = 21, 10000回の平均, z:10^10")
persp(y$covinv, zlim = c(-0.5, 10^5), xlab = "x", main = "M = 20, N = 21, 500回の平均, z:10^5")
persp(z$covinv, zlim = c(-0.5, 10^3), xlab = "x", main = "M = 20, N = 22, 500回の平均, z:10^3")
persp(w$covinv, zlim = c(-0.5, 10), xlab = "x", main = "M = 20, N = 100, 500回の平均, z:10")
par(mfrow = c(1, 1))
検証していただきありがとうございます。
そうですか、計算結果は同じようになりましたか。
> 繰り返し回数を増やしたら各成分が小さくなると思っていたのですが、逆に大きくなりました。
有限長の小数で演算を繰り返すほどに、誤差が膨らんでいくという可能性が考えられますね。
私が用いた言語と異なる言語で検証していただき、たいへん参考になりました。ありがとうございました。
No.3
- 回答日時:
M列N行だと思いますが、
NがMに対して少ないので、
独立したガウスノイズになっているか疑問です。
なので逆行列が正確に求められるかも疑問です。
・・・がinvでも似たような結果になるのでは・・・・。
> M列N行だと思いますが、
> NがMに対して少ないので、
> 独立したガウスノイズになっているか疑問です。
N行M列でした。失礼しました。独立したガウスノイズになっていることも確認しました。
ご指摘ありがとうございます。
No.1
- 回答日時:
>行列X(M行N列,成分の値は分散1のガウスノイズ)
これは、N行M列の間違いですか?
で、発生させたガウスノイズ自体は独立(共分散=0)ということなんですかね。
ということは、求めた分散共分散行列は、ほとんど単位行列に近いものと思ってよいですか?
多分、その現象は、数値計算上の不安定性(誤差)によるものだと思います。S行列は、基本的には逆行列が存在するんだと思いますが、pinvではなくて、普通にinvを使えばどうでしょう。
ご回答ありがとうございます。
> N行M列の間違いですか?
はい、間違いでした。N行M列ですね。失礼しました。
> 発生させたガウスノイズ自体は独立(共分散=0)ということなんですかね。
> ということは、求めた分散共分散行列は、ほとんど単位行列に近いものと思ってよいですか?
はい、仰るとおりです。
> 多分、その現象は、数値計算上の不安定性(誤差)によるものだと思います。S行列は、基本的には逆行列が存在す> るんだと思いますが、pinvではなくて、普通にinvを使えばどうでしょう。
pinvの場合とinvの場合で計算結果はほぼ同じでした。
ただ、現象の原因としては、rabbit_catさんの仰るように、数値の誤差が膨らんだことは十分に考えられます。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- C言語・C++・C# このプログラミング誰か教えてくれませんか 1 2022/06/02 15:27
- Excel(エクセル) 別シートに毎回異なるデータをコピーする 7 2022/06/24 09:02
- Excel(エクセル) ExcelのIF関数について 4 2023/05/24 12:54
- C言語・C++・C# このプログラミングの問題を教えてほしいです。 キーボードからデータ数nとn個のデータを入力し、平均値 3 2022/12/19 22:51
- Excel(エクセル) エクセル 自動計算 1 2023/01/30 13:28
- 数学 確率について 事象Aが起こる確率が0.25である独立行列において、試行回数を5回とした時Aの起こった 2 2022/06/06 19:46
- その他(プログラミング・Web制作) パイソンのプログラミングについての質問です 2 2023/05/22 12:39
- Excel(エクセル) エクセルで2つの表を比較して、文字列が同じだが、その行のある値が違うものを抽出したい 1 2022/10/06 21:48
- Visual Basic(VBA) vba 重複データ合算 5 2023/07/05 18:55
- Excel(エクセル) エクセルで最初に値が入っているセルを見つける方法はありますか? 2 2023/07/18 14:58
このQ&Aを見た人はこんなQ&Aも見ています
おすすめ情報
このQ&Aを見た人がよく見るQ&A
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
3行3列の行列の和と積の計算...
-
線形代数です。 正方行列A,BがA...
-
行列式 証明
-
基本行列の積
-
[☆急いでます!!☆] 基本変形の解...
-
単因子の計算問題
-
核空間と像空間の求め方
-
直交補空間の問題が分かりませ...
-
にゃんこ先生の自作問題、ヴァ...
-
掃き出し法は、行の基本変形と...
-
行列を分割するメリットを教え...
-
diag(-1,1)
-
Vandermondeの逆行...
-
det(-A) = (-1)^n detA
-
AとBは同じサイズの正方行列と...
-
行列の階数(rank)を求める
-
一般逆行列の求め方
-
【数値解析】行列の可約、既約...
-
零因子の問題です・・
-
n次の正方行列の逆行列の行列式...
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
3行3列の行列の和と積の計算...
-
線形代数です。 正方行列A,BがA...
-
基本行列の積
-
行列式を帰納てきに求めるにあ...
-
数学「行列」の実生活への応用
-
行列と行列式の違いは?
-
matlabで条件をみたしたデータ...
-
逆行列(AB)^-1について
-
Aはn次正方行列で、どんなn次...
-
diag(-1,1)
-
直交補空間の問題が分かりませ...
-
行列の消去法のコツなど教えて...
-
大学数学を忘れました。3×3行列...
-
行列の階数(rank)を求める
-
行列の平方根?のようなもの
-
det(-A) = (-1)^n detA
-
[☆急いでます!!☆] 基本変形の解...
-
線形代数学のユニタリ行列の質問
-
行列式 証明
-
NPU付きのPCを買ったのですが、...
おすすめ情報