Commit 0442b931 authored by Philipp Arras's avatar Philipp Arras
Browse files

Fixes

parent ee9f5673
Pipeline #70589 failed with stages
in 11 minutes and 13 seconds
......@@ -67,7 +67,7 @@ if __name__ == '__main__':
# Compute likelihood and Hamiltonian
position = ift.from_random('normal', harmonic_space)
likelihood = ift.BernoulliEnergy(data)(p)
likelihood = ift.BernoulliEnergy(data) @ p
ic_newton = ift.DeltaEnergyController(
name='Newton', iteration_limit=100, tol_rel_deltaE=1e-8)
minimizer = ift.NewtonCG(ic_newton)
......
......@@ -96,7 +96,7 @@ if __name__ == '__main__':
data = lamb(mock_position)
data = np.random.poisson(data.val.astype(np.float64))
data = ift.Field.from_raw(d_space, data)
likelihood = ift.PoissonianEnergy(data)(lamb)
likelihood = ift.PoissonianEnergy(data) @ lamb
# Settings for minimization
ic_newton = ift.DeltaEnergyController(
......
......@@ -122,8 +122,7 @@ if __name__ == '__main__':
N_samples = 20
# Set up likelihood and information Hamiltonian
likelihood = ift.GaussianEnergy(
mean=data, inverse_covariance=N.inverse)(signal_response)
likelihood = ift.GaussianEnergy(mean=data, inverse_covariance=N.inverse) @ signal_response
H = ift.StandardHamiltonian(likelihood, ic_sampling)
# Begin minimization
......
......@@ -81,63 +81,63 @@ class PolynomialResponse(ift.LinearOperator):
return ift.makeField(self._tgt(mode), out)
# Generate some mock data
N_params = 10
N_samples = 100
size = (12,)
x = np.random.random(size) * 10
y = np.sin(x**2) * x**3
var = np.full_like(y, y.var() / 10)
var[-2] *= 4
var[5] /= 2
y[5] -= 0
# Set up minimization problem
p_space = ift.UnstructuredDomain(N_params)
params = ift.full(p_space, 0.)
R = PolynomialResponse(p_space, x)
ift.extra.consistency_check(R)
d_space = R.target
d = ift.makeField(d_space, y)
N = ift.DiagonalOperator(ift.makeField(d_space, var))
IC = ift.DeltaEnergyController(tol_rel_deltaE=1e-12, iteration_limit=200)
likelihood = ift.GaussianEnergy(d, N)(R)
Ham = ift.StandardHamiltonian(likelihood, IC)
H = ift.EnergyAdapter(params, Ham, want_metric=True)
# Minimize
minimizer = ift.NewtonCG(IC)
H, _ = minimizer(H)
# Draw posterior samples
metric = Ham(ift.Linearization.make_var(H.position, want_metric=True)).metric
samples = [metric.draw_sample(from_inverse=True) + H.position
for _ in range(N_samples)]
# Plotting
plt.errorbar(x, y, np.sqrt(var), fmt='ko', label='Data with error bars')
xmin, xmax = x.min(), x.max()
xs = np.linspace(xmin, xmax, 100)
sc = ift.StatCalculator()
for ii in range(len(samples)):
sc.add(samples[ii])
ys = polynomial(samples[ii], xs)
if ii == 0:
plt.plot(xs, ys, 'k', alpha=.05, label='Posterior samples')
continue
plt.plot(xs, ys, 'k', alpha=.05)
ys = polynomial(H.position, xs)
plt.plot(xs, ys, 'r', linewidth=2., label='Interpolation')
plt.legend()
plt.savefig('fit.png')
plt.close()
# Print parameters
mean = sc.mean.val
sigma = np.sqrt(sc.var.val)
for ii in range(len(mean)):
print('Coefficient x**{}: {:.2E} +/- {:.2E}'.format(ii, mean[ii],
sigma[ii]))
if __name__ == '__main__':
# Generate some mock data
N_params = 10
N_samples = 100
size = (12,)
x = np.random.random(size) * 10
y = np.sin(x**2) * x**3
var = np.full_like(y, y.var() / 10)
var[-2] *= 4
var[5] /= 2
y[5] -= 0
# Set up minimization problem
p_space = ift.UnstructuredDomain(N_params)
params = ift.full(p_space, 0.)
R = PolynomialResponse(p_space, x)
ift.extra.consistency_check(R)
d_space = R.target
d = ift.makeField(d_space, y)
N = ift.DiagonalOperator(ift.makeField(d_space, var))
IC = ift.DeltaEnergyController(tol_rel_deltaE=1e-12, iteration_limit=200)
likelihood = ift.GaussianEnergy(d, N) @ R
Ham = ift.StandardHamiltonian(likelihood, IC)
H = ift.EnergyAdapter(params, Ham, want_metric=True)
# Minimize
minimizer = ift.NewtonCG(IC)
H, _ = minimizer(H)
# Draw posterior samples
metric = Ham(ift.Linearization.make_var(H.position, want_metric=True)).metric
samples = [metric.draw_sample(from_inverse=True) + H.position
for _ in range(N_samples)]
# Plotting
plt.errorbar(x, y, np.sqrt(var), fmt='ko', label='Data with error bars')
xmin, xmax = x.min(), x.max()
xs = np.linspace(xmin, xmax, 100)
sc = ift.StatCalculator()
for ii in range(len(samples)):
sc.add(samples[ii])
ys = polynomial(samples[ii], xs)
if ii == 0:
plt.plot(xs, ys, 'k', alpha=.05, label='Posterior samples')
continue
plt.plot(xs, ys, 'k', alpha=.05)
ys = polynomial(H.position, xs)
plt.plot(xs, ys, 'r', linewidth=2., label='Interpolation')
plt.legend()
plt.savefig('fit.png')
plt.close()
# Print parameters
mean = sc.mean.val
sigma = np.sqrt(sc.var.val)
for ii in range(len(mean)):
print('Coefficient x**{}: {:.2E} +/- {:.2E}'.format(ii, mean[ii], sigma[ii]))
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment