Skip to content
Snippets Groups Projects
Commit 4457a718 authored by Wuttke, Joachim's avatar Wuttke, Joachim
Browse files

bayesian example now reading data as 1D table

parent 53be1e0b
No related branches found
No related tags found
1 merge request!1980cleanup DataUtil; improve table reading
Pipeline #113432 passed
......@@ -57,9 +57,12 @@ def run_simulation(points, ni_thickness, ti_thickness):
if __name__ == '__main__':
filepath = os.path.join(datadir, "specular/genx_alternating_layers.dat.gz")
data = ba.readData2D(filepath).npArray()
data[:, 0] *= np.pi/360 # convert incident angles from deg to rad
data[:, 2] = data[:, 1]*0.1 # arbitrary uncertainties of 10%
flags = ba.ImportSettings1D("2alpha (deg)", "#", "", 1, 2)
data = ba.readData1D(filepath, ba.csv1D, flags)
q = data.npXcenters()
y = data.npArray()
dy = y * 0.1 # arbitrary uncertainties
def log_likelihood(P):
"""
......@@ -70,9 +73,6 @@ if __name__ == '__main__':
:array yerr: the ordinate uncertainty (dR-values)
:return: log-likelihood
"""
q = data[:, 0]
y = data[:, 1]
dy = data[:, 2]
y_sim = run_simulation(q, *P)
sigma2 = dy**2 + y_sim**2
return -0.5*np.sum((y - y_sim)**2/sigma2 + np.log(sigma2))
......@@ -104,14 +104,10 @@ if __name__ == '__main__':
plt.show()
# Plot and show MLE and data of reflectivity
plt.errorbar(data[:, 0],
data[:, 1],
data[:, 2],
marker='.',
ls='')
plt.errorbar(q, y, dy, marker='.', ls='')
plt.plot(
data[:, 0],
run_simulation(data[:, 0], *flat_samples.mean(axis=0)),
q,
run_simulation(q, *flat_samples.mean(axis=0)),
'-')
plt.xlabel('$\\alpha$/rad')
plt.ylabel('$R$')
......
......@@ -57,9 +57,12 @@ def run_simulation(points, ni_thickness, ti_thickness):
if __name__ == '__main__':
filepath = os.path.join(datadir, "specular/genx_alternating_layers.dat.gz")
data = ba.readData2D(filepath).npArray()
data[:, 0] *= np.pi/360 # convert incident angles from deg to rad
data[:, 2] = data[:, 1]*0.1 # arbitrary uncertainties of 10%
flags = ba.ImportSettings1D("2alpha (deg)", "#", "", 1, 2)
data = ba.readData1D(filepath, ba.csv1D, flags)
q = data.npXcenters()
y = data.npArray()
dy = y * 0.1 # arbitrary uncertainties
def log_likelihood(P):
"""
......@@ -70,9 +73,6 @@ if __name__ == '__main__':
:array yerr: the ordinate uncertainty (dR-values)
:return: log-likelihood
"""
q = data[:, 0]
y = data[:, 1]
dy = data[:, 2]
y_sim = run_simulation(q, *P)
sigma2 = dy**2 + y_sim**2
return -0.5*np.sum((y - y_sim)**2/sigma2 + np.log(sigma2))
......@@ -104,14 +104,10 @@ if __name__ == '__main__':
# plt.show()
# Plot and show MLE and data of reflectivity
plt.errorbar(data[:, 0],
data[:, 1],
data[:, 2],
marker='.',
ls='')
plt.errorbar(q, y, dy, marker='.', ls='')
plt.plot(
data[:, 0],
run_simulation(data[:, 0], *flat_samples.mean(axis=0)),
q,
run_simulation(q, *flat_samples.mean(axis=0)),
'-')
plt.xlabel('$\\alpha$/rad')
plt.ylabel('$R$')
......
......@@ -57,9 +57,12 @@ def run_simulation(points, ni_thickness, ti_thickness):
if __name__ == '__main__':
filepath = os.path.join(datadir, "specular/genx_alternating_layers.dat.gz")
data = ba.readData2D(filepath).npArray()
data[:, 0] *= np.pi/360 # convert incident angles from deg to rad
data[:, 2] = data[:, 1]*0.1 # arbitrary uncertainties of 10%
flags = ba.ImportSettings1D("2alpha (deg)", "#", "", 1, 2)
data = ba.readData1D(filepath, ba.csv1D, flags)
q = data.npXcenters()
y = data.npArray()
dy = y * 0.1 # arbitrary uncertainties
def log_likelihood(P):
"""
......@@ -70,9 +73,6 @@ if __name__ == '__main__':
:array yerr: the ordinate uncertainty (dR-values)
:return: log-likelihood
"""
q = data[:, 0]
y = data[:, 1]
dy = data[:, 2]
y_sim = run_simulation(q, *P)
sigma2 = dy**2 + y_sim**2
return -0.5*np.sum((y - y_sim)**2/sigma2 + np.log(sigma2))
......@@ -104,14 +104,10 @@ if __name__ == '__main__':
<%= sm ? "# " : "" %>plt.show()
# Plot and show MLE and data of reflectivity
plt.errorbar(data[:, 0],
data[:, 1],
data[:, 2],
marker='.',
ls='')
plt.errorbar(q, y, dy, marker='.', ls='')
plt.plot(
data[:, 0],
run_simulation(data[:, 0], *flat_samples.mean(axis=0)),
q,
run_simulation(q, *flat_samples.mean(axis=0)),
'-')
plt.xlabel('$\\alpha$/rad')
plt.ylabel('$R$')
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment