global-warming slides

### Global Warming¶

#### Example for Bayesian Linear Regression¶

Assumption:

• Factors of influence for the global warming are greenhouse gases and solar irridiance.
• Super simple linear model.

### Data for atmospheric CO$_2$ concentration¶

As greenhouse gas we consider only CO$_2$. But the increase of CO$_2$ in the atmosphere should correlate strongly with the emission of other greenhouse gases.

In [3]:
# data from http://www.esrl.noaa.gov/gmd/ccgg/trends/
# since 1959 direct measurements!
# set year as index
co2_concentration = co2_concentration.set_index("Year")
co2_concentration.tail(3)

Out[3]:
EPICA Dome C and Vostok Station, Antarctica Law Dome, Antarctica (75-year smoothed) Siple Station, Antarctica Mauna Loa, Hawaii Barrow, Alaska Cape Matatula, American Samoa South Pole, Antarctica Cape Grim, Australia Lampedusa Island, Italy Shetland Islands, Scotland
Year
2012 NaN NaN NaN 393.82 394.940000 391.568333 390.060000 NaN NaN NaN
2013 NaN NaN NaN 396.48 397.928333 394.257500 392.779167 NaN NaN NaN
2014 NaN NaN NaN 398.55 NaN NaN NaN NaN NaN NaN

Use Pandas for computing the mean of each row. NaN values are ignored:

In [5]:
# set new mean series
co2_concentration['mean'] = co2_concentration.mean(axis=1)

In [6]:
co2_concentration.tail(3)

Out[6]:
EPICA Dome C and Vostok Station, Antarctica Law Dome, Antarctica (75-year smoothed) Siple Station, Antarctica Mauna Loa, Hawaii Barrow, Alaska Cape Matatula, American Samoa South Pole, Antarctica Cape Grim, Australia Lampedusa Island, Italy Shetland Islands, Scotland mean
Year
2012 NaN NaN NaN 393.82 394.940000 391.568333 390.060000 NaN NaN NaN 392.597083
2013 NaN NaN NaN 396.48 397.928333 394.257500 392.779167 NaN NaN NaN 395.361250
2014 NaN NaN NaN 398.55 NaN NaN NaN NaN NaN NaN 398.550000
In [7]:
plt.plot(co2_concentration.index, co2_concentration['mean'])
plt.xlim((1400, 2020))
plt.xlabel("Year")
plt.ylabel("CO$_2$ concentration / ppm")

Out[7]:
<matplotlib.text.Text at 0x10b3b9e50>
In [12]:
# convert c02 concentation to numpy
co2_year = co2_concentration.index.get_values()
co2_mean = co2_concentration['mean'].get_values()


In [ ]:
lean_df = pd.read_csv("./lean2000_irradiance_only_data.csv", sep=',')
#lean_df

In [15]:
# solar irridiance
year_mask = (lean_df["YEAR"]- 0.5 >= df_air.Year.min())
scy = "11yrCYCLE+BKGRND"
#scy = "11yrCYCLE"

smo_ = np.ndarray(len(lean_df.YEAR))
import pandas.stats.moments
smoothed = pandas.stats.moments.ewma(lean_df[scy], halflife=3.)

# according to a graph in http://www.climate4you.com/Sun.htm
# the irridience goes sligtly down sice 2000, to be more conservative we leave it constant
# from 2000-2014
plt.xlabel("Year")

/Users/christian/.virtualenvs/alternative_1/lib/python2.7/site-packages/ipykernel/__main__.py:10: FutureWarning: pd.ewm_mean is deprecated for Series and will be removed in a future version, replace with

Out[15]:
<matplotlib.text.Text at 0x10b59ee90>

#### Global air temperature¶

In [16]:
#mean = (df_air["J-D"]*0.01).mean()
temp_dev_air = -0.2 + ((df_air["J-D"]-df_air["J-D"].min())*0.01).get_values()
#temp_dev_lo = ((df_land_and_ocean["J-D"] - (df_land_and_ocean["J-D"]).min())*0.01).get_values()

In [18]:
plot_tempe()

In [26]:
plt.plot(co2_original, alo_temp, 'bo')
plt.xlabel("CO$_2$ / ppm")
plt.ylabel("$\Delta$ temperature / ${}^\circ$C")

Out[26]:
<matplotlib.text.Text at 0x10b5f1dd0>
In [27]:
l = (year_mask).sum()
plt.xlabel("solar irridiance / ?")
plt.ylabel("$\Delta$ temperature / ${}^\circ$C")

Out[27]:
<matplotlib.text.Text at 0x10ba22ed0>
In [28]:
import pymc3
model = pymc3.Model()
with model:
co_2_weight = pymc3.Normal("co_2_weight", mu=0, sd=1.)
irridiation_weight = pymc3.Normal("irridiaton_weight", mu=0, sd=1.)

offset = pymc3.Normal("offset", mu=0, sd=1.)

i_year = year - year[0]
i_co2_year = year - co2_year[0] # + co_2_delay
temp_mu = offset + irridiation_weight * smo[i_year] + co_2_weight * co2_
temp_sigma = pymc3.HalfNormal('sigma', sd=0.5)

temp_dev_ = pymc3.Normal("temp_dev_", mu=temp_mu, sd=temp_sigma, observed = alo_temp)
temp_pred = pymc3.Normal("temp_pred", mu=temp_mu, sd=temp_sigma, shape = alo_temp.shape )

Couldn't import dot_parser, loading of dot files will not be possible.
Applied log-transform to sigma and added transformed sigma_log to model.

In [30]:
with model:
# obtain starting values via MAP
start = pymc3.find_MAP()#fmin=optimize.fmin_powell)

#step = pymc3.Slice()
tr = pymc3.sample(5000, start=start)

Assigned NUTS to co_2_weight
Assigned NUTS to irridiaton_weight
Assigned NUTS to offset
Assigned NUTS to sigma_log
Assigned NUTS to temp_pred
[-----------------100%-----------------] 5000 of 5000 complete in 7.2 sec
In [32]:
plt.figure(figsize=(7, 7))
pymc3.traceplot(tr)
plt.tight_layout();

<matplotlib.figure.Figure at 0x1143897d0>

Compare the influence of irridiation and CO$_2$(Normalized Data):

In [34]:
irridiation_weight_.mean()

Out[34]:
0.045895907084405878
In [35]:
co_2_weight_.mean()

Out[35]:
0.33630017388761918
In [36]:
plt.figure(figsize=(12,10))
temp_pred_ = tr.get_values(temp_pred)[4000:]
for t in temp_pred_:
plt.plot(year, t, 'b.', alpha=0.02)
plt.xlabel("Year")
plt.ylabel("$\Delta T$ / ${}^\circ$C")
plt.plot(year, alo_temp, 'k.')

Out[36]:
[<matplotlib.lines.Line2D at 0x113cc7b10>]
In [37]:
plt.figure(figsize=(12,10))
for t in temp_.T:
plt.plot(year, t, 'b-', alpha=0.04)
plt.xlabel("Year")
plt.ylabel("$\Delta T$ / ${}^\circ$C")
plt.plot(year, alo_temp, 'k.')

Out[37]:
[<matplotlib.lines.Line2D at 0x113cc7b90>]