Bimodal eigenvectors from simple multivariate normal model

[Mike Lawrence]

I generated some fake 3-variable multivariate data via MASS::mvrnorm, and ran it through the following standard-ish multivariate model:

  data {

    # N: number of cases

    int<lower=1> N ;

    # K: outcomes per case

    int<lower=2> K ;

    # Y: matrix of outcomes for each case & variable

    matrix[N,K] Y ;

  }

  parameters {

    # means: vector of population-level means for intercept & slope

    vector[K] means ;

    # sds: vector of population-level sds for intercept & slope

    vector<lower=0>[K] sds ;

    # cor: population-level correlation between intercept & slope

    corr_matrix[K] cor ;

  }

  transformed parameters{

    # cov: covariance matrix (derived from cor & sds)

    matrix[K,K] cov ; #used in getting correlation to work

    #compute covariance matrix from cor & sds

    cov = quad_form_diag(cor,sds) ;

  }

  model {

    #priors on population parameters; presumes data has been scale()ed

    means ~ normal(0,1) ;

    sds ~ weibull(2,1) ;

    cor ~ lkj_corr(1) ;

    #assert sampling of case-level outcomes on each variable given

    #  multivariate normal population

    for(s in 1:N){

      Y[s,] ~ multi_normal(means,cov) ;

    }

  }

I was interested in doing something akin to principal component analysis, so in R I computed the eigenvectors given each iteration's covariance matrix:

  #get the samples of the covariance matrix

  cov = rstan::extract(post,pars='cov')[[1]]

  #define a function to compute eigenvectors on each sample

  # covariance matrix

  get_evec = function(i){

    e = eigen(cov[i,,])

    v = tibble::as_tibble(eigen(z[i,,])$vectors)

    v$variable = 1:nrow(v)

    v %>% tidyr::gather(

      key = 'component'

      , value = 'value'

      , -variable

    ) 

  }

  #compute all the eigenvectors

  evec = purrr::map_df(

    .x = 1:dim(z)[1]

    , .f = get_evec

  )

  #visualize the eigenvectors

  evec %>%

    ggplot(

      mapping = aes(

        x = value

      )

    )+

    geom_density(

      colour = 'transparent'

      , fill = 'black'

      , alpha = .5

    )+

    facet_wrap(

      ~ component * variable

      , scale = 'free'

      , labeller = labeller(component = label_both,variable=label_both)

    )

Which yields this plot that has clear bimodality in a few of the eigenvectors:

Is this expected behaviour? 

Mike

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

[Ben Goodrich]

On Wednesday, November 9, 2016 at 8:40:37 PM UTC-5, Mike Lawrence wrote:

Is this expected behaviour? 

Yes. If X * v = lambda * v, you can let w = -v and X * w = lambda * w. If necessary, one way to proceed is to multiply the j-th eigenvector by -1 iff V[1,j] < 0.

Ben

[Mike Lawrence]

Ah, would that be achieved by doing the following for each sample eigenvector matrix v?

  e = eigen(cov[i,,]) #compute on ith sample for the covariance matrix v = e$vectors #get the eigenvector matrix

  for(j in 1:ncol(v)){

  if(v[1,j]<0){

  v[,j] = -v[,j]

  }

  }

When I do this, I till get a little bit of bimodality ("component: v2; variable: 2" & "component: V2; variable: 3" plots:

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Ben Goodrich]

On Friday, November 11, 2016 at 11:13:52 AM UTC-5, Mike Lawrence wrote:

Ah, would that be achieved by doing the following for each sample eigenvector matrix v?

Yes.
 

When I do this, I till get a little bit of bimodality ("component: v2; variable: 2" & "component: V2; variable: 3" plots:

I don't think genuine bimodality can be ruled out because the eigenvectors are a complicated function of the covariance matrix. Also, weird numerical things can happen if either the eigenvalues get too close together or too close to zero.

Ben