<jupyter_start><jupyter_text>Setting up the model$$\cos (\Delta I) = \cos i_\mathrm{disk} \cos i_\star + \sin i_\mathrm{disk} \sin i_\star \cos(\Omega_\mathrm{disk} - \Omega_\star)$$The prior probability distributions are$$p(i_\mathrm{disk}, \Omega_\mathrm{disk}) = \frac{\sin i_\mathrm{disk}}{2} $$<jupyter_code># Define the likelihood by conditioning on the observations 
data = np.array([np.cos(33.5 * deg), 256.0 * deg, np.cos(33.42 * deg)])
    
# converting the uncertainty requires a change of variables 
# sigma_cos_i = sin(i) * sigma_i (where sigma_i is in radians)
sigma = np.array([ np.sin(33.5 * deg) * 1.4 * deg, 1.0 * deg, np.sin(33.42 * deg) * 0.01 * deg])
print(data)
print(sigma)
# instantiate a PyMC3 model class 
with pm.Model() as model:
    
    # priors on i_disk, Omega_disk, i_star, Omega_star
    # everything's assumed to be in radians
    # a uniform prior on cos(i) means p(i) = sin(i)
    cos_i_disk = pm.Uniform("cosIdisk", lower=-1.0, upper=1.0)
    Omega_disk = pm.Uniform("OmegaDisk", lower=0.0, upper=2*np.pi)
     
    cos_i_star = pm.Uniform("cosIstar", lower=-1.0, upper=1.0)
    Omega_star = pm.Uniform("OmegaStar", lower=0.0, upper=2*np.pi)
    
    # calculate the intermediate quantities for determining mutual inclination
    # to go from cos(i) -> sin(i), do 
    sin_i_disk = tt.sqrt(1.0 - cos_i_disk**2)
    sin_i_star = tt.sqrt(1.0 - cos_i_star**2)
        
    # store Delta I as a deterministic byproduct
    cos_Delta_I = pm.Deterministic("cosDeltaI", cos_i_disk * cos_i_star + \
                               sin_i_disk * sin_i_star * tt.cos(Omega_disk - Omega_star))

    # assemble a vector to compare to the observations
    mu = [cos_i_disk, Omega_disk, cos_i_star]
    obs = pm.Normal("obs", mu=mu, sd=sigma, observed=data)
    
model.free_RVs
model.deterministics
# We can test the model likelihood if we don't let PyMC3 automatically transform the variables
# transform=None in the prior definitions. 
# obs.logp({"cos_i_disk":0.833, "Omega_disk":4.46, "cos_i_star":0.83, "Omega_star":4.0})
# sample the model!
with model:
    trace = pm.sample(draws=1000, tune=2000, chains=2)
# visualize output
pm.summary(trace)
pm.traceplot(trace)
pm.model_to_graphviz(model)
# get the samples on cosDeltaI and map these into Delta I. Large cos(Delta I) values are actually small Delta I 
# values.
pm.__version__
import corner  

samples = pm.trace_to_dataframe(trace, varnames=["cosIdisk", "OmegaDisk", "cosIstar", "OmegaStar", "cosDeltaI"])
# convert to actual parameter space


# corner.corner(samples) #, truths=[true_m, true_b, true_logs]);
samples.head()
asamples = samples.copy()
asamples.rename(index=str, columns={"cosIdisk":"iDisk", "cosIstar":"iStar", "cosDeltaI":"deltaI"}, inplace=True);
asamples.head()
icos2deg = lambda x: np.arccos(x)/deg

asamples["iDisk"] = asamples["iDisk"].apply(icos2deg)
asamples["OmegaDisk"] = asamples["OmegaDisk"] / deg
asamples["iStar"] = asamples["iStar"].apply(icos2deg)
asamples["OmegaStar"] = asamples["OmegaStar"] / deg
asamples["deltaI"] = asamples["deltaI"].apply(icos2deg)
asamples.head()
corner.corner(asamples) #, truths=[true_m, true_b, true_logs]);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14721610569361435, 'fraction_numerical': 0.022963195973576597, 'mean_word_length': 4.435897435897436, 'pattern_counts': {'":': 7, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '282792', 'license': 'MIT', 'n_tokens_mistral': 1177, 'n_tokens_neox': 1124, 'n_words': 276, 'path': 'PyMC3 sample one disk.ipynb', 'repo_name': 'iancze/hierarchical-mutual-inclinations', 'size': '3179'}