Actually, your code will not generate values from a t-distribution (not quite, at least).
This is because t.test
by default uses the Welch test.
There are two things that prevent your results having a t-distribution:
i) the Welch-statistic only has approximately a t-distribution -- the denominator is only approximately chi-square with the nominal df of the Welch test.
ii) the df are computed as a function the sample values, and thus vary each time you generate a new sample. Thus you'll actually get some mixture of these not-quite-t-distributions-that-are-conditional-on-sample-specific-variance-ratios.
Here's how to actually get a t-distribution; if you want a two-sample t-test, you'll need to use an equal-variance t-test:
n <- 10 # this will get you a t with 18 df
res <- replicate(5000,{x=rnorm(n);y=rnorm(n);(mean(x)-mean(y))/sqrt((var(x)+var(y))/n)})
This statistic is a simplified version (because of equal sample size) of an ordinary equal variance two-sample t-statistic.
Alternatively, if you want to use t.test
:
n <- 10 # this will get you a t with 18 df
res <- replicate(5000,t.test(rnorm(n),rnorm(n),var.equal=TRUE)$statistic)
You can check it gives a t-distribution with the supposed d.f. by transforming the output using the cdf and comparing the result with a uniform:
plot(ecdf(pt(res,2*n-2)))
abline(0,1,col=2,lty=3)
giving
However, it's probably simpler to start with a one-sample test.
Who exactly made the tables of critical values for t distributions found in statistical text books?
Each set of tables - unless it credits someone else as a source - will likely have been constructed from scratch, in a variety of ways depending on when they were published.
Critical values are effectively tables of the inverse cdf, but if you can compute the cdf you can find the value that corresponds to a certain point of the cdf relatively quickly (e.g. via interpolation and say a Newton or secant-method step); there were also approaches like the use of finite difference methods that could also be used to achieve smooth interpolated results in use before people used computers for everything.
Note that for moderate to low d.f. it's quite reasonable to compute integrals directly (you can use integration by parts to write it in terms of an integral with d.f. that is two lower than the original, and the two lowest d.f. yield to standard integration methods. There's also numerical integration, which can be made quite accurate. There are series approximations, and there are highly accurate functional approximations (if I recall correctly, some are found in Abramowitz and Stegun for example). There are particular (highly accurate) functional approximations that are suitable for computer implementation, which is how tables are constructed now.
A typical set of t-tables might contain 100 to 200 values to about 3 or 4 figure accuracy; these could be computed by hand in a fairly short period of time - probably a couple of days at most if you have several people checking the computations (people who did such computations - often on mechanical calculators - were called computers). Perhaps 3 or 4 days with a single person doing the calculations and the checking.
The t-distribution itself dates back a long way; to Helmert in the mid 1870s (who was also responsible for the chi-square). Its use for the t-test is due to Gosset ("Student"), but I think the first tables of the inverse cdf for that purpose would be due to Fisher. For example, he gave t-tables in his 1925 book Statistical Methods for Research Workers (see the bottom of the page).
However, as mentioned by Fisher here (1922)$^{[1]}$, Gosset produced tables of the probability integral (i.e. the cdf rather than the inverse cdf) in 1917$^{[2]}$ (up to 30df); it would be practically possible to do the t-test from those. Indeed, Gosset's 1908 paper $^{[3]}$ has* a table for an integral related to the cdf of the $t$ for a set of small d.f. from which (with a little additional work), it's actually possible to in effect do t-tests (as he does in the paper, though keep in mind the hypothesis testing framework we're used to isn't entirely in place at that time).
*(see p19 here; a link to a scan of the original paper is on the Wikipedia page for the Student's t distribution, see note 9)
[1] Fisher, R.A. (1922),
"The goodness of fit of regression formulae, and the distribution of regression coefficients".
Journal of the Royal Statistical Society 85 (4): 597–612.
[2] Student, (1917),
"Tables for estimating the Probability that the Mean of a unique Sample of Observations lies between $-\infty$ and any given Distance of the Mean of the Population from which the Sample is drawn",
Biometrika 11 (4): 414-417
[3] Student, (1908),
"The probable error of a mean",
Biometrika 6 (1): 1–25.