-
Notifications
You must be signed in to change notification settings - Fork 2
/
emergent_constraint.py
103 lines (74 loc) · 3.39 KB
/
emergent_constraint.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
#!/usr/bin/env python
# coding: utf-8
#------------------------------------------------------------------------------------#
# #
# SCRIPT for basic functions to calculate emergent constraint method #
# (based on the Cox et al., 2013 paper) #
# #
# BY : KM.Noh #
# DATE : 2023.03.15 #
# #
#------------------------------------------------------------------------------------#
## Modules for Calculate netCDF
import numpy as np
import xarray as xr
##########################################################################################
###### Emergent Constraint Method (Cox et al 2013)
##########################################################################################
### Calculate Gaussian PDF using mean and std
def calc_GAUSSIAN_PDF(mu,sigma,x):
PDF = 1/np.sqrt(2*np.pi*sigma**2) * np.exp(-1/2*((x-mu)/sigma)**2)
return PDF
### Calculate PDF of Constrained Projections
def calc_PDF_EC(tmp_x,tmp_y,x,y,PDF_x):
dx = x[1]-x[0]
### Inter-model Spread in CMIP Data
xn = tmp_x.values
yn = tmp_y.values
### Linear Regression for inter-model Diversity
N = len(xn)
sigma_x = np.sqrt(np.var(xn))
sigma_xy = np.sqrt(np.cov(xn,yn))[0,1]
b = (sigma_xy/sigma_x)**2
a = -1/N*(b*xn-yn).sum()
fn = a + b*xn
### Prediction Error
s = np.sqrt(1/(N-2)*((yn-fn)**2).sum())
sigma_b = (s/sigma_x) * np.sqrt(N)
### Calculate PDF(y) [posteriori] by given PDF(y|x) [Priori] and PDF(x) [OBS]
PDF_y = np.zeros(len(y))
for ind_y in range(len(y)):
for ind_x in range(len(x)):
sigma_fx = s * np.sqrt(1 + 1/N + (x[ind_x]-xn.mean())**2/(N*sigma_x**2))
PDF_yx = 1/np.sqrt(2*np.pi*sigma_fx**2) \
* np.exp(-(y[ind_y]-(a+b*x[ind_x]))**2/(2*sigma_fx**2))
PDF_y[ind_y] += PDF_yx * PDF_x[ind_x] * dx
thres = 0.341 ###
sigma_y = find_std_from_PDF(thres,y,PDF_y)
mean_y = y[PDF_y.argmax()]
return PDF_y, sigma_y, mean_y
### Calculate standard-deviation of Gaussian PDF
def find_std_from_PDF(thres,y,PDF):
ind_max = PDF.argmax()
for ind_point in range(len(y)):
PDF_sum = PDF[ind_point:ind_max+1].sum()/PDF.sum()
if (PDF_sum < thres):
sigma = y[ind_max]-y[ind_point]
break
return sigma
### Calculate the Prior of Probability in CMIP projections
def calc_PDF_EC_PRIOR(tmp_x,tmp_y,x,y):
xn = tmp_x.values
yn = tmp_y.values
N = len(xn)
sigma_x = np.sqrt(np.var(xn))
sigma_xy = np.sqrt(np.cov(xn,yn))[0,1]
b = (sigma_xy/sigma_x)**2
a = -1/N*(b*xn-yn).sum()
fn = a + b*xn
## prediction error
s = np.sqrt(1/(N-2)*((yn-fn)**2).sum())
sigma_b = (s/sigma_x) * np.sqrt(N)
sigma_fx = s * np.sqrt(1 + 1/N + (x-xn.mean())**2/(N*sigma_x**2))
PDF_yx = 1/np.sqrt(2*np.pi*sigma_fx**2) * np.exp(-(y-(a+b*x))**2/(2*sigma_fx**2))
return PDF_yx,sigma_fx, a+b*x