Lecture 14: Introduction to Age Models

1. The importance of knowing time
2. Building an age model
   • Markov chain Monte Carlo approaches
   • constant sedimentation rates
   • varying sedimentation rates

We acknowledge and respect the lək̓ʷəŋən peoples on whose traditional territory the university stands and the Songhees, Esquimalt and W̱SÁNEĆ peoples whose historical relationships with the land continue to this day.

Importance of knowing time

• we have talked before about the difficulty of recovering signals from stratigraphic sections
• in the transport module, the challenge was primarily to correlate surfaces across the basin

Importance of knowing time

• transgressive regressive cycles are way clearer on a wheeler diagram
• can see the migration of the shoreline
  • basinward during progradation / sea level fall
  • landward during transgression
  • shoreline can be followed by just tracking intertidal deposits
• in addition to being able to qualitatively view the record more clearly, a time model also makes this data more amenable to quantitative analysis

• this is a time series of the wheeler diagram
• non-zero numbers are where net deposition occurred, and the thickness of that unit of deposition
• zeros are time periods of non-deposition
• very similar to columns dataset that we looked at last week (i.e., the presence of these zeros)

• this is the sea level history that forced the model

• what feature in the data does the peak correspond to?
• can also do this as simple binary - presence or absence of sediment?
  • in other words, just identify the gaps
• but signal is totally shredded if we don't have the gaps!

• an FFT of the sea level history reveals the same peak
• some notes:
  • in the real world, it would be hard to detect this signal, partly because the time series is not really long enough
  • forcing may not have to be periodic
  • the "binary" view maybe could be a poisson process rather than a periodic sea level forcing!
• the point is: you need time to see the same signal in data and the forcing

Age models are important: how do we get them?

1. Cyclostratigraphy
2. Biostratigraphy
3. Absolute ages
   • U-Pb (volcanics), Ar-Ar (volcanics), Re-Os (sediments)
4. Signal matching
   • magnetostratigraphy
   • chemostratigraphy

Biostratigraphy

• based on the unique, sequential, nonrepeating appearance of fossils through time

Biostratigraphy

• based on the unique, sequential, nonrepeating appearance of fossils through time

• most important are first appearance datums and last appearance datums
• in principle, a very simple concept

Biostratigraphy

• what is wrong with this picture?

• Tangled-fence diagram of the observed contradictions between 7 sections that preserve 62 taxa of the Cambrian Riley Formation of Texas (data from Palmer, 1954; Shaw, 1964).
• attempted correlations when FADs and LADs are used AS IS in different sections
• lines are meant to represent time lines, so equal time

Contradictory ranges

• what could be causing this?

• problems arise because relationship between taxa are not same in each section
• in this example, in one these taxa coexist, and in another, they do not
• ranges need to be MODIFIED to be consistent
  • assuming coexistence is real, then at section 6, this range needs to be extended

Resolving contradictory ranges

• ranges of original data need modifcation

• in actuality, there are 6 different possibilities that could be considered
• time is older on the left, and younger on the right
• we can rule out some of these, because we have direct evidence of coexistence
• this is pretty tractable, and we could pretty quickly come up with a solution
  • a guiding principle is, do the MINIMUM amount of modification needed to make the ranges consistent
• OK, but what if we are dealing with THREE TAXA
  • number of possibilities grows to 90

Number of possible sequences

• requires constrained optimization (CONOP9; Sadler and Cooper, 2008)

• quickly becomes a Herculean task, if we wanted to actually look at EVERY possible combination!
• number of atoms in universe= 10^82
• need algorithms that only look at the likliest solutions, and avoid those that can be easily discounted (called constrained opitimization)
• familiar to anyone who has played chess - 10^120 possible games

Best-fit solution

• what are some features of this that look familiar to our model outputs?

Building an age model

Building an age model

fig=plt.figure(1,figsize=(15,6))
ax=fig.add_subplot(121)
ax,trace,boxes,liths=rando_strat(ax,num_box=5,height=2,boxes=[],liths=[]) #stratigraphy for fun
#ashes in the strat column
for a in ashes: 
    tmp_w=trace[(trace[:,0]<=ashes[a]['height']) & (trace[:,1]>ashes[a]['height']),2]
    ax.plot([0,tmp_w],[ashes[a]['height'],ashes[a]['height']],'r--')
#KT boundary
tmp_w=trace[(trace[:,0]<=KT) & (trace[:,1]>KT),2]
ax.plot([0,tmp_w],[KT,KT],'k--')
ax.text(tmp_w,KT,'  KT boundary',verticalalignment='center',fontsize=20)
#ash ages
ax=fig.add_subplot(122)
for a in ashes:
    ax.plot(ashes[a]['age'],ashes[a]['height'],'rs')
ax.plot([66.2,65.95],[KT,KT],'k--')
ax.set_xlim([66.2,65.95]); ax.set_ylim([-0.1,2])
ax.set_ylabel('meters'); _=ax.set_xlabel('age (Ma)')

• so how do we get an age for the KT boundary?

Building an age model

fig=plt.figure(1,figsize=(15,6))
ax=fig.add_subplot(121)
ax,trace,boxes,liths=rando_strat(ax,num_box=5,height=2,boxes=boxes,liths=liths) #stratigraphy for fun
#ashes in the strat column
for a in ashes: 
    tmp_w=trace[(trace[:,0]<=ashes[a]['height']) & (trace[:,1]>ashes[a]['height']),2]
    ax.plot([0,tmp_w],[ashes[a]['height'],ashes[a]['height']],'r--')
#KT boundary
tmp_w=trace[(trace[:,0]<=KT) & (trace[:,1]>KT),2]
ax.plot([0,tmp_w],[KT,KT],'k--')
#ash ages
ax=fig.add_subplot(122)
ash_table=np.array([(ashes[a]['age'],ashes[a]['height']) for a in ashes])
ax.plot(ash_table[:,0],ash_table[:,1],'rs-',lw=0.5)
#calculate KT age
KT_age=calc_age(ash_table,KT)
ax.plot([66.2,KT_age],[KT,KT],'k--')
ax.plot([KT_age,KT_age],[KT,KT],'w*',mec='k',markersize=20)
ax.text(KT_age,KT,'   %2.2f Ma.' % (KT_age),horizontalalignment='left',verticalalignment='center',fontsize=20)
ax.set_xlim([66.2,65.95]); ax.set_ylim([-0.1,2]); ax.set_ylabel('meters'); _=ax.set_xlabel('age (Ma)')

• so is that it? are we done?

Building an age model

fig=plt.figure(1,figsize=(15,6))
ax=fig.add_subplot(121)
ax,trace,boxes,liths=rando_strat(ax,num_box=5,height=2,boxes=boxes,liths=liths) #stratigraphy for fun
#ashes in the strat column
for a in ashes: 
    tmp_w=trace[(trace[:,0]<=ashes[a]['height']) & (trace[:,1]>ashes[a]['height']),2]
    ax.plot([0,tmp_w],[ashes[a]['height'],ashes[a]['height']],'r--')
#KT boundary
tmp_w=trace[(trace[:,0]<=KT) & (trace[:,1]>KT),2]
ax.plot([0,tmp_w],[KT,KT],'k--')
#ash ages
ax=fig.add_subplot(122)
for a in ashes:
    ax.plot(ashes[a]['age'],ashes[a]['height'],'rs')
    ax.plot([ashes[a]['age']+ashes[a]['2sd'],ashes[a]['age']-ashes[a]['2sd']],
            [ashes[a]['height'],ashes[a]['height']],'r-')
ax.plot([66.2,65.95],[KT,KT],'k--')
ax.set_xlim([66.2,65.95]); ax.set_ylim([-0.1,2]); ax.set_ylabel('meters');_=ax.set_xlabel('age (Ma)')

• bars represent uncertainty in the age determination, 2-sigma level or 95% confidence interval
• often always drawn as lines, but why is this misleading?

Building an age model

fig=plt.figure(1,figsize=(15,6))
ax=fig.add_subplot(121)
ax,trace,boxes,liths=rando_strat(ax,num_box=5,height=2,boxes=boxes,liths=liths) #stratigraphy for fun
#ashes in the strat column
for a in ashes: 
    tmp_w=trace[(trace[:,0]<=ashes[a]['height']) & (trace[:,1]>ashes[a]['height']),2]
    ax.plot([0,tmp_w],[ashes[a]['height'],ashes[a]['height']],'r--')
#ash ages as distributions
ax=fig.add_subplot(122)
for a in ashes:
    #normal distribution
    mu=ashes[a]['age']
    sigma=ashes[a]['2sd']/2
    x = np.linspace(mu - 3*sigma, mu + 3*sigma, 300)
    #scaled for display
    distro=stats.norm.pdf(x, mu, sigma)
    distro=distro/max(distro)*0.3
    ax.plot(x, distro+ashes[a]['height'],'r')
    #ash strat height
    ax.plot([ashes[a]['age']+ashes[a]['2sd']*3/2,ashes[a]['age']-ashes[a]['2sd']*3/2],
            [ashes[a]['height'],ashes[a]['height']],'r--',lw=0.5)
ax.plot([66.2,65.95],[KT,KT],'k--')
ax.set_xlim([66.2,65.95]); ax.set_ylim([-0.1,2]); ax.set_ylabel('meters'); _=ax.set_xlabel('age (Ma)')

• distributions, here assumed to be Gaussian (a reasonale assumption for this kind of data)
• how can we incorporate this uncertainty into our KT age?
• there is an analytical solution, but I will show you a computational-based approached (Monte Carlo)

Building an age model

fig=plt.figure(1,figsize=(15,6))
ax=fig.add_subplot(121)
ax,trace,boxes,liths=rando_strat(ax,num_box=5,height=2,boxes=boxes,liths=liths) #stratigraphy for fun
for a in ashes: 
    tmp_w=trace[(trace[:,0]<=ashes[a]['height']) & (trace[:,1]>ashes[a]['height']),2]
    ax.plot([0,tmp_w],[ashes[a]['height'],ashes[a]['height']],'r--')
#picking ages from distributions
ax=fig.add_subplot(122)
for a in ashes:
    mu=ashes[a]['age']
    sigma=ashes[a]['2sd']/2
    x = np.linspace(mu - 3*sigma, mu + 3*sigma, 300)
    distro=stats.norm.pdf(x, mu, sigma)
    pick=np.random.normal(mu,sigma,1)
    pick_distro=stats.norm.pdf(pick, mu, sigma)
    pick_distro=pick_distro/max(distro)*0.3
    distro=distro/max(distro)*0.3
    ax.plot(x, distro+ashes[a]['height'],'r')
    ax.plot([ashes[a]['age']+ashes[a]['2sd']*3/2,ashes[a]['age']-ashes[a]['2sd']*3/2],
            [ashes[a]['height'],ashes[a]['height']],'r--',lw=0.5)
    ax.plot(pick, pick_distro+ashes[a]['height'],'r*',mec='k',markersize=20)
ax.plot([66.2,65.95],[KT,KT],'k--')
ax.set_xlim([66.2,65.95]); ax.set_ylim([-0.1,2]); ax.set_ylabel('meters'); _=ax.set_xlabel('age (Ma)')

• if I collect all my picks, what will the mean be of all those picks?
• what about the standard deviation of my picks

Building an age model

fig=plt.figure(1,figsize=(15,6))
ax=fig.add_subplot(121)
ax,trace,boxes,liths=rando_strat(ax,num_box=5,height=2,boxes=boxes,liths=liths) #stratigraphy for fun
#ashes in the strat column
for a in ashes: 
    tmp_w=trace[(trace[:,0]<=ashes[a]['height']) & (trace[:,1]>ashes[a]['height']),2]
    ax.plot([0,tmp_w],[ashes[a]['height'],ashes[a]['height']],'r--')
ax=fig.add_subplot(122)
for a in ashes:
    ax.plot([ashes[a]['age']+ashes[a]['2sd'],ashes[a]['age']-ashes[a]['2sd']],
            [ashes[a]['height'],ashes[a]['height']],'r-')
#plot ten paths
np.random.shuffle(paths) #shuffle the deck
for p in paths[0:100]:
    ax.plot(p,[ashes[a]['height'] for a in ashes],'r-',lw=0.5)
#plot MCMC KT age
ax.plot([np.mean(KT_ages)+np.std(KT_ages)*2,np.mean(KT_ages)-np.std(KT_ages)*2],
        [KT,KT],'k-')

ax.plot([np.mean(KT_ages),np.mean(KT_ages)],[KT,KT],'w*',mec='k',markersize=20)
ax.text(np.mean(KT_ages)-np.std(KT_ages)*2,KT,r'   %2.2f$\pm$%2.2f Ma.' % (np.mean(KT_ages),np.std(KT_ages)*2),
        horizontalalignment='left',verticalalignment='center',fontsize=15)
ax.set_xlim([66.2,65.95]); ax.set_ylim([-0.1,2]); ax.set_ylabel('meters'); _=ax.set_xlabel('age (Ma)')

• should there by any rules for my paths? in other words, should I reject any?
  • Markov Chain Monte Carlo

Power of Markov Chain Monte Carlo approaches

• prior distributions need not be Gaussian (figure from Haslett and Parnell, 2008)

• MCMC is especially helpful when prior distributions are not normally distributed
  • common for carbon 14 ages, because calibration curve is non-linear
• MCMC is GENERAL
  • allows to explore uncertainty for any type of model
  • x = A + BxC - D (normal, uniform, bimodal, complex)
• what is a major assumption of this model?
  • constant sedimentation rate between tie points

How many samples is enough?

fig=plt.figure(1,figsize=(15,6))
ax=fig.add_subplot(111)
idx=list(range(1000)) +list(range(1000,len(KT_ages)+1,100))
#statistics of KT age as function of number of trials
KT_evolve=[]
for i in idx:
    KT_evolve.append((np.mean(KT_ages[0:i+1]),np.std(KT_ages[0:i+1]),len(KT_ages[0:i+1])))
KT_evolve=np.array(KT_evolve)
#plot the results
ax.plot(KT_evolve[:,2],KT_evolve[:,0],'r-',label='mean age')
ax.fill_between(KT_evolve[:,2],KT_evolve[:,0]+KT_evolve[:,1],KT_evolve[:,0]-KT_evolve[:,1],
                color='k',alpha=0.5,zorder=0,label='mean +/- 2s.d.')
ax.set_ylabel('age (Ma.)'); ax.set_xlabel('number of trials');_=ax.legend(loc='best')

• spend rest of class developing a more complex age modeling approach
• first I want to address the question from Grayson - how many trials are enough?
• 1st look at evolving mean KT age from 1000 trials
• 2nd look at evolving mean KT age from 10000 trials
• 3rd consider in context of 2 sigma uncertainty

Constant sedimentation rate assumption

def make_path(pts,g_shape,g_scale,delta):
    #separate a number of points (X) along a [0,1] path
    if pts!=0:
        #X numbers are drawn from a gamma distribution (G)
        #--> these represent the distances between successive points
        path=np.cumsum(np.random.gamma(g_shape,g_scale, pts))
    else:
        path=[1]#if X = 0, then path is simply [0,1]
    path=np.hstack((0,path))#add 0 as the path beginning
    #scale first to be between 0 and 1, and then 0 to delta
    path=path/path[-1]*delta
    return(path)

fig=plt.figure(1,figsize=(6,6)); ax=fig.add_subplot(111)
#assumes constant sedimentation rate between anchors
ax.plot([0,1],[0,1],'-s');ax.set_xlabel('change in time'); _=ax.set_ylabel('change in height')

• goal: how can we get away from the assumption of constant sedimentation rate?

Constant sedimentation rate assumption

fig=plt.figure(1,figsize=(6,6)); ax=fig.add_subplot(111); 
ax.set_xlabel('change in time'); _=ax.set_ylabel('change in height')
num_paths=200 #add knots to the path
for i in range(num_paths):
    pts=np.random.poisson(5, 1)[0] #discrete integer
    t_path=make_path(pts,1,1,1) #age change between knots
    h_path=make_path(pts,1,1,1) #position change between knots
    ax.plot([0]+list(np.cumsum([1/pts]*pts)),[0]+list(np.cumsum([1/pts]*pts)),'-s')
#     ax.plot(t_path,h_path,'-s') #plot the path
    #plot many paths
    ax.plot(t_path,h_path,'k-',lw=0.5,alpha=0.25) 

Constant sedimentation rate assumption

...
for i in range(num_paths):
    pts=np.random.poisson(5, 1)[0] #discrete integer
    ...

Poisson distribution is a discrete probability distribution that expresses the probability of a given number of events occurring in a fixed interval of time or space if these events occur with a known constant mean rate and independently of the time since the last event.

Exponential distribution is the probability distribution of the time between events in a Poisson point process

• approach: between our two anchors, randomly insert N points
• now its a line with N+1 segments
• 1st still assume constant sedimentation rate
  • but we can allow sed rate to vary ...
• 2nd allow sed rate to vary
  • one path, two paths, 5 paths
• 3rd look at 1000 paths

import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
fig=plt.figure(figsize=(10,10))
a = np.random.uniform(0,1,100000)
b = a<0.1  
sns.kdeplot(np.diff(np.where(b)[0]),color=sns.color_palette('deep')[2],fill=True,lw=3)
c = np.random.exponential(10,100000)
sns.kdeplot(c,color=sns.color_palette('deep')[1],fill=True,lw=3)
sns.despine()

Deciding on number of sedimentation rate changes

#number of timesteps
t=200000
#probability of change
prob=0.25
#start with sediment on
state=1
#switch state
sed=np.array([state*-1 if np.random.random()<prob else state for i in range(t)])
#plot the sedimentation rate history
fig=plt.figure(1,figsize=(25,7.5));ax=fig.add_subplot(111)
#plot as timeseries
ax.plot(sed[0:100]); ax.set_ylabel('sediment switch'); _=ax.set_xlabel('time')

• How do we pick number of see rate changes?
• Draw from a Poisson distribution
• For a not quite exact analogy, we can relate to a Poisson process
• What is frequency of sedimentation event duration?

Deciding on number of sedimentation rate changes

sns.histplot(np.random.poisson(1,10000),
             binwidth=1)
_=plt.gca().set_xlabel("discrete events over an interval of time")

sns.histplot(np.random.exponential(1,10000))
_=plt.gca().set_xlabel('interval of time between events')
# plt.gca().set_yscale('log')

• 1st show on linear scale: exponential
  • Shortest events most common
  • Longer are rare
• 2nd show on a logarithmic one

Deciding on number of sedimentation rate changes

fig=plt.figure(1,figsize=(7.5,7.5)); ax=fig.add_subplot(111); 
ax.set_xlabel('change in time'); _=ax.set_ylabel('change in height')
num_paths=1 #add knots to the path
for i in range(num_paths):
    pts=np.random.poisson(5, 1)[0]
    #age change between knots
    t_path=make_path(pts,1,1,1)
    #position change between knots
    h_path=make_path(pts,1,1,1)
#     #plot the knot number
#     ax.plot([0]+list(np.cumsum([1/pts]*pts)),[0]+list(np.cumsum([1/pts]*pts)),'-s')
    #plot the path
    ax.plot(t_path,h_path,'-s')

• so the number of sedimentation rate changes determined from a Poisson distribution
• now we will turn our attention to how we estimate sedimentation rates between these points
• before we do, just using what we know so far, does anyone know how we can force a constant sed rate assumption?
  • 1st is more obvious, make number of points = 0
  • 2nd is make the numbr of points ENORMOUS, lambda = 5000

Picking sedimentation rates

fig=plt.figure(1, figsize=(12,6)); ax=fig.add_subplot(111)
g_shape=1.5; g_loc=0; g_scale=1
gam_hist=np.random.gamma(g_shape,g_scale, 10000) #10K draws with shape = 1.5 and scale = 1
x=np.linspace(0,max(gam_hist),1000) 
gam=stats.gamma.pdf(x,g_shape,g_loc,g_scale) #continuous function
ax.hist(gam_hist,density=True,alpha=0.75,color='#6495ED', 
        label=r'$\lambda$ = %2.1f; $\mu$ = %2.3f' % (g_shape*g_scale,np.mean(gam_hist))); ax.plot(x,gam,'k--')
g_shape=10.5; g_loc=0; g_scale=1 #10K draws with shape = 10.5 and scale = 1
gam_hist=np.random.gamma(g_shape,g_scale, 10000)
x=np.linspace(0,max(gam_hist),1000); gam=stats.gamma.pdf(x,g_shape,g_loc,g_scale) #continuous function
ax.hist(gam_hist,density=True,alpha=0.75,color='#B22222',
        label=r'$\lambda$ = %2.1f; $\mu$ = %2.3f' % (g_shape*g_scale,np.mean(gam_hist))); ax.plot(x,gam,'k--')
#plot labels
ax.set_title('examples of gamma distributions')
ax.legend(); _=ax.text(0.99,0.725,', '.join(['%2.3f' %(g) for g in gam_hist[0:6].tolist()]) + ', ... ',
                       transform=ax.transAxes,horizontalalignment='right',fontsize=15)

• now we need to pick the sedimentation rates along our path
• to do this, for each segment, we need a dh and a dt
• draw both from gamma distributions
• very similar to Poisson distributions, except they are continous
  • blue and red here should look very similar to the blue and red Poisson distros
• 2nd change this parameter to 1, get an exponential distribution (just like the Poisson process experiment)
  • again, red can come from the same math, and we will gain some insight into how this is possible in a sec

gamma distributions in sedimentology

• gamma distributions are very common in sedimentology and stratigraphy
• this paper showed that the distribution of thicknesses of limestone mud beds followed a gamma distribution

gamma distributions in sedimentology

• logic is the same as what we are after
• have some path, with a start and a stop
  • randomly chop up that path
  • this mimics a Poisson process of building sediment in random increments
  • length of segments = gamma distribtion (exponential distribution)
• used as an argument AGAINST cyclcity (without using time series analysis)

gamma distributions in sedimentology

• what happens when we bundle couplets in another Poisson process?
• and bundling is itself a Poisson process
• get towards our red distributions - look more and more Gaussian

compound Poisson gamma age models

fig=plt.figure(1,figsize=(7.5,7.5)); ax=fig.add_subplot(111); 
ax.set_xlabel('change in time'); _=ax.set_ylabel('change in height')
num_paths=200 #add knots to the path
for i in range(num_paths):
    pts=np.random.poisson(5, 1)[0]
    #age change between knots
    t_path=make_path(pts,1,1,1)
    #position change between knots
    h_path=make_path(pts,1,1,1)
#     #plot the knot number
#     ax.plot([0]+list(np.cumsum([1/pts]*pts)),[0]+list(np.cumsum([1/pts]*pts)),'-s')
    #plot the path
    ax.plot(t_path,h_path,'k-',lw=0.5,alpha=0.25)
_=ax.plot(t_path,h_path,'-s')

• so throughout the lecture yesterday, we built another type of age model that does not need to assume constant sedimentation rate between two anchors
• called "compound Poisson gamma age models", which is a mouthful
• all it is saying is number of changes in sedimentation rate = drawn from a Poisson distribution
• and those sedimentation rates (defined by dt and dh) = drawn from gamma distributions
• both distrubutions look very similar
  • one discrete, one continuous
  • low mean --> exponential
  • high mean --> Gaussian
• NB: connection back to your problem set:
  • a Poisson process can lead to a gamma distribution
  • like sedimentation on/off --> thickness of couplets
  • sorry names are confusing!!

compound Poisson gamma age models

import pandas as pd
import time
#local way to load data
sluggan=[{'position': 44.5, 'age': 985.0, '2sd': 45.0}, {'position': 49.5, 'age': 1225.0, '2sd': 65.0},
        {'position': 69.0, 'age': 1635.0, '2sd': 75.0}, {'position': 102.0, 'age': 2130.0, '2sd': 45.0},
        {'position': 125.0, 'age': 2930.0, '2sd': 85.0}, {'position': 165.0, 'age': 3945.0, '2sd': 85.0},
        {'position': 185.0, 'age': 4180.0, '2sd': 90.0}, {'position': 232.5, 'age': 4556.666667, '2sd': 76.66666667},
        {'position': 239.0, 'age': 4965.0, '2sd': 75.0}, {'position': 272.5, 'age': 5320.0, '2sd': 65.0},
        {'position': 297.5, 'age': 6760.0, '2sd': 90.0}, {'position': 332.0, 'age': 7855.0, '2sd': 115.0},
        {'position': 367.5, 'age': 8176.666667, '2sd': 65.0},{'position': 407.0, 'age': 8540.0, '2sd': 120.0},
        {'position': 427.0, 'age': 9360.0, '2sd': 150.0},{'position': 447.5, 'age': 9475.0, '2sd': 145.0},
        {'position': 461.0, 'age': 9610.0, '2sd': 130.0},{'position': 484.5, 'age': 10805.0, '2sd': 125.0},
        {'position': 499.0, 'age': 10995.0, '2sd': 160.0},{'position': 509.0, 'age': 11625.0, '2sd': 160.0},
        {'position': 518.0, 'age': 12265.0, '2sd': 125.0}]
print('number of ages: %s' %(len(sluggan)))
#convert list of dictionaries into a DataFrame
sluggan=pd.DataFrame(sluggan)
sluggan['position']=sluggan['position']*-1
sluggan=sluggan.sort_values(by='position')
sluggan=sluggan.reset_index()
sluggan.head()

number of ages: 21

	index	position	age	2sd
0	20	-518.0	12265.0	125.0
1	19	-509.0	11625.0	160.0
2	18	-499.0	10995.0	160.0
3	17	-484.5	10805.0	125.0
4	16	-461.0	9610.0	130.0

compound Poisson gamma age models

compound Poisson gamma age models

fig=plt.figure(1,figsize=(12,6)); ax=fig.add_subplot(111)
#randomly pick ~100 age models (could be less)
idx=np.random.randint(len(true_paths))
for a in true_paths[idx:idx+100]:
    ax.plot(a[:,0],a[:,1]/100,'-',lw=0.5)
#statistics on interpolated ages
# ax.fill_betweenx(path_stats['position']/100,path_stats['2.5'],path_stats['97.5'],
#                 color='k',edgecolor='none',alpha=0.3,zorder=0)
#plot the median age (constant sed rate)
# ax.plot(path_stats['50'],path_stats['position']/100,'k--',lw=1)
#plot age constraints 
# for i,a in sluggan.iterrows():
#     ax.plot([a['age']+a['2sd'],a['age']-a['2sd']],[a['position']/100,a['position']/100],'k-',lw=2)
#plot age of demo horizon
# ax.plot([demo['age']+demo['2sd'],demo['age']-demo['2sd']],[demo['position']/100,demo['position']/100],'r-',lw=2)
#plot labels
ax.set_xlim([13000,0]); ax.set_xlabel('age (years BP)');_=ax.set_ylabel('depth (m)')

• 1st plot 100 paths and full range of all paths
• 2nd 95% of all paths and an interpolated age
• 3rd plot the median of all paths
• what does this represent?

compound Poisson gamma age models

#RUN AGE MODEL FOR KT U-Pb DATA
tic=time.time()
#pick self-consistent ages (based on Markov Chain Monte Carlo) 
num_trials=1000
KT_picks=age_pick(KTashes,num_trials)
#define lambda for poisson draws
#NB: if lambda=0, "constant sed rate model" is result
p_lam=0
p_lam=10.5
#define shape and scale for gamma draws
#NB: age model only sensitive to shape
g_shape=1.5
g_scale=1\
#number of interpolation points per age anchor
num_interp=1000
#run compound Poission-Gamma age model
true_paths,interp_paths,path_stats,age_models = make_age_model(KTashes,KT_picks,p_lam,g_shape,g_scale,num_interp)
#age of KT Boundary
KT={'position':0.6,
       'age': np.mean([a(0.6) for a in age_models]),
        '2sd': 2*np.std([a(0.6) for a in age_models])}
toc=time.time()      
print(r'number of trials: %s' %(num_trials))
print(r'model run time: %2.2f seconds' % (toc-tic))
print(r'age of KT boundary produce by model: %2.2f +/- %2.2f Ma.' %(KT['age'],KT['2sd']))

number of trials: 1000
model run time: 0.47 seconds
age of KT boundary produce by model: 66.08 +/- 0.04 Ma.

compound Poisson gamma age models

fig=plt.figure(1,figsize=(12,6)); gs = fig.add_gridspec(1,4); ax=fig.add_subplot(gs[0])
ax,trace,boxes,liths=rando_strat(ax,num_box=5,height=2,boxes=boxes,liths=liths) #stratigraphy for fun
ax=fig.add_subplot(gs[1:])
#statistics on interpolated ages
ax.fill_betweenx(path_stats['position'],path_stats['2.5'],path_stats['97.5'],color='k',edgecolor='none',alpha=0.3,zorder=0)
#plot age constraints
for i,a in KTashes.iterrows():
    ax.plot([a['age']+a['2sd'],a['age']-a['2sd']],
            [a['position'],a['position']],'k-',lw=2)
#     ax.plot(np.mean(KT_picks[:,i]),a['position'],'r+',alpha=0.5,markersize=10)
#     ax.plot([np.mean(KT_picks[:,i])+2*np.std(KT_picks[:,i]),np.mean(KT_picks[:,i])-2*np.std(KT_picks[:,i])],
#            [a['position'],a['position']],'r-',lw=3,alpha=0.5)
#plot and label KT age
ax.text(KT['age']-KT['2sd'],KT['position'],'   %2.2f$\pm$%2.2f Ma.' % (KT['age'],KT['2sd']),fontsize=15,horizontalalignment='left',verticalalignment='center')
ax.plot([KT['age']+KT['2sd'],KT['age']-KT['2sd']],
        [KT['position'],KT['position']],'k--',lw=2); ax.plot(KT['age'],KT['position'],'w*',mec='k',markersize=20)
ax.set_xlim([66.2,65.95]);ax.set_xlabel('age (Ma.)');_=ax.set_ylabel('height (m)')

• arrive at the same mean age, but uncertainty has been doubled
• with age models, how well we know the answer is just as important as the answer!
• also, more information could be brought into this framework
  • for example, what if this is a shallowing upward cycle, and behaves like our model?

compound Poisson gamma age models

• can apply the same concept to a much larger scale
• no longer single sections, but an entire geologic Period
• compared to linear model, mean answer is similar, but uncertainty is larger

Appendix: more on gamma distributions in sedimentology

#timesteps; probability of change; starting state
t=20000; prob=0.25; state=1
#minimum beds in a parasequence
min_couplet_num=3

#switch state
sed=np.array([state*-1 if np.random.random()<prob else state for i in range(t)])
#find the "on" periods
on=np.where(sed==1)[0]

#find indices of consecutive "on" periods
on_breaks=list(np.where(np.diff(on)!=1)[0]+1)

#group into pairs
on_breaks=[0]+on_breaks+[len(on)]
on_breaks=list(zip(on_breaks[0:-1],on_breaks[1:]))

#calculate lengths of "sediment on" period
on_time=[]
for o in on_breaks:
    on_time.append(len(on[o[0]:o[1]]))
    
#here, the "on" periods are used to bundle couplets
#--> value = number of couplets
on_time=np.cumsum(np.array(on_time)+min_couplet_num-1)
couplet_idx=list(zip(list(on_time[0:-1]),list(on_time[1:])))
couplet_idx=[tuple((0, couplet_idx[0][0]))]+couplet_idx
#generate thickess of those couplets
couplets=np.random.gamma(1,1, couplet_idx[-1][1])
#package them up
bundles=[np.sum(couplets[i[0]:i[1]]) for i in couplet_idx]

Appendix: more on gamma distributions in sedimentology

#first hundred couplets for plotting
demo_couplets=list(np.cumsum(couplets)[0:100])
demo_couplets=[0]+demo_couplets
demo_couplets=list(zip(demo_couplets[0:-1],demo_couplets[1:]))
#make thickness vs. height boxes
couplet_boxes=[]
for d in demo_couplets:
    xy=np.array([(0,d[0]),(d[1]-d[0],d[0]),(d[1]-d[0],d[1]),(0,d[1])])
    rect = Polygon(xy,closed=True,facecolor="#"+''.join([np.random.choice([a for a in '0123456789ABCDEF']) for j in range(6)]))
    couplet_boxes.append(rect)
#bundles of couplets
demo_bundles=np.cumsum(bundles)[0:100]
demo_bundles=demo_bundles[demo_bundles<demo_couplets[-1][1]]

Appendix: more on gamma distributions in sedimentology

fig=plt.figure(1,figsize=(20,11)); ax=fig.add_subplot(121); ax.plot(sed[0:100],range(100)) #plot sedimentation history
ax.set_ylabel('time'); ax.set_xlabel('sediment switch'); ax.set_xticks([-1.0,1.0]); ax.set_xticklabels(['off','on'])
ax=fig.add_subplot(122); _=[ax.add_patch(r) for r in couplet_boxes] #plot bed thickness vs. height
ax.set_xlim([0,np.ceil(max(couplets[0:100]))]); ax.set_ylim([0,np.ceil(demo_couplets[-1][1])])
ax.set_ylabel('height (meters)');ax.set_xlabel('bed thickness (meters)')
_=[ax.plot(ax.get_xlim(),[d,d],'k--') for d in demo_bundles] #plot the parasequence boundaries
_=ax.text(ax.get_xlim()[1],demo_bundles[0],'parasequence top',horizontalalignment='right',verticalalignment='bottom',fontsize=15)

Appendix: more on gamma distributions in sedimentology

#histogram of bed thickness
fig=plt.figure(1,figsize=(20,11)); ax=fig.add_subplot(221)
ax.hist(couplets,bins=30,density=True); ax.set_xlabel('bed thickness (meters)')
#histogram of number of beds in a parasequence
ax=fig.add_subplot(222); ax.hist([i[1]-i[0] for i in couplet_idx],bins=30,density=True); ax.set_xlabel('beds per parasequence')
#histogram of parasequence thickness
ax=fig.add_subplot(212); ax.hist(bundles,bins=30,density=True);_=ax.set_xlabel('parasequence thickness (meters)')