#  Stochastic Calculus, Courant Institute, NYU, Fall 2014
#  http://www.math.nyu.edu/faculty/goodman/teaching/StochCalc2014/index.html

#  (your name and email here)

#  Partial code for assignment 1

#  Explore the relationship between Gaussian and t densities
#  Compare the variances, with the same sigma, plot the density functions

#-------------  procedures to calculate the un-normalized densities ---------------------

#  return the un-normalized Gaussian pdf, e^{-(x-mu)^2/sig^2}

g <- function( x, mu, sig){
  z = (x - mu)/sig                  # distance of x from mean, in units of sigma
  un_normalized_pdf = exp( -z*z/2)  # you can use a long variable name for things not used in formulas.
  return( un_normalized_pdf )
 }

#   return the un-normalized t pdf, with n degrees of freedom

f <- function( x, mu, sig, n){
  z = (x - mu)/sig
  pow = .5*(n+1)
  un_normalized_pdf = 1/( ( 1 + z*z/n)^pow )
  return( un_normalized_pdf )
 }

#------ The main program, compute the normalization constants, the variance of t, and make plots


cat("Welcome to R world\n")        # output to the terminal

mu    = 0.                         # use these mean and "variance" parameters for ...
sig   = 2.                         # ... the whole computation
n     = 5.                         # number of degrees of freedom for t

output_string = sprintf( "Computing with mu = %8.3f, sigma = %8.3f, and n = %8.3f\n", mu, sig, n)
cat( output_string)

#   Compute integrals numerically (approximately) using the trapezoid rule and equally spaced points

x_min_int = mu - 8.*sig                    #  endpoints for the integration interval.
x_max_int = mu + 8.*sig                    #  should be large enough to include "all" the mass in the integrals.
nx        = 1000                           #  the number of integration points
dx        = (x_max_int - x_min_int)/(nx-1) #  the distance between integration points

#    The Gaussian integral, for the normalization constant

g_int = 0.                               # initialize the variable that will hold the sum g(x_i)*dx
for (i in (1:nx)){                       # this makes i go from 1 to nx
  x     = x_min_int + (i-1)* dx          # the location of the integration point x_i.  i=1 makes x_i = x_min
  g_int = g_int + g( x, mu, sig)*dx      # add the contribution for integration point x_i
 }

#    The Gaussian variance integral, using the normalized PDF

g_var = 0.                               # will be the integral that defines the variance
for (i in (1:nx)){
  x     = x_min_int + (i-1)* dx
  g_var = g_var + ( g( x, mu, sig)/g_int)*dx # add the contribution for integration point x_i
 }

#     Print the numerical answers and the errors, which you know for the Gaussian

twopi         = 2*3.14159
g_int_exact   = sqrt( twopi*sig*sig)   # the exact answer, by formula
g_int_error   = g_int - g_int_exact
g_var_exact   = sig*sig                # the exact variance, by formula
g_var_error   = g_var - g_var_exact
output_string = sprintf( "numerical Gaussian integral = %8.4f, error = %9.3e\n", g_int, g_int_error)
cat( output_string)
output_string = sprintf( "numerical Gaussian variance = %8.4f, error = %9.3e\n", g_var, g_var_error)
cat( output_string)

#    The t integral, for the normalization constant, more or less the same as the gaussian integral

t_int = 0.
for (i in (1:nx)){
  x     = x_min_int + (i-1)* dx
  t_int = t_int + f( x, mu, sig, n)*dx
 }
cat( sprintf("t integral = %8.3f\n", t_int))

#    The t variance integral, using the normalized PDF

t_var = 0.                               # will be the integral that defines the variance
for (i in (1:nx)){
  x     = x_min_int + (i-1)* dx
  t_var = t_var + ( f( x, mu, sig, n)/t_int)*dx # add the contribution for integration point x_i
 }

output_string = sprintf( "Gaussian variance = %8.3f, t variance = %8.3f, t var - Gaussian var = %9.3e\n", g_var_exact, t_var, t_var - g_var_exact)
cat( output_string)

nd = ( t_int - g_int)*n
cat( sprintf("normalized difference is = %8.3f\n", nd))

x_plot_max = mu + 5.*sig                        # upper and lower limits of the plot interval, smaller than the compute interval
x_plot_min = mu - 5.*sig
nx_plot    = 200                                # number of points to plot, for each curve
dx_plot = (x_plot_max-x_plot_min)/(nx_plot-1)
g_vals = rep(0, nx_plot)                        # the gaussian values, etc.
t_vals = rep(0, nx_plot)
x_vals = rep(0, nx_plot)

for ( i in (1:nx_plot)){
  x = x_plot_min + dx_plot*(i-1)
  g_vals[i] = g( x, mu, sig)/g_int
  t_vals[i] = f( x, mu, sig, n)/t_int
  x_vals[i] = x
 }

#-------- create the plot and write it to a .pdf file -----------------

# pdf( "PDF_plots.pdf")                           # PDF is for probability density function

# title = sprintf( title and parameter values, mu, sig, n)

# plot(  the gaussian plot, with some kind of line )


# lines( the t density plot, with some other kind of lines )   #  a blue line

#  grid()    so you can read the plot

#  labels = c(two labels, for the legend)

#  legend(... ) #  characters.  NA means "none"

#  dev.off( which = dev.cur() )