{"cell_type":"code","exec_count":0,"id":"39a91d","input":"Nsample <- 500\nlibrary(boot)\t\t\t\t\t# boot is one of ~30 packages pre-installed with R\n\t\t\t\t\t \t\t# library() shows all CRAN packages currently installed in your computer\nmedian.stat <- function(x,i) median(x[i]) # you can bootstrap a simple or complicated statistic of the data\nyboot <- boot(y, median.stat, R=500)\t # calculate the bootstrap resamples\nstr(yboot)\t\t\t\t\t\t# lots of information in the bootstrap output! See help(boot) for details.\nyboot\t\t\t\t\t\t\t# +/- 0.8 (bootstrap estimate of standard error of the median)\n\nboot.ci(yboot, conf=0.95, type=c('basic', 'norm', 'perc')) \t# three estimates of the 95% confidence interval\n\nhist(yboot$t, breaks=50)\t\t# note that the distribution of medians is not smooth, even for large Nsample","pos":9,"state":"done","type":"cell"}
{"cell_type":"code","exec_count":0,"id":"5d22f5","input":"# We calculate the ages of the universe for three different scenarios:\n# the benchmark model, a matter-dominated universe, and a cosmological constant dominated universe\nbenchmark = timeUniverse(z, omega_m, omega_lambda)\nmatterDom = timeUniverse(z, 1, 0)\ncosmoConstantDom = timeUniverse(z, 0.01, 0.99)","pos":5,"state":"done","type":"cell"}
{"cell_type":"code","exec_count":0,"id":"5df721","input":"# We define a function for the equation described above\n# The function has three variables: redshift, the matter density parameter, \n# and the cosmological constant density parameter\n\ntimeUniverse = function(redshift, omegaM, omegaLambda) {\n\n # We define the integrand within the equation\n integrand = function(redshift) {1/((1+redshift) * sqrt(omegaM*(1+redshift)^3 + omegaLambda))}\n\n # Using the integrand, we can integrate this function for each redshift given\n # We then define a function that fully solves the equation above\n # We also add in an extra factor in order to put the time in units of Gigayears\n ageUniverse = function(redshift) { \n age = (1/H0_convert) * integrate(integrand, redshift, Inf)$value / (10^9 * 31557600)\n return(age) }\n \n # We calculate the ages associated with each redshift and return this vector of values\n ageUni = Vectorize(ageUniverse)(z)\n \n return(ageUni)\n}","pos":4,"state":"done","type":"cell"}
{"cell_type":"code","exec_count":0,"id":"6d42ca","input":"# We plot the three different scenarios\nplot(z, benchmark, type='l', lty=1, lwd=2, xlab='z (redshift)', ylab='time (Gyr)', ylim=c(0,30))\nlines(z, matterDom, lty=2, lwd=2)\nlines(z, cosmoConstantDom, lty=3, lwd=2)\n\n# We create a legend with labels for each different value of the density parameters\nlegend(5.0, 30, lty=c(1,2,3), lwd=c(2,2,2), legend=c(expression(paste(Omega[m], \n \"= 0.31 \", Omega[Lambda], \"= 0.69 \")), expression(paste(Omega[m], \"= 1 \", \n Omega[Lambda], \"= 0 \")), expression(paste(Omega[m], \"= 0.01 \", \n Omega[Lambda], \"= 0.99 \"))))\n\ntitle('Age of the Universe')","pos":6,"state":"done","type":"cell"}
{"cell_type":"code","exec_count":0,"id":"aacb4e","input":"x <- seq(1,10,by=0.01) \t\t# 500 evenly-spaced numbers between 1 and 10\ny <- x^2 - 4*x*sin(x) \t\t\t# sample with weird distribution\nplot(x,y, type='l') \t\t\t# look at the sample … I guess the median is around 30\nmedian(y) \t\t\t\t# the median is actually 33.0\n1.48 * mad(y) / sqrt((length(y) -1))\t# +/- 0.9 (MAD estimate of standard error of the median)","pos":8,"state":"done","type":"cell"}
{"cell_type":"code","exec_count":0,"id":"ba1caa","input":"# First, we generate a set of redshifts from 0 to 10\nz = seq(0.0, 10, 0.1)\n\n# We define some helpful constants\nH0 = 68 # km/s/Mpc\nkmtoMpc = 3.08559*10^19 # km/Mpc\nH0_convert = H0/kmtoMpc # s^-1\n\n# We define the benchmark values that are the accepted current values of the matter density parameter, \n# and the cosmological constant density parameter\nomega_m = 0.31\nomega_lambda = 0.69","pos":3,"state":"done","type":"cell"}
{"cell_type":"markdown","id":"3df1bd","input":"> **Exercise 1.** Practice with some elementary R functions: arithmetic and algebra, 2-dimensional array manipulation, producing multi-element lists. Write a brief program exercising program flow control (if, ifelse, when, repeat, break, stop, etc.)","pos":1,"state":"done","type":"cell"}
{"cell_type":"markdown","id":"4835e7","input":"> **Exercise 2.**\nEstimate the age of the Universe as a function of redshift for a standard $\\Lambda$CDM universe model: \n$$ t(z) = H_0^{-1} \\int_{z}^{\\infty}\\frac{dz'}{(1+z')h(z')} $$ \nwhere $ h(z) = \\sqrt{(1-\\Omega_{total})(1+z)^2 + \\Omega_m(1+z)^3 + \\Omega_{\\Lambda} } $, $\\Omega_m$ is the matter density parameter, and $\\Omega_{\\Lambda}$ is the dark energy density parameter. Plot the age of the Universe vs. redshift ($z=0$ to 10) for three hypothetical universes: matter-dominated ($\\Omega_m=1.0$ and $\\Omega_{\\Lambda}=0.0$), dark-energy-dominated ($\\Omega_m=0.01$ and $\\Omega_{\\Lambda}=0.99$), and a realistic universe ($\\Omega_m=0.31$ and $\\Omega_{\\Lambda}=0.69$). This problem and solution is courtesy of graduate student Phoebe Sandhaus, Penn State.","pos":2,"state":"done","type":"cell"}
{"cell_type":"markdown","id":"643c93","input":"> **Exercise 3.** Use bootstrap resampling (random sampling with replacement) to estimate uncertainties of a statistic. Create a univariate sample with a weird distribution … maybe sampling from a polynomial or nonlinear function over some interval. First, calculate the median and a robust measure of its standard error: $1.48*(MAD$ where MAD is the median absolute deviation and the 1.48 scales it to the standard deviation for a Gaussian distribution. Second, estimate the uncertainty of the median from a bootstrap resampling. Give the standard error, 95% confidence interval, and plot a histogram of the bootstrap medians.","pos":7,"state":"done","type":"cell"}
{"cell_type":"markdown","id":"bcecd9","input":"## Exercises for Tutorial 1 Jupyter notebook\n### Summer School for Astronomers 2021\n### Eric Feigelson","pos":0,"state":"done","type":"cell"}
