Notes from discussion with H, K, S

[bpbond]

• Joint parameters 😄
• Note this is atmospheric CO2-C, not all C
• Figure 3.2 - perhaps rework as histogram focusing on 2100?
• Carbon cycle filter re results
• LUC as a source of uncertainty
• H: might be a GMD/pedagogical paper - but 3.10 goes beyond
• S: start runs every e.g. 50 years and ask, 100 years from emissions, how much C is in atmosphere based on when emitted?
• I.e., does it matter when C is emitted?
• H: 1T FFI C ≠ 1T LUC C - "one ton is not created equal"
• One goal: "if we want to constrain the C cycle further, here's what we have to learn/measure"
• Present at IAMC?

[bpbond]

L-B meeting 4/12:

• [x] Make a histogram of Ca for 1970 and 2020 and add observed lines
• [ ] Implement a joint PDF toy example. A ~ N(0, 1) and B ~ N(0, 1). Let's say their covariance is N(1, 1).

COV_SLOPE <- 1

# Uncorrelated A and B from their distributions
df <- data.frame(A = rnorm(1000), B = rnorm(1000))

# Given A value, compute how far B is from 1:1
# (we're assuming they're positively corelated 1:1)
df$ABdiff <- abs(df$A - df$B) # <---- add in COV_SLOPE
hist(df$ABdiff)

# Generate a random 'acceptance' value: mean of 1, s.d. of 0.5
df$random <- rnorm(1000, 1, 0.5)
df$ACCEPT <- df$ABdiff < df$random
ggplot(df, aes(A, B, color = ACCEPT)) + geom_point() + geom_abline()
ggplot(df, aes(A, B, color = ACCEPT)) + geom_point() + geom_abline(slope = COV_SLOPE) + facet_wrap(~ACCEPT)

# Now compute cov(A, B (accepted values) ) and confirm it's close to COV_SLOPE above

☝️ we want to generalize to an arbitrary correlation (not just 1:1, i.e. slope of 1)

[bpbond]

@leeyap for further documentation:

Here are the A and B draws before we do anything: df %>% pivot_longer(cols = c("A","B", "ABdiff")) %>% ggplot(aes(x = value)) + facet_wrap(~name, scales="free") + geom_histogram(bins = 20)

Here's the plot showing accept/reject:

Here's the post-reject plot showing the distribution of A, B, and their difference:

[bpbond]

Eventually I think we'd like to start with a correlation coefficient and base our rejection criterion on that.

[bpbond]

Note to myself for later

diffusivity <- function(n) rnorm(n)
climate_sensitivity <- function(n) rnorm(n)

# Generate a joint distribution for A and B using rejections sampling
# https://en.wikipedia.org/wiki/Rejection_sampling
# Parameters:
#   fA    Function that generates n values of A (marginal distribution, ignoring B)
#   fB    Function that generates n values of B (marginal distribution, ignoring A) 
#   n     Total number of values we want
#   pcor  Pearson's correlation between A and B, from -1 to 1
joint_pdf <- function(fA, fB, n, pcor) {

    A <- rep(NA_real_, n)
    B <- rep(NA_real_, n)
    accept <- rep(FALSE, n)
    #randoms <- rnorm(n)
    i <- 0
    while (any(!accept)) {
        i <- i+1
        not_accepteds <- which(!accept)
        n_not_acc <- length(not_accepteds)

        A[not_accepteds] <- fA(n_not_acc)
        B[not_accepteds] <- fB(n_not_acc)

        # test which ones are not accepted, updating accept
        accept[not_accepteds] <- sample(c(TRUE, FALSE), n_not_acc, replace = TRUE)

        message(i, " n_accept = ", sum(accept))
    }
    data.frame(A = A, B = B)  # better: get caller function names and use those
}

x <- joint_pdf(diffusivity, climate_sensitivity, 1000, 1.0)

[bpbond]

MASS::mvrnorm()

[bpbond]

From kalyn:

re the prior I think there is a pretty well established body of literature about it kunti & forest have done a lot of work on it