R : Copyright 2004, The R Foundation for Statistical Computing
Version 1.9.1  (2004-06-21), ISBN 3-900051-00-3

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

[Previously saved workspace restored]

> nseq1
  [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
 [18]  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34
 [35]  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51
 [52]  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68
 [69]  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85
 [86]  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100

## this was assigned earlier, using nseq <- 1:100

> stirling
function(n){
sqrt(2*pi)*n^(n+.5)*exp(-n)}

## stirling was also made earlier; crude approximation to n!

> plot(nseq1,log(factorial(nseq1)))

## I plotted the log because the function blows up too fast

> lines(nseq1,log(stirling(nseq1)),col="red",lwd=2)

## this adds lines to an open plot, 'lwd' stands for 'line width'

> nseq <- 100: 1000  # now try really large n

> plot(nseq,log(factorial(nseq)))
> lines(nseq,log(stirling(nseq)),col="red",lwd=2)

## this plot bombs our around n=200

> logstirling
function(n){
(n+.5)*log(n)-n+.5*sqrt(2*pi)}


## now I'll compute it on a more sensible scale

> plot(nseq,logstirling(nseq))
> lines(nseq,log(factorial(nseq)),col="red",lwd=2)

## Note that using the built in function 'factorial' is still blowing up

## In fact 'factorial(n)' is actually computed using 'gamma(n+1)' 
## and there is also a function called 'lgamma' which returns the log

> lines(nseq,lgamma(nseq+1),col="green")  # this is perfectly sensible

> save.image(file="/Users/nancy/teaching/410/jan6/stirling.rda")

> q()
Save workspace image? [y/n/c]: n