Lab 2.2

Lab 2.2#

Part 1: Correlation between two independent random variables.#

In sample statistics we talked about how one can get spurious correlations between two independent random variables.

Let’s say we have two set of $N$ observations from two random variables $X_1,X_2$ and we want to test the null hypothesis that they are independent using squared correlation ($\rho^2$) where $\rho$ is the Pearson correlation coefficient.

Question 1.1#

For $N=20$ simulate numerically then null distribution of $\rho^2$ under the nully hypothesis that $X_1$ and $X_2$ are independent Gaussian random variables with the same variance. Plot the distribution and the critical value for rejecting the null hypothesis with 95% confidence.

# Solution

Question 1.2#

Plot the critical value as a function of $N$ for values of $N$ between 5 and 50.

# Solution

Part 2: Obliquity Pacing of Pleistocene Ice Ages#

Required reading: Huybers and Wunsch (2004).

In this lab we will apply hypothesis testing to evaluate whether obliquity might be pacing the Pleistocene Ice Ages. This lab broadly follows Huybers and Wunsch (2004).

Since the data from that paper is not available, we will use a slightly different dataset and technique. Ideally, the result of our hypothesis test should not be too sensitive to small methodological differences.

Data: Pleistocene $\delta^{18}$O and Obliquity#

We’ll use the Lisiecky and Raymo stack of sediment $\delta^{18}O$ measurements, which is representative of ice volume during the Pleistocene Ice Ages.

We’ll use the Berger ‘91 Orbital solution to look at variations in Obliquity over the last 700kyrs.

# Let's load the data and take a look at it
url='https://raw.githubusercontent.com/cdds-uiuc/simles-book/main/content/Module02/data/LR04.csv'
df_LR04   =pd.read_csv(url);
age=df_LR04['age']
d18=df_LR04['d18O']

# Load orbital solutions
url='https://github.com/cdds-uiuc/simles-book/raw/main/content/Module02/data/orbital.csv'
df_orbital=pd.read_csv(url);
orb_age=-df_orbital['age']
obliq  =df_orbital['obliquity']

../../_images/99fe85ed9dd5d4a7abc9dfa5a434be79d2870dde6dd2d8a8650fce6bcdac01c7.png

Following hallowed paleoclimate tradition we will simply identify the mid-points of terminations by hand. We’ll then visualize the data

terminations=np.asarray([622, 534, 423, 336, 244, 220, 132, 13])

../../_images/6974e6a92f398a040967443a08ea6159ea446a87169a7ceabfdfa3cd421d66a8.png

Phase of terminations#

Next, we’ll need to compute the phase of obliquity at the moment of terminations. Huybers Wunsch 2004 mention a “stroboscopic” sampling of obliquity, but do not give much detail. So instead, we’ll use a very rough methodology illustrated below, where we take advantage of the fact that obliquity varies with a period of aproximately 41kyr.

Illustration of phase calculation:#

We’ll take at a point in time $t_0$ and we’ll look at the previous $41$ kyr before $t_0$ (one full period of obliquity). We’ll normalize obliquity to have zero mean, and a range of -1 to 1, and call the normalized obliquity over the 41ka preceeding $t_0$ as $y(t_0)$.

\[y(t_0)=\frac{obliq(t_0)-\overline{obliq(t_0)}}{\left[\text{max}(obliq(t_0))-\text{min}(obliq(t_0))\right]/2}\]

We’ll compare this normalized obliquity with a cosine function

$$\hat y(t_0,\phi) =\cos\left(2\pi\frac{t-t_0}{T} +\phi\right)$$,

where$\phi$ is the phase difference.

If the best fit between normalized obliquity, $y(t_0)$ and the cosine $\hat y(t_0,\phi)$ is obtained for $\phi=0$ that means that obliquity peaks right at $t_0$. In the example below, for $t_0$ the best fit is for $\phi=60$, indicating that the obliquity peaked before $t_0$. Since the entire cycle of 360 degrees takes 41 years, it means obliquity peaked at a time $t=60/360*41\approx 9$ ka.

# Some functions we need to compute phase
t0=42
period=41
time=t0+np.arange(41)

def normalize(obliq):
    y=(obliq-np.mean(obliq))/(np.max(obliq)-np.min(obliq))*2
    return y

def y_hat(phi):
    period=41
    y_hat = np.cos(2*np.pi*np.arange(41)/41+phi*np.pi/180)
    return y_hat

(42.0, 84.0)

../../_images/35a7b2098a3c131f44ad62178bc0a9f687eb64813b81471692ad8bbea1242424.png

Phase estimation#

At any point $t_0$ we compute the phase of obliquity by finding the phase difference $\phi$ that minimizes the misfit over the 41 ka preceeding $t_0$ between the normalized obliquity and a shifted cosine

\[\phi =\arg \text{min}_\phi \left[J(\phi,t_0)\right]\]

\[J(\phi, t_0)=\sum_{t=t_0}^{t_0+41}\left[ y(t)-\cos\left(2\pi\frac{t-t_0}{T} +\phi\right) \right]^2\]

\[J(\phi, t_0)=\sum_{t=t_0}^{t_0+41}\left[ y(t)-\hat y(\phi,t,t_0) \right]^2\]

# define a cost function that computes the phase misfit 
def J_cost(phi,y):
    res=np.sum((normalize(y)-y_hat(phi))**2)
    return res

obliq_phase=np.zeros(700)
t0=np.arange(700)

# for each time point over the last 700kyr, compute the phase as the phase that minimizes the misfit:
for j in np.arange(700):
    time=t0[j]+np.arange(41)
    y=obliq[j:j+41]
    out=optimize.least_squares(J_cost,0,args=[y]);
    obliq_phase[j]=out.x[0]

We can plot the phase of obqliquity (normalized such that $-180^\circ:180^\circ$ corresponds to $-1:1$). Our estimation seems to have worked pretty well. Notice that the obliquity minima correspond to $\phi=-1\equiv 180^\circ$ degrees of phase difference, and the obliquity maxima correspond to a phase of $\phi=0\equiv 0^\circ$

../../_images/4a3ec88a5938e771b0d78a892e51e90f26ddfaf5f1681f29d6179ec7511a6769.png

Next we’ll evaluate the phase of obliquity at glacial terminations, plot the phasors and compute Rayleigh’s R (Fig.1 in Huybers & Wunsch 2004).

Rayleigh’s R is the test statistic for our null hypothesis test. It is measure of how tight are obliquity and terminations coupled. If R is large, than terminations always happen near the same phase of obliquity. If R is small, terminations can occur at any point in an an obliquity cycle.

# compute the phase of obliquity at terminations and transform the phase in radians
phi_obliq_terms    =obliq_phase[np.asarray(terminations)]*np.pi/180

# compute Raileigh's R
R_obliq=1/8*np.abs(np.sum(np.cos(phi_obliq_terms)+1j*np.sin(phi_obliq_terms)))

../../_images/ae66a561e369ab800911b14f96042938542690617ca89e7bad2a1d704c78f1df.png

In the figure above the circle represents the phase of obliquity. Each red dot represents the phase of obliquity at each of the terminations (1-7). These would be called phasors.

Rayleigh’s $R$ represents a vector average of the phasors for all terminations.

The Null Hypothesis#

Rayleigh’s R is the test statistic for our null hypothesis test. It is measure of how tight are obliquity and terminations coupled.

We want to know whether such tight coupling could arise by chance. To do so will generate a null hypothesis wherein glacial terminations are not affected by obliquity. We will then compute the $R$ statistic between obliquity and this model for termination. The model is described in Huybers & Wunsch 2004:

\[ \text{if } V_t<T_0:\;V_{t+1}=V_t+\eta_t\]

where $\eta_t \sim \mathcal N (\mu=1,\sigma=2)$. The non-zero mean will cause an average increase in ice volume. When the ice volume crosses the threshold $V_t>T_0$, it will deglaciate completely over a period of 10 kyr.

from scipy import stats

def H0(obliq_phase):
    V=np.zeros(700)   #initialize ice volume
    t=699             #starting point
    T0=90             #threshold for deglaciation
    terminations=[]   #list to store time of termination
    
    while t>0:    
        if V[t]<T0:    # below threshold
            V[t-1]=V[t]+stats.norm.rvs(loc=1.1,scale=2) 
            t=t-1;        
        elif t>10:     #crossing threshold at >10ky before present
            terminations.append(t-5)            #note termination
            V[t-10:t]=np.linspace(0,V[t],10)  #deglaciate in 10kye
            t=t-10                            
        else:          
            #we have an exception to catch.
            #if we cross the threshold with less than 10ky left before present day, 
            # we can't just skip 10ky. 
            
            V[0:t]=np.linspace(0,V[t],10)[-t:]  #only keep data until present day
            if t>=5:    
                terminations.append(t-5)
            t=0

    np.size(terminations)
    phi_obliq_terms    =obliq_phase[np.asarray(terminations)]*np.pi/180
    R_obliq=1/8*np.abs(np.sum(np.cos(phi_obliq_terms)+1j*np.sin(phi_obliq_terms)))

        
    return R_obliq, terminations,phi_obliq_terms,V

from scipy import stats

def H0(obliq_phase):
    V=np.zeros(700)   #initialize ice volume
    t=0             #starting point
    T0=90             #threshold for deglaciation
    terminations=[]   #list to store time of termination
    
    while t<699:    
        if V[t]<T0:    # below threshold
            V[t+1]=V[t]+stats.norm.rvs(loc=1.1,scale=2) 
            t=t+1;        
        elif t<689:     #crossing threshold at >10ky before present
            terminations.append(t+5)            #note termination
            V[t:t+10]=np.linspace(V[t],0,10)  #deglaciate in 10kye
            t=t+10  
        else:
            t=700;
            #we have an exception to catch.
            #if we cross the threshold with less than 10ky left before present day, 
            # we can't just skip 10ky. 
            
            #V[]=np.linspace(0,V[t],10)[t:698]  #only keep data until present day
            #if t>=5:    
            #    terminations.append(t+5)
            #t=0

    np.size(terminations)
    phi_obliq_terms    =obliq_phase[np.asarray(terminations)]*np.pi/180
    R_obliq=1/8*np.abs(np.sum(np.cos(phi_obliq_terms)+1j*np.sin(phi_obliq_terms)))

        
    return R_obliq, terminations,phi_obliq_terms,V

Let’s generate two realizations of the null nypothesis and plot them. We’ll plot obliquity two realizations of the d18O stack under the null hypothesis, and the phasor R.

../../_images/aea13e086a6fd602c1ac24651266cd804bc77660c0398d7dc0c8fb89b05456a5.png

Question 2.1: Evaluate $H_0$#

Generate a distribution of the null hypothesis, $H_0$, that is a distribution of the R-statistic for samples of the noise-model. Evaluate the p-value for the R-statistics of the glacial terminations.

To do this, follow the example in the Null Hypothesis notebook. Generate 5000 realizations of $R_{H_0}$.Plot the histogram of $R_{H_0}$, the critical value $R_{crit}$ and the observed value $R_{obs}$.

Can we still reject the null hypothesis, despite our simpler approach?

# Solution:

Question 2.2: Alternate hypothesis#

Generate an alternate hypothesis as described in Eqn (2) in Huybers & Wunsch 2004. You will need to take the function H0 (modelling Eqn 1) and make a very simialr funcion H1 (modelling Eqn 2) to model the alternate hypothesis. In this situation, the threshold needs to be dependent on normalized obliquity.

Use a value of a=20 for Eqn 2. (rather than a=15 suggested in the paper).

Are the observations consistent with the alternate hypothesis?