Sum of Gaussian Distributions + noise.ipynb

Kernel: Unknown Kernel

import pymc as pm
import numpy as np
import pylab as pl


def GaussFunc(x, amplitude, centroid, sigma):
    return amplitude * np.exp(-0.5 * ((x - centroid) / sigma) ** 2)

wavelength = np.arange(5000, 5050, 0.02)

# Profile 1
centroid_one = 5025.0
sigma_one = 2.2
height_one = 0.8
profile1 = GaussFunc(wavelength, height_one, centroid_one, sigma_one, )

# Profile 2
centroid_two = 5027.0
sigma_two = 1.2
height_two = 0.5
profile2 = GaussFunc(wavelength, height_two, centroid_two, sigma_two, )

# Measured values
noise = np.random.normal(0.0, 0.02, len(wavelength))
combined = profile1 + profile2 + noise

# If you want to plot what this looks like
pl.plot(wavelength, combined, label="Measured")
pl.plot(wavelength, profile1, color='red', linestyle='dashed', label="1")
pl.plot(wavelength, profile2, color='green', linestyle='dashed', label="2")
pl.title("Feature One and Two")
pl.legend()

<matplotlib.legend.Legend at 0x6075af0>

# Unknowns

est_centroid_one = pm.Uniform("est_centroid_one", 5000, 5050)
est_centroid_two = pm.Uniform("est_centroid_two", 5000, 5050)

est_sigma_one = pm.Uniform("est_sigma_one", 0, 5)
est_sigma_two = pm.Uniform("est_sigma_two", 0, 5)

est_height_one = pm.Uniform("est_height_one", 0, 5)
est_height_two = pm.Uniform("est_height_two", 0, 5)

std_deviation = 1. / mc.Uniform("std", 0, 1) ** 2

# Set up the inference


@pm.deterministic(trace=False)
def est_profile_1(x=wavelength, centroid=est_centroid_one, sigma=est_sigma_one, height=est_height_one):
    return GaussFunc(x, height, centroid, sigma)


@pm.deterministic(trace=False)
def est_profile_2(x=wavelength, centroid=est_centroid_two, sigma=est_sigma_two, height=est_height_two):
    return GaussFunc(x, height, centroid, sigma)


@pm.deterministic(trace=False)
def mean(profile_1=est_profile_1, profile_2=est_profile_2):
    return profile_1 + profile_2


observations = pm.Normal("obs", mean, std_deviation, value=combined, observed=True)

model = pm.Model([est_centroid_one,
                  est_centroid_two,
                  est_height_one,
                  est_height_two,
                  est_sigma_one,
                  est_sigma_two,
                  std_deviation])

map_ = pm.MAP(model)
map_.fit()

mcmc = pm.MCMC(model)
mcmc.sample(70000, 60000)

[****************100%******************]  70000 of 70000 complete

mcplot(mcmc)

Plotting est_height_one
Plotting est_sigma_two
Plotting est_height_two
Plotting est_centroid_two
Plotting est_centroid_one
Plotting est_sigma_one
Plotting std

mcplot?