Significance/confidence intervals for PCA or factor loadings - how can such be defined? Current discussions here in SSE made me to reconsider the PCA and FA models and procedures. I got curious how one would determine confidence intervals for the components/factor loadings by assuming the structure in your samples correlations is a predictor for the populations correlations and even more:        
if you assume that in the population with some measured items there is a true factor structure latent which reflects varimax-, pc-, paf-models (or any theoretical structure in the loadings) and in your sample you try to reproduce that models/structure - can one assume some confidence intervals for the estimates for the population?                 
What I've tried was so far to find such confidence intervals by bootstrapping; for instance assume some synthetical data of four items for the population and a model of pc/factor loadings, then draw samples with smaller N and look how the pc/factor loadings reflect the "true" (population's) loadings and getting from this confidence intervals.              
I do not yet know whether such an approach is meaningful at all, but just produced some results for models of pc, varimax-rotated pc, "little jiffy" (see end of posting), and paf with 2 factors. I took a set of S=256 samples each of N=40 and another set of N=160 each. Here are a couple of scatter plots of the pc-loadings of the first two pc's the same for the varimax and for the paf models.     
Q: is such an approach meaningful at all? And if, how could I proceed? (The confidence intervals seem to be either elliptical or circular or twovariate normal or ???) 
Remark: while I'm writing this I recall, that an answer likely exists in SEM-methodology, but until now I didn't look in this deeper because of the heavy machinery in some papers which I've tried to read...  So I tried for the beginning some intuition by my described experiments. 

Here are some pictures from my website where I have some more images and more  explanations of the procedure and of observations.
The colours of the clouds indicate the four different items; the white circles in the clouds indicate the factor loadings in the population (according to the referred modeling) 

model pc , axes (ax1,ax2), n=40

model pc , axes (ax3,ax4), n=40

model pc , axes (ax1,ax2), n=160

model varimax , axes (ax1,ax2), n=160

model paf (2 factor) , axes (ax1,ax2), n=40

model paf (2 factors) , axes (ax3,ax4) = (pc1,pc2) of the residual covariance, n=40
 

Here is an excerpt of an article of H.Kaiser describing how the term "little jiffy" was coined:

 A: I am not sure what exactly the question here is. I think the approach is valid, and I recall recently upvoting your answer in How to identify variables with significant loadings in PCA?, where you essentially suggested this bootstrapping scheme.

The confidence intervals seem to be either elliptical or circular or twovariate normal or ???

Perhaps you are asking about the "weird" shape of the bootstrapped clouds that you get. I repeated your analysis on the Fisher Iris dataset: I do PCA on the full data (once on covariance, once on correlation), and then bootstrap it 1000 times (note that one should take care that all bootstrapped principal axes are aligned similarly and not flipped). In each case, I plot the standardized PC1-2 scores of the full data as black dots, full loadings as coloured lines, and 1000 resampled loadings as coloured dots:

The clouds on the left look okay, but the clouds on the right look obviously circular, similar to the ones on your plots. This is easy to understand. Recall that the length of the loading vector approximates the variance of the original variable. When I do the analysis on the covariance matrix, the variance can be anything; but when I do the analysis on the correlation matrix I essentially fix the individual variances to $1$. Hence the length of the loading vectors is also approximately $1$, and the jitter happens only along its direction, not the length.
Nevertheless, I think the procedure makes sense, and you can perfectly well compute the boostrapped confidence intervals e.g. on the first principal axis. For covariance, I get  $$(0.36\pm 0.016, -0.085\pm  0.017, 0.86\pm    0.0045, 0.36\pm   0.0081).$$

Matlab code to produce this figure
clear all
close all
load fisheriris

n = size(meas,1); %// n = 150
varNames = {'Sepal length', 'Sepal width', 'Petal length', 'Petal width'};
col = [0 0 1; 1 0 1; 1 0 0; 0 0.6 0];
modes = {'PCA on covaraince', 'PCA on correlation'};

figure('Position', [100 100 1400 500])

%// covariance or correlation matrix
for mode = 1:2

    %// PCA of the full dataset
    X = bsxfun(@minus, meas, mean(meas));   %// centering
    if mode == 2
        X = bsxfun(@times, X, 1./std(X));   %// scaling
    end
    [U0,S0,V0] = svd(X,0);      %// SVD
    L0 = V0*S0/sqrt(n-1);       %// loadings

    %// boostrapping (1000 repetitions)
    for rep=1:1000
        ind = randi(n, n, 1);

        %// generating a bootstrapped dataset and SVDing it
        X = bsxfun(@minus, meas(ind,:), mean(meas(ind,:)));
        if mode == 2
            X = bsxfun(@times, X, 1./std(X));
        end
        [U,S,V] = svd(X,0);

        %// adjust the signs of components to align them with SVD of the full
        %// dataset (otherwise signs can flip)
        for i=1:size(U,2)
            if sum((V(:,i)-V0(:,i)).^2) > sum((V(:,i)+V0(:,i)).^2)
                V(:,i) = -V(:,i);
                U(:,i) = -U(:,i);
            end
        end

        %// saving loadings and principal axes of the first two PCs
        loadings(rep,:) = [V(:,1)'*S(1,1) V(:,2)'*S(2,2)] / sqrt(n-1);
        axes(rep,:) = [V(:,1)' V(:,2)'];
    end

    %// display principal axes of the full data, mean of the bootstrapped
    %// axes, standard deviations
    num2str([V0(:,1) mean(axes(:,1:4))' std(axes(:,1:4))'], 2)

    subplot(1,2,mode)
    hold on
    %// plotting standardized scores
    scatter(U0(:,1)*sqrt(n-1), U0(:,2)*sqrt(n-1), [], 'k')
    %// plotting loadings
    for i=1:4
        scatter(loadings(:,i), loadings(:,i+4), [], col(i,:))

        plot([0 L0(i,1)], [0 L0(i,2)], 'Color', col(i,:))
        pos = L0(i,1:2) + L0(i,1:2)/norm(L0(i,1:2)) * 0.4;
        text(pos(1), pos(2), varNames{i}, 'Color', col(i,:));
    end
    rectangle('Position', [-1 -1 2 2], 'Curvature', [1 1])
    axis([-1 1 -1 1] * 3)
    axis square
    title(modes{mode})
    xlabel('PC1')
    ylabel('PC2')
end

