Hi !
I’m trying to work on MLE and MAP with pyro for my internship.
Here is my issue :
I’ve got 3 parameter (I0, beta, T) that I would like to maximize.
However my program doesn’t run as it should the optimizer : loss doesn’t reduce over iterations.
import math
import numpy as np
import matplotlib.pyplot as plt
import os
import torch
import torch.distributions.constraints as constraints
import pyro
from pyro.optim import Adam
from pyro.infer import SVI, Trace_ELBO
import pyro.distributions as dist
this is for running the notebook in our testing framework
smoke_test = “CI” in os.environ
n_steps = 2 if smoke_test else 2000
assert pyro.version.startswith(“1.8.1”)
clear the param store in case we’re in a REPL
pyro.clear_param_store()
Lambda_250 = torch.tensor(250e-6)
Lambda_360 = torch.tensor(360e-6)
Lambda_520 = torch.tensor(520e-6)
I0_vrai = torch.tensor(1.0)
beta_vrai = torch.tensor(1.0)
T_vrai = torch.tensor(18.0)
def f_beta(Lambda, T):
k = np.float64(1.38e-23) # constante de Boltztmann
h = np.float64(6.626e-34) # constante de Planck
c = np.float64(3e8) # célérité
return (
2
* h
* pow(c, 2)
/ pow(Lambda, 5)
* 1
/ (torch.exp(h * c / (Lambda * k * T)) - 1)
)
def I(I0, beta, T, Lambda):
return I0 * torch.pow(Lambda, -beta) * f_beta(Lambda, T)
def H(I0, beta, T):
return torch.tensor(
[
I(I0, beta, T, Lambda_250),
I(I0, beta, T, Lambda_360),
I(I0, beta, T, Lambda_520),
]
)
data = torch.tensor(
[
I(I0_vrai, beta_vrai, T_vrai, Lambda_250),
I(I0_vrai, beta_vrai, T_vrai, Lambda_360),
I(I0_vrai, beta_vrai, T_vrai, Lambda_520),
]
)
noise = torch.div(data, 100, rounding_mode=“trunc”)
data bruitées
data = torch.tensor(
[
dist.Normal(I(I0_vrai, beta_vrai, T_vrai, Lambda_250), noise[0]).sample(),
dist.Normal(I(I0_vrai, beta_vrai, T_vrai, Lambda_360), noise[1]).sample(),
dist.Normal(I(I0_vrai, beta_vrai, T_vrai, Lambda_520), noise[2]).sample(),
]
)
écart type des lois normales
std = torch.tensor(100.0)
Déinition du modèle
def model(data):
# sample f from the Uniform prior
I0 = pyro.sample("I0", dist.Uniform(0, 5))
beta = pyro.sample("beta", dist.Uniform(0, 2))
T = pyro.sample("T", dist.Uniform(5, 100))
# loop over the observed data
with pyro.plate("data", data.size(0)):
pyro.sample("obs", dist.Normal(H(I0, beta, T), std), obs=data)
def guide(data):
# register the three variational parameters with Pyro.
I0_q = pyro.param(“I0_q”, torch.tensor(1), constraint=constraints.positive)
beta_q = pyro.param(
“beta_q”, torch.tensor(1), constraint=constraints.interval(0.0, 2.0)
)
T_q = pyro.param(
“T_q”, torch.tensor(18), constraint=constraints.interval(5, 100)
)
# sample latent_intesity from the distribution Normale(I(I0_q, beta_q, T_q)
pyro.sample(“latent_intensity”, dist.Normal(H(I0_q, beta_q, T_q), std))
autoguide = pyro.infer.autoguide.AutoDelta(model)
Set up de la stochastic variational inference (svi)
Set up de la stochastic variational inference (svi)
def train(model, guide, lr=0.005, n_steps=201):
pyro.clear_param_store()
adam_params = {“lr”: lr}
adam = pyro.optim.Adam(adam_params)
svi = SVI(model, guide, adam, loss=Trace_ELBO())
for step in range(n_steps):
loss = svi.step(data)
if step % 50 == 0:
print("[iter {}] loss: {:.4f}".format(step, loss))
train(model, autoguide)
grab the learned variational parameters
I0_q = pyro.param(“I0_q”).item()
beta_q = pyro.param(“beta_q”).item()
T_q = pyro.param(“T_q”).item()
print(I0_q, beta_q, T_q)
print(
“Our estimated parameter I0 is {:.3f}”.format(autoguide.median(data)[“I0”].item())
)
print(
“Our estimated parameter beta is {:.3f}”.format(
autoguide.median(data)[“beta”].item()
)
)
print(“Our estimated parameter T is {:.3f}”.format(autoguide.median(data)[“T”].item()))
I0_est = torch.tensor(autoguide.median(data)[“I0”].item())
beta_est = torch.tensor(autoguide.median(data)[“beta”].item())
T_est = torch.tensor(autoguide.median(data)[“T”].item())
def model_mle(data):
# sample f from the Uniform prior
# I0_mle = pyro.param("I0_mle", I0_est, constraint=constraints.positive,)
# beta_mle = pyro.param("beta_mle", beta_est, constraint=constraints.positive,)
# T_mle = pyro.param("T_mle", T_est, constraint=constraints.positive,)
I0_mle = pyro.param("I0_mle", torch.tensor(2.0), constraint=constraints.positive,)
beta_mle = pyro.param(
"beta_mle", torch.tensor(0.8), constraint=constraints.positive,
)
T_mle = pyro.param("T_mle", torch.tensor(35.0), constraint=constraints.positive,)
# loop over the observed data
with pyro.plate("data", data.size(0)):
pyro.sample("obs", dist.Normal(H(I0_mle, beta_mle, T_mle), std), obs=data)
def guide_mle(data):
pass
train(model_mle, guide_mle)
mle_estimate_I0 = pyro.param(“I0_mle”).item()
mle_estimate_beta = pyro.param(“beta_mle”).item()
mle_estimate_T = pyro.param(“T_mle”).item()
print(“Our MLE estimate of the I0 parameter is {:.3f}”.format(mle_estimate_I0))
print(“Our MLE estimate of the beta parameter is {:.3f}”.format(mle_estimate_beta))
print(“Our MLE estimate of the T parameter is {:.3f}”.format(mle_estimate_T))
Have you any idea ?