1 Circular data ぱらぱらめくる『Directional Statistics』
# Generate 100 observations from a von Mises distribution. # with mean direction 0 and concentration 3. data.vm <- rvonmises(n=1000, mu=circular(3), kappa=0.1) # Shrink the plot so that all points fit. plot(data.vm, stack=TRUE, bins=150, shrink=1.5)
- von Mises Fisher distributionは
という確率密度分布なので、(こちらでRcppを練習していることもあり)、Rcppを使って重みづけサンプリングをして、von Mises Fisher distribution乱数を作ってみよう
#include <RcppArmadilloExtensions/sample.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp ;
// [[Rcpp::export]]
List RvonMisesFisher(int n,NumericVector mu,double k) {
int d = mu.size();
NumericMatrix ret(n,d);
int cnt = 0;
RNGScope scope;
while(cnt<n){
NumericVector tmp(d);
double rsq = 0;
for(int i=0;i<d;++i){
tmp[i] = ::Rf_rnorm(0,1);
rsq += tmp[i]*tmp[i];
}
rsq = sqrt(rsq);
if(rsq==0){
for(int i=0;i<d;++i){
tmp[i] = 0;
}
}else{
for(int i=0;i<d;++i){
tmp[i] /= rsq;
}
}
double tmp2 = 0;
for(int i=0;i<d;++i){
tmp2 += mu[i]*tmp[i];
}
NumericVector tmpr = runif(1);
if(exp(tmp2*k)/exp(k) > tmpr[0]){
for(int i=0;i<d;++i){
ret(cnt,i) = tmp[i];
}
cnt ++;
}
}
return List::create(
Named("") = ret
);
}d <-2 mu <- c(1,rep(0,d-1)) k <- 10 n <- 1000 out <- RvonMisesFisher(n,mu,k) tmp <- out[[1]] tmp <- rbind(tmp,rep(1,d)) tmp <- rbind(tmp,rep(-1,d)) plot(tmp) d <-3 mu <- c(1,rep(0,d-1)) k <- 10 n <- 1000 out <- RvonMisesFisher(n,mu,k) library(rgl) tmp <- out[[1]] tmp <- rbind(tmp,rep(1,d)) tmp <- rbind(tmp,rep(-1,d)) plot3d(tmp)
*[方向統計学][Directional Statistics][R][circular][ぱらぱらめくるシリーズ]2 要約統計量 ぱらぱらめくる『Directional Statistics』 -中心からの距離と方向ベクトルに分解する >|cpp| #include <RcppArmadilloExtensions/sample.h> // [[Rcpp::depends(RcppArmadillo)]] using namespace Rcpp ; // [[Rcpp::export]] List circCoords(NumericMatrix x) { int n = x.nrow(), k = x.ncol(); NumericVector r(n); NumericMatrix v(n,k); for(int i=0;i<n;++i){ for(int j=0;j<k;++j){ r[i] += x(i,j)*x(i,j); } r[i] = sqrt(r[i]); } for(int i=0;i<n;++i){ if(r[i] > 0){ for(int j=0;j<n;j++){ v(i,j) = x(i,j)/r[i]; } }else{ for(int j=0;j<n;j++){ v(i,j) = 0; } r[i] = 0; } } return List::create( Named("r") = r, Named("v") = v ); }
x <- matrix(rnorm(20),5,4) circCoords(x) */
> circCoords(x) $r [1] 2.912546 3.375059 1.299661 1.943983 1.778478 $v [,1] [,2] [,3] [,4] [1,] -0.795953396 -0.2943604 0.52891325 -0.007805327 [2,] 0.003528010 0.9706041 -0.12037580 0.208386263 [3,] -0.007654404 -0.7251910 0.38246179 -0.572505330 [4,] 0.986174603 0.1533123 0.03668821 -0.051078015 [5,] 0.157526272 0.2688395 -0.58064677 -0.752170286 > circCoords(x) -> cC > apply(cC$v^2,1,sum) [1] 1 1 1 1 1
- 位置平均
- 距離の平均、方向ベクトルの平均
#include <RcppArmadilloExtensions/sample.h> // [[Rcpp::depends(RcppArmadillo)]] using namespace Rcpp ; // [[Rcpp::export]] List meanCirc(NumericMatrix x) { int n = x.nrow(), k = x.ncol(); NumericVector m(k); for(int i=0;i<n;++i){ for(int j=0;j<k;++j){ m[j] += x(i,j)/n; } } double r=0; for(int i=0;i<k;++i){ r += m[i]*m[i]; } r = sqrt(r); return List::create( Named("m") = m, Named("r") = r ); }
d <-3 mu <- c(1,rep(0,d-1)) k <- 10 n <- 1000 out <- RvonMisesFisher(n,mu,k) meanCirc(out[[1]])
> meanCirc(out[[1]]) $mean.vector [1] 0.896681731 -0.006023112 0.001850635 $mean.resultant.length [1] 0.8967039 $circular.variance [1] 0.1032961
- ばらつきの評価
- (任意次元)球面上のベクトルの重心ベクトルのノルム
について
をcircular variance, Vと呼ぶ
- このVは次のようにしても計算できる
// [[Rcpp::export]] List varianceCirc(NumericMatrix x) { int n = x.nrow(); int k = x.ncol(); NumericVector m(k); for(int i=0;i<n;++i){ for(int j=0;j<k;++j){ m[j] += x(i,j)/n; } } double r=0; for(int i=0;i<k;++i){ r += m[i]*m[i]; } r = sqrt(r); double v = 1-r; double v2 = 0; if(r!=0){ for(int i=0;i<n;++i){ double tmp2 = 0; for(int j=0;j<k;++j){ tmp2 += m[j]*x(i,j)/r; } v2 += (1-tmp2)/n; } } double v3 = 0; for(int i=0;i<n;++i){ for(int j=0;j<n;++j){ double tmp2 = 0; for(int j2=0;j2<k;++j2){ tmp2 += x(i,j2)*x(j,j2); } v3 += (1-tmp2)/n; } } return List::create( Named("mean.vector") = m, Named("mean.resultant.length") = r, Named("circular.variance") = v, Named("circular.variance2") = v2, Named("circular.variance3") = v3 ); }
-
- von Mises-Fisherの係数kとばらつきの指標との関係を作成した関数を用いて調べてみる
library(Rcpp) library(RcppArmadillo) sourceCpp("circStat.cpp") n.iter <- 100 n <- 1000 d <- 3 mu <- c(1,rep(0,d-1)) ks <- seq(from=0,to=100,length=n.iter) vs <- rep(0,n.iter) vs2 <- vs3 <- vs rs <- vs for(i in 1:n.iter){ k <- ks[i] tmp <- RvonMisesFisher(n,mu,k) tmp.m <- varianceCirc(tmp[[1]]) rs[i] <- tmp.m[[2]] vs[i] <- tmp.m[[3]] vs2[i] <- tmp.m[[4]] vs3[i] <- tmp.m[[5]] } plot(ks,vs) plot(vs,vs2) plot(vs2,vs3) plot(1-rs^2,vs3)
- Variance matrix
#include <RcppArmadilloExtensions/sample.h> // [[Rcpp::depends(RcppArmadillo)]] using namespace Rcpp ; // [[Rcpp::export]] List varianceMatrix(NumericMatrix x) { int n = x.nrow(); int k = x.ncol(); NumericVector m(k); for(int i=0;i<n;++i){ for(int j=0;j<k;++j){ m[j] += x(i,j)/n; } } double r=0; for(int i=0;i<k;++i){ r += m[i]*m[i]; } r = sqrt(r); NumericMatrix s(k,k); for(int i=0;i<n;++i){ for(int i2=0;i2<k;++i2){ for(int i3=0;i3<k;++i3){ s(i2,i3) += (x(i,i2)-m[i2])*(x(i,i3)-m[i3])/n; } } } double trs = 0; for(int i=0;i<k;++i){ trs += s(i,i); } double R = 1-trs; return List::create( Named("mean.vector") = m, Named("mean.resultant.length") = r, Named("R-squared") = R, Named("traceS") = trs, Named("variance.matrix") = s ); }
d <-3 mu <- c(1,rep(0,d-1)) k <- 10 n <- 1000 out <- RvonMisesFisher(n,mu,k) vM.out <- varianceMatrix(out[[1]]) vM.out vM.out[[2]]^2 vM.out[[3]]
> d <-3 > mu <- c(1,rep(0,d-1)) > k <- 10 > n <- 1000 > out <- RvonMisesFisher(n,mu,k) > > vM.out <- varianceMatrix(out[[1]]) > > vM.out $mean.vector [1] 0.898436769 -0.009032775 0.009338527 $mean.resultant.length [1] 0.8985307 $`R-squared` [1] 0.8073574 $traceS [1] 0.1926426 $variance.matrix [,1] [,2] [,3] [1,] 0.009088391 -0.001069172 -0.001199073 [2,] -0.001069172 0.089571139 -0.002792975 [3,] -0.001199073 -0.002792975 0.093983042 > vM.out[[2]]^2 [1] 0.8073574 > vM.out[[3]] [1] 0.8073574
