---
title: "Numeric Integration and the Expanding Cosmos"
author: Jed Rembold
date: April 17, 2025
slideNumber: true
theme: tokyo-night-light
highlightjs-theme: tokyo-night-light
width: 1920
height: 1080
transition: slide
p5js: true
---


## Announcements
- Debriefing for HW5 is available until midnight tonight! Get it done!
- Check in for HW6 is available this weekend. At least have touched base with your partner.
- Poll link for partner preferences on the final project went out (or available [here]()). Complete it by tonight if you want your preferences expressed!
- Homework 6
  - Partner assignments at the end of class
  - Some edits being made to second problem (for the better), but should come out today
- Quiz 2 handed back!
  - You can add 1.5 pts to the total written on your last page

## Recap
- Rotation or velocity curves measure the speed of orbiting objects as you move outwards from a central body
  - Commonly assume circular or largely circular orbits
  - Generally measured through Doppler shifts
- The shape of a rotation curve tells you about the distribution of material throughout the region
- The rotation curves of galaxies do not match what we would expect, thus leading to conjectures about dark matter

## Today
- Understanding numeric integration
  - Why do we care?
  - How can we leverage existing tools?
- Hubble's law and implications


# Numeric Integration
## The Case for Integration
- Models can often require a calculation of all the area beneath a curve
  - Generally this equates to summing up all the contributions over the independent axis
- Examples:
  - The overall brightness of a star is the area under the Planck curve: the brightness at each wavelength summed up
  - The overall mass of something is the area under the mass curve: the mass at each distance summed up
- If one of these things is constant, this calculation becomes easy, as it is just the area of a rectangle
  - But they usually aren't constant!


## Analytic Integration
::::::cols
::::col
- Mathematically, this is the domain of calculus
- Can use integration to directly compute these areas ...sometimes
- The real world is often messy in a way that nice and neat functions almost never are 
- So often an analytic solution to the integration simply doesn't exist, or it is simply so hard to calculate that it may not be worth the time
::::

::::col
$$\begin{align}
f(x) &= 4x^2 + 3 \\
\int_0^5 f(x)\,dx &= \int_0^5 4x^2 + 3\,dx \\
&= \left.\frac{4}{3}x^3 + 3x\right|_0^5 \\
&= \left(\frac{4}{3}5^3 + 3*5\right) - \left(\frac{4}{3}0^3 + 3*0\right) \\
&= 181\tfrac{2}{3}
\end{align}$$
::::
::::::

## Numeric Integration
- Numeric integration is the task of computationally computing the area under these curves
- We actually already visited some problems it has when we introduced MCMC, but for lower dimensional integration problems, we will be ok with classic methods
- These methods generally involve dividing the data into small bins, approximating the area under each bin, and then summing them all up
- While you **could** write a function to do this yourself, for this class you should stick with existing defined tools to accomplish this, as they will be much more accurate


## In Python: `quad`
- The [`quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) function from `scipy.integrate` is built to handle this sort of numeric integration
- Takes 3 main arguments:
  - The function to integrate, where the dimension of integration must be the first parameter of the function
  - A lower bound to the integration
  - An upper bound to the integration
- It will return to you a list of values, of which the first is the integration total value
  - The second is an estimate in the error of that computed value


## In R: `integrate`
- The aptly named [`integrate`](https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/integrate) function in R functions in basically an identical way
- Still takes 3 main arguments:
  - The function to integrate, where the dimension of integration must be the first parameter of the function
  - A lower bound to the integration
  - An upper bound to the integration
- Returns an associative array, so if you just want the integration total, access it with `$value`


## Integration Example
- Suppose we wanted to determine the overall brightness of a Plank curve with emissivity of 1 and a temperature of 8000K
- Reminder:
  $$ B(\lambda, T) = \varepsilon\cdot\frac{2hc^2}{\lambda^5} \frac{1}{\exp\left(\frac{hc}{\lambda k_B T}\right) - 1} $$

<!--
## Emcee
- While it isn't difficult to create our own MCMC sampler, there can be some benefits to using existing libraries for more complicated use cases
  - Tend to give more intuitive interfaces through which to accomplish the basic tasks that we would want
- In the Python landscape, the main option used is the `emcee` package, which you will probably need to install through `pip`
- `Emcee` operates as an abstract data type, wherein you create a sampling object and then can interact with it and run samples using defined methods
- In the R landscape, there is the `mcmc` package, which also seems reasonably strong, though it doesn't seem to have all of the flexibility of `emcee`


## Emcee Creation
- When you create the object initially, you need to provide:
  - The number of walkers you'd like to generate
  - The number of dimensions you'll be walking (number of parameters)
  - The function you'll be walking (your log posterior)
  - Any extra arguments that your function will need beyond the parameters
  - A pool should you be using multiprocessing
```python
sampler = emcee.EnsembleSampler(
                num_walkers, num_dims,
                log_function, args=[extra arguments]
                )
```

## Emcee Running
- You can then start a sampling run by telling the sampler where all the walkers should begin and how many steps they should take
- Generating starting points usually done with some variation of a random gaussian near a starting point:
  ```python
  starts = np.random.multivariate_normal(
              mean = [0,1,10],
              cov = [[1,0,0],[0,0.5,0], [0,0,5]],
              size = num_dims
              )
  ```
- Then you can just run the sampler:
  ```python
  sampler.run_mcmc(starts, num_iterations)
  ```

## Emcee Retrieving Chains
- You can get the iteration chains back from the sampler after a run using `.get_chain()`
  - This will usually return a 3D array, indexing over the parameter, walker, and iteration
  - Can visualize a particular parameter over all walkers using
    ```python
    plt.plot(sampler.get_chain()[:,:,0], 'k', alpha=0.3)
    ```
- After examining, will commonly want to discard the burn in and flatten all the individual walkers:
  ```python
  flat_samples = sampler.get_chain(discard=num_dis, 
                                   flat=True)
  ```

## Interpreting Results
- Commonly several things you'll want to look at after flattening the chains:
  - Histograms of the individual parameter distributions
  - 2D Histograms/KDE plots of pair-wise parameter combinations
  - Visualizing the best fit with errors back on the original dataset
- The first two can be done individually, or they can be streamlined using the `corner` package
  - `corner` will automatically generate both individual parameter distributions and all pair-wise 2D histograms
-->


# Going Further
## The Distance Ladder
\begin{tikzpicture}%%width=75%
  \foreach \x in {-1.0,-.5,...,8}{
    \draw (\x,0) --+(0,5);
  }
  \node[above] at (3.5,5.5) {Distance (in pc)};
  \begin{scope}[font=\scriptsize]
    \node[above] at (-1.0,5) {$10^{-6}$};
    \node[above] at (0.5,5) {$10^{-3}$};
    \node[above] at (2.0,5) {$1$};
    \node[above] at (3.5,5) {$10^3$};
    \node[above] at (5.0,5) {$10^6$};
    \node[above] at (6.5,5) {$10^9$};
    \node[above] at (8.0,5) {$10^{12}$};
    \node[black,fill=Cyan, text=white, anchor=west, minimum width=1cm] at (-1,1) {Radar};
    \node[black,fill=Yellow,text=white,  anchor=west, minimum width=2.5cm] at (-.25,2) {E Parallax};
    \node[black,fill=Red,text=white,  anchor=west, minimum width=1.75cm] at (2.0,3) {S Parallax};
    \node[black,fill=Purple,text=white,  anchor=west, minimum width=1.75cm] at (3.75,4) {Cepheids};
    \node[black,fill=Blue,text=white,  anchor=west, minimum width=2.0cm] at (2.5,1) {MS Fits};
    \node[black,fill=Green,text=white,  anchor=west, minimum width=1.5cm] at (5.25,2) {Type I};
    \node[draw, very thick, dashed, Red, anchor=west, minimum width=2.5cm] at (6,3) {?};
  \end{scope}

\end{tikzpicture}

## A Red Tale
- Vesto Melvin Slipher - 1912
  - While observing spiral galaxies, found that they **all** seemed to be redshifted by some amount!
  - This would imply that they are **all** moving away from us, which seemed somewhat odd
- Edwin Hubble - 1929
  - Used Type Ia supernova to estimate the distances to distant galaxies
  - Found that the more distant galaxies were more redshifted

:::{.block name='The Hubble Relation' .fragment .alert}
The more distant an object is, the faster it is moving away from us!
:::

## Hubble's Law

![](../images/python/hubble.png){width=80%}

## The Hubble Constant
::::::cols
::::col
- Has varied greatly throughout its history
- Still empirically (experimentally) determined
- Current estimates range between 65-78
- In a bit of a crisis atm, as competing methods give very different answers
  - From stars: $\approx 73$ km / s / Mpc
  - From CMB: $\approx 68$ km / s / Mpc
::::

::::col

![](../images/python/hubble_evolution.png){width=100%}

::::
::::::


## Puffing Up
- So all galaxies are moving away from us, but surely we aren't actually in the center?
  - **Nope!**
  - But then again, neither is anyone else!
- The _Cosmological Principle_: at a given cosmic time, the universe looks basically the same to all observers.
- Everyone must see everything moving away because the entire universe is actually expanding!


## They are a crusty bunch...
::::::cols
::::col
- The Raisin Bread Analogy
  - Raisins are the galaxies (or stars)
  - The dough is space
  - At it rises and cooks, all the raisins move away from one another
- Fails the cosmological principle though!
  - The creatures of the crust
::::

::::col
![](../images/ch16_raisinbread.jpg){width=100%}

::::
::::::

## Snakes on a @#$*&@ Plane!
- Imagine yourself a smooshed interstellar snake, living on a flat sheet of paper
- You can move around on the paper, but can not lift off or dig into the paper

\begin{tikzpicture}%%width=30%
  \draw[ultra thick, fill=Yellow!50] (0,0) --++(2,4) --++(3,0) --++(-2,-4) --cycle;
  \node[circle,rotate=0, xslant=.9, inner sep=0pt, outer sep=0pt] (sm) at (3,3) {\includegraphics[width=1cm]{ch26_snake.pdf}};
  \draw[very thick, -latex, Green] (sm.40)--+(40:.5);
  \draw[very thick, -latex, Green] (sm.-40)--+(-40:.5);
  \draw[very thick, -latex, Red] (.5,.5) --+(0,1);
  \draw[very thick, -latex, Red] (.5,0) --+(0,-.7);
\end{tikzpicture}

- This is just the raisin loaf so far


## Snakes on a @#%*^% Sphere?
::::::cols
::::col
![](../images/snake_sphere.svg)

::::

::::{.col style='font-size:.9em'}
- Suppose now we connected the ends of the paper to make it into a perfect sphere
- Now your "universe" has:
  - No center
  - No edge
  - Looks the same regardless of where you are at!
- Inflating the sphere will increase all the distances
- The Hubble constant measures how quickly the sphere inflates
::::
::::::


## Expanding Space
:::{style='font-size:.9em'}
- Galaxies appear to move only because the space between is expanding
  - Galaxies are just conveniently glowing marked points in space
- Space was expanding long before there were galaxies though
- Galaxies themselves remain largely the same size
  - Gravity holds them together and determines size
  - We only observe the effects of expanding space when looking over immense regions of space
- Light is redshifted because the space expands, stretching the wavelength
:::

\begin{tikzpicture}%%width=80%
[scale=.4]
  \draw[cyan, domain=0:9*pi, smooth, samples=100, ultra thick] plot (\x, {sin(deg(\x))});
  \draw[red, domain=0:9*pi, smooth, samples=100, ultra thick] plot (\x, {sin(deg(0.5*\x))-4});
\end{tikzpicture}


## Einstein's Demon
- Einstein's general relativity + a homogeneous universe predicts either an expansion or contraction of space
- Einstein hated this, and was convinced it couldn't be true
- Originally added an extra term, a _"cosmological constant"_ to his equations to allow for a static, unchanging universe

::::::cols
::::col
![](../images/ch16_fieldeqns.jpg){width=80%}
::::

::::col

$$R_{ab}-\frac{1}{2}Rg_{ab} = -8\pi T_{ab} + \Lambda g_{ab}$$

::::
::::::

## Expansion and Age
- If everything is expanding, we can reverse it to figure out how old the universe is
- "Hubble Time"
  $$ t_H \approx \frac{1}{H} $$
- Comes out to about 14 billion years with current estimates
- Note that this is assuming the rate of the universe is constant!


## Homework 6 Partners
- Below are your partners for the coming HW6. Check in with them!

::::::cols
::::col
- Conor & Jared
- Lucca & Gabby
- M & Pearson
- Maddie & Ema
- Sawyer & Aurora
- Sergio & Evyn
::::

::::col
- Evan & Mamadou
- Sage & Felicity
- Salem & Luca 
- Izzy & Tegan
- Sadie & Greg
- Oscar, Elliott, and Luna
::::
::::::