多変量正規分布(その2)
前回は2変量に相関がないケースを見たが、今回は相関があるケースを見てみる。
各変量の平均と分散は前回と同じ。相関係数だけ変えてみる。
分散共分散行列は以下のとおりとなる。
> Sigma.positive_cor [,1] [,2] [1,] 2 1 [2,] 1 4 set.seed(123); data <- rmvnorm(n=N, mean=c(x_mu, y_mu), sigma=Sigma.positive_cor) [,1] [,2] [1,] 0.15661631 2.3779756 [2,] 3.17597129 3.6033091 [3,] 1.68913945 6.4304202 (snip) [9998,] -0.8679178275 1.856718741 [9999,] 1.1757277988 4.769811627 [10000,] 3.8768372227 5.398123198
チェック。
# xの標本平均 mean(data[,1]) [1] 0.9964789 # yの標本平均 mean(data[,2]) [1] 2.977462 # xの標本分散(自作関数を使用) variance(data[,1]) [1] 2.010994 # yの標本分散(自作関数を使用) variance(data[,2]) [1] 3.980254 # xとyの標本共分散(自作関数を使用) covariance(data[,1], data[,2]) [1] 1.001641 # xとyの相関係数 > cor(data[,1], data[,2]) [1] 0.3540394
2次元散布図
3次元散布図
変量xと変量yの平均値での確率密度を求めてみる。
> x_mu <- 1 > y_mu <- 3 > dmvnorm(matrix(c(1, 3), ncol=2), mean=c(x_mu, y_mu), sigma=Sigma.positive_cor) [1] 0.06015491
相関がないケースでの確率密度は「0.05626977」であった。
改めて確率密度の式を見てみる。
分母の行列式の部分を比べてみる。
# 相関係数が0.35 > sqrt(det(Sigma.positive_cor)) [1] 2.645751 # 相関係数が0 > sqrt(det(Sigma.zero)) [1] 2.828427
相関がある方が、行列式の値が小さくなっている。
これは、確率密度が描く山を上に押し上げている要素だと思われる。
比率で見てみる。
# 確率密度での比較(相関係数が0の場合を基準にして何倍か) > 0.06015491 / 0.05626977 [1] 1.069045 # 行列式の平方根での比較(相関係数が0の場合を基準にして何倍か) > 2.828427 / 2.645751 [1] 1.069045
また、指数関数の中の、分散共分散行列の逆行列を比べてみる。
# 相関係数が0.35 > solve(Sigma.positive_cor) [,1] [,2] [1,] 0.5714286 -0.1428571 [2,] -0.1428571 0.2857143 # 相関係数が0 > solve(Sigma.zero) [,1] [,2] [1,] 0.5 0.00 [2,] 0.0 0.25
これは、確率密度が描く山を楕円から扁平に押しつぶしている要素だと思われる。
視点を変えてみる。変量yに着目した形状。
左図が相関あり、右図が相関なし
上から見下ろした形状。
左図が相関あり、右図が相関なし
各変量が互いに関連していれば、同時に分析することで、データの特徴が分かる。
一方、各変量に相関がない場合は、前回見たように、2変量の確率密度は1変量正規分布の確率密度の積になる。
これは同時に分析せずとも、それぞれの変量を別個に調べればよいことを示している。
多変量正規分布
多変量混合正規分布のパラメータ推定に入る前に、どのようなモデルを扱おうとしているのか整理しておく。
以下の視点で考えた時、
- 変量が「1」or「多」
- 「混合分布である」or「混合分布ではない」
以下の4パターンに分類される。
今回は上記3の「多変量正規分布」について見てみる。
変量の数が「2」、「変量x」と「変量y」で表されるとすると、確率密度関数は以下のとおり。
ここで、∑は「変量x」と「変量y」の分散共分散行列を表す。
以下のサイトを参考にさせていただいた。
http://cse.naro.affrc.go.jp/minaka/R/R-binormal.html
具体的に値を求めてみる。
であるとすると、分散共分散行列は以下のとおりとなる。
このデータはどんな形をしているか。
Rの「mvtnorm」パッケージで、多変量正規分布の乱数を作れる「rmvnorm」関数があるので、これを使う。
# 各変量の平均 x_mu <- 1 y_mu <- 3 # 各変量の分散と共分散 x_sigma2 <- 2 y_sigma2 <- 4 xy_cov <- 0 # 分散共分散行列 Sigma.zero <- matrix(c(x_sigma2, xy_cov, xy_cov, y_sigma2), nrow=2, byrow=TRUE) # 乱数生成 N <- 10000 set.seed(123); data <- rmvnorm(n=N, mean=c(x_mu, y_mu), sigma=Sigma.zero) [,1] [,2] [1,] 0.2073677 2.539645 [2,] 3.2043464 3.141017 [3,] 1.1828405 6.430130 (snip) [9998,] -0.7927958106 2.2253545108 [9999,] 0.9042444874 4.8101121211 [10000,] 3.6597604578 4.8591241077
チェック。
標本分散、標本共分散については自作関数を使用。
work-in-progress.hatenablog.com
# xの標本平均 mean(data[,1]) [1] 0.9998627 # yの標本平均 mean(data[,2]) [1] 2.977238 # xの標本分散(自作関数を使用) variance(data[,1]) [1] 2.011234 # yの標本分散(自作関数を使用) variance(data[,2]) [1] 3.978216 # xとyの標本共分散(自作関数を使用) covariance(data[,1], data[,2]) [1] 0.002540545 # xとyの相関係数 > cor(data[,1], data[,2]) [1] 0.0008981552
2次元散布図
3次元散布図
確率密度の計算。
先ほどの式をRで実装してみる。
f_2var_norm <- function(x, y, xmean, ymean, Sigma){ # x と y をセンタリングしたcolumn vector(縦ベクトル) A <- matrix(c((x-xmean), (y-ymean)), nrow=2, byrow=TRUE) # 指数関数のべき乗部分 power <- -0.5 * (t(A) %*% solve(Sigma) %*% A) denominator <- 2 * pi * sqrt(det(Sigma)) return ( 1 / denominator * exp(power) ) }
x=0.5、y=2.5の確率密度を計算してみる。
x_mu <- 1 y_mu <- 3 f_2var_norm(0.5, 2.5, x_mu, y_mu, Sigma.zero) [,1] [1,] 0.05123421
検算。
「mvtnorm」パッケージの「dmvnorm」関数の値。
x_mu <- 1 y_mu <- 3 dmvnorm(matrix(c(0.5, 2.5), ncol=2), mean=c(x_mu, y_mu), sigma=Sigma.zero) [1] 0.05123421
一致していることを確認。
各変量(xとy)ごとの形状を見るため、横からの角度でプロットしてみる。
左図は変量xに着目した分布、右図は変量yに着目した分布
比較用として、1変量の確率密度を出力してみる。
1変量正規分布における、平均値での確率密度は
# N(1, 2)における、平均値(=1)での確率密度 > dnorm(1, mean=1, sd=sqrt(2)) [1] 0.2820948 # N(3, 4)における、平均値(=3)での確率密度 > dnorm(3, mean=3, sd=sqrt(4)) [1] 0.1994711
一方、2変量正規分布における、変量xと変量yの平均値での確率密度は
> x_mu <- 1 > y_mu <- 3 > dmvnorm(matrix(c(1, 3), ncol=2), mean=c(x_mu, y_mu), sigma=Sigma.zero) [1] 0.05626977
これは、1変量正規分布における平均値での確率密度を掛け合わせたものと等しい。
> dnorm(1, mean=1, sd=sqrt(2)) * dnorm(3, mean=3, sd=sqrt(4)) [1] 0.05626977
パラメータ推定(混合正規分布)
混合分布を使うことで、より柔軟に分布をモデル化できる。
例えば、ピークが2つあるような分布。
これは以下2つの分布を混合したもの。
分布1 : 正規分布N(2, 5)、混合比「0.6」
分布2 : 正規分布N(10, 5)、混合比「0.4」
この混合分布は音声認識において使われていて、隠れマルコフモデル(HMM)の各状態の出力確率が混合正規分布(混合ガウス分布)に従うとしたモデルを「GMM-HMM」と呼ぶ。
具体的にどのようにモデル化しているかは以下の記事に書いた。
work-in-progress.hatenablog.com
混合数をKとおく。
混合パラメータ
観測データがN個のときの対数尤度関数
xが与えられたとき、それに対応する確率密度を各分布ごとに分解した割合を仮置きするために、潜在変数zを導入する。
潜在変数zは、K個の要素からなり、どれか1つだけ1を取り、他は全て0となる。
xについてn番目の値が与えられたとき、k番目のzが1となる確率を負担率と呼び、以下の式で計算する。
分母は、混合分布における、n番目のxの確率密度。
分子は、k番目の分布における、n番目のxの確率密度となるので、負担率は
この辺は以下のサイトを参考にさせていただいた。
混合ガウスモデルとEM - 人工知能に関する断創録
負担率はN個のxのそれぞれについて求まる。
この負担率を使って、パラメータ(平均、分散、混合比)の最尤値を計算する。
負担率を求めるためにはパラメータが必要で、パラメータの最尤値を計算するためには負担率が必要と、循環してしまうので、以下の2ステップに分けて計算する。
- Eステップ:現在の混合比、平均、分散に対応する負担率を計算する
- Mステップ:Eステップで計算した負担率をもとに、混合比、平均、分散の最尤値を求める
混合ガウス分布の数が「2」、各分布が以下であるサンプルデータ「500」個用意し、推定の過程を追ってみる。
分布1 : 正規分布、平均「2」、分散「5」、混合比「0.6」
分布2 : 正規分布、平均「10」、分散「5」、混合比「0.4」
rnorm関数に与えるseedは「123」を指定。
[1] 0.7467384 1.4853075 5.4853777 2.1576616 2.2890962 9.7872438 [7] 3.0306400 -0.8287629 0.4641503 1.0034695 9.7822757 2.8045682 [13] 2.8961522 10.4597300 8.6216076 5.9956592 9.7052800 -2.3974896 (snip) [486] 3.741583972 4.045495218 11.760084269 5.637816986 13.701085639 [491] 11.511050698 5.139551842 4.893659583 9.523030728 10.700400827 [496] 9.798809343 12.393746629 2.368595386 8.831393539 3.234662183
密度プロット(赤い線はcurve関数で書いた真の確率密度関数)
「R」には「mixtool」ライブラリで「normalmixEM」関数が用意されているので、これを使えば最尤値を求められる。
検算の意味で最初に値を求めておく。
mixtool_val <- normalmixEM(data, k=2)
最尤値での対数尤度
mixtool_val$all.loglik[length(mixtool_val$all.loglik)] -1388.84
最尤値(混合比、平均、分散)
mixtool_val$lambda [1] 0.589741 0.410259 # 混合比 mixtool_val$mu [1] 1.969816 10.129432 # 平均 mixtool_val$sigma^2 [1] 4.498699 4.566088 # 分散
先ほどの式で最尤値を100回、計算してみた結果。
初期値は、以下の値を与えた。
分布1 : 混合比0.2、平均1、分散2
分布2 : 混合比0.8、平均5、分散2
負担率の推移。
パラメータ推定(ニュートン・ラフソン法)
前回の続き。
「ニュートン・ラフソン法」を使って、2変数パラメータの最尤推定値を求めてみる。
以下のサイトを参考にさせていただいた。
http://www012.upp.so-net.ne.jp/doi/sas/numerical/Newton_Raphson/Newton_Raphson.pdf
正規分布の対数尤度関数を再掲。
データ
# 平均2、分散5の正規分布に従うサンプルを200個 set.seed(123); data <- rnorm(200, mean=2, sd=sqrt(5)) [1] 0.746738355 1.485307487 5.485377748 2.157661556 2.289096164 [6] 5.835001897 3.030639969 -0.828762916 0.464150333 1.003469540 [11] 4.737130109 2.804568177 2.896152207 2.247494077 0.757101438 (snip) [186] 1.559101297 4.481857217 2.189478346 3.686115522 0.883549109 [191] 2.479514290 1.273980231 2.211495199 -0.002093333 -0.931041334 [196] 6.465904894 3.343225764 -0.797927823 0.633391465 -0.650814055
初期パラメータは、平均1、分散1とおいて、先ほどの式を計算した結果。
最尤値となるパラメータは、平均「1.980836」、分散「4.425514」
検算。Rの関数で計算した結果。
# 準ニュートン法による最適解 > optim(c(1, 1), fn=log_likelihood_norm, obs=data, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
結果
$par [1] 1.980839 4.425508 $value [1] -432.5263 $counts function gradient 31 10 $convergence [1] 0 $message NULL $hessian [,1] [,2] [1,] -4.519255e+01 3.119993e-05 [2,] 3.119993e-05 -5.105931e+00
小数点第五位まで一致していることが確認できた。
また、k回目の最尤値でのヘッセ行列を固有値分解してみる。
H_0 <- matrix(c(-200.000000 , -196.167181 , -196.167181 , -977.510617 ), nrow=2) H_1 <- matrix(c(-152.232794 , -35.661895 , -35.661895 , -340.743870 ), nrow=2) H_2 <- matrix(c(-108.715893 , -7.280127 , -7.280127 , -113.101423 ), nrow=2) H_3 <- matrix(c( -79.549506 , -1.429134 , -1.429134 , -39.900363 ), nrow=2) H_4 <- matrix(c( -61.130669 , -0.254284 , -0.254284 , -15.933080 ), nrow=2) H_5 <- matrix(c( -50.656386 , -0.036104 , -0.036104 , -7.966423 ), nrow=2) H_6 <- matrix(c( -46.162098 , -0.002919 , -0.002919 , -5.555946 ), nrow=2) H_7 <- matrix(c( -45.231581 , -0.000058 , -0.000058 , -5.123588 ), nrow=2) H_8 <- matrix(c( -45.192558 , -0.000000 , -0.000000 , -5.105933 ), nrow=2) H_9 <- matrix(c( -45.192490 , -0.000000 , -0.000000 , -5.105903 ), nrow=2)
固有値分解の結果
最尤値をとるパラメータを与えたとき、ヘッセ行列の固有ベクトルは標準基底になっている。
また、各パラメータの2階偏微分の値が固有値になっている。
推定するパラメータが2つの場合、フィッシャー情報行列は以下として定義される。
今回は1セットのデータのみ用いており、パラメータ推定値も1セットのみであるから、フィッシャー情報行列は、対数尤度関数の2階偏微分にマイナスを掛けたものになる。
今回、各繰り返しにおけるフィッシャー情報量と対数尤度は以下のとおりであった。
平均を固定して、分散を変えてみると以下のとおりであった。
フィッシャー情報行列は、対数尤度関数のスコア関数(1階偏微分)の分散として定義されるから、以下のことがいえる。
- スコア関数の分散が大きい(=フィッシャー情報量が大きい)ほど、相対的に尤度は下がる
- スコア関数の分散が小さい(=フィッシャー情報量が小さい)ほど、相対的に尤度は上がる
パラメータ推定(最尤推定法)
今、手元にn個の観測値(標本)が得られているとする。
統計モデルとして表現するには、
が必要となる。
1番目の確率分布の形が分かったとして、2番目のパラメータを推定する方法の1つである「最尤推定法」に関するメモ。
ある1つの観測値が、他の観測値に影響を及ぼさない(確率変数は独立している、という)場合、それぞれの観測値の確率密度を掛けあわせたもの(同時確率)が尤度関数となる。
観測値つまり確率変数をx、未知パラメータをθとおくと、尤度関数は以下の形で表される。
\[
L(\theta) = \displaystyle \prod_{ i = 1 }^n f(x_i | \theta)
\]
尤度関数すなわち同時確率を最大にするパラメータθの値が「最尤推定値」になる。
小数の乗算はアンダーフローを引き起こすので、尤度関数の対数をとった「対数尤度関数」を最大化するパラメータθを求める。
具体的に見てみる。
確率分布の形が「正規分布」だとすると、推定したいパラメータθ は、
となる。
(μは平均、σは標準偏差を表す)
を対数尤度関数の式にあてはめて、式変形する。
導出過程は以下のサイトを参照。
統計学的画像再構成法である
これをRの関数として定義する。
第1項は定数なので、最大化を考えるにあたっては無視できるが、今回は確認用に残しておく。
対数尤度関数(正規分布)
# arg1 : vector(size 2), [1]mean, [2]variance # arg2 : vector(size N), observation log_likelihood_norm <- function(param, obs){ mu <- param[1]; sigma2 <- param[2]; term1 <- (-1) * length(obs) * 0.5 * log(2*pi) term2 <- (-1) * length(obs) * 0.5 * log(sigma2); term3 <- (-1) * 0.5 / sigma2 * sum((obs-mu)^2); return ( term1 + term2 + term3 ); }
検算。
平均「2」、分散「5」の正規分布に従う乱数を10個用意する。
> data <- rnorm(10, mean=2, sd=sqrt(5)) > data [1] -0.01731340 1.79861872 2.99598784 1.78515021 2.42412228 -1.61575147 [7] 1.74058535 3.84369888 1.93563703 0.03983286
Rのdnorm関数を使って求めた結果
value <- 0 for (i in 1:length(data)) { value <- value + log(dnorm(data[i], mean=2, sd=sqrt(5))) } value # -19.80805 value <- 1 for (i in 1:length(data)) { value <- value * dnorm(data[i], mean=2, sd=sqrt(5)) } log(value) # -19.80805
作成した関数での結果
log_likelihood_norm(c(2, 5), data) # -19.80805
一致していることを確認。
値が最大になるパラメータをRのoptim関数を使って探してみる。
最適解(サンプル数10)
optim(c(0.5, 1), fn=log_likelihood_norm, obs=data, control=list(fnscale=-1), hessian=FALSE) # $par # [1] 1.493105 2.313280 # $value # [1] -18.38531 # $counts # function gradient 73 NA # $convergence # [1] 0
出てきた最適値は、平均が「1.493105」、分散が「2.313280」
統計量(サンプル数10)
mean(data) # [1] 1.493057 # 不偏分散 var(data) # [1] 2.571645 # 標本分散 variance <- function(x) { var(x) * (length(x)-1) / length(x) } variance(data) # [1] 2.31448
サンプル数を100に増やしてみる。
最適解(サンプル数100)
data <- rnorm(100, mean=2, sd=sqrt(5)) optim(c(0.5, 1), fn=log_likelihood_norm, obs=data, control=list(fnscale=-1), hessian=FALSE) # $par # [1] 1.918190 4.851344 # $value # [1] -220.8375 # $counts # function gradient 59 NA # $convergence # [1] 0
出てきた最適値は、平均が「1.918190」、分散が「4.851344」
統計量(サンプル数100)
mean(data) # [1] 1.917897 var(data) # [1] 4.898471 variance(data) # [1] 4.849486
さらにサンプル数を増やしてみる。
最適解(サンプル数1000)
data <- rnorm(1000, mean=2, sd=sqrt(5)) optim(c(0.5, 1), fn=log_likelihood_norm, obs=data, control=list(fnscale=-1), hessian=FALSE) # $par # [1] 1.891001 4.752660 # $value # [1] -2198.32 # $counts # function gradient 63 NA # $convergence # [1] 0
統計量(サンプル数1000)
mean(data) # [1] 1.891252 var(data) # [1] 4.7577 variance(data) # [1] 4.752943
サンプル数が多くなるほど、標本分散と不偏分散の差が少なくなっている。
対数尤度関数が返す値の3次元散布図を書いてみる。
meanが「2の前後3、-1から5まで1刻みの7通り」、sigma2は「5の前後3、2から8まで1刻みの7通り」の49パターンとした。
一番高くなるところのパラメータが最尤推定値となる。
最尤推定値を求めるにあたって、Rのoptim関数では以下の最適化アルゴリズムが使える模様。
method="Nelder-Mead" # シンプレックス法、関数が微分できなくても可 method="BFGS" # 準ニュートン法 method="L-BFGS-B" # 準ニュートン法の変種(推定するパラメータのlower/upperを指定) method="CG" # 共役勾配法 method="SANN" # シミュレーテッド・アニーリング法、関数が微分できなくても可
methodパラメータの省略時デフォルトは「Nelder-Mead」(シンプレックス法)
optim function | R Documentation
The default method is an implementation of that of Nelder and Mead (1965), that uses only function values and is robust but relatively slow.
It will work reasonably well for non-differentiable functions.
(参考)多変量正規分布のパラメータ推定で、optim関数のmethodを変えて比較してみた結果
パラメータ推定(多変量正規分布) - ichou1のブログ
ヘッセ行列の固有値による凸性の判定
最適化を考える時、どの方向に向かって修正していくかを決める必要がある。
扱う対象が関数化できるなら、その最大値(あるいは最小値)を求める問題に帰着するので、関数の凸性が大事な概念になってくる。
今回は2変数関数を対象に、2階偏微分であるヘッセ行列を使って凸性を調べる方法を見てみる。
以下のサイトを参考にさせていただいた。
http://www2.kaiyodai.ac.jp/~yoshi-s/Lectures/Optimization/2013/lecture_1.pdf
偏微分の順番で微分結果は変わらないので、ヘッセ行列は対称行列になる。
3パターンの関数で見てみる。
[パターン1] 狭義凸関数
ヘッセ行列が正定値(固有値が全て0より大きい)
グラフが全ての方向で丸く突き出ている
[パターン2] 凸関数(狭義凸関数ではない)
ヘッセ行列が半正定値(固有値が全て0以上)
2点を結ぶ線分とグラフが重なる部分がある
勾配ベクトル
ヘッセ行列
ヘッセ行列を固有値分解する。
固有値は「2」であることを確認
> eigen(matrix(c(2, 0, 0, 2), nrow=2, ncol=2, byrow=TRUE)) eigen() decomposition $values [1] 2 2 $vectors [,1] [,2] [1,] 0 -1 [2,] 1 0
角度を変えてみる
y軸が奥行きと重なるように見たところ
x軸が奥行きと重なるように見たところ
全ての方向で丸く突き出た形をしている。
[パターン2] 凸関数(狭義凸関数ではない)
勾配ベクトル
ヘッセ行列
ヘッセ行列を固有値分解する。
固有値は「2」と「0」であることを確認
> eigen(matrix(c(0, 0, 0, 2), nrow=2, ncol=2, byrow=TRUE)) eigen() decomposition $values [1] 2 0 $vectors [,1] [,2] [1,] 0 -1 [2,] 1 0
角度を変えてみる。
y軸が奥行きと重なるように見たところ
グラフの底が、線分の形をしている。
x軸が奥行きと重なるように見たところ
[パターン3] 凸面と凹面が混在する関数
勾配ベクトル
行列の低ランク近似と近似度
特異値分解を使った「行列の低ランク近似」に関するメモ。
まずは、視覚的に確認しやすい画像データで試してみる。
高さ200px、幅320px、インデックスカラーのRAWデータを用意する。
1pxを1バイトで表現、0〜255の範囲の値をとる200行 x 320列の行列ができる。
以下は左上を拡大した画像。
先頭6バイト分のみ値を添えてみた。値が小さいほど黒に近く、値が大きいほど白に近くなる。
CSVに落としたものをRに読み込み、
A <- data.matrix(read.table("200_320_gray_scale.csv", header=F, sep=","))
ランク数を「r」とおいて、特異値分解する。
A.svd <- svd(A) u <- A.svd$u[,1:r, drop=F] d <- diag(A.svdA$d)[1:r, 1:r, drop=F] v <- A.svd$v[,1:r, drop=F]
そして、再度、csvに書き戻す。
write.table(round(u %*% d %*% t(v)), str_c("200_320_gray_scale_low_rank_approximate_", r), ".csv"), sep=",", row.names=FALSE, col.names=FALSE)
r = 1
r = 1 (拡大)
r = 10
r = 10 (拡大)
r = 20
フロベニウスノルム(元々の行列)
norm(A, "F") # 26607.34
続けて、近似した行列のフロベニウスノルム
フロベニウスノルム(ランク1、5、10、15、20、25、50、75、100で近似した行列)
norm(A_rank1, "F") # 25696.45 norm(A_rank5, "F") # 26186.34 norm(A_rank10, "F") # 26333.71 norm(A_rank15, "F") # 26412.39 norm(A_rank20, "F") # 26460.39 norm(A_rank25, "F") # 26493.76 norm(A_rank50, "F") # 26570.16 norm(A_rank75, "F") # 26594.95 norm(A_rank100, "F") # 26603.22
近似に使うランク数を大きくするほど、フロベニウスノルムも大きくなっていることが確認できる。
尚、特異値の上位20は以下のとおりであった。
> A.svd$d[1:20] [1] 25695.8853 3754.6842 2374.1272 1782.2953 1588.0870 1466.8724 [7] 1308.7928 1248.0294 1112.3381 1042.8967 983.9771 963.9730 [13] 916.1656 841.3105 837.3214 794.9816 735.2430 722.9311 [19] 685.1815 634.1618
下位20は以下のとおり。
> A.svd$d[180:200] [1] 15.898742 15.417982 14.987853 14.435300 14.174012 13.721447 13.289164 [8] 12.747505 12.436280 11.677728 11.473563 11.297093 10.792612 10.623457 [15] 9.714841 9.523776 9.005887 8.843195 8.200154 8.091183 7.484766
各ランクの特異値(ランク1〜200)
各ランクの特異値(ランク11〜200)
(参考)特異値を全部足すと以下のとおり。
> sum(A.svd$d) [1] 75484.84