simulation-d20-k5-misspecification.ipynb

{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"3.6.3"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This R environment comes with many helpful analytics packages installed\n# It is defined by the kaggle/rstats Docker image: https://github.com/kaggle/docker-rstats\n# For example, here's a helpful package to load\n\nlibrary(tidyverse) # metapackage of all tidyverse packages\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nlist.files(path = \"../input\")\n","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2021-06-09T10:54:14.424025Z","iopub.execute_input":"2021-06-09T10:54:14.427746Z","iopub.status.idle":"2021-06-09T10:54:16.209462Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"install.packages(\"tailDepFun\")","metadata":{"execution":{"iopub.status.busy":"2021-06-09T10:54:16.219526Z","iopub.execute_input":"2021-06-09T10:54:16.346023Z","iopub.status.idle":"2021-06-09T10:54:43.961293Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Code for simulations in Janssen & Wan (2019)\n# Setting: d=20, k=2, model correctly specified\n\nlibrary(skmeans)\nrequire(matrixStats)\nrequire(gtools)\nlibrary(tailDepFun)\nlibrary(transport)\nlibrary(ramify)","metadata":{"execution":{"iopub.status.busy":"2021-06-09T10:54:43.965579Z","iopub.execute_input":"2021-06-09T10:54:43.970101Z","iopub.status.idle":"2021-06-09T10:54:47.462728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# below is the relevant content for the CRPS estimation from Yuen and Stoev, except that all \"set.seed's\" are deleted, since that interfered with the overall simulation\n############# Begin of the CRPS part\n\n## function maxlinear: same as A %*% Z except replaces plus\n## operation by max\nmaxlinear <- function(A,Z){\n  d <- nrow(A)\n  p <- ncol(A)\n  if(class(Z)==\"numeric\"){\n    Z <- matrix(Z,length(Z),1)\n  }\n  if(p != nrow(Z)) stop(\"non-conformable arguments\")\n  n <- ncol(Z)\n  X <- matrix(rowMaxs(A[rep(1:d,n),]*t(Z[,rep(1:n,each = d)])),d,n)\n  return(X)\n}\n############# End of maxlinear","metadata":{"execution":{"iopub.status.busy":"2021-06-09T10:54:47.466245Z","iopub.execute_input":"2021-06-09T10:54:47.468484Z","iopub.status.idle":"2021-06-09T10:54:47.796734Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"################# Start of the k PC method ##########\n#######################\n#single iteration\n#centroids is a k*d matrix with current proposals\nclusterPC_iter=function(data, centroids){\n  k=length(centroids[,1])\n  n=length(data[,1])\n  d=length(data[1,])\n#     data-size(n*d); centroids-size(k*d);  so we transpose centroids\n  M=data%*%t(centroids)\n  #find current value, this value is our criteria corresponds to minimizing the expected c2 dissimilarity:\n  v=mean(apply(M,1,max))\n#     find the index gr where the maximum is obtained\n  gr=argmax(M,rows=T)\n  for (i in 1:k){\n    seldata=data[gr==i,]\n    if (length(seldata)==d) #interpretation problem when just one vector\n      seldata=t(seldata)\n#       sig is a d * d matrix of seldata^2\n    Sig=t(seldata)%*%seldata/n\n    res=eigen(Sig)\n    centroids[i,]=abs(res$vectors[,1]) #use the first eigenvector, the entries of which are necessarily positive\n#      the first eigenvector for each cluster as the new centroid_i\n  }\n  list(centroids,v,gr)\n}\n################\n#pick randomly the initial centers\nclusterPCOnce=function(data,k,tol,startFromMeans=FALSE){\n  val=0\n  n=length(data[,1])\n  if (startFromMeans)\n    centroids=clusterMeans(data,k)  \n  else{\n    centroids=data[sample(1:n,k),]\n    if (k==1)   #make sure it is a matrix\n      centroids=t(as.matrix(centroids))\n  }\n  niter=0\n  repeat{\n    niter=niter+1\n    res=clusterPC_iter(data, centroids)\n    centroids=res[[1]]\n    diff=res[[2]]-val\n    val=res[[2]]\n    cluster=res[[3]]\n    if(diff<tol)\n      break\n  }\n#   print(niter)\n  list(centroids,val,cluster)\n}\n###############\n#iterate nrep times and pick the best\nclusterPC=function(data,k,tol=10^(-5),nrep=100,startFromMeans=FALSE){\n  maxval=0\n  for( i in 1:nrep){\n    res=clusterPCOnce(data,k,tol,startFromMeans && (i==1))\n    if (res[[2]]>maxval){\n      cluster=res[[3]]\n      maxval=res[[2]]\n      centroids=res[[1]]\n    }\n  }\n#   centroids\nlist(centroids,maxval,cluster)\n}","metadata":{"execution":{"iopub.status.busy":"2021-06-09T10:54:47.80347Z","iopub.execute_input":"2021-06-09T10:54:47.808441Z","iopub.status.idle":"2021-06-09T10:54:47.841692Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"############# Start of additional setup, necessary functions\n\n# L2-norm is used below:\nnorm_vec <- function(x) sqrt(sum(x^2))\n\n\n# Bmatrix transforms theta-parameters into a matrix A, makes sure that standardization is met\nBmatrix <- function(m,theta){\n  B<- cbind(m,1-theta)\n  return(B)\n}\n\n############# End of additional setup, necessary functions","metadata":{"execution":{"iopub.status.busy":"2021-06-09T10:54:47.848197Z","iopub.execute_input":"2021-06-09T10:54:47.853322Z","iopub.status.idle":"2021-06-09T10:54:47.883035Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"############ Start of core simulation part\nm<-100\n# m=number of simulations \n# Initialization if evaluation characteristics:\n# sk relates to our method, spherical means\n# dist* relates to d_s (Table 1), wasser* relates to W_1 (Tables 2 and 3)\nwasserkPC2<-rep(0,m)\nwassersk2<-rep(0,m)\nduration_sk2<-rep(0,m)\nduration_kPC2<-rep(0,m)\n\nwasserkPC8<-rep(0,m)\nwassersk8<-rep(0,m)\nduration_sk8<-rep(0,m)\nduration_kPC8<-rep(0,m)\n\nset.seed(10)\n","metadata":{"execution":{"iopub.status.busy":"2021-06-09T10:54:47.889567Z","iopub.execute_input":"2021-06-09T10:54:47.894455Z","iopub.status.idle":"2021-06-09T10:54:47.946624Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# misspecified to k = 2 / 8 when correct k=5\n\nfor(i in 1:m){\nprint(i)\n# Random model initialization, see paper for model description\ntheta1 <- c(runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,\n           runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,\n           runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,\n           runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5)\ntheta2 <- c(runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,\n           0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)\ntheta3 <- c(0,0,0,0,0,\n           runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,\n           0,0,0,0,0,0,0,0,0,0)\ntheta4 <- c(0,0,0,0,0,0,0,0,0,0,\n           runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,runif(1)*0.5,\n           0,0,0,0,0)\n    \ntemp_sum <- theta1 + theta2 + theta3 + theta4\nfirst_m <- cbind(theta1, theta2, theta3, theta4)\nB <- Bmatrix(first_m,temp_sum)\n\n# C are the true points of support (as rows) of the spectral measure and trueprobs the corresponding probabilities\nC<-t(cbind(B[,1]/norm_vec(B[,1]),B[,2]/norm_vec(B[,2]),B[,3]/norm_vec(B[,3]),B[,4]/norm_vec(B[,4]),B[,5]/norm_vec(B[,5])))\ntrueprobs <- apply(B, 2, norm_vec)/sum(apply(B, 2, norm_vec))\n\n\nn<-1000\n# n=sample size per simulation\n# Generate Z vectors and sample data\nZ <- matrix(1/rexp(5*n),5,n)\ndata <- t(maxlinear(B,Z))\n# Transform to unit Pareto marginals via rank transformation for sk-means procedure \nx <- apply(data, 2, function(i) n/(n + 1 - rank(i)))\n\n# For k-means estimator, chose transformed observations with largest L2-norms\nl2norm<-apply(x, 1, norm_vec)\n# we use the largest 100 observations:\nu<-0.90\nex<-x[l2norm>quantile(l2norm,u),]/l2norm[l2norm>quantile(l2norm,u)]\n\n# apply spherical 2-means:\nstart_time_sk <- Sys.time()\nskout<-skmeans(ex,2,control=list(nruns = 100))\nend_time_sk <- Sys.time()\n# estimated points of mass (as rows) and corresponding probabilities\nCsk<-skout$prototypes\nskprobs<-table(skout$cluster)/sum(table(skout$cluster))\n\n# The characteristic W_1 is derived below, with the help of the package 'transport'\nwassersk2[i]<-wasserstein(wpp(C,trueprobs),wpp(Csk,skprobs),p=1)\n           \n# record running time for each specification\nduration_sk2[i]<- round(end_time_sk - start_time_sk,4)\nprint(\"k-means misspecification k = 2\")\nprint(c(wassersk2[i],duration_sk2[i]))\n\n# run the current simulated data with spherical k-PC method          \n# apply kPC 2-means:\nstart_time_pc <- Sys.time()\nkPCout<-clusterPC(ex,2)\nend_time_pc <- Sys.time()\n# estimated points of mass (as rows) and corresponding probabilities\nCkPC<-kPCout[[1]]\nkPCprobs<-table(kPCout[[3]])/sum(table(kPCout[[3]])) \n\n# The characteristic W_1 is derived below, with the help of the package 'transport'\nwasserkPC2[i]<-wasserstein(wpp(C,trueprobs),wpp(CkPC,kPCprobs),p=1)\n\n#record the running time of kPC\nduration_kPC2[i]<- round(end_time_pc - start_time_pc,4)\n\nprint(\"k-Principle Component misspecification k = 2\")\nprint(c(wasserkPC2[i],duration_kPC2[i]))\n           \n# apply spherical 8-means:\nstart_time_sk <- Sys.time()\nskout<-skmeans(ex,8,control=list(nruns = 100))\nend_time_sk <- Sys.time()\n# estimated points of mass (as rows) and corresponding probabilities\nCsk<-skout$prototypes\nskprobs<-table(skout$cluster)/sum(table(skout$cluster))\n\n# The characteristic W_1 is derived below, with the help of the package 'transport'\nwassersk8[i]<-wasserstein(wpp(C,trueprobs),wpp(Csk,skprobs),p=1)\n           \n# record running time for each specification\nduration_sk8[i]<- round(end_time_sk - start_time_sk,4)\nprint(\"k-means misspecification k = 8\")\nprint(c(wassersk8[i],duration_sk8[i]))\n\n# run the current simulated data with spherical k-PC method          \n# apply kPC 2-means:\nstart_time_pc <- Sys.time()\nkPCout<-clusterPC(ex,8)\nend_time_pc <- Sys.time()\n# estimated points of mass (as rows) and corresponding probabilities\nCkPC<-kPCout[[1]]\nkPCprobs<-table(kPCout[[3]])/sum(table(kPCout[[3]])) \n\n# The characteristic W_1 is derived below, with the help of the package 'transport'\nwasserkPC8[i]<-wasserstein(wpp(C,trueprobs),wpp(CkPC,kPCprobs),p=1)\n\n#record the running time of kPC\nduration_kPC8[i]<- round(end_time_pc - start_time_pc,4)\n\nprint(\"k-Principle Component misspecification k = 8\")\nprint(c(wasserkPC8[i],duration_kPC8[i]))\n}","metadata":{"execution":{"iopub.status.busy":"2021-06-09T10:54:47.956269Z","iopub.execute_input":"2021-06-09T10:54:47.965458Z","iopub.status.idle":"2021-06-09T11:15:25.211583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"round(mean(wassersk2),4)\nround(sd(wassersk2),4)\nround(mean(wassersk8),4)\nround(sd(wassersk8),4)","metadata":{"execution":{"iopub.status.busy":"2021-06-09T11:15:25.215291Z","iopub.execute_input":"2021-06-09T11:15:25.217577Z","iopub.status.idle":"2021-06-09T11:15:25.26576Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"round(mean(wasserkPC2),4)\nround(sd(wasserkPC2),4)\nround(mean(wasserkPC8),4)\nround(sd(wasserkPC8),4)","metadata":{"execution":{"iopub.status.busy":"2021-06-09T11:15:25.269464Z","iopub.execute_input":"2021-06-09T11:15:25.271684Z","iopub.status.idle":"2021-06-09T11:15:25.318967Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"round(mean(duration_sk2),4)\nround(sd(duration_sk2),4)\nround(mean(duration_sk8),4)\nround(sd(duration_sk8),4)","metadata":{"execution":{"iopub.status.busy":"2021-06-09T11:15:25.321268Z","iopub.execute_input":"2021-06-09T11:15:25.322617Z","iopub.status.idle":"2021-06-09T11:15:25.350375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"round(mean(duration_kPC2),4)\nround(sd(duration_kPC2),4)\nround(mean(duration_kPC8),4)\nround(sd(duration_kPC8),4)","metadata":{"execution":{"iopub.status.busy":"2021-06-09T11:15:25.352506Z","iopub.execute_input":"2021-06-09T11:15:25.353834Z","iopub.status.idle":"2021-06-09T11:15:25.381915Z"},"trusted":true},"execution_count":null,"outputs":[]}]}