Questions on Interpreting Factor Analysis Results and Scores
Is my approach in generating a distribution based on the samples' factor scores flawed? How do they do it in practice?


I found this somewhat difficult to follow. But in general, you should be able to approximate a set of test scores using a multivariate normal distribution where the covariance matrix implies positive correlations between all tests. Some might be larger and some smaller, but the idea is that all ability tests are correlated. And general mental ability can be estimated as the first unrotated factor that results from such tests.


  Is the log-likelihood of a sample even an appropriate score to use? (If not, what alternative ways are there to score a sample?)


This sounds more like how you evaluate a model. E.g., how you evaluate a factor analytic solutions. 
In general, factor saved scores will be a weighted composite of the scores on the  component tests. 

In R, you can use factanal

factanal(x, factors, data = NULL, covmat = NULL, n.obs = NA,
         subset, na.action, start = NULL,
         scores = c("none", "regression", "Bartlett"),
         rotation = "varimax", control = NULL, ...)


See the scores argument. There are a few different methods.


  I have gone ahead and generated the scores using the score_samples() method for all the samples, but they range between -4 and -49. Is there a reason they would be negative?


I don't know Python. But in general, factor saved scores are typically quantified in such a way that they are z-scores (e.g., mean = 0, sd = 1).


  If you are only looking for 1 latent factor, is it good practice to set the number of factors to 1 or should you leave it unspecified anyways?


You need to either extract only one factor or ensure that you apply no rotation to the extract factors. Without a rotation, the first factor will be equivalent to just one factor. If you rotate, variation will be partitioned across the extracted factors.