Commit c76a3039 authored by Maximilian Schanner's avatar Maximilian Schanner
Browse files

Lat, lon was mixed up in evaluation, reduced the sample size for faster examples.

parent 784d5977
...@@ -148,8 +148,7 @@ def location(posterior, mu_dips, cov_dips, n_points=10000, ...@@ -148,8 +148,7 @@ def location(posterior, mu_dips, cov_dips, n_points=10000,
The dipole location pdf, evaluated at the gridpoints The dipole location pdf, evaluated at the gridpoints
""" """
# set up grid to evaluate on # set up grid to evaluate on
(ti, tf), (pi, pf) = bounds theta, phi = zip(*equi_sph(n_points, bounds=bounds))
theta, phi = zip(*equi_sph(n_points, ti, tf, pi, pf))
phi = np.asarray(phi) phi = np.asarray(phi)
theta = np.pi/2. - np.asarray(theta) theta = np.pi/2. - np.asarray(theta)
# allocate array for the density # allocate array for the density
...@@ -176,8 +175,8 @@ def location(posterior, mu_dips, cov_dips, n_points=10000, ...@@ -176,8 +175,8 @@ def location(posterior, mu_dips, cov_dips, n_points=10000,
dens += posterior[np.unravel_index(it, posterior.shape)] \ dens += posterior[np.unravel_index(it, posterior.shape)] \
/ np.sqrt((2*np.pi)**3 * np.linalg.det(cov_dips[it])) \ / np.sqrt((2*np.pi)**3 * np.linalg.det(cov_dips[it])) \
* np.cos(theta) * zk / np.sqrt(zalpha**3) * dummy * np.cos(theta) * zk / np.sqrt(zalpha**3) * dummy
lat = np.rad2deg(phi) lat = np.rad2deg(theta)
lon = np.rad2deg(theta) lon = np.rad2deg(phi)
return lat, lon, dens return lat, lon, dens
...@@ -286,8 +285,7 @@ def coeffs(posterior, mu_coeffs, cov_coeffs, r_ref, r_at=REARTH): ...@@ -286,8 +285,7 @@ def coeffs(posterior, mu_coeffs, cov_coeffs, r_ref, r_at=REARTH):
cov_coeffs) cov_coeffs)
gm_weights = posterior / posterior.sum() gm_weights = posterior / posterior.sum()
ens = sample_GM(gm_weights.flatten(), mu_coeffs, cov_coeffs, ens = sample_GM(gm_weights.flatten(), mu_coeffs, cov_coeffs)
n_samps=50000)
err_16, err_84 = np.percentile(ens, (16, 84), axis=1) err_16, err_84 = np.percentile(ens, (16, 84), axis=1)
mean = (mu_coeffs.T * gm_weights.flatten()).sum(axis=1) mean = (mu_coeffs.T * gm_weights.flatten()).sum(axis=1)
...@@ -338,8 +336,7 @@ def spectrum(posterior, mu_coeffs, cov_coeffs, r_ref, r_at=REARTH): ...@@ -338,8 +336,7 @@ def spectrum(posterior, mu_coeffs, cov_coeffs, r_ref, r_at=REARTH):
gm_weights = posterior / posterior.sum() gm_weights = posterior / posterior.sum()
ens = sample_GM(gm_weights.flatten(), mu_coeffs, cov_coeffs, ens = sample_GM(gm_weights.flatten(), mu_coeffs, cov_coeffs)
n_samps=50000)
# 1st two moments of mixture # 1st two moments of mixture
mean = np.sum(mu_coeffs * gm_weights.flatten()[:, None], axis=0) mean = np.sum(mu_coeffs * gm_weights.flatten()[:, None], axis=0)
......
Markdown is supported
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