# Simulate Godfrey Thomson's "sampling model" of mental abilities, and perform
# factor analysis on the resulting test scores.


# MASS package used to make random multivariate normal vectors
require(MASS)
# stats package used to do factor analysis, but if it's missing you've got big
# problems!
require(stats)

# Set the mean and standard deviation of each elementary ability to that of an
# IQ test - not actually important!
ability.mean = 100
ability.sd = 15

thomson.model <- function(n,d,a,q,output="variance") {
    # Using incomprehensible parameter names is bad
    # number of testees = n
    # number of tests = d
    # number of shared abilities = a
    # max. number of specific abilities per test = q

    # assign abilities to tests
    general.per.test = floor(runif(d,min=0,max=a)) + 1
    specifics.per.test = floor(runif(d,min=0,max=q)) + 1

    # Define the matrix assigning abilities to tests
    general.to.tests = matrix(0,a,d)
    # The use of a for loop here is maybe a little slower than some sort
    # of vectorization, but it's sanity-preserving; so is using the temporary
    # variable "abilties" to hold the sample.
    for (i in 1:d) {
        abilities = sample(1:a,size=general.per.test[i],replace=FALSE)
        general.to.tests[abilities,i] = 1
    }

    # Covariance matrix of the general abilities
    sigma = matrix(0,a,a)
    diag(sigma) = (ability.sd)^2
    mu=rep(ability.mean,a)
    x = mvrnorm(n,mu,sigma) # person-by-abilities matrix of abilities

    # The "general" part of the tests
    general.tests = x %*% general.to.tests

    specific.tests = matrix(0,n,d)
    noisy.tests = matrix(0,n,d)
    # Each test gets its own specific abilities, which are independent for each
    # person
    # Again, I could probably vectorize, but d is small and this is saner
    for (i in 1:d) {
        # Each test has noises.per.test disturbances, each of which has s.d. 15,
        # since these are all independent their variances add
        j = specifics.per.test[i]
        specifics = rnorm(n,mean=100*j,sd=ability.sd*sqrt(j))
        specific.tests[,i] = general.tests[,i] + specifics
        # Finally, for extra realism, some mean-zero trial-to-trial noise, so
        # that if we re-use this combination of general and specific ability
        # scores, we won't get the exact same test scores twice
        noises = rnorm(n,mean=0,sd=ability.sd)
        noisy.tests[,i] = specific.tests[,i] + noises
    }

    # do the factor analysis
    my.fa = factanal(noisy.tests,1)
    # extract the loadings
    v = as.vector(my.fa$loadings)
    # compute how much of the variance (really, correlation) the factor
    # describes
    varexp = sum(v^2)/d
    # output: should do a return() instead to be completely pukka
    switch(output,
        variance = varexp,
        everything = list(general.ability.pattern = general.to.tests,
                          numbers.of.specifics = specifics.per.test,
                          ability.matrix = x, specific.tests = specific.tests,
                          data = noisy.tests, loadings = v,
                          varexp = varexp),
        "I can't tell what you want me to output")
}