bondscell_results $daa349a6-545b-4dfc-a3b0-92884da98e5equeued¤logsrunning¦outputbody G
Luminosity distance from Distance modulus
The standard definition of the distance modulus is $μ = 5 \log_{10}\left(\frac{d}{\mathrm{pc}}\right) - 5$. We'll work in megaparsecs (Mpc), so we get $μ = 5 \log_{10}\left(\frac{d}{\mathrm{Mpc}}\right) + 25$. Solving for $d$ gives $d = 10^{(\mu-25)/5}$ Since we're measuring distances based on the luminosity of supernovae, we're using a luminosity distance here. Hence I've included the subscript L.)
Computes the maximum likelihood estimator for θ for the linear model y = A θ where measurements errors of y_obs are normally distributed with covariance matrix covar about the true values of y.
mimetext/htmlrootassigneelast_run_timestampAGpersist_js_state·has_pluto_hook_features§cell_id$556d8f17-8a03-42c8-82c7-52c8125b604ddepends_on_disabled_cells§runtime õpublished_object_keysdepends_on_skipped_cells§errored$77b17593-6a28-4294-9b17-41303d57d5a9queued¤logsrunning¦outputbody%loss (generic function with 1 method)mimetext/plainrootassigneelast_run_timestampA8 Xpersist_js_state·has_pluto_hook_features§cell_id$77b17593-6a28-4294-9b17-41303d57d5a9depends_on_disabled_cells§runtime #published_object_keysdepends_on_skipped_cells§errored$2482544e-ad05-495c-9787-0797939445bequeued¤logsrunning¦outputbody٦
It incorporated data from 833 supernovae (SNe), drawn from 19 datasets. Of these, 580 SNe passed usability cuts. We present new data from the HST Cluster Survey. All SNe were fit using a single lightcurve fitter (SALT2-1) and uniformly analyzed. All analyses and cuts were developed in a blind manner, i.e. with the cosmology hidden. The Union2.1 paper is available at Suzuki et al. (2012).
The data are stored in a simple text file. There are many ways to read them. Julia has a fast reader for CSV files, so we'll demonstrate that below. For now, you don't need to pay attention to the syntax, but at some point you'll likely want to refer back to examples like this to see how to accomplish common tasks.
A diagonal matrix containing the squared measurement uncertainties, $\Sigma = \mathrm{diag}(\sigma_{d_{L}}^{2}) = \begin{bmatrix} \sigma_{d,1}^{-2} & & \\
& \sigma_{d,2}^{-2} & & \\
& & \ddots & \\
& & & \sigma_{d,N_{\mathrm{obs}}}^{-2} \end{bmatrix},$ and the design matrix contains the velocities $A = \left( \begin{matrix} v_1 \\
v_2 \\
... \\
v_{N_{\mathrm{obs}}} \end{matrix} \right).$ In this model, we're presuming that the velocities are known exactly. While this isn't technically true, the velocities are typically measured much more precisely than the distances, so it can be a reasonable approximation.
After we've specified the general form of the statistcal model, we create a specific version of the model that is specialized for the valuse of the observations.
mimetext/htmlrootassigneelast_run_timestampA-rw)persist_js_state·has_pluto_hook_features§cell_id$92b2c9f6-5102-4e1a-aff1-05e9ce8c19a8depends_on_disabled_cells§runtime Ipublished_object_keysdepends_on_skipped_cells§errored$6c86d0f1-5c6f-4204-9414-60577164f4b8queued¤logsrunning¦outputbodyprefixٹDynamicPPL.Model{typeof(model_hubble_law_lognormal_prior), (:v, :d_obs, :σ_d_obs), (), (), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}elementsfBmodel_hubble_law_lognormal_prior (generic function with 2 methods)text/plainargselementsvprefixFloat64elements8418.88text/plain14627.6text/plain15447.5text/plain20276.6text/plain18199.8text/plainmore57893.3text/plaintypeArrayprefix_shortobjectidcfe55747b0413d3e!application/vnd.pluto.tree+objectd_obsprefixFloat64elements117.305text/plain217.007text/plain230.961text/plain308.565text/plain313.821text/plainmore1291.84text/plaintypeArrayprefix_shortobjectid14301e3138d6cedc!application/vnd.pluto.tree+objectσ_d_obsprefixFloat64elements12.0956text/plain16.6721text/plain16.5664text/plain22.5181text/plain22.5595text/plainmore145.219text/plaintypeArrayprefix_shortobjectid6af18a9b6659a5c3!application/vnd.pluto.tree+objecttypeNamedTupleobjectid59d0421f6398110f!application/vnd.pluto.tree+objectdefaultselementstypeNamedTupleobjectidffffffff412bd36e!application/vnd.pluto.tree+objectcontextprefixDynamicPPL.DefaultContextelementstypestructprefix_shortDefaultContextobjectidffffffff9dc0ec1d!application/vnd.pluto.tree+objecttypestructprefix_shortModelobjectidd30f1722f056fac0mime!application/vnd.pluto.tree+objectrootassignee+model_hubble_law_lognormal_prior_given_datalast_run_timestampAF@!persist_js_state·has_pluto_hook_features§cell_id$6c86d0f1-5c6f-4204-9414-60577164f4b8depends_on_disabled_cells§runtime LUpublished_object_keysdepends_on_skipped_cells§errored$a7981271-0b12-4b75-b5f0-5e91af761419queued¤logsrunning¦outputbody
In this case, we had a very simple model and the linear algebra was straight forward. In more complicated cases, there could be many input variables (here velocity), multiple observables being predicted (here distance), correlations between the different observables, and many possible models to try. Setting up all the matrices and getting all the algebra right for each model is error prone. Therefore, most modern high-level programming languages have packages that make it easy to fit linear models such as this one. Below, we'll repeat the analysis, but using Julia's GLM.jl package for general linear models.
If one puts in a fair bit of algebra, then one can find an analytical expression for the value of the model parameters ($\theta$) that minimizes the likelihood for a linear model.
You're likely to see the derivation of the solution in class on staistics (or maybe linear algebra), but we'll just use the result:
where $y_{\mathrm{obs}}$ is the vector of measurement values, $\Sigma^{-1}$ is the covariance matrix, and $A$ is the design matrix. For our problem, these simplify to:
Notice that this code is very different from a tradditional program where we'd specify how to compute each variable in order (known as the imperative programming model).
Q1e: How does the best-fit value computed above compare to your previous estimate based on the plot of $\chi^2$ vs $H_0$? What is likely to explain any differences?
Does Hubble's law provide a good model for these observations? If not, explain your concern.
The Hubble "law" ($v = H_0 \times d$) for the expansion of the universe is based on velocity versus distance. That's not what this data file provided, but we can compute those using standard expressions that you may have seen in an introduction to astronomy class. You learn or refresh your memory about each by clicking below.
mimetext/htmlrootassigneelast_run_timestampA-\persist_js_state·has_pluto_hook_features§cell_id$92ca099c-b81b-4cfe-a3fa-60833657f4efdepends_on_disabled_cells§runtime ;published_object_keysdepends_on_skipped_cells§errored$03fc47ea-9246-4cb8-8398-612fc2fb09faqueued¤logsrunning¦outputbodymissingmimetext/plainrootassigneeresponse_1blast_run_timestampA8^Wpersist_js_state·has_pluto_hook_features§cell_id$03fc47ea-9246-4cb8-8398-612fc2fb09fadepends_on_disabled_cells§runtime0published_object_keysdepends_on_skipped_cells§errored$d3db2bc4-a26b-4e01-a193-b155a92dab48queued¤logsrunning¦outputbodyBmodel_hubble_law_lognormal_prior (generic function with 2 methods)mimetext/plainrootassigneelast_run_timestampA:>)persist_js_state·has_pluto_hook_features§cell_id$d3db2bc4-a26b-4e01-a193-b155a92dab48depends_on_disabled_cells§runtimeRʐpublished_object_keysdepends_on_skipped_cells§errored$9db3ff24-6f54-4a37-b787-de4b14694d95queued¤logsrunning¦outputbody
calc_fisher_matrix_linear_model
calc_fisher_matrix_linear_model(A, covar)
Computes the Fisher information matrix for θ for the linear model y = A θ where measurements errors of y_obs are normally distributed and have covariance covar. Note that y_obs is not passed to this function, because it does not affect results.
mimetext/htmlrootassigneelast_run_timestampA8lpersist_js_state·has_pluto_hook_features§cell_id$0c0312da-558a-4bc1-9e4a-b02ecd1c0572depends_on_disabled_cells§runtimegRpublished_object_keysdepends_on_skipped_cells§errored$b15b1913-d265-4a35-b43e-67df53c93329queued¤logsrunning¦outputbodyLHow to access metadata
You can access the metadata from a DataFrame using function calls like metadata(df_in) or colmetadata(df_in,:z).
As we develop more complex models, it will often not be practical to find the maximum likelihood estimator values analyticly or using a general linear modeling package. Or there may be significant prior knowledge that we'd like to incorporate into our analysis. Or we might want to perform a more rigorous analysis of the uncertainties on the model parameters. In any of these cases, it can be very useful to use a probabilistic modeling language (PPL), such as Turing or Stan.
Q2a: How do the mean and standard deviation for the Hubble constant estimated with our Bayesian model and Turing compare to the analytic estimates using linear model (and the same dataset)? Are any differences large or small compare to the reported uncertainties?
This value is nothing like what we got above. Why? Remember that to write this in a form ammenable to the analytic solution, we predicted distances based on velocities $d = v / H_0$, so the best-fit constant of proportionality is the recipocal of $H_0$. Let's transform that back to a best-fit value for the Hubble constant, $\hat{H_0}_{,bf} = 1/\hat{\theta}_{bf}$.
You likely noticed that Hubble's "Law" wasn't a great model for this dataset. That's because it's an incomplete model that works well in the local universe, but breaks down once astronomers look far enough away that we can see the effects of dark energy causing the rate of expansion to accelerate. If we want to measure the local Hubble constant, then we could fit a subset of our dataset, picking those SN that have velocities less than some threshold.
Minimizing $\mathcal{L}$ is equivalent to minimizing $\ell = \log(\mathcal{L})$. Since the first factor doesn't depend on the measured values of distance, we can minimize $\mathcal{L}$ and $\ell$ by minimizing $\chi^2$.
mimetext/htmlrootassigneelast_run_timestampA-rTpersist_js_state·has_pluto_hook_features§cell_id$a56871c9-affa-4612-b3f5-66363720f274depends_on_disabled_cells§runtime ۵published_object_keysdepends_on_skipped_cells§errored$4f27d8d1-0d14-49e5-abee-9f277369a78bqueued¤logsrunning¦outputbody@model_hubble_law_uniform_prior (generic function with 2 methods)mimetext/plainrootassigneelast_run_timestampA:VUpersist_js_state·has_pluto_hook_features§cell_id$4f27d8d1-0d14-49e5-abee-9f277369a78bdepends_on_disabled_cells§runtime S4published_object_keysdepends_on_skipped_cells§errored$0144843f-b649-4ea2-a546-e8c14a928670queued¤logsrunning¦outputbody[
Q2b: How do the estimates of the posterior mean and standard deviation using two different choices for the prior compare to each other? Are any differences large or small compare to the reported uncertainties? What does this imply about the sensitivity of our analysis to the choice of prior?
Q1b: Do you notice anything potentially concerning about the plot? If so, describe your concerns. Optionally, use the cell below to perform any tests to address your concerns.
Traditionally, programmers provide explicit instructions for what steps the program should perform in what order (known as the imperative programming model). But there are alternatives. For example, probabilistic programming languages (PPLs) allow programmers to specify target probability distributions to be computed (or sampled from), without specifying how to perform inference with that target distribution. In this exercise, we'll demonstrate how to use the Turing.jl PPL. Turing is implemented entirely in Julia and is highly composable, so that one can perform inference with complex models written in Julia. In this notebook, we'll demonstrate applying it to the simple model for the Hubble expansion.
Usually, we don't know the true value of all our model's parameters (in this case the Hubble constant). Therefore, we want to use astronomical data to measure (or more generally to improve our knowldge of) the model parameters. For any value of the model parameter(s), we can compute the residuals between the observations and the model predicitons, $\Delta_i = d_i - v_{\mathrm{obs},i}/H_0$.
If the model and model parameters were the perfect truth, then $\Delta_i \sim \mathrm{Normal}(0,\sigma_i^2).$
Even if the model isn't perfect, we can compute a goodness-of-fit statistic,
since $E[\Delta_i^2/\sigma_i^2] = 1$ for each $i$.
If we evaluate $\chi^2$ using model parameters that differ from the truth, then the expected value of the value of $\chi^2$ statistic increases. We can find the value of the model parameters than minimize $\chi^2$ to provide us with the best-fit model.
This model for observations is so common that statistican have given a name to the probabilitiy distribution that the $\chi^2$ statistic follows... the $\chi^2$ distribution.
Q1h: Which dataset do you think provides a better estimate of the Hubble constant? Why? For which dataset is the formal measurement uncertainty larger? Explain your why that is.
Use the slider below to pick a threshold velocity below which the model of a linear releationship between $v$ and $d$ appears to be accurate. The figure will automatically update to show only the data you've selected and the predictions with the best-fit $H_0$ for that subset of data.
In this lab, we'll measure the Hubble constant ($H_0$), the rate of expansion of the local universe using observations of Type 1a Supernovae. First, we'll download and read data from the Supernovae Cosmology Project and do a few transformations, so that we can plot velocity versus distance. Second, we'll develop a statistical model to describe our data based on our understanding of the astrophysics and measurement process. We'll compute the maximum likelihood estimate for the model parameters using conventional linear algebra techniques. This is a very useful and common method. (I'd say its analogous the harmonic oscillator in physics!) Next, we'll preview multiple other ways we could have fit the same model, but using different techniques. You'll think about some of the results as we go and the advantages/disadvantages/limitations of each approach.
Q1d: Based on the figure above, what would you estiamte for the best-fit value of $H_0$? How did your previous estimate of the best-fit value for $H_0$ compare to the minimum in the figure above?
Q1i: Is the difference between the best-fit values for the two datasets larger or smaller than the formal uncertainty estimates? What does that imply for the true uncertainty in the Hubble constant?
In the function model_hubble_law_uniform_prior, we specified the probability distribution for each random variable using the syntax variable ~ distribution. For example, we specified the prior for $H_0$ and the likelihood, $\mathcal{L} = p(d_{\mathrm{obs}} | d_{\mathrm{true}}, σ_d)$. We specify how the true distance can be computed from the velocity and Hubble constant, but we . The input arguements for the model are the data to be analyzed, including $d_\mathrm{obs}$. But we won't specify a value for $H_0$ or the values of $d_\mathrm{true}$. The computer will figure out the distributions implied for both of these variables given the data we provide.
Before we can use Julia's CSV reader, first we need to load the CSV.jl package. That's done at the bottom of the notebook with the command using CSV. Then we can call the function CSV.read(...). The first two function arguements specify the filename and that we should read the data into a DataFrame. The remaining arguements are option and allow us to fine-tune the behavior of the CSV reader. This data file used tabs (written using the escaple code \t) to separate the columns rather than the standard commas, so we had to specify that with the optional delim parameter. Additionally, the file included some comments/metadata in the first five rows, so we had to tell the CSV reader to skip to row 6 using the optional skipto parameter. Finally, the file didn't include column heading, so we provided those manually.
Packages like GLM.jl allow us to specify a model with a modeling language. The compiler figures out how to build the matrices and what calculations to peform in order to compute the MLEs. The simplest possible model formula is d ~ v meaning the distance is proportional to velocity (the constant to be fit is implicit). We then call the lm function to return the best fit model for a given model and dataset.
In the figure above, the vertical line indicates the best-fit value for each dataset and the the horizontal bar shows the the measurement uncertainty for each dataset. We could estimate the uncertainties by calculating the second derivative of the log likelihood.
The output table provides the best-fit value for each of the parameters in the linear model. In this case there is the the coefficient to v ($1/H_0$) as well as an intercept. That's because it fit the model $d = c_v v + c_{\mathrm{intercept}}$, rather than just $d = c_v v$. We can tell it that we don't want to include an intercept using the model formula dₗ ~ 0 + v.
Q2c: Given your findings above, which of the following choices is likely to have the most significant effect on your measurement of the Hubble constant?
Choice of supernovae to be included in analysis
Choice of whether to estimate $H_0$ with an MLE or MAP estimator
When using a PPL like Turing, we need to specify a prior distribution for each model parameter and the likelihood function describing how the observations are generated from the model. For example, we could adopt a Normal prior for $H_0$:
Q2d: Can you think of another improvement we could make to improve the quality of inferences about the expansion rate of the universe? Explain your thinking.
One nice feature of using a PPL like Turing is that they make it easy to repeat analyses using slightly different assumptions. For example, we could switch to using a LogNormal prior distribution for $H_0$.
Computes the uncertainty in the maximum likelihood estimate of θ for the linear model y = A θ where measurements errors of y_obs are normally distributed and have covariance covar. Note that y_obs is not passed to this function, because it does not affect results.
Hmm... This value is significantly different that what we found above. Why? In the calls above, we used the method of Ordinary Least Squares which ignores the measurement uncertainties and weights each point equally. We can instead request it use the method of Weighted Least Squares and tell it to use inverse variance weighting as shown below.
Q2e: How could using a probabilistic programming language like Turing make it relatively easy to improve the quality of our inferences about the expansion rate of our universe?
The data file provide measurement uncertainties for the distance modulus. We'd like measurement uncertainties for the distance. Unders several assumptions (e.g., small magnitude of uncertainties, uncorrelated uncertainties), we can estimate the measurement uncertainty in the distance by applying the Propagation of uncertainties method. If you haven't seen this before, then I'd suggest that you just trust me on this for now rather than worrying about the derivation of the equation used.
In this case, it may be tempting to rush to fitting a line to the data. But we want to use this example as a chance to think through the process of building models. It's often useful to divide the full model for astronomical observations into two parts: a physical model and a measurement model.
Physical model
In this case, our initial physical model is that the local universe is expanding away from us at a velocity ($v$) proportion to the distance ($d$). The rate of proportionality is known as the Hubble constant ($H_0$). If this model were perfectly accurate and there were no measurement errors, then we'd have $d = H_0^{-1} v.$
Measurement model
Astronomical measurements almost always have some measurement errors. It can be useful to be explicit about when we're referring to the true value and when we're referring to the observed value. For example, if we were measuring $v$, then we might write: $d_{\rm obs} = d_{\rm true} + \epsilon,$ where $\epsilon$ is the measurement error. By definition the true value of the measurement error ($\epsilon$) is almost always unknown.
Typically, astronomers design the measurements such that the expected value of the measurement is equal the true value ($E[d_{\rm obs}] = d_{\rm true}$), and characterize the distribution of the measurement errors. One very common model/approximation is that the measurement errors follow a Normal (or Gaussian) distribution with a scale parameter, $\sigma$. We denote that with
$$\epsilon \sim \mathrm{Normal}(0,\sigma^2).$$
Based on the properties of the Normal distribution, $E[\epsilon]=0$, and $E[\epsilon^2] = \sigma^2$. For this lab, we'll assume that's a good approximation here.
When we have many measurements, we could write the errors explicitly like $d_{i,\rm obs} = d_{i,\rm true} + \epsilon_i = v_i/H_0 + \epsilon_i,$ where $\epsilon_i \sim \mathrm{Normal}(0,\sigma_i^2)$ for each obsevation $i\in 1...N_{\mathrm obs}$
Evaluating the posterior PDF direclty is hard (often due to the noramlization constant). To get around the challenge, Turing allows us to generate samples from the posterior PDF. There are many interesting algorithms for sampling from a complex posterior distribution. Here, we'll adopt the "NUTS" sampler which based on Hamiltonian Monte Carlo. The algorithm is quite complex, so we won't go into the details. But it's good to konw that this is often a good choice for models with a few to dozens of model parameters. We'll specify that we want it to draw num_mcmc_samples samples and repeat the calculation num_chains different times (running in parallel).
One disadvantage of using simple text files to share data is that they don't have a uniform way to store metadata about the file. For example, it would be nice to include a longer description of the of each column. So after we've read in the datafile, we can add that metadata ourselves.
Thinking a little more broadly, sometimes the goodness-of-fit statistic that we want to minimize is more complicated. So often people refer to minimizing an objective function or a loss function. We'll implement a $\chi^2$ loss function below.
mimetext/htmlrootassigneelast_run_timestampA-mHpersist_js_state·has_pluto_hook_features§cell_id$8ab41174-d7d9-42c7-8129-381c54906280depends_on_disabled_cells§runtime jpublished_object_keysdepends_on_skipped_cells§errored$7f61c098-1c54-4c82-be40-f9fd154156a1queued¤logsrunning¦outputbodyChoice of Priors
The choice of prior probability distributions is often subjective. For example, one astronomer might say that they want to assume as little as practical about the Hubble constant. This has some philosophical appeal, but it can be tricky to put into practice. For example, should we assume that the Hubble constant is positive? Should we assume that there's some maximum velocity (maybe the speed of light)? Should we assume that very large Hubble constants are less likely that small values? If so, what counts as "large"?
Another astronomer might take a different approach and want to approximation to our prior knowledge about the value of the Hubble constant based on previous studies (critically, excluding the data that is currently being analyzing). This also comes with follow-up questions. Should we approximate our prior knowledge with a simple, named distribution? Or should we use a more complex prior to better describe the conclusions of previous research?
There are almost always some very bad choices for priors. But often there's not a single best choice. Rather than trying to find "the correct" prior, it's often wise to perform multiple analyses with multiple reasonable choices for priors. You may find the the difference in outcome are so small that they don't affect your conclusions. Or you might discover than even small differences in the choice of the priors can have a substantive impact on the conclusions. Knowing which regime you're in is often more a valuable result than computing the result under one single choice of priors.
Here, if we assumed a broad uniform prior for $H_0$, then the maximum a posteriori estimate (MAP) would match the MLE. Altneratively, we might adopt as Log-Normal distribution, since it has non-zero probability for all positive values and we're able to pick a location for the peak and scale that are qualitatively similar to what we had already learned from previous studies.
Now we'll call a function (implemented at the bottom of the notebook) to compute the best-fit model parameters for a linear model using linear algebra.
Q1f: What velocity threshold did you pick? What is the results best-fit Hubble constant? How does this compare to the best-fit Hubble constant when fitting to the full data set?
For a linear model, we can estimate the measurement uncertainty associated with our maximum likelihood estimators. It's most intuitive to plot the $\Delta\chi^2$ as a function of $H_0$. By construction, the slope of $\chi^2(H_0)$ is zero at the best-fit value. The magnitude of the measurement uncertainty is inversely proportional to the curvature (absolute value of the second derivative) of $\chi^2(H_0)$.
We can plot histograms for each of the calculations. We can compute an estimtae of the posterior density ($p(H_0|\bf{v},\bf{d},\bf{\sigma_d})$) using the samples from each chain.
and propagate the uncertainties for $H_0$. We'll skip the tedious derivation, but you can see the implementation in the function calc_sigma_mle_linear_model at the bottom of this notebook.
mimetext/htmlrootassigneelast_run_timestampA-xpersist_js_state·has_pluto_hook_features§cell_id$aab05339-dea2-4c7f-8b38-bc7d7971f1d9depends_on_disabled_cells§runtime Ϳpublished_object_keysdepends_on_skipped_cells§errored$46eb1088-0265-46a6-91ec-44dfa3bd862dqueued¤logsrunning¦outputbody
mimeimage/svg+xmlrootassigneelast_run_timestampA:persist_js_state·has_pluto_hook_features§cell_id$46eb1088-0265-46a6-91ec-44dfa3bd862ddepends_on_disabled_cells§runtime0published_object_keysdepends_on_skipped_cells§errored$25ef65bd-8b86-4f88-b0b7-32f05ce82adaqueued¤logsrunning¦outputbodymissingmimetext/plainrootassigneeresponse_1alast_run_timestampA8persist_js_state·has_pluto_hook_features§cell_id$25ef65bd-8b86-4f88-b0b7-32f05ce82adadepends_on_disabled_cells§runtime+Եpublished_object_keysdepends_on_skipped_cells§errored±cell_dependencies $daa349a6-545b-4dfc-a3b0-92884da98e5eprecedence_heuristic cell_id$daa349a6-545b-4dfc-a3b0-92884da98e5edownstream_cells_mapupstream_cells_map@md_strdetailsgetindex$556d8f17-8a03-42c8-82c7-52c8125b604dprecedence_heuristic cell_id$556d8f17-8a03-42c8-82c7-52c8125b604ddownstream_cells_mapcalc_mle_linear_model$2f3d492c-4a5e-4a56-a773-e8078e14476a$4dff84e9-3ca6-4ffb-893b-ef7ec03939b8upstream_cells_map @docDocs.docAbstractVectorAbstractMatrix@assertBase.throwDocs#___this_pluto_module_namesizeadjointPDiagMatlengthBasePDMat\>=*Union==conjBase.AssertionError$77b17593-6a28-4294-9b17-41303d57d5a9precedence_heuristic cell_id$77b17593-6a28-4294-9b17-41303d57d5a9downstream_cells_maploss$0d0917aa-674e-47c4-a8dc-4d1585ebb70e$ce45a056-a949-4a9c-850e-9b5a997ecb0a$67e9aa94-2d3a-4888-97b7-c2cb019eb04a$ccdf5418-12e4-49ee-bdea-90efc202583aupstream_cells_mapAbstractDataFramesumabs2-/Real$2482544e-ad05-495c-9787-0797939445beprecedence_heuristic cell_id$2482544e-ad05-495c-9787-0797939445bedownstream_cells_mapmax_v_plt$0991f15c-b734-480b-8b12-6d2beda547db$b54f8319-dce6-4009-873e-bd47dc2a4619upstream_cells_map @md_strminimumPlutoUI.Slider:maxCorePlutoUI$56d8e6c4-1f19-4bb4-a3db-7bf00675efd9islessBase.get@bindceilBasePlutoRunnerPlutoRunner.create_bond/Core.applicabledf_all$da9eb594-e0fe-4e9c-be69-4150c9a44fb9*maximumgetindex$cd4fc36e-74f1-4b33-b92c-f0b7913a7203precedence_heuristic cell_id$cd4fc36e-74f1-4b33-b92c-f0b7913a7203downstream_cells_mapupstream_cells_map@md_strdetailsgetindex$0735d5d6-9328-47f7-a764-24b4b8d92ca3precedence_heuristic cell_id$0735d5d6-9328-47f7-a764-24b4b8d92ca3downstream_cells_mapupstream_cells_map@md_strgetindex$beeeb261-4dcb-4093-b1cd-08682785eadeprecedence_heuristic cell_id$beeeb261-4dcb-4093-b1cd-08682785eadedownstream_cells_mapupstream_cells_map@md_strasidetipgetindex$da9eb594-e0fe-4e9c-be69-4150c9a44fb9precedence_heuristic cell_id$da9eb594-e0fe-4e9c-be69-4150c9a44fb9downstream_cells_mapdf_all$46eb1088-0265-46a6-91ec-44dfa3bd862d$316f65a1-9062-4cb0-a3fb-af54149afb4d$0d0917aa-674e-47c4-a8dc-4d1585ebb70e$d9058d92-b1c9-44f2-af5c-3540932dec48$fa745b36-7465-443e-9f50-c786fe418c6a$2f3d492c-4a5e-4a56-a773-e8078e14476a$d95d2e2d-c6ff-4cab-8b49-432ef61065bc$466dacda-1835-42fd-b8f4-cac5c22f6756$8ddbc921-4541-4c45-90c5-5c51ec2e391a$b4a68a89-bd1f-49d5-adb2-6239317a8b82$0991f15c-b734-480b-8b12-6d2beda547db$2482544e-ad05-495c-9787-0797939445be$ce45a056-a949-4a9c-850e-9b5a997ecb0a$67e9aa94-2d3a-4888-97b7-c2cb019eb04aupstream_cells_mapselect@__dot__colmetadata!-logspeed_of_light_mps$a7d8f517-4117-4585-8089-117510c944c5/^+df_in$4518f1f7-a2a3-4783-86ae-4413fc4291b5added_metadata$d81347d3-9246-42d9-9a85-42c773ef2a1d*$d19e7f11-d1ff-4624-8603-38ffec0c7beaprecedence_heuristic cell_id$d19e7f11-d1ff-4624-8603-38ffec0c7beadownstream_cells_mapupstream_cells_map@md_strgetindex$c6f03d5d-18e2-4491-90fe-fe0008f774c4precedence_heuristic cell_id$c6f03d5d-18e2-4491-90fe-fe0008f774c4downstream_cells_mapupstream_cells_map@md_strgetindex$629dae3a-46ec-40dc-b9ed-92de9eb21adbprecedence_heuristic cell_id$629dae3a-46ec-40dc-b9ed-92de9eb21adbdownstream_cells_mapupstream_cells_mapresponse_1e$bd1c72d5-1d0d-4e36-aba4-95c9c9e4fab1ismissingstill_missing$4518f1f7-a2a3-4783-86ae-4413fc4291b5precedence_heuristic cell_id$4518f1f7-a2a3-4783-86ae-4413fc4291b5downstream_cells_mapdf_in$d81347d3-9246-42d9-9a85-42c773ef2a1d$da9eb594-e0fe-4e9c-be69-4150c9a44fb9upstream_cells_mapCSV$56d8e6c4-1f19-4bb4-a3db-7bf00675efd9filename$3f2ef930-d367-11ef-3b5c-95b4989c2e21DataFrameCSV.read$e4597b96-6e06-44c2-a845-3d9967d2423cprecedence_heuristic cell_id$e4597b96-6e06-44c2-a845-3d9967d2423cdownstream_cells_mapupstream_cells_mapismissingresponse_1i$751eff6b-2ed2-4073-b8b5-02ad9981edb6still_missing$92b2c9f6-5102-4e1a-aff1-05e9ce8c19a8precedence_heuristic cell_id$92b2c9f6-5102-4e1a-aff1-05e9ce8c19a8downstream_cells_mapupstream_cells_map@md_strgetindex$6c86d0f1-5c6f-4204-9414-60577164f4b8precedence_heuristic cell_id$6c86d0f1-5c6f-4204-9414-60577164f4b8downstream_cells_map+model_hubble_law_lognormal_prior_given_data$cf6fbe84-9a51-4ddb-9def-c22893a17ed3upstream_cells_mapdf_fit$0991f15c-b734-480b-8b12-6d2beda547db model_hubble_law_lognormal_prior$d3db2bc4-a26b-4e01-a193-b155a92dab48$a7981271-0b12-4b75-b5f0-5e91af761419precedence_heuristic cell_id$a7981271-0b12-4b75-b5f0-5e91af761419downstream_cells_mapupstream_cells_map@md_strgetindex$276a850b-09da-41f0-8661-5d63bb2afdf8precedence_heuristic cell_id$276a850b-09da-41f0-8661-5d63bb2afdf8downstream_cells_mapchain_uniform_prior$1eb1face-a795-4cb3-9c51-f995f73e2509$0a4aaa8a-20da-486f-ace5-fbadcf646609$2ad2563a-8a72-47e8-a36d-f6d2a49cb517$915fc532-7d8d-4ff5-81fd-af8f6e9e64daupstream_cells_map)model_hubble_law_uniform_prior_given_data$340dc0ff-dd1c-44e0-907c-f357adb122d5num_mcmc_samples$2c1c7d8a-eda6-4900-9be7-e4b908fed81aNUTSnum_chains$2c1c7d8a-eda6-4900-9be7-e4b908fed81aMCMCThreadssampleready_to_run_mcmc$c433e938-31f7-4417-b3c0-a2f49fea8640nothing$bd560d5d-3bc7-4b4d-aab7-67b1884a974bprecedence_heuristic cell_id$bd560d5d-3bc7-4b4d-aab7-67b1884a974bdownstream_cells_mapupstream_cells_mapscatterloss_vs_H₀_all$0d0917aa-674e-47c4-a8dc-4d1585ebb70eH₀s_plt$0d0917aa-674e-47c4-a8dc-4d1585ebb70eplot!xlabel!ylabel!$0236d055-0d6d-440d-a44b-0d5cb0ab2273precedence_heuristic cell_id$0236d055-0d6d-440d-a44b-0d5cb0ab2273downstream_cells_mapresponse_2b$2d97268e-560b-4398-982e-1cf5cd9ab7ddupstream_cells_mapmissing$efde0b1e-ec47-4180-aba7-ed8843395b50precedence_heuristic cell_id$efde0b1e-ec47-4180-aba7-ed8843395b50downstream_cells_mapupstream_cells_map@md_strgetindex$48e246bc-4937-40d6-b686-73524f3ce057precedence_heuristic cell_id$48e246bc-4937-40d6-b686-73524f3ce057downstream_cells_mapupstream_cells_map@md_strgetindex$2786b37a-207f-4831-983e-81a380c37045precedence_heuristic cell_id$2786b37a-207f-4831-983e-81a380c37045downstream_cells_mapupstream_cells_map@md_strgetindex$0d0917aa-674e-47c4-a8dc-4d1585ebb70eprecedence_heuristic cell_id$0d0917aa-674e-47c4-a8dc-4d1585ebb70edownstream_cells_maploss_vs_H₀_all$bd560d5d-3bc7-4b4d-aab7-67b1884a974bH₀s_plt$bd560d5d-3bc7-4b4d-aab7-67b1884a974bupstream_cells_mapmaprangedf_all$da9eb594-e0fe-4e9c-be69-4150c9a44fb9loss$77b17593-6a28-4294-9b17-41303d57d5a9$d3e82a29-c1b5-4f5a-ba96-a64d22636493precedence_heuristic cell_id$d3e82a29-c1b5-4f5a-ba96-a64d22636493downstream_cells_mapupstream_cells_map@md_strgetindex$49090825-2843-4b3e-8ee6-658187232f96precedence_heuristic cell_id$49090825-2843-4b3e-8ee6-658187232f96downstream_cells_mapupstream_cells_map@md_strgetindex$e0e610d7-a237-42ea-8706-8e09858fba03precedence_heuristic cell_id$e0e610d7-a237-42ea-8706-8e09858fba03downstream_cells_mapupstream_cells_mapismissingresponse_2e$adf117fe-0f54-4f2a-b107-0635f56263f2still_missing$ce45a056-a949-4a9c-850e-9b5a997ecb0aprecedence_heuristic cell_id$ce45a056-a949-4a9c-850e-9b5a997ecb0adownstream_cells_mapupstream_cells_map H₀_bf_all$f76f4c3c-74f1-44b0-9b3a-5895aec22e70!σ_H₀_fit$4dff84e9-3ca6-4ffb-893b-ef7ec03939b8df_fit$0991f15c-b734-480b-8b12-6d2beda547dbmaprangeplot!fillscatter!response_1f$69ca76b8-8766-47fe-89d6-31fcc7a60e88H₀_wls_fit$0991f15c-b734-480b-8b12-6d2beda547dbsizeσ_H₀_all$1d70572f-214e-4161-b804-2c3b7c077c56ismissingloss$77b17593-6a28-4294-9b17-41303d57d5a9-ylims!plotdf_all$da9eb594-e0fe-4e9c-be69-4150c9a44fb9+*H₀_wls_all$d67f0b9e-10be-4511-93ab-edceb3aef49c$1d70572f-214e-4161-b804-2c3b7c077c56precedence_heuristic cell_id$1d70572f-214e-4161-b804-2c3b7c077c56downstream_cells_mapσ_H₀inv_allσ_H₀_all$ce45a056-a949-4a9c-850e-9b5a997ecb0aupstream_cells_mapH₀_bf_all$f76f4c3c-74f1-44b0-9b3a-5895aec22e70/Σ_all$fa745b36-7465-443e-9f50-c786fe418c6acalc_sigma_mle_linear_model$c086c05f-ad73-44c7-bf8b-a304d3533690A_all$d9058d92-b1c9-44f2-af5c-3540932dec48*H₀inv_bf_all$2f3d492c-4a5e-4a56-a773-e8078e14476aabs$23ee2544-2a2f-4302-8a62-2dd868cbb5d1precedence_heuristic cell_id$23ee2544-2a2f-4302-8a62-2dd868cbb5d1downstream_cells_mapresponse_2d$a01c7aad-429d-447b-8943-a4e89dbac6f6upstream_cells_mapmissing$69ca76b8-8766-47fe-89d6-31fcc7a60e88precedence_heuristic cell_id$69ca76b8-8766-47fe-89d6-31fcc7a60e88downstream_cells_mapresponse_1f$0adf35a1-d7a8-4e7b-ad3f-ac83bd44ef18$ce45a056-a949-4a9c-850e-9b5a997ecb0a$67e9aa94-2d3a-4888-97b7-c2cb019eb04a$ccdf5418-12e4-49ee-bdea-90efc202583a$4dff84e9-3ca6-4ffb-893b-ef7ec03939b8upstream_cells_map$18106638-c96c-47d6-990d-b35b6b0ba218precedence_heuristic cell_id$18106638-c96c-47d6-990d-b35b6b0ba218downstream_cells_maphistogram_or_density$1eb1face-a795-4cb3-9c51-f995f73e2509$915fc532-7d8d-4ff5-81fd-af8f6e9e64daupstream_cells_map@md_strCorehistogramBase.get@bind=>BasePlutoRunnerPlutoRunner.create_bonddensityCore.applicableSelectgetindex$f90594d8-92a1-4b0b-93d0-0e41e76b3f55precedence_heuristic cell_id$f90594d8-92a1-4b0b-93d0-0e41e76b3f55downstream_cells_mapupstream_cells_map@md_strgetindex$aaa80da1-f7de-42ed-9691-c367e6bd9a28precedence_heuristic cell_id$aaa80da1-f7de-42ed-9691-c367e6bd9a28downstream_cells_mapupstream_cells_mapTableOfContents$0adf35a1-d7a8-4e7b-ad3f-ac83bd44ef18precedence_heuristic cell_id$0adf35a1-d7a8-4e7b-ad3f-ac83bd44ef18downstream_cells_mapupstream_cells_mapresponse_1f$69ca76b8-8766-47fe-89d6-31fcc7a60e88ismissingstill_missing$92ca099c-b81b-4cfe-a3fa-60833657f4efprecedence_heuristic cell_id$92ca099c-b81b-4cfe-a3fa-60833657f4efdownstream_cells_mapupstream_cells_map@md_strgetindex$03fc47ea-9246-4cb8-8398-612fc2fb09faprecedence_heuristic cell_id$03fc47ea-9246-4cb8-8398-612fc2fb09fadownstream_cells_mapresponse_1b$0c0312da-558a-4bc1-9e4a-b02ecd1c0572upstream_cells_mapmissing$d3db2bc4-a26b-4e01-a193-b155a92dab48precedence_heuristic cell_id$d3db2bc4-a26b-4e01-a193-b155a92dab48downstream_cells_map model_hubble_law_lognormal_prior$6c86d0f1-5c6f-4204-9414-60577164f4b8upstream_cells_map LogNormalmissingBase.throwMvNormal@model@assert!H₀sizeAnyNamedTupleBase===log/==Base.AssertionError$9db3ff24-6f54-4a37-b787-de4b14694d95precedence_heuristic cell_id$9db3ff24-6f54-4a37-b787-de4b14694d95downstream_cells_mapcalc_fisher_matrix_linear_model$c086c05f-ad73-44c7-bf8b-a304d3533690upstream_cells_map @docBase.throwAbstractMatrix@assertTuple#___this_pluto_module_namesizeadjointBase\>=*Union==conjBase.AssertionError$c6af6d4a-ff58-4012-b4b2-4e3b09d53638precedence_heuristic cell_id$c6af6d4a-ff58-4012-b4b2-4e3b09d53638downstream_cells_mapupstream_cells_mapresponse_1h$22da4ae4-6aea-4d45-ab2c-6fe935783e41ismissingstill_missing$0c0312da-558a-4bc1-9e4a-b02ecd1c0572precedence_heuristic cell_id$0c0312da-558a-4bc1-9e4a-b02ecd1c0572downstream_cells_mapupstream_cells_mapresponse_1b$03fc47ea-9246-4cb8-8398-612fc2fb09faismissingstill_missing$b15b1913-d265-4a35-b43e-67df53c93329precedence_heuristic cell_id$b15b1913-d265-4a35-b43e-67df53c93329downstream_cells_mapupstream_cells_map@md_strdetailsgetindex$4655178c-a633-4b50-8f92-6ee028e21d0eprecedence_heuristic cell_id$4655178c-a633-4b50-8f92-6ee028e21d0edownstream_cells_mapupstream_cells_map@md_strgetindex$3bdb27da-451c-4260-bf43-5568985d24f2precedence_heuristic cell_id$3bdb27da-451c-4260-bf43-5568985d24f2downstream_cells_mapupstream_cells_map@md_strgetindex$2f7c11e8-d429-41eb-873e-ae554e028195precedence_heuristic cell_id$2f7c11e8-d429-41eb-873e-ae554e028195downstream_cells_mapupstream_cells_map@md_strgetindex$8764ccfe-81cd-4bff-985a-c112fc5c902fprecedence_heuristic cell_id$8764ccfe-81cd-4bff-985a-c112fc5c902fdownstream_cells_mapupstream_cells_map@md_strgetindex$b7e2903b-be56-4e36-885f-e16a9ddc94c6precedence_heuristic cell_id$b7e2903b-be56-4e36-885f-e16a9ddc94c6downstream_cells_mapupstream_cells_map@md_strgetindex$4eb123a1-f976-4fad-9a96-e3a61bdbf264precedence_heuristic cell_id$4eb123a1-f976-4fad-9a96-e3a61bdbf264downstream_cells_mapresponse_2a$33045766-a7c5-42ff-bdf8-bde1fce9c56bupstream_cells_mapmissing$089c50b6-5494-4f27-af3e-a8af972e2869precedence_heuristic cell_id$089c50b6-5494-4f27-af3e-a8af972e2869downstream_cells_mapresponse_2c$fb2545e1-d624-4669-9e07-c94a598ba8edupstream_cells_mapmissing$34ccd3d7-e914-4f6f-a8eb-f421ed17e709precedence_heuristic cell_id$34ccd3d7-e914-4f6f-a8eb-f421ed17e709downstream_cells_mapupstream_cells_map@md_strgetindex$e02dd163-59bc-4d48-ae43-6fc559536aa6precedence_heuristic cell_id$e02dd163-59bc-4d48-ae43-6fc559536aa6downstream_cells_mapupstream_cells_map@md_strgetindex$a56871c9-affa-4612-b3f5-66363720f274precedence_heuristic cell_id$a56871c9-affa-4612-b3f5-66363720f274downstream_cells_mapupstream_cells_map@md_strgetindex$4f27d8d1-0d14-49e5-abee-9f277369a78bprecedence_heuristic cell_id$4f27d8d1-0d14-49e5-abee-9f277369a78bdownstream_cells_mapmodel_hubble_law_uniform_prior$340dc0ff-dd1c-44e0-907c-f357adb122d5upstream_cells_map missingBase.throwMvNormal@model@assert!UniformH₀sizeAnyNamedTupleBase-speed_of_light_mps$a7d8f517-4117-4585-8089-117510c944c5===/==Base.AssertionError$0144843f-b649-4ea2-a546-e8c14a928670precedence_heuristic cell_id$0144843f-b649-4ea2-a546-e8c14a928670downstream_cells_mapupstream_cells_map@md_strgetindex$a7d8f517-4117-4585-8089-117510c944c5precedence_heuristic cell_id$a7d8f517-4117-4585-8089-117510c944c5downstream_cells_mapspeed_of_light_mps$da9eb594-e0fe-4e9c-be69-4150c9a44fb9$4f27d8d1-0d14-49e5-abee-9f277369a78bupstream_cells_map$18dd7cac-fac3-4c40-9230-f6387e41751bprecedence_heuristic cell_id$18dd7cac-fac3-4c40-9230-f6387e41751bdownstream_cells_mapupstream_cells_map@md_strgetindex$fa745b36-7465-443e-9f50-c786fe418c6aprecedence_heuristic cell_id$fa745b36-7465-443e-9f50-c786fe418c6adownstream_cells_mapΣ_all$2f3d492c-4a5e-4a56-a773-e8078e14476a$1d70572f-214e-4161-b804-2c3b7c077c56upstream_cells_map^df_all$da9eb594-e0fe-4e9c-be69-4150c9a44fb9PDiagMat$b4a68a89-bd1f-49d5-adb2-6239317a8b82precedence_heuristic cell_id$b4a68a89-bd1f-49d5-adb2-6239317a8b82downstream_cells_mapwls_all$d67f0b9e-10be-4511-93ab-edceb3aef49cupstream_cells_mapweightsformula_wo_intercept$2656d886-60dc-44a7-b318-100335dcf05b^df_all$da9eb594-e0fe-4e9c-be69-4150c9a44fb9lm$a8d72fee-302c-44f6-8ca3-9e0b1a130ec6precedence_heuristic cell_id$a8d72fee-302c-44f6-8ca3-9e0b1a130ec6downstream_cells_mapupstream_cells_map@md_strdetailsgetindex$c433e938-31f7-4417-b3c0-a2f49fea8640precedence_heuristic cell_id$c433e938-31f7-4417-b3c0-a2f49fea8640downstream_cells_mapready_to_run_mcmc$276a850b-09da-41f0-8661-5d63bb2afdf8$cf6fbe84-9a51-4ddb-9def-c22893a17ed3upstream_cells_mapCore@md_strBasePlutoRunner.create_bondPlutoRunnerCheckBoxCore.applicable@bindBase.getgetindex$3aba1134-7c60-41f6-a1ac-73306efb7824precedence_heuristic cell_id$3aba1134-7c60-41f6-a1ac-73306efb7824downstream_cells_mapupstream_cells_map@md_strgetindex$5e635a20-e205-48f0-908d-cb6d0c55b43bprecedence_heuristic cell_id$5e635a20-e205-48f0-908d-cb6d0c55b43bdownstream_cells_mapupstream_cells_map@md_strgetindex$4dff84e9-3ca6-4ffb-893b-ef7ec03939b8precedence_heuristic cell_id$4dff84e9-3ca6-4ffb-893b-ef7ec03939b8downstream_cells_mapH₀inv_fitΣ_fitσ_H₀_fit$ce45a056-a949-4a9c-850e-9b5a997ecb0aA_fitH₀_fit$ccdf5418-12e4-49ee-bdea-90efc202583aσ_H₀inv_fitupstream_cells_map!response_1f$69ca76b8-8766-47fe-89d6-31fcc7a60e88ismissingdf_fit$0991f15c-b734-480b-8b12-6d2beda547dbsizePDiagMatcalc_mle_linear_model$556d8f17-8a03-42c8-82c7-52c8125b604dreshape/^calc_sigma_mle_linear_model$c086c05f-ad73-44c7-bf8b-a304d3533690*abs$7be54420-b2ec-4ef4-adc8-848817305d07precedence_heuristic cell_id$7be54420-b2ec-4ef4-adc8-848817305d07downstream_cells_mapupstream_cells_map@md_strgetindex$16b0a2c1-5803-4310-9ee7-884b1333fb74precedence_heuristic cell_id$16b0a2c1-5803-4310-9ee7-884b1333fb74downstream_cells_mapupstream_cells_map@md_strgetindex$0c20b9a7-69a9-49c8-9ba4-45b8770393d0precedence_heuristic cell_id$0c20b9a7-69a9-49c8-9ba4-45b8770393d0downstream_cells_mapupstream_cells_map@md_strgetindex$56d8e6c4-1f19-4bb4-a3db-7bf00675efd9precedence_heuristiccell_id$56d8e6c4-1f19-4bb4-a3db-7bf00675efd9downstream_cells_mapCSV$4518f1f7-a2a3-4783-86ae-4413fc4291b5Downloads$3f2ef930-d367-11ef-3b5c-95b4989c2e21StatsBaseMCMCChainsPDMatsDistributionsPlotsPlutoUI$d70d95e9-f0e7-4225-8e06-80735278cbca$2482544e-ad05-495c-9787-0797939445beStatsPlotsPlutoTeachingToolsLinearAlgebraDataFramesTuringLaTeXStringsGLMupstream_cells_map$d70d95e9-f0e7-4225-8e06-80735278cbcaprecedence_heuristic cell_id$d70d95e9-f0e7-4225-8e06-80735278cbcadownstream_cells_mapplt_resid$316f65a1-9062-4cb0-a3fb-af54149afb4dplt_yerr$316f65a1-9062-4cb0-a3fb-af54149afb4dH₀_plt_1$316f65a1-9062-4cb0-a3fb-af54149afb4dupstream_cells_map@md_strPlutoUI.Slider:CorePlutoUI$56d8e6c4-1f19-4bb4-a3db-7bf00675efd9Base.get@bindPlutoUI.CheckBoxBasePlutoRunnerPlutoRunner.create_bondCore.applicablegetindex$9f4d8983-d2d7-4bc9-bb31-5460a18dac5bprecedence_heuristic cell_id$9f4d8983-d2d7-4bc9-bb31-5460a18dac5bdownstream_cells_mapupstream_cells_map@md_strdetailsgetindex$22da4ae4-6aea-4d45-ab2c-6fe935783e41precedence_heuristic cell_id$22da4ae4-6aea-4d45-ab2c-6fe935783e41downstream_cells_mapresponse_1h$c6af6d4a-ff58-4012-b4b2-4e3b09d53638upstream_cells_mapmissing$a8f89fe8-110b-452e-9ecb-53769aa02a2eprecedence_heuristic cell_id$a8f89fe8-110b-452e-9ecb-53769aa02a2edownstream_cells_mapupstream_cells_map@md_strgetindex$a01c7aad-429d-447b-8943-a4e89dbac6f6precedence_heuristic cell_id$a01c7aad-429d-447b-8943-a4e89dbac6f6downstream_cells_mapupstream_cells_mapresponse_2d$23ee2544-2a2f-4302-8a62-2dd868cbb5d1ismissingstill_missing$cf6fbe84-9a51-4ddb-9def-c22893a17ed3precedence_heuristic cell_id$cf6fbe84-9a51-4ddb-9def-c22893a17ed3downstream_cells_mapchain_lognormal_prior$18edb530-8472-4932-946c-917d8ab9dbde$915fc532-7d8d-4ff5-81fd-af8f6e9e64daupstream_cells_mapnum_mcmc_samples$2c1c7d8a-eda6-4900-9be7-e4b908fed81aNUTSnum_chains$2c1c7d8a-eda6-4900-9be7-e4b908fed81aMCMCThreadssampleready_to_run_mcmc$c433e938-31f7-4417-b3c0-a2f49fea8640+model_hubble_law_lognormal_prior_given_data$6c86d0f1-5c6f-4204-9414-60577164f4b8nothing$098afd51-2c39-42ef-91bf-cfc326827ea1precedence_heuristic cell_id$098afd51-2c39-42ef-91bf-cfc326827ea1downstream_cells_mapupstream_cells_map@md_strgetindex$d95d2e2d-c6ff-4cab-8b49-432ef61065bcprecedence_heuristic cell_id$d95d2e2d-c6ff-4cab-8b49-432ef61065bcdownstream_cells_mapupstream_cells_mapplot_hubble_diagram$2b3729fc-62b2-478d-b35c-e69e6b0efc93H₀_bf_all$f76f4c3c-74f1-44b0-9b3a-5895aec22e70title!df_all$da9eb594-e0fe-4e9c-be69-4150c9a44fb9$a1b5939b-b1ed-4d7c-907d-5bf3e7dd4479precedence_heuristic cell_id$a1b5939b-b1ed-4d7c-907d-5bf3e7dd4479downstream_cells_mapupstream_cells_map@md_strgetindex$78eb6aff-3752-41d8-a4d9-0cc2b31fd56dprecedence_heuristic cell_id$78eb6aff-3752-41d8-a4d9-0cc2b31fd56ddownstream_cells_mapupstream_cells_map@md_strgetindex$ac47f631-787d-4a3c-96cf-5fee86a3b2f4precedence_heuristic cell_id$ac47f631-787d-4a3c-96cf-5fee86a3b2f4downstream_cells_mapH₀_ols_allupstream_cells_map/coefols_all$8ddbc921-4541-4c45-90c5-5c51ec2e391a$42366f09-d437-4ba2-93ec-f7b8278eb09dprecedence_heuristic cell_id$42366f09-d437-4ba2-93ec-f7b8278eb09ddownstream_cells_mapupstream_cells_map@md_strgetindex$d0f003fa-f72d-4a6f-ba60-1cad3844ef6aprecedence_heuristic cell_id$d0f003fa-f72d-4a6f-ba60-1cad3844ef6adownstream_cells_mapupstream_cells_map@md_strgetindex$f0220da0-6f3a-438c-b533-acb346382418precedence_heuristic cell_id$f0220da0-6f3a-438c-b533-acb346382418downstream_cells_mapupstream_cells_map@md_strdetailsgetindex$8bfe0413-782b-4244-8624-a9f95e393d17precedence_heuristic cell_id$8bfe0413-782b-4244-8624-a9f95e393d17downstream_cells_mapupstream_cells_map@md_strgetindex$3ec3e49c-04ce-4ffa-b1ca-8c58df0307ddprecedence_heuristic cell_id$3ec3e49c-04ce-4ffa-b1ca-8c58df0307dddownstream_cells_mapupstream_cells_map@md_strgetindex$92238689-d9df-46c2-bd72-3fcc41a22594precedence_heuristic cell_id$92238689-d9df-46c2-bd72-3fcc41a22594downstream_cells_mapupstream_cells_map@md_strgetindex$d5c3923b-8515-4002-9d93-6768f087aef6precedence_heuristic cell_id$d5c3923b-8515-4002-9d93-6768f087aef6downstream_cells_mapupstream_cells_map@md_strgetindex$4f11d0e2-7126-47f1-865a-3474298597dcprecedence_heuristic cell_id$4f11d0e2-7126-47f1-865a-3474298597dcdownstream_cells_mapupstream_cells_map@md_strgetindex$337a560b-49f3-4fc5-926f-1689911eacf2precedence_heuristic cell_id$337a560b-49f3-4fc5-926f-1689911eacf2downstream_cells_mapupstream_cells_map@md_strgetindex$7b5b3e8a-61ac-43cf-8b85-4a8a8e4846f6precedence_heuristic cell_id$7b5b3e8a-61ac-43cf-8b85-4a8a8e4846f6downstream_cells_mapupstream_cells_map@md_strgetindex$d81347d3-9246-42d9-9a85-42c773ef2a1dprecedence_heuristic cell_id$d81347d3-9246-42d9-9a85-42c773ef2a1ddownstream_cells_mapadded_metadata$da9eb594-e0fe-4e9c-be69-4150c9a44fb9upstream_cells_mapmetadata!df_in$4518f1f7-a2a3-4783-86ae-4413fc4291b5url$3f2ef930-d367-11ef-3b5c-95b4989c2e21colmetadata!$adf117fe-0f54-4f2a-b107-0635f56263f2precedence_heuristic cell_id$adf117fe-0f54-4f2a-b107-0635f56263f2downstream_cells_mapresponse_2e$e0e610d7-a237-42ea-8706-8e09858fba03upstream_cells_mapmissing$2ad2563a-8a72-47e8-a36d-f6d2a49cb517precedence_heuristic cell_id$2ad2563a-8a72-47e8-a36d-f6d2a49cb517downstream_cells_mapH₀_posterior_meanH₀_posterior_stdupstream_cells_mapvecchain_uniform_prior$276a850b-09da-41f0-8661-5d63bb2afdf8mean_and_std==!=nothing$c086c05f-ad73-44c7-bf8b-a304d3533690precedence_heuristic cell_id$c086c05f-ad73-44c7-bf8b-a304d3533690downstream_cells_mapcalc_sigma_mle_linear_model$1d70572f-214e-4161-b804-2c3b7c077c56$4dff84e9-3ca6-4ffb-893b-ef7ec03939b8upstream_cells_map @docinvdiagBase.throwAbstractMatrix@assertTuplesqrt#___this_pluto_module_namesizecalc_fisher_matrix_linear_model$9db3ff24-6f54-4a37-b787-de4b14694d95BasePDMat>=Union==Base.AssertionError$0a4aaa8a-20da-486f-ace5-fbadcf646609precedence_heuristic cell_id$0a4aaa8a-20da-486f-ace5-fbadcf646609downstream_cells_mapupstream_cells_mapdescribechain_uniform_prior$276a850b-09da-41f0-8661-5d63bb2afdf8!===nothing$eb72dda6-69f0-4538-a75a-25967653759fprecedence_heuristic cell_id$eb72dda6-69f0-4538-a75a-25967653759fdownstream_cells_mapupstream_cells_map@md_strgetindex$2d97268e-560b-4398-982e-1cf5cd9ab7ddprecedence_heuristic cell_id$2d97268e-560b-4398-982e-1cf5cd9ab7dddownstream_cells_mapupstream_cells_mapismissingresponse_2b$0236d055-0d6d-440d-a44b-0d5cb0ab2273still_missing$28b844a8-650a-47d8-b9f6-586f22d6f603precedence_heuristic cell_id$28b844a8-650a-47d8-b9f6-586f22d6f603downstream_cells_mapupstream_cells_map@md_strgetindex$c78f0ad9-03c0-4066-a390-193e7bfa84b1precedence_heuristic cell_id$c78f0ad9-03c0-4066-a390-193e7bfa84b1downstream_cells_mapupstream_cells_map@md_strgetindex$1c35257a-9b98-4af6-bde0-a945665c91cdprecedence_heuristic cell_id$1c35257a-9b98-4af6-bde0-a945665c91cddownstream_cells_mapupstream_cells_map@md_strdetailsgetindex$751eff6b-2ed2-4073-b8b5-02ad9981edb6precedence_heuristic cell_id$751eff6b-2ed2-4073-b8b5-02ad9981edb6downstream_cells_mapresponse_1i$e4597b96-6e06-44c2-a845-3d9967d2423cupstream_cells_mapmissing$5974dfc3-8416-4bae-91de-d3f0a4f0c7dbprecedence_heuristic cell_id$5974dfc3-8416-4bae-91de-d3f0a4f0c7dbdownstream_cells_mapupstream_cells_mapresponse_1a$25ef65bd-8b86-4f88-b0b7-32f05ce82adaismissingstill_missing$20f739ff-cc5b-460e-851a-141003d0c9d0precedence_heuristic cell_id$20f739ff-cc5b-460e-851a-141003d0c9d0downstream_cells_mapupstream_cells_map@md_strgetindex$f76f4c3c-74f1-44b0-9b3a-5895aec22e70precedence_heuristic cell_id$f76f4c3c-74f1-44b0-9b3a-5895aec22e70downstream_cells_mapH₀_bf_all$d95d2e2d-c6ff-4cab-8b49-432ef61065bc$ce45a056-a949-4a9c-850e-9b5a997ecb0a$67e9aa94-2d3a-4888-97b7-c2cb019eb04a$1d70572f-214e-4161-b804-2c3b7c077c56upstream_cells_map/H₀inv_bf_all$2f3d492c-4a5e-4a56-a773-e8078e14476a$a98ace81-9afd-4618-88a6-fadaf04eb364precedence_heuristic cell_id$a98ace81-9afd-4618-88a6-fadaf04eb364downstream_cells_mapupstream_cells_mapresponse_1d$0c25a4cc-9f53-4729-985b-de692278c433ismissingstill_missing$bd1c72d5-1d0d-4e36-aba4-95c9c9e4fab1precedence_heuristic cell_id$bd1c72d5-1d0d-4e36-aba4-95c9c9e4fab1downstream_cells_mapresponse_1e$629dae3a-46ec-40dc-b9ed-92de9eb21adbupstream_cells_mapmissing$95bc4b33-d619-4a95-af19-5771d209336dprecedence_heuristic cell_id$95bc4b33-d619-4a95-af19-5771d209336ddownstream_cells_maptopic$e09e30f7-838f-4571-bb9c-6d0b425807a8$bf7ca6ea-4c07-46e4-baff-b205e1663505title$e09e30f7-838f-4571-bb9c-6d0b425807a8$bf7ca6ea-4c07-46e4-baff-b205e1663505week$e09e30f7-838f-4571-bb9c-6d0b425807a8upstream_cells_map$4bdd7b2d-9a67-4161-a0a0-93073d9edcf0precedence_heuristic cell_id$4bdd7b2d-9a67-4161-a0a0-93073d9edcf0downstream_cells_mapupstream_cells_map@md_strgetindex$2532e9de-8fa7-44ce-8502-c2473928aa8eprecedence_heuristic cell_id$2532e9de-8fa7-44ce-8502-c2473928aa8edownstream_cells_mapupstream_cells_mapWidthOverDocs$4a2c109d-8e18-4561-913b-769f7eaeb6e9precedence_heuristic cell_id$4a2c109d-8e18-4561-913b-769f7eaeb6e9downstream_cells_mapupstream_cells_map@md_strgetindex$2f3d492c-4a5e-4a56-a773-e8078e14476aprecedence_heuristic cell_id$2f3d492c-4a5e-4a56-a773-e8078e14476adownstream_cells_mapH₀inv_bf_all$f76f4c3c-74f1-44b0-9b3a-5895aec22e70$1d70572f-214e-4161-b804-2c3b7c077c56upstream_cells_mapcalc_mle_linear_model$556d8f17-8a03-42c8-82c7-52c8125b604dA_all$d9058d92-b1c9-44f2-af5c-3540932dec48df_all$da9eb594-e0fe-4e9c-be69-4150c9a44fb9Σ_all$fa745b36-7465-443e-9f50-c786fe418c6a$c43f2686-da2d-4b67-b87e-7c22bd09257cprecedence_heuristic cell_id$c43f2686-da2d-4b67-b87e-7c22bd09257cdownstream_cells_mapresponse_1g$3a06141a-a92d-4c31-9950-13f0bf2c023aupstream_cells_mapmissing$b54f8319-dce6-4009-873e-bd47dc2a4619precedence_heuristic cell_id$b54f8319-dce6-4009-873e-bd47dc2a4619downstream_cells_mapupstream_cells_mapstringplot_hubble_diagram$2b3729fc-62b2-478d-b35c-e69e6b0efc93floortitle!Int*H₀_wls_fit$0991f15c-b734-480b-8b12-6d2beda547dbdf_fit$0991f15c-b734-480b-8b12-6d2beda547dbmax_v_plt$2482544e-ad05-495c-9787-0797939445be$ddf74f60-42eb-4eab-b77b-b6e78f0e90a7precedence_heuristic cell_id$ddf74f60-42eb-4eab-b77b-b6e78f0e90a7downstream_cells_mapupstream_cells_map@md_strgetindex$a6cf5e52-8692-4e3c-95aa-11167b118bd4precedence_heuristic cell_id$a6cf5e52-8692-4e3c-95aa-11167b118bd4downstream_cells_mapupstream_cells_map@md_strgetindex$0f6aa1eb-2f89-400f-b1e8-bfc07012179dprecedence_heuristic cell_id$0f6aa1eb-2f89-400f-b1e8-bfc07012179ddownstream_cells_mapupstream_cells_map@md_strgetindex$0991f15c-b734-480b-8b12-6d2beda547dbprecedence_heuristic cell_id$0991f15c-b734-480b-8b12-6d2beda547dbdownstream_cells_mapwls_fitH₀_wls_fit$b54f8319-dce6-4009-873e-bd47dc2a4619$ce45a056-a949-4a9c-850e-9b5a997ecb0adf_fit$b54f8319-dce6-4009-873e-bd47dc2a4619$ce45a056-a949-4a9c-850e-9b5a997ecb0a$ccdf5418-12e4-49ee-bdea-90efc202583a$4dff84e9-3ca6-4ffb-893b-ef7ec03939b8$340dc0ff-dd1c-44e0-907c-f357adb122d5$6c86d0f1-5c6f-4204-9414-60577164f4b8upstream_cells_map<=weightsformula_wo_intercept$2656d886-60dc-44a7-b318-100335dcf05b^df_all$da9eb594-e0fe-4e9c-be69-4150c9a44fb9/lmcoeffiltermax_v_plt$2482544e-ad05-495c-9787-0797939445be$5c9b0395-4d2e-4679-8d86-b9581fa0f66cprecedence_heuristic cell_id$5c9b0395-4d2e-4679-8d86-b9581fa0f66cdownstream_cells_mapupstream_cells_map@md_strgetindex$340dc0ff-dd1c-44e0-907c-f357adb122d5precedence_heuristic cell_id$340dc0ff-dd1c-44e0-907c-f357adb122d5downstream_cells_map)model_hubble_law_uniform_prior_given_data$276a850b-09da-41f0-8661-5d63bb2afdf8upstream_cells_mapmodel_hubble_law_uniform_prior$4f27d8d1-0d14-49e5-abee-9f277369a78bdf_fit$0991f15c-b734-480b-8b12-6d2beda547db$0c25a4cc-9f53-4729-985b-de692278c433precedence_heuristic cell_id$0c25a4cc-9f53-4729-985b-de692278c433downstream_cells_mapresponse_1d$a98ace81-9afd-4618-88a6-fadaf04eb364upstream_cells_mapmissing$ad758d43-b135-418e-a29c-36cea3a24457precedence_heuristic cell_id$ad758d43-b135-418e-a29c-36cea3a24457downstream_cells_mapupstream_cells_map@md_strasidetipgetindex$18edb530-8472-4932-946c-917d8ab9dbdeprecedence_heuristic cell_id$18edb530-8472-4932-946c-917d8ab9dbdedownstream_cells_mapupstream_cells_mapchain_lognormal_prior$cf6fbe84-9a51-4ddb-9def-c22893a17ed3describe!===nothing$2b3729fc-62b2-478d-b35c-e69e6b0efc93precedence_heuristic cell_id$2b3729fc-62b2-478d-b35c-e69e6b0efc93downstream_cells_mapplot_hubble_diagram$316f65a1-9062-4cb0-a3fb-af54149afb4d$d95d2e2d-c6ff-4cab-8b49-432ef61065bc$b54f8319-dce6-4009-873e-bd47dc2a4619upstream_cells_map string@assertTupleReallengthnamesrange/∈inroundscatter!Base.AssertionError@docBase.throwBool#___this_pluto_module_namesizexlabel!predict_linear_model$a1eb8ff3-ff5b-4e50-a75c-480cf2a7d861AbstractDataFrameBase-reshape>=plot*Unionmaximumylabel!$92af3764-6c53-4e03-a338-81b1de202a04precedence_heuristic cell_id$92af3764-6c53-4e03-a338-81b1de202a04downstream_cells_mapupstream_cells_map@md_strgetindex$0a998eb0-4bf5-4d4e-b193-54964b54f1abprecedence_heuristic cell_id$0a998eb0-4bf5-4d4e-b193-54964b54f1abdownstream_cells_mapupstream_cells_map@md_strgetindex$8ab41174-d7d9-42c7-8129-381c54906280precedence_heuristic cell_id$8ab41174-d7d9-42c7-8129-381c54906280downstream_cells_mapupstream_cells_map@md_strgetindex$7f61c098-1c54-4c82-be40-f9fd154156a1precedence_heuristic cell_id$7f61c098-1c54-4c82-be40-f9fd154156a1downstream_cells_mapupstream_cells_map@md_strdetailsgetindex$ef78506f-0326-4bd0-8fb3-af6f34a06841precedence_heuristic cell_id$ef78506f-0326-4bd0-8fb3-af6f34a06841downstream_cells_mapupstream_cells_map@md_strasidetipgetindex$316f65a1-9062-4cb0-a3fb-af54149afb4dprecedence_heuristic cell_id$316f65a1-9062-4cb0-a3fb-af54149afb4ddownstream_cells_mapupstream_cells_mapplt_resid$d70d95e9-f0e7-4225-8e06-80735278cbcaplot_hubble_diagram$2b3729fc-62b2-478d-b35c-e69e6b0efc93df_all$da9eb594-e0fe-4e9c-be69-4150c9a44fb9plt_yerr$d70d95e9-f0e7-4225-8e06-80735278cbcaH₀_plt_1$d70d95e9-f0e7-4225-8e06-80735278cbca$8ddbc921-4541-4c45-90c5-5c51ec2e391aprecedence_heuristic cell_id$8ddbc921-4541-4c45-90c5-5c51ec2e391adownstream_cells_mapols_all$ac47f631-787d-4a3c-96cf-5fee86a3b2f4upstream_cells_mapStatsModels~df_all$da9eb594-e0fe-4e9c-be69-4150c9a44fb9+lmStatsModels.ConstantTerm@formulaStatsModels.Term$d67f0b9e-10be-4511-93ab-edceb3aef49cprecedence_heuristic cell_id$d67f0b9e-10be-4511-93ab-edceb3aef49cdownstream_cells_mapH₀_wls_all$ce45a056-a949-4a9c-850e-9b5a997ecb0aupstream_cells_mapwls_all$b4a68a89-bd1f-49d5-adb2-6239317a8b82/coef$1eb1face-a795-4cb3-9c51-f995f73e2509precedence_heuristic cell_id$1eb1face-a795-4cb3-9c51-f995f73e2509downstream_cells_mapupstream_cells_mapminimumchain_uniform_prior$276a850b-09da-41f0-8661-5d63bb2afdf8ceilxlabel!nothingfloorrangetitle!!===histogram_or_density$18106638-c96c-47d6-990d-b35b6b0ba218maximum$ccdf5418-12e4-49ee-bdea-90efc202583aprecedence_heuristic cell_id$ccdf5418-12e4-49ee-bdea-90efc202583adownstream_cells_mapupstream_cells_map!sqrtresponse_1f$69ca76b8-8766-47fe-89d6-31fcc7a60e88H₀_fit$4dff84e9-3ca6-4ffb-893b-ef7ec03939b8ismissingdf_fit$0991f15c-b734-480b-8b12-6d2beda547dbloss$77b17593-6a28-4294-9b17-41303d57d5a9-/^+*$2efb7fb2-0402-4eec-b58a-a52fd530e3f5precedence_heuristic cell_id$2efb7fb2-0402-4eec-b58a-a52fd530e3f5downstream_cells_mapupstream_cells_map@md_strgetindex$3a06141a-a92d-4c31-9950-13f0bf2c023aprecedence_heuristic cell_id$3a06141a-a92d-4c31-9950-13f0bf2c023adownstream_cells_mapupstream_cells_mapstill_missingismissingresponse_1g$c43f2686-da2d-4b67-b87e-7c22bd09257c$bf7ca6ea-4c07-46e4-baff-b205e1663505precedence_heuristic cell_id$bf7ca6ea-4c07-46e4-baff-b205e1663505downstream_cells_mapupstream_cells_map@md_strtopic$95bc4b33-d619-4a95-af19-5771d209336dtitle$95bc4b33-d619-4a95-af19-5771d209336dgetindex$2656d886-60dc-44a7-b318-100335dcf05bprecedence_heuristic cell_id$2656d886-60dc-44a7-b318-100335dcf05bdownstream_cells_mapformula_wo_intercept$b4a68a89-bd1f-49d5-adb2-6239317a8b82$0991f15c-b734-480b-8b12-6d2beda547dbupstream_cells_mapStatsModels~+StatsModels.ConstantTermStatsModels.Term@formula$2c1c7d8a-eda6-4900-9be7-e4b908fed81aprecedence_heuristic cell_id$2c1c7d8a-eda6-4900-9be7-e4b908fed81adownstream_cells_mapnum_mcmc_samples$276a850b-09da-41f0-8661-5d63bb2afdf8$cf6fbe84-9a51-4ddb-9def-c22893a17ed3num_chains$276a850b-09da-41f0-8661-5d63bb2afdf8$cf6fbe84-9a51-4ddb-9def-c22893a17ed3upstream_cells_map$a1eb8ff3-ff5b-4e50-a75c-480cf2a7d861precedence_heuristic cell_id$a1eb8ff3-ff5b-4e50-a75c-480cf2a7d861downstream_cells_mappredict_linear_model$2b3729fc-62b2-478d-b35c-e69e6b0efc93upstream_cells_map@docBase.throwAbstractVectorAbstractMatrix@assertTuple#___this_pluto_module_namesizelengthBase*Union==Base.AssertionError$9ef67704-5063-4f75-853b-078261e50684precedence_heuristic cell_id$9ef67704-5063-4f75-853b-078261e50684downstream_cells_mapupstream_cells_map@md_strgetindex$095bdd00-ffa9-434d-9fc5-0978019b060bprecedence_heuristic cell_id$095bdd00-ffa9-434d-9fc5-0978019b060bdownstream_cells_mapupstream_cells_map$3f2ef930-d367-11ef-3b5c-95b4989c2e21precedence_heuristic cell_id$3f2ef930-d367-11ef-3b5c-95b4989c2e21downstream_cells_mapfilename$4518f1f7-a2a3-4783-86ae-4413fc4291b5url$d81347d3-9246-42d9-9a85-42c773ef2a1dupstream_cells_map Downloads$56d8e6c4-1f19-4bb4-a3db-7bf00675efd9Downloads.download@assertfilesizejoinpathisless>!endtitle$95bc4b33-d619-4a95-af19-5771d209336dweek$95bc4b33-d619-4a95-af19-5771d209336dgetindex$a2b56e1c-9735-49b0-91da-40f489823390precedence_heuristic cell_id$a2b56e1c-9735-49b0-91da-40f489823390downstream_cells_mapupstream_cells_map@md_strasidetipgetindex$4dcd9f29-8dd5-40ff-9ce9-50bd38f33348precedence_heuristic cell_id$4dcd9f29-8dd5-40ff-9ce9-50bd38f33348downstream_cells_mapupstream_cells_map@md_strgetindex$e4ca4972-42c8-4426-94ae-e316662513c3precedence_heuristic cell_id$e4ca4972-42c8-4426-94ae-e316662513c3downstream_cells_mapupstream_cells_map@md_strgetindex$6a75dd68-5034-4b16-89a0-7388cfae0cbdprecedence_heuristic cell_id$6a75dd68-5034-4b16-89a0-7388cfae0cbddownstream_cells_mapupstream_cells_map@md_strgetindex$d68dd4b4-c091-4435-a286-3de589b19009precedence_heuristic cell_id$d68dd4b4-c091-4435-a286-3de589b19009downstream_cells_mapupstream_cells_map@md_strgetindex$33045766-a7c5-42ff-bdf8-bde1fce9c56bprecedence_heuristic cell_id$33045766-a7c5-42ff-bdf8-bde1fce9c56bdownstream_cells_mapupstream_cells_mapresponse_2a$4eb123a1-f976-4fad-9a96-e3a61bdbf264ismissingstill_missing$9555ae09-be3d-48e8-9ee2-face2bef8cceprecedence_heuristic cell_id$9555ae09-be3d-48e8-9ee2-face2bef8ccedownstream_cells_mapupstream_cells_map@md_strgetindex$466dacda-1835-42fd-b8f4-cac5c22f6756precedence_heuristic cell_id$466dacda-1835-42fd-b8f4-cac5c22f6756downstream_cells_mapols_try1upstream_cells_mapStatsModels~df_all$da9eb594-e0fe-4e9c-be69-4150c9a44fb9lmStatsModels.Term@formula$915fc532-7d8d-4ff5-81fd-af8f6e9e64daprecedence_heuristic cell_id$915fc532-7d8d-4ff5-81fd-af8f6e9e64dadownstream_cells_mapupstream_cells_mapminimumchain_lognormal_prior$cf6fbe84-9a51-4ddb-9def-c22893a17ed3vcatchain_uniform_prior$276a850b-09da-41f0-8661-5d63bb2afdf8ceilxlabel!nothingfloorrangeplottitle!!===histogram_or_density$18106638-c96c-47d6-990d-b35b6b0ba218maximum$fb2545e1-d624-4669-9e07-c94a598ba8edprecedence_heuristic cell_id$fb2545e1-d624-4669-9e07-c94a598ba8eddownstream_cells_mapupstream_cells_mapresponse_2c$089c50b6-5494-4f27-af3e-a8af972e2869ismissingstill_missing$67e9aa94-2d3a-4888-97b7-c2cb019eb04aprecedence_heuristic cell_id$67e9aa94-2d3a-4888-97b7-c2cb019eb04adownstream_cells_mapupstream_cells_mapH₀_bf_all$f76f4c3c-74f1-44b0-9b3a-5895aec22e70!sqrtresponse_1f$69ca76b8-8766-47fe-89d6-31fcc7a60e88ismissingloss$77b17593-6a28-4294-9b17-41303d57d5a9-/df_all$da9eb594-e0fe-4e9c-be69-4150c9a44fb9+^*$389e8239-6f62-46fb-8b13-d336516e5e2fprecedence_heuristic cell_id$389e8239-6f62-46fb-8b13-d336516e5e2fdownstream_cells_mapupstream_cells_mapLogNormallogtitle!plotxlabel!ylabel!$54b51239-59b0-4767-9bc8-018d77ee5e44precedence_heuristic cell_id$54b51239-59b0-4767-9bc8-018d77ee5e44downstream_cells_mapupstream_cells_map@md_strgetindex$d62cf934-a592-482b-b0b3-7fd55f9e0617precedence_heuristic cell_id$d62cf934-a592-482b-b0b3-7fd55f9e0617downstream_cells_mapupstream_cells_map@md_strgetindex$aab05339-dea2-4c7f-8b38-bc7d7971f1d9precedence_heuristic cell_id$aab05339-dea2-4c7f-8b38-bc7d7971f1d9downstream_cells_mapupstream_cells_map@md_strgetindex$46eb1088-0265-46a6-91ec-44dfa3bd862dprecedence_heuristic cell_id$46eb1088-0265-46a6-91ec-44dfa3bd862ddownstream_cells_mapupstream_cells_mapscatterdf_all$da9eb594-e0fe-4e9c-be69-4150c9a44fb9$25ef65bd-8b86-4f88-b0b7-32f05ce82adaprecedence_heuristic cell_id$25ef65bd-8b86-4f88-b0b7-32f05ce82adadownstream_cells_mapresponse_1a$5974dfc3-8416-4bae-91de-d3f0a4f0c7dbupstream_cells_mapmissingcell_execution_order $56d8e6c4-1f19-4bb4-a3db-7bf00675efd9$95bc4b33-d619-4a95-af19-5771d209336d$e09e30f7-838f-4571-bb9c-6d0b425807a8$bf7ca6ea-4c07-46e4-baff-b205e1663505$aaa80da1-f7de-42ed-9691-c367e6bd9a28$a8f89fe8-110b-452e-9ecb-53769aa02a2e$7be54420-b2ec-4ef4-adc8-848817305d07$d19e7f11-d1ff-4624-8603-38ffec0c7bea$3f2ef930-d367-11ef-3b5c-95b4989c2e21$cd4fc36e-74f1-4b33-b92c-f0b7913a7203$2efb7fb2-0402-4eec-b58a-a52fd530e3f5$0735d5d6-9328-47f7-a764-24b4b8d92ca3$4518f1f7-a2a3-4783-86ae-4413fc4291b5$f0220da0-6f3a-438c-b533-acb346382418$ef78506f-0326-4bd0-8fb3-af6f34a06841$0a998eb0-4bf5-4d4e-b193-54964b54f1ab$d81347d3-9246-42d9-9a85-42c773ef2a1d$b15b1913-d265-4a35-b43e-67df53c93329$4655178c-a633-4b50-8f92-6ee028e21d0e$92ca099c-b81b-4cfe-a3fa-60833657f4ef$daa349a6-545b-4dfc-a3b0-92884da98e5e$9f4d8983-d2d7-4bc9-bb31-5460a18dac5b$1c35257a-9b98-4af6-bde0-a945665c91cd$4bdd7b2d-9a67-4161-a0a0-93073d9edcf0$2f7c11e8-d429-41eb-873e-ae554e028195$25ef65bd-8b86-4f88-b0b7-32f05ce82ada$5974dfc3-8416-4bae-91de-d3f0a4f0c7db$18dd7cac-fac3-4c40-9230-f6387e41751b$095bdd00-ffa9-434d-9fc5-0978019b060b$03fc47ea-9246-4cb8-8398-612fc2fb09fa$0c0312da-558a-4bc1-9e4a-b02ecd1c0572$4a2c109d-8e18-4561-913b-769f7eaeb6e9$3aba1134-7c60-41f6-a1ac-73306efb7824$efde0b1e-ec47-4180-aba7-ed8843395b50$299e9d1f-bbb2-4fc3-8682-51c3daab0ba3$0e42628b-95a9-45a0-b375-51f62dc0ce0d$d70d95e9-f0e7-4225-8e06-80735278cbca$0f6aa1eb-2f89-400f-b1e8-bfc07012179d$8ab41174-d7d9-42c7-8129-381c54906280$77b17593-6a28-4294-9b17-41303d57d5a9$beeeb261-4dcb-4093-b1cd-08682785eade$a1b5939b-b1ed-4d7c-907d-5bf3e7dd4479$098afd51-2c39-42ef-91bf-cfc326827ea1$0c25a4cc-9f53-4729-985b-de692278c433$a98ace81-9afd-4618-88a6-fadaf04eb364$5c9b0395-4d2e-4679-8d86-b9581fa0f66c$e02dd163-59bc-4d48-ae43-6fc559536aa6$92af3764-6c53-4e03-a338-81b1de202a04$2786b37a-207f-4831-983e-81a380c37045$ac1a64e8-7a5c-43b5-923a-633dac4fb664$c6f03d5d-18e2-4491-90fe-fe0008f774c4$6a75dd68-5034-4b16-89a0-7388cfae0cbd$a2b56e1c-9735-49b0-91da-40f489823390$9ef67704-5063-4f75-853b-078261e50684$b7e2903b-be56-4e36-885f-e16a9ddc94c6$f90594d8-92a1-4b0b-93d0-0e41e76b3f55$bd1c72d5-1d0d-4e36-aba4-95c9c9e4fab1$629dae3a-46ec-40dc-b9ed-92de9eb21adb$e4ca4972-42c8-4426-94ae-e316662513c3$a7981271-0b12-4b75-b5f0-5e91af761419$8bfe0413-782b-4244-8624-a9f95e393d17$92238689-d9df-46c2-bd72-3fcc41a22594$28b844a8-650a-47d8-b9f6-586f22d6f603$2656d886-60dc-44a7-b318-100335dcf05b$20f739ff-cc5b-460e-851a-141003d0c9d0$9555ae09-be3d-48e8-9ee2-face2bef8cce$34ccd3d7-e914-4f6f-a8eb-f421ed17e709$0c20b9a7-69a9-49c8-9ba4-45b8770393d0$4dcd9f29-8dd5-40ff-9ce9-50bd38f33348$69ca76b8-8766-47fe-89d6-31fcc7a60e88$0adf35a1-d7a8-4e7b-ad3f-ac83bd44ef18$ddf74f60-42eb-4eab-b77b-b6e78f0e90a7$c43f2686-da2d-4b67-b87e-7c22bd09257c$3a06141a-a92d-4c31-9950-13f0bf2c023a$78eb6aff-3752-41d8-a4d9-0cc2b31fd56d$d68dd4b4-c091-4435-a286-3de589b19009$3ec3e49c-04ce-4ffa-b1ca-8c58df0307dd$d62cf934-a592-482b-b0b3-7fd55f9e0617$16b0a2c1-5803-4310-9ee7-884b1333fb74$22da4ae4-6aea-4d45-ab2c-6fe935783e41$c6af6d4a-ff58-4012-b4b2-4e3b09d53638$42366f09-d437-4ba2-93ec-f7b8278eb09d$751eff6b-2ed2-4073-b8b5-02ad9981edb6$e4597b96-6e06-44c2-a845-3d9967d2423c$5e635a20-e205-48f0-908d-cb6d0c55b43b$3bdb27da-451c-4260-bf43-5568985d24f2$4f11d0e2-7126-47f1-865a-3474298597dc$ad758d43-b135-418e-a29c-36cea3a24457$d3e82a29-c1b5-4f5a-ba96-a64d22636493$a8d72fee-302c-44f6-8ca3-9e0b1a130ec6$d0f003fa-f72d-4a6f-ba60-1cad3844ef6a$7f61c098-1c54-4c82-be40-f9fd154156a1$92b2c9f6-5102-4e1a-aff1-05e9ce8c19a8$a6cf5e52-8692-4e3c-95aa-11167b118bd4$c433e938-31f7-4417-b3c0-a2f49fea8640$54b51239-59b0-4767-9bc8-018d77ee5e44$18106638-c96c-47d6-990d-b35b6b0ba218$a56871c9-affa-4612-b3f5-66363720f274$8764ccfe-81cd-4bff-985a-c112fc5c902f$4eb123a1-f976-4fad-9a96-e3a61bdbf264$33045766-a7c5-42ff-bdf8-bde1fce9c56b$49090825-2843-4b3e-8ee6-658187232f96$7b5b3e8a-61ac-43cf-8b85-4a8a8e4846f6$389e8239-6f62-46fb-8b13-d336516e5e2f$d3db2bc4-a26b-4e01-a193-b155a92dab48$eb72dda6-69f0-4538-a75a-25967653759f$0144843f-b649-4ea2-a546-e8c14a928670$0236d055-0d6d-440d-a44b-0d5cb0ab2273$2d97268e-560b-4398-982e-1cf5cd9ab7dd$d5c3923b-8515-4002-9d93-6768f087aef6$089c50b6-5494-4f27-af3e-a8af972e2869$fb2545e1-d624-4669-9e07-c94a598ba8ed$337a560b-49f3-4fc5-926f-1689911eacf2$23ee2544-2a2f-4302-8a62-2dd868cbb5d1$a01c7aad-429d-447b-8943-a4e89dbac6f6$c78f0ad9-03c0-4066-a390-193e7bfa84b1$adf117fe-0f54-4f2a-b107-0635f56263f2$e0e610d7-a237-42ea-8706-8e09858fba03$48e246bc-4937-40d6-b686-73524f3ce057$2532e9de-8fa7-44ce-8502-c2473928aa8e$2c1c7d8a-eda6-4900-9be7-e4b908fed81a$a7d8f517-4117-4585-8089-117510c944c5$da9eb594-e0fe-4e9c-be69-4150c9a44fb9$46eb1088-0265-46a6-91ec-44dfa3bd862d$0d0917aa-674e-47c4-a8dc-4d1585ebb70e$bd560d5d-3bc7-4b4d-aab7-67b1884a974b$d9058d92-b1c9-44f2-af5c-3540932dec48$fa745b36-7465-443e-9f50-c786fe418c6a$466dacda-1835-42fd-b8f4-cac5c22f6756$8ddbc921-4541-4c45-90c5-5c51ec2e391a$ac47f631-787d-4a3c-96cf-5fee86a3b2f4$b4a68a89-bd1f-49d5-adb2-6239317a8b82$d67f0b9e-10be-4511-93ab-edceb3aef49c$2482544e-ad05-495c-9787-0797939445be$0991f15c-b734-480b-8b12-6d2beda547db$6c86d0f1-5c6f-4204-9414-60577164f4b8$cf6fbe84-9a51-4ddb-9def-c22893a17ed3$18edb530-8472-4932-946c-917d8ab9dbde$4f27d8d1-0d14-49e5-abee-9f277369a78b$340dc0ff-dd1c-44e0-907c-f357adb122d5$276a850b-09da-41f0-8661-5d63bb2afdf8$1eb1face-a795-4cb3-9c51-f995f73e2509$0a4aaa8a-20da-486f-ace5-fbadcf646609$2ad2563a-8a72-47e8-a36d-f6d2a49cb517$915fc532-7d8d-4ff5-81fd-af8f6e9e64da$aab05339-dea2-4c7f-8b38-bc7d7971f1d9$a1eb8ff3-ff5b-4e50-a75c-480cf2a7d861$556d8f17-8a03-42c8-82c7-52c8125b604d$2f3d492c-4a5e-4a56-a773-e8078e14476a$f76f4c3c-74f1-44b0-9b3a-5895aec22e70$67e9aa94-2d3a-4888-97b7-c2cb019eb04a$9db3ff24-6f54-4a37-b787-de4b14694d95$c086c05f-ad73-44c7-bf8b-a304d3533690$1d70572f-214e-4161-b804-2c3b7c077c56$4dff84e9-3ca6-4ffb-893b-ef7ec03939b8$ce45a056-a949-4a9c-850e-9b5a997ecb0a$ccdf5418-12e4-49ee-bdea-90efc202583a$2b3729fc-62b2-478d-b35c-e69e6b0efc93$316f65a1-9062-4cb0-a3fb-af54149afb4d$d95d2e2d-c6ff-4cab-8b49-432ef61065bc$b54f8319-dce6-4009-873e-bd47dc2a4619last_hot_reload_time shortpathex1.jlprocess_statusreadypath"/home/runner/work/lab3/lab3/ex1.jlpluto_versionv0.19.47last_save_timeA, \cell_order $95bc4b33-d619-4a95-af19-5771d209336d$e09e30f7-838f-4571-bb9c-6d0b425807a8$bf7ca6ea-4c07-46e4-baff-b205e1663505$aaa80da1-f7de-42ed-9691-c367e6bd9a28$a8f89fe8-110b-452e-9ecb-53769aa02a2e$7be54420-b2ec-4ef4-adc8-848817305d07$d19e7f11-d1ff-4624-8603-38ffec0c7bea$3f2ef930-d367-11ef-3b5c-95b4989c2e21$cd4fc36e-74f1-4b33-b92c-f0b7913a7203$2efb7fb2-0402-4eec-b58a-a52fd530e3f5$0735d5d6-9328-47f7-a764-24b4b8d92ca3$4518f1f7-a2a3-4783-86ae-4413fc4291b5$f0220da0-6f3a-438c-b533-acb346382418$ef78506f-0326-4bd0-8fb3-af6f34a06841$0a998eb0-4bf5-4d4e-b193-54964b54f1ab$d81347d3-9246-42d9-9a85-42c773ef2a1d$b15b1913-d265-4a35-b43e-67df53c93329$4655178c-a633-4b50-8f92-6ee028e21d0e$92ca099c-b81b-4cfe-a3fa-60833657f4ef$daa349a6-545b-4dfc-a3b0-92884da98e5e$9f4d8983-d2d7-4bc9-bb31-5460a18dac5b$1c35257a-9b98-4af6-bde0-a945665c91cd$da9eb594-e0fe-4e9c-be69-4150c9a44fb9$4bdd7b2d-9a67-4161-a0a0-93073d9edcf0$46eb1088-0265-46a6-91ec-44dfa3bd862d$2f7c11e8-d429-41eb-873e-ae554e028195$25ef65bd-8b86-4f88-b0b7-32f05ce82ada$5974dfc3-8416-4bae-91de-d3f0a4f0c7db$18dd7cac-fac3-4c40-9230-f6387e41751b$095bdd00-ffa9-434d-9fc5-0978019b060b$03fc47ea-9246-4cb8-8398-612fc2fb09fa$0c0312da-558a-4bc1-9e4a-b02ecd1c0572$4a2c109d-8e18-4561-913b-769f7eaeb6e9$3aba1134-7c60-41f6-a1ac-73306efb7824$efde0b1e-ec47-4180-aba7-ed8843395b50$299e9d1f-bbb2-4fc3-8682-51c3daab0ba3$0e42628b-95a9-45a0-b375-51f62dc0ce0d$316f65a1-9062-4cb0-a3fb-af54149afb4d$d70d95e9-f0e7-4225-8e06-80735278cbca$0f6aa1eb-2f89-400f-b1e8-bfc07012179d$8ab41174-d7d9-42c7-8129-381c54906280$77b17593-6a28-4294-9b17-41303d57d5a9$beeeb261-4dcb-4093-b1cd-08682785eade$a1b5939b-b1ed-4d7c-907d-5bf3e7dd4479$0d0917aa-674e-47c4-a8dc-4d1585ebb70e$bd560d5d-3bc7-4b4d-aab7-67b1884a974b$098afd51-2c39-42ef-91bf-cfc326827ea1$0c25a4cc-9f53-4729-985b-de692278c433$a98ace81-9afd-4618-88a6-fadaf04eb364$5c9b0395-4d2e-4679-8d86-b9581fa0f66c$e02dd163-59bc-4d48-ae43-6fc559536aa6$92af3764-6c53-4e03-a338-81b1de202a04$2786b37a-207f-4831-983e-81a380c37045$ac1a64e8-7a5c-43b5-923a-633dac4fb664$c6f03d5d-18e2-4491-90fe-fe0008f774c4$6a75dd68-5034-4b16-89a0-7388cfae0cbd$d9058d92-b1c9-44f2-af5c-3540932dec48$fa745b36-7465-443e-9f50-c786fe418c6a$a2b56e1c-9735-49b0-91da-40f489823390$9ef67704-5063-4f75-853b-078261e50684$2f3d492c-4a5e-4a56-a773-e8078e14476a$b7e2903b-be56-4e36-885f-e16a9ddc94c6$f76f4c3c-74f1-44b0-9b3a-5895aec22e70$d95d2e2d-c6ff-4cab-8b49-432ef61065bc$f90594d8-92a1-4b0b-93d0-0e41e76b3f55$bd1c72d5-1d0d-4e36-aba4-95c9c9e4fab1$629dae3a-46ec-40dc-b9ed-92de9eb21adb$e4ca4972-42c8-4426-94ae-e316662513c3$a7981271-0b12-4b75-b5f0-5e91af761419$8bfe0413-782b-4244-8624-a9f95e393d17$466dacda-1835-42fd-b8f4-cac5c22f6756$92238689-d9df-46c2-bd72-3fcc41a22594$8ddbc921-4541-4c45-90c5-5c51ec2e391a$ac47f631-787d-4a3c-96cf-5fee86a3b2f4$28b844a8-650a-47d8-b9f6-586f22d6f603$2656d886-60dc-44a7-b318-100335dcf05b$b4a68a89-bd1f-49d5-adb2-6239317a8b82$d67f0b9e-10be-4511-93ab-edceb3aef49c$20f739ff-cc5b-460e-851a-141003d0c9d0$9555ae09-be3d-48e8-9ee2-face2bef8cce$34ccd3d7-e914-4f6f-a8eb-f421ed17e709$0991f15c-b734-480b-8b12-6d2beda547db$0c20b9a7-69a9-49c8-9ba4-45b8770393d0$b54f8319-dce6-4009-873e-bd47dc2a4619$2482544e-ad05-495c-9787-0797939445be$4dcd9f29-8dd5-40ff-9ce9-50bd38f33348$69ca76b8-8766-47fe-89d6-31fcc7a60e88$0adf35a1-d7a8-4e7b-ad3f-ac83bd44ef18$ddf74f60-42eb-4eab-b77b-b6e78f0e90a7$c43f2686-da2d-4b67-b87e-7c22bd09257c$3a06141a-a92d-4c31-9950-13f0bf2c023a$78eb6aff-3752-41d8-a4d9-0cc2b31fd56d$d68dd4b4-c091-4435-a286-3de589b19009$ce45a056-a949-4a9c-850e-9b5a997ecb0a$3ec3e49c-04ce-4ffa-b1ca-8c58df0307dd$67e9aa94-2d3a-4888-97b7-c2cb019eb04a$ccdf5418-12e4-49ee-bdea-90efc202583a$d62cf934-a592-482b-b0b3-7fd55f9e0617$1d70572f-214e-4161-b804-2c3b7c077c56$4dff84e9-3ca6-4ffb-893b-ef7ec03939b8$16b0a2c1-5803-4310-9ee7-884b1333fb74$22da4ae4-6aea-4d45-ab2c-6fe935783e41$c6af6d4a-ff58-4012-b4b2-4e3b09d53638$42366f09-d437-4ba2-93ec-f7b8278eb09d$751eff6b-2ed2-4073-b8b5-02ad9981edb6$e4597b96-6e06-44c2-a845-3d9967d2423c$5e635a20-e205-48f0-908d-cb6d0c55b43b$3bdb27da-451c-4260-bf43-5568985d24f2$4f11d0e2-7126-47f1-865a-3474298597dc$ad758d43-b135-418e-a29c-36cea3a24457$4f27d8d1-0d14-49e5-abee-9f277369a78b$d3e82a29-c1b5-4f5a-ba96-a64d22636493$a8d72fee-302c-44f6-8ca3-9e0b1a130ec6$d0f003fa-f72d-4a6f-ba60-1cad3844ef6a$7f61c098-1c54-4c82-be40-f9fd154156a1$92b2c9f6-5102-4e1a-aff1-05e9ce8c19a8$340dc0ff-dd1c-44e0-907c-f357adb122d5$a6cf5e52-8692-4e3c-95aa-11167b118bd4$c433e938-31f7-4417-b3c0-a2f49fea8640$276a850b-09da-41f0-8661-5d63bb2afdf8$54b51239-59b0-4767-9bc8-018d77ee5e44$18106638-c96c-47d6-990d-b35b6b0ba218$1eb1face-a795-4cb3-9c51-f995f73e2509$a56871c9-affa-4612-b3f5-66363720f274$0a4aaa8a-20da-486f-ace5-fbadcf646609$8764ccfe-81cd-4bff-985a-c112fc5c902f$4eb123a1-f976-4fad-9a96-e3a61bdbf264$33045766-a7c5-42ff-bdf8-bde1fce9c56b$49090825-2843-4b3e-8ee6-658187232f96$2ad2563a-8a72-47e8-a36d-f6d2a49cb517$7b5b3e8a-61ac-43cf-8b85-4a8a8e4846f6$389e8239-6f62-46fb-8b13-d336516e5e2f$d3db2bc4-a26b-4e01-a193-b155a92dab48$6c86d0f1-5c6f-4204-9414-60577164f4b8$cf6fbe84-9a51-4ddb-9def-c22893a17ed3$18edb530-8472-4932-946c-917d8ab9dbde$eb72dda6-69f0-4538-a75a-25967653759f$915fc532-7d8d-4ff5-81fd-af8f6e9e64da$0144843f-b649-4ea2-a546-e8c14a928670$0236d055-0d6d-440d-a44b-0d5cb0ab2273$2d97268e-560b-4398-982e-1cf5cd9ab7dd$d5c3923b-8515-4002-9d93-6768f087aef6$089c50b6-5494-4f27-af3e-a8af972e2869$fb2545e1-d624-4669-9e07-c94a598ba8ed$337a560b-49f3-4fc5-926f-1689911eacf2$23ee2544-2a2f-4302-8a62-2dd868cbb5d1$a01c7aad-429d-447b-8943-a4e89dbac6f6$c78f0ad9-03c0-4066-a390-193e7bfa84b1$adf117fe-0f54-4f2a-b107-0635f56263f2$e0e610d7-a237-42ea-8706-8e09858fba03$48e246bc-4937-40d6-b686-73524f3ce057$2532e9de-8fa7-44ce-8502-c2473928aa8e$56d8e6c4-1f19-4bb4-a3db-7bf00675efd9$2c1c7d8a-eda6-4900-9be7-e4b908fed81a$a7d8f517-4117-4585-8089-117510c944c5$aab05339-dea2-4c7f-8b38-bc7d7971f1d9$a1eb8ff3-ff5b-4e50-a75c-480cf2a7d861$556d8f17-8a03-42c8-82c7-52c8125b604d$9db3ff24-6f54-4a37-b787-de4b14694d95$c086c05f-ad73-44c7-bf8b-a304d3533690$2b3729fc-62b2-478d-b35c-e69e6b0efc93published_objectsnbpkginstall_time_ns t0instantiatedòinstalled_versionsCSV0.10.15DownloadsstdlibMCMCChains6.0.7StatsBase0.34.4PDMats0.11.31Distributions0.25.116PlutoUI0.7.60Plots1.40.9StatsPlots0.15.7PlutoTeachingTools0.3.1LinearAlgebrastdlibDataFrames1.7.0Turing0.35.5LaTeXStrings1.4.0GLM1.9.0terminal_outputs CSV@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`Downloads@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`MCMCChains@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`PDMats@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`Distributions@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`PlutoUI@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`StatsPlots@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`PlutoTeachingTools@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`LinearAlgebra@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`LaTeXStrings@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`StatsBase@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`Plots@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`nbpkg_sync@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`DataFrames@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`Turing@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`GLM@
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Project.toml`
[32m[1m No Changes[22m[39m to `/tmp/jl_6itfGm/Manifest.toml`
[0m[1mPrecompiling...[22m
[90m===[39m
[32m[1m Activating[22m[39m project at `/tmp/jl_6itfGm`enabled÷restart_recommended_msgrestart_required_msgbusy_packageswaiting_for_permission,waiting_for_permission_but_probably_disabled«cell_inputs $daa349a6-545b-4dfc-a3b0-92884da98e5ecell_id$daa349a6-545b-4dfc-a3b0-92884da98e5ecodeBdetails(md"""Luminosity distance from Distance modulus""",
md"""
The standard definition of the [distance modulus](https://en.wikipedia.org/wiki/Distance_modulus) is $μ = 5 \log_{10}\left(\frac{d}{\mathrm{pc}}\right) - 5$.
We'll work in megaparsecs (Mpc), so we get
$μ = 5 \log_{10}\left(\frac{d}{\mathrm{Mpc}}\right) + 25$.
Solving for $d$ gives
$d = 10^{(\mu-25)/5}$
Since we're measuring distances based on the luminosity of supernovae, we're using a [luminosity distance](https://en.wikipedia.org/wiki/Luminosity_distance) here. Hence I've included the subscript L.)
""")metadatashow_logsèdisabled®skip_as_script«code_folded$556d8f17-8a03-42c8-82c7-52c8125b604dcell_id$556d8f17-8a03-42c8-82c7-52c8125b604dcode[begin
"""
`calc_mle_linear_model(A, y_obs, covar)`
Computes the maximum likelihood estimator for θ for the linear model
`y = A θ`
where measurements errors of `y_obs` are normally distributed with covariance matrix `covar` about the true values of `y`.
"""
function calc_mle_linear_model end
function calc_mle_linear_model(A::AbstractMatrix, y_obs::AbstractVector, covar::PDMat)
@assert size(A,1) == length(y_obs) == size(covar,1) == size(covar,2)
@assert size(A,2) >= 1
(A' * (covar \ A)) \ (A' * (covar \ y_obs) )
end
function calc_mle_linear_model(A::AbstractMatrix, y_obs::AbstractVector, covar::AbstractMatrix)
calc_mle_linear_model(A,y_obs,PDMat(covar))
end
function calc_mle_linear_model(A::AbstractMatrix, y_obs::AbstractVector, covar::AbstractVector)
calc_mle_linear_model(A,y_obs,PDiagMat(covar))
end
Docs.doc(calc_mle_linear_model)
endmetadatashow_logsèdisabled®skip_as_script«code_folded$77b17593-6a28-4294-9b17-41303d57d5a9cell_id$77b17593-6a28-4294-9b17-41303d57d5a9codejfunction loss(H::Real, data::AbstractDataFrame)
sum(abs2.((data.dₗ - data.v / H ) ./data.σ_dₗ ))
endmetadatashow_logsèdisabled®skip_as_script«code_folded$2482544e-ad05-495c-9787-0797939445becell_id$2482544e-ad05-495c-9787-0797939445becodeٛmd"""
Maximum v to plot: $(@bind max_v_plt PlutoUI.Slider(max(0,ceil(minimum(df_all.v)/200)*200):200:maximum(df_all.v), default=maximum(df_all.v)/2))
"""metadatashow_logsèdisabled®skip_as_script«code_folded$cd4fc36e-74f1-4b33-b92c-f0b7913a7203cell_id$cd4fc36e-74f1-4b33-b92c-f0b7913a7203codedetails(
md"""
This data is from the [Supernovae Cosmology Project](https://supernova.lbl.gov/)'s ["Union2.1" SN Ia compilation](https://supernova.lbl.gov/Union/). Click if you're curious for more details.""",
md"""It incorporated data from 833 supernovae (SNe), drawn from 19 datasets. Of these, 580 SNe passed usability cuts. We present new data from the HST Cluster Survey.
All SNe were fit using a single lightcurve fitter (SALT2-1) and uniformly analyzed.
All analyses and cuts were developed in a blind manner, i.e. with the cosmology hidden.
The Union2.1 paper is available at [Suzuki et al. (2012)](https://ui.adsabs.harvard.edu/abs/2012ApJ...746...85S/abstract).""")metadatashow_logsèdisabled®skip_as_script«code_folded$0735d5d6-9328-47f7-a764-24b4b8d92ca3cell_id$0735d5d6-9328-47f7-a764-24b4b8d92ca3codeJmd"""
The data are stored in a simple text file. There are many ways to read them. Julia has a fast reader for CSV files, so we'll demonstrate that below. For now, you don't need to pay attention to the syntax, but at some point you'll likely want to refer back to examples like this to see how to accomplish common tasks.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$beeeb261-4dcb-4093-b1cd-08682785eadecell_id$beeeb261-4dcb-4093-b1cd-08682785eadecode7aside(tip(md"""What does `abs2.(x)` mean?
First, `abs2(x)` would return the squared absolute value of x.
Second, when we're working with arrays rather than scalars, then we can specify element-wise arithmetic by adding a dot after a function name.
Similarly, we can add a dot before arithmetic operator (e.g., `./` instead of `/`) to specify element-wise arithmetic. (The difference is particularly important for multiplication when we need to distinguish between matrix-vector or matrix-matrix multiplication and element-wise multiplication.)
"""), v_offset=-100)metadatashow_logsèdisabled®skip_as_script«code_folded$da9eb594-e0fe-4e9c-be69-4150c9a44fb9cell_id$da9eb594-e0fe-4e9c-be69-4150c9a44fb9codeSbegin
added_metadata
df_all = select(df_in,[:id,:z]) # select only the columns we need
# Perform variabile transofrmations informed by astrophysics
local c = speed_of_light_mps/1000 # convention is to use km/s
df_all.v = @. c-((2*c)/((df_in.z+1)^2+1)) # see equations above
df_all.dₗ = @. 10.0^((df_in.μ-25)/5) # see equations above
df_all.σ_dₗ = @. log(10) * df_all.dₗ * df_in.σ_μ/5 # propagation of uncertainty
# Now that columns have dimensions, add those as metadata
colmetadata!(df_all, :dₗ, "label", "velocity")
colmetadata!(df_all, :v, "units", "km/s")
colmetadata!(df_all, :dₗ, "label", "luminosity distance")
colmetadata!(df_all, :σ_dₗ, "label", "uncertainty in luminosity distance")
colmetadata!(df_all, :dₗ, "units", "Mpc")
colmetadata!(df_all, :σ_dₗ, "units", "Mpc")
df_all
endmetadatashow_logsèdisabled®skip_as_script«code_folded$d19e7f11-d1ff-4624-8603-38ffec0c7beacell_id$d19e7f11-d1ff-4624-8603-38ffec0c7beacodemd"""
## Download the Data
"""metadatashow_logsèdisabled®skip_as_script«code_folded$c6f03d5d-18e2-4491-90fe-fe0008f774c4cell_id$c6f03d5d-18e2-4491-90fe-fe0008f774c4codetmd"""
A diagonal matrix containing the squared measurement uncertainties,
$\Sigma = \mathrm{diag}(\sigma_{d_{L}}^{2}) =
\begin{bmatrix}
\sigma_{d,1}^{-2} & & \\
& \sigma_{d,2}^{-2} & & \\
& & \ddots & \\
& & & \sigma_{d,N_{\mathrm{obs}}}^{-2}
\end{bmatrix},$
and the design matrix contains the velocities
$A = \left( \begin{matrix}
v_1 \\
v_2 \\
... \\
v_{N_{\mathrm{obs}}}
\end{matrix}
\right).$
In this model, we're presuming that the velocities are known exactly. While this isn't technically true, the velocities are typically measured much more precisely than the distances, so it can be a reasonable approximation.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$629dae3a-46ec-40dc-b9ed-92de9eb21adbcell_id$629dae3a-46ec-40dc-b9ed-92de9eb21adbcode.if ismissing(response_1e)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$4518f1f7-a2a3-4783-86ae-4413fc4291b5cell_id$4518f1f7-a2a3-4783-86ae-4413fc4291b5codeubegin
df_in = CSV.read(filename, DataFrame, delim="\t", skipto=6, header=[:id,:z,:μ,:σ_μ,:prob_low_mass_gal])
endmetadatashow_logsèdisabled®skip_as_script«code_folded$e4597b96-6e06-44c2-a845-3d9967d2423ccell_id$e4597b96-6e06-44c2-a845-3d9967d2423ccode.if ismissing(response_1i)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$92b2c9f6-5102-4e1a-aff1-05e9ce8c19a8cell_id$92b2c9f6-5102-4e1a-aff1-05e9ce8c19a8code٫md"""
After we've specified the general form of the statistcal model, we create a specific version of the model that is specialized for the valuse of the observations.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$6c86d0f1-5c6f-4204-9414-60577164f4b8cell_id$6c86d0f1-5c6f-4204-9414-60577164f4b8codeumodel_hubble_law_lognormal_prior_given_data = model_hubble_law_lognormal_prior(df_fit.v, df_fit.dₗ, df_fit.σ_dₗ)metadatashow_logsèdisabled®skip_as_script«code_folded$a7981271-0b12-4b75-b5f0-5e91af761419cell_id$a7981271-0b12-4b75-b5f0-5e91af761419codemd"""
In this case, we had a very simple model and the linear algebra was straight forward. In more complicated cases, there could be many input variables (here velocity), multiple observables being predicted (here distance), correlations between the different observables, and many possible models to try. Setting up all the matrices and getting all the algebra right for each model is error prone. Therefore, most modern high-level programming languages have packages that make it easy to fit linear models such as this one. Below, we'll repeat the analysis, but using Julia's [GLM.jl](https://juliastats.org/GLM.jl/stable/) package for general linear models."""metadatashow_logsèdisabled®skip_as_script«code_folded$276a850b-09da-41f0-8661-5d63bb2afdf8cell_id$276a850b-09da-41f0-8661-5d63bb2afdf8codeٿif ready_to_run_mcmc
chain_uniform_prior = sample(model_hubble_law_uniform_prior_given_data, NUTS(0.65), MCMCThreads(), num_mcmc_samples, num_chains)
else
chain_uniform_prior = nothing
end;metadatashow_logsèdisabled®skip_as_script«code_folded$bd560d5d-3bc7-4b4d-aab7-67b1884a974bcell_id$bd560d5d-3bc7-4b4d-aab7-67b1884a974bcodebegin
scatter(H₀s_plt, loss_vs_H₀_all, markercolor=:blue, label=:none)
plot!(H₀s_plt, loss_vs_H₀_all, linecolor=:blue, label=:none)
xlabel!("H₀ (km/s/Mpc)")
ylabel!("loss = χ²(H₀)")
endmetadatashow_logsèdisabled®skip_as_script«code_folded$0236d055-0d6d-440d-a44b-0d5cb0ab2273cell_id$0236d055-0d6d-440d-a44b-0d5cb0ab2273coderesponse_2b = missingmetadatashow_logsèdisabled®skip_as_script«code_folded$efde0b1e-ec47-4180-aba7-ed8843395b50cell_id$efde0b1e-ec47-4180-aba7-ed8843395b50codeقmd"""**Q1c:** Use the slider below to try different values of $H_0$. Approximately what value looks like the best-fit to you?"""metadatashow_logsèdisabled®skip_as_script«code_folded$48e246bc-4937-40d6-b686-73524f3ce057cell_id$48e246bc-4937-40d6-b686-73524f3ce057code*# hideall
md"""
# Setup & Helper Code
"""metadatashow_logsèdisabled®skip_as_script«code_folded$2786b37a-207f-4831-983e-81a380c37045cell_id$2786b37a-207f-4831-983e-81a380c37045codeomd"""
If one puts in a fair bit of algebra, then one can find an analytical expression for the value of the model parameters ($\theta$) that minimizes the likelihood for a linear model.
You're likely to see the derivation of the solution in class on staistics (or maybe linear algebra), but we'll just use the result:
$\hat{\theta} = (A' \Sigma^{-1} A)^{-1} (A' \Sigma^{-1} \ y_\mathrm{obs} ),$
where $y_{\mathrm{obs}}$ is the vector of measurement values, $\Sigma^{-1}$ is the covariance matrix, and $A$ is the **[design matrix](https://en.wikipedia.org/wiki/Design_matrix)**.
For our problem, these simplify to:
"""metadatashow_logsèdisabled®skip_as_script«code_folded$0d0917aa-674e-47c4-a8dc-4d1585ebb70ecell_id$0d0917aa-674e-47c4-a8dc-4d1585ebb70ecodebegin
# Specify a range of values to try
H₀s_plt = range(40, stop=100, step=2)
# Apply loss function for each value in the range
loss_vs_H₀_all = map(H₀->loss(H₀,df_all),H₀s_plt)
end;metadatashow_logsèdisabled®skip_as_script«code_folded$d3e82a29-c1b5-4f5a-ba96-a64d22636493cell_id$d3e82a29-c1b5-4f5a-ba96-a64d22636493codemd"""
Notice that this code is very different from a tradditional program where we'd specify how to compute each variable in order (known as the [imperative programming](https://en.wikipedia.org/wiki/Imperative_programming) model).
"""metadatashow_logsèdisabled®skip_as_script«code_folded$49090825-2843-4b3e-8ee6-658187232f96cell_id$49090825-2843-4b3e-8ee6-658187232f96code,md"""
## Sensitivity to Choice of Priors
"""metadatashow_logsèdisabled®skip_as_script«code_folded$e0e610d7-a237-42ea-8706-8e09858fba03cell_id$e0e610d7-a237-42ea-8706-8e09858fba03code.if ismissing(response_2e)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$ce45a056-a949-4a9c-850e-9b5a997ecb0acell_id$ce45a056-a949-4a9c-850e-9b5a997ecb0acode|if !ismissing(response_1f)
let
H₀s_plt_zoom = range(50, stop=75, step=0.5)
# Apply loss function for each value in the range
loss_vs_H₀_all = map(H₀->loss(H₀,df_all),H₀s_plt_zoom)
# Plot χ² surface for all data in blue
loss_H₀_bf_all = loss(H₀_bf_all,df_all)
dof_all = size(df_all,1)-1
plt = plot(xlabel = "H₀ (km/s/Mpc)",ylabel = "Δχ²", legend=:topleft)
ylims!(0,size(df_all,1)*0.4)
scatter!(plt, H₀s_plt_zoom, loss_vs_H₀_all .- loss_H₀_bf_all, markercolor=:blue, label=:none)
plot!(plt, H₀s_plt_zoom, loss_vs_H₀_all .- loss_H₀_bf_all, linecolor=:blue, label=:none)
plot!(plt,fill(H₀_wls_all,2), [0, 1000], linecolor=:blue, label="Best-fit for All" )
plot!(plt,[H₀_wls_all-σ_H₀_all,H₀_wls_all+σ_H₀_all], fill(100,2), linecolor=:blue, lw=3, label=:none )
# Plot χ² surface for selected data in red
loss_vs_H₀ = map(H₀->loss(H₀,df_fit),H₀s_plt_zoom)
loss_H₀_bf_fit = loss(H₀_wls_fit,df_fit)
dof_fit = size(df_fit,1)-1
scatter!(plt, H₀s_plt_zoom, loss_vs_H₀ .- loss_H₀_bf_fit, markercolor=:red, label=:none)
plot!(plt, H₀s_plt_zoom, loss_vs_H₀ .- loss_H₀_bf_fit, linecolor=:red, label=:none)
plot!(plt,fill(H₀_wls_fit,2), [0, 1000], linecolor=:red, label="Best-fit for Subset" )
plot!(plt,[H₀_wls_fit-σ_H₀_fit,H₀_wls_fit+σ_H₀_fit], fill(100,2), linecolor=:red, lw=3, label=:none )
end
endmetadatashow_logsèdisabled®skip_as_script«code_folded$1d70572f-214e-4161-b804-2c3b7c077c56cell_id$1d70572f-214e-4161-b804-2c3b7c077c56codeَbegin
σ_H₀inv_all = calc_sigma_mle_linear_model(A_all, Σ_all)[1]
σ_H₀_all = abs(H₀_bf_all * σ_H₀inv_all / H₀inv_bf_all[1])
endmetadatashow_logsèdisabled®skip_as_script«code_folded$23ee2544-2a2f-4302-8a62-2dd868cbb5d1cell_id$23ee2544-2a2f-4302-8a62-2dd868cbb5d1coderesponse_2d = missingmetadatashow_logsèdisabled®skip_as_script«code_folded$69ca76b8-8766-47fe-89d6-31fcc7a60e88cell_id$69ca76b8-8766-47fe-89d6-31fcc7a60e88coderesponse_1f = 0 # missingmetadatashow_logsèdisabled®skip_as_script«code_folded$18106638-c96c-47d6-990d-b35b6b0ba218cell_id$18106638-c96c-47d6-990d-b35b6b0ba218code}md"""
Select type of plot to show: $(@bind histogram_or_density Select([histogram => "histogram", density => "density"]))
"""metadatashow_logsèdisabled®skip_as_script«code_folded$f90594d8-92a1-4b0b-93d0-0e41e76b3f55cell_id$f90594d8-92a1-4b0b-93d0-0e41e76b3f55codemd"""**Q1e:** How does the best-fit value computed above compare to your previous estimate based on the plot of $\chi^2$ vs $H_0$? What is likely to explain any differences?
Does Hubble's law provide a good model for these observations? If not, explain your concern."""metadatashow_logsèdisabled®skip_as_script«code_folded$aaa80da1-f7de-42ed-9691-c367e6bd9a28cell_id$aaa80da1-f7de-42ed-9691-c367e6bd9a28code# hideall
TableOfContents()metadatashow_logsèdisabled®skip_as_script«code_folded$0adf35a1-d7a8-4e7b-ad3f-ac83bd44ef18cell_id$0adf35a1-d7a8-4e7b-ad3f-ac83bd44ef18code.if ismissing(response_1f)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$92ca099c-b81b-4cfe-a3fa-60833657f4efcell_id$92ca099c-b81b-4cfe-a3fa-60833657f4efcodeRmd"""The Hubble "law" ($v = H_0 \times d$) for the expansion of the universe is based on velocity versus distance. That's not what this data file provided, but we can compute those using standard expressions that you may have seen in an introduction to astronomy class. You learn or refresh your memory about each by clicking below."""metadatashow_logsèdisabled®skip_as_script«code_folded$03fc47ea-9246-4cb8-8398-612fc2fb09facell_id$03fc47ea-9246-4cb8-8398-612fc2fb09facoderesponse_1b = missingmetadatashow_logsèdisabled®skip_as_script«code_folded$d3db2bc4-a26b-4e01-a193-b155a92dab48cell_id$d3db2bc4-a26b-4e01-a193-b155a92dab48code@model function model_hubble_law_lognormal_prior(v, d_obs, σ_d_obs)
# Verify that size of each data vector matches
@assert size(v) == size(d_obs) == size(σ_d_obs)
# Specify priors for model parameters
H₀ ~ LogNormal(log(100), 1.0)
# Specify how to compute the model's predictions given parameter values
d_true = v / H₀
# Specify distribution for observations given the true values and measurement uncertainties
d_obs ~ MvNormal(d_true, σ_d_obs)
endmetadatashow_logsèdisabled®skip_as_script«code_folded$9db3ff24-6f54-4a37-b787-de4b14694d95cell_id$9db3ff24-6f54-4a37-b787-de4b14694d95code"""
`calc_fisher_matrix_linear_model(A, covar)`
Computes the Fisher information matrix for θ for the linear model
`y = A θ`
where measurements errors of `y_obs` are normally distributed and have covariance `covar`.
Note that `y_obs` is not passed to this function, because it does not affect results.
"""
function calc_fisher_matrix_linear_model(A::AbstractMatrix, covar::AbstractMatrix)
@assert size(A,1) == size(covar,1) == size(covar,2)
@assert size(A,2) >= 1
(A' * (covar \ A))
endmetadatashow_logsèdisabled®skip_as_script«code_folded$c6af6d4a-ff58-4012-b4b2-4e3b09d53638cell_id$c6af6d4a-ff58-4012-b4b2-4e3b09d53638code.if ismissing(response_1h)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$0c0312da-558a-4bc1-9e4a-b02ecd1c0572cell_id$0c0312da-558a-4bc1-9e4a-b02ecd1c0572code.if ismissing(response_1b)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$b15b1913-d265-4a35-b43e-67df53c93329cell_id$b15b1913-d265-4a35-b43e-67df53c93329code٤details("How to access metadata",
md"""
You can access the metadata from a `DataFrame` using function calls like `metadata(df_in)` or `colmetadata(df_in,:z)`.
""")metadatashow_logsèdisabled®skip_as_script«code_folded$4655178c-a633-4b50-8f92-6ee028e21d0ecell_id$4655178c-a633-4b50-8f92-6ee028e21d0ecode,md"""
## Apply Variable Transformations
"""metadatashow_logsèdisabled®skip_as_script«code_folded$3bdb27da-451c-4260-bf43-5568985d24f2cell_id$3bdb27da-451c-4260-bf43-5568985d24f2code md"""
As we develop more complex models, it will often not be practical to find the maximum likelihood estimator values analyticly or using a general linear modeling package.
Or there may be significant prior knowledge that we'd like to incorporate into our analysis.
Or we might want to perform a more rigorous analysis of the uncertainties on the model parameters.
In any of these cases, it can be very useful to use a probabilistic modeling language (PPL), such as [Turing](https://turinglang.org/) or [Stan](https://mc-stan.org/).
"""metadatashow_logsèdisabled®skip_as_script«code_folded$2f7c11e8-d429-41eb-873e-ae554e028195cell_id$2f7c11e8-d429-41eb-873e-ae554e028195code{md"""**Q1a:** Is this plot consistent with your expectations for velocity versus distance? If not, what surprised you?"""metadatashow_logsèdisabled®skip_as_script«code_folded$8764ccfe-81cd-4bff-985a-c112fc5c902fcell_id$8764ccfe-81cd-4bff-985a-c112fc5c902fcodemd"""**Q2a:** How do the mean and standard deviation for the Hubble constant estimated with our Bayesian model and Turing compare to the analytic estimates using linear model (and the same dataset)? Are any differences large or small compare to the reported uncertainties?"""metadatashow_logsèdisabled®skip_as_script«code_folded$b7e2903b-be56-4e36-885f-e16a9ddc94c6cell_id$b7e2903b-be56-4e36-885f-e16a9ddc94c6codezmd"""
This value is nothing like what we got above. Why? Remember that to write this in a form ammenable to the analytic solution, we predicted distances based on velocities $d = v / H_0$, so the best-fit constant of proportionality is the recipocal of $H_0$. Let's transform that back to a best-fit value for the Hubble constant, $\hat{H_0}_{,bf} = 1/\hat{\theta}_{bf}$.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$4eb123a1-f976-4fad-9a96-e3a61bdbf264cell_id$4eb123a1-f976-4fad-9a96-e3a61bdbf264coderesponse_2a = missingmetadatashow_logsèdisabled®skip_as_script«code_folded$089c50b6-5494-4f27-af3e-a8af972e2869cell_id$089c50b6-5494-4f27-af3e-a8af972e2869coderesponse_2c = missingmetadatashow_logsèdisabled®skip_as_script«code_folded$34ccd3d7-e914-4f6f-a8eb-f421ed17e709cell_id$34ccd3d7-e914-4f6f-a8eb-f421ed17e709codemd"""
You likely noticed that Hubble's "Law" wasn't a great model for this dataset.
That's because it's an incomplete model that works well in the local universe, but breaks down once astronomers look far enough away that we can see the effects of dark energy causing the rate of expansion to accelerate.
If we want to measure the local Hubble constant, then we could fit a subset of our dataset, picking those SN that have velocities less than some threshold.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$e02dd163-59bc-4d48-ae43-6fc559536aa6cell_id$e02dd163-59bc-4d48-ae43-6fc559536aa6codemd"""
For the above model, then likelihood is given by
$\mathcal{L} = \frac{1}{\sqrt{(2\pi)^{N_{\mathrm{obs}}} \prod_{i=1}^{N_{\mathrm{obs}}} \sigma_{d,i}^2} } \; e^{-\chi^2/2}$
Minimizing $\mathcal{L}$ is equivalent to minimizing $\ell = \log(\mathcal{L})$.
Since the first factor doesn't depend on the measured values of distance, we can minimize $\mathcal{L}$ and $\ell$ by minimizing $\chi^2$.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$a56871c9-affa-4612-b3f5-66363720f274cell_id$a56871c9-affa-4612-b3f5-66363720f274code5md"""
and we can compute some summary statistics.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$4f27d8d1-0d14-49e5-abee-9f277369a78bcell_id$4f27d8d1-0d14-49e5-abee-9f277369a78bcode@model function model_hubble_law_uniform_prior(v, d_obs, σ_d_obs)
# Verify that size of each data vector matches
@assert size(v) == size(d_obs) == size(σ_d_obs)
# Specify priors for model parameters
H₀ ~ Uniform(-speed_of_light_mps/1000, speed_of_light_mps/1000)
# Specify how to compute the model's predictions given parameter values
d_true = v / H₀
# Specify distribution for observations given the true values and measurement uncertainties
d_obs ~ MvNormal(d_true, σ_d_obs)
endmetadatashow_logsèdisabled®skip_as_script«code_folded$0144843f-b649-4ea2-a546-e8c14a928670cell_id$0144843f-b649-4ea2-a546-e8c14a928670code2md"""**Q2b:** How do the estimates of the posterior mean and standard deviation using two different choices for the prior compare to each other?
Are any differences large or small compare to the reported uncertainties? What does this imply about the sensitivity of our analysis to the choice of prior?"""metadatashow_logsèdisabled®skip_as_script«code_folded$a7d8f517-4117-4585-8089-117510c944c5cell_id$a7d8f517-4117-4585-8089-117510c944c5code!speed_of_light_mps = 299_792_458;metadatashow_logsèdisabled®skip_as_script«code_folded$18dd7cac-fac3-4c40-9230-f6387e41751bcell_id$18dd7cac-fac3-4c40-9230-f6387e41751bcodeپmd"""**Q1b:** Do you notice anything potentially concerning about the plot? If so, describe your concerns. Optionally, use the cell below to perform any tests to address your concerns."""metadatashow_logsèdisabled®skip_as_script«code_folded$fa745b36-7465-443e-9f50-c786fe418c6acell_id$fa745b36-7465-443e-9f50-c786fe418c6acode$Σ_all = PDiagMat(df_all.σ_dₗ.^2)metadatashow_logsèdisabled®skip_as_script«code_folded$b4a68a89-bd1f-49d5-adb2-6239317a8b82cell_id$b4a68a89-bd1f-49d5-adb2-6239317a8b82codeMwls_all = lm(formula_wo_intercept, df_all, wts=weights(df_all.σ_dₗ.^(-2)))metadatashow_logsèdisabled®skip_as_script«code_folded$a8d72fee-302c-44f6-8ca3-9e0b1a130ec6cell_id$a8d72fee-302c-44f6-8ca3-9e0b1a130ec6codedetails("Probabilistic Programming Languages",
md"""
Traditionally, programmers provide explicit instructions for what steps the program should perform in what order (known as the [imperative programming](https://en.wikipedia.org/wiki/Imperative_programming) model).
But there are alternatives.
For example, [probabilistic programming languages](https://en.wikipedia.org/wiki/Probabilistic_programming) (PPLs) allow programmers to specify target probability distributions to be computed (or sampled from), without specifying how to perform inference with that target distribution.
In this exercise, we'll demonstrate how to use the [Turing.jl](https://turinglang.org/) PPL.
Turing is implemented entirely in Julia and is highly composable, so that one can perform inference with complex models written in Julia.
In this notebook, we'll demonstrate applying it to the simple model for the Hubble expansion.
""")metadatashow_logsèdisabled®skip_as_script«code_folded$c433e938-31f7-4417-b3c0-a2f49fea8640cell_id$c433e938-31f7-4417-b3c0-a2f49fea8640codeqmd"""
Check once you're ready to run the Bayesian models: $(@bind ready_to_run_mcmc CheckBox(; default=true))
"""metadatashow_logsèdisabled®skip_as_script«code_folded$3aba1134-7c60-41f6-a1ac-73306efb7824cell_id$3aba1134-7c60-41f6-a1ac-73306efb7824codemd"""
## Goodness-of-fit statistic
Usually, we don't know the true value of all our model's parameters (in this case the Hubble constant).
Therefore, we want to use astronomical data to measure (or more generally to improve our knowldge of) the model parameters.
For any value of the model parameter(s), we can compute the residuals between the observations and the model predicitons, $\Delta_i = d_i - v_{\mathrm{obs},i}/H_0$.
**If** the model and model parameters were the perfect truth, then
$\Delta_i \sim \mathrm{Normal}(0,\sigma_i^2).$
Even if the model isn't perfect, we can compute a goodness-of-fit [statistic](https://en.wikipedia.org/wiki/Statistic),
$\chi^2 = \sum_{i=1}^{N_\mathrm{obs}} \Delta_i^2.$
If the model is perfect, then
$\chi^2 = \sum_{i=1}^{N_\mathrm{obs}} \Delta_i^2/\sigma_i^2$
and its expected value becomes
$E\left[\sum_{i=1}^{N_\mathrm{obs}} \Delta_i^2/\sigma_i^2 \right] = N_{\mathrm{obs}},$
since $E[\Delta_i^2/\sigma_i^2] = 1$ for each $i$.
If we evaluate $\chi^2$ using model parameters that differ from the truth, then the expected value of the value of $\chi^2$ statistic increases.
We can find the value of the model parameters than minimize $\chi^2$ to provide us with the **best-fit** model.
This model for observations is so common that statistican have given a name to the probabilitiy distribution that the $\chi^2$ statistic follows... the [$\chi^2$ distribution](https://en.wikipedia.org/wiki/Chi-squared_distribution).
"""metadatashow_logsèdisabled®skip_as_script«code_folded$5e635a20-e205-48f0-908d-cb6d0c55b43bcell_id$5e635a20-e205-48f0-908d-cb6d0c55b43bcode2md"""
# Bayesian Inference with a Linear Model
"""metadatashow_logsèdisabled®skip_as_script«code_folded$4dff84e9-3ca6-4ffb-893b-ef7ec03939b8cell_id$4dff84e9-3ca6-4ffb-893b-ef7ec03939b8codeLif !ismissing(response_1f)
A_fit = reshape(df_fit.v, size(df_fit,1), 1);
Σ_fit = PDiagMat(df_fit.σ_dₗ.^2);
H₀inv_fit = calc_mle_linear_model(A_fit, df_fit.dₗ, Σ_fit)[1]
σ_H₀inv_fit = calc_sigma_mle_linear_model(A_fit, Σ_fit)[1]
H₀_fit = 1/H₀inv_fit
σ_H₀_fit = abs(H₀_fit * σ_H₀inv_fit/H₀inv_fit)
endmetadatashow_logsèdisabled®skip_as_script«code_folded$7be54420-b2ec-4ef4-adc8-848817305d07cell_id$7be54420-b2ec-4ef4-adc8-848817305d07code&md"""
# Injests & Prepare the Data
"""metadatashow_logsèdisabled®skip_as_script«code_folded$16b0a2c1-5803-4310-9ee7-884b1333fb74cell_id$16b0a2c1-5803-4310-9ee7-884b1333fb74codemd"""**Q1h:** Which dataset do you think provides a better estimate of the Hubble constant? Why?
For which dataset is the formal measurement uncertainty larger? Explain your why that is."""metadatashow_logsèdisabled®skip_as_script«code_folded$0c20b9a7-69a9-49c8-9ba4-45b8770393d0cell_id$0c20b9a7-69a9-49c8-9ba4-45b8770393d0code(md"""
Use the slider below to pick a threshold velocity below which the model of a linear releationship between $v$ and $d$ appears to be accurate. The figure will automatically update to show only the data you've selected and the predictions with the best-fit $H_0$ for that subset of data.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$56d8e6c4-1f19-4bb4-a3db-7bf00675efd9cell_id$56d8e6c4-1f19-4bb4-a3db-7bf00675efd9code# hideall
begin
using Downloads
using CSV, DataFrames
using Plots, LaTeXStrings
using StatsBase, GLM, PDMats, LinearAlgebra
using Distributions, Turing, MCMCChains, StatsPlots
using PlutoUI, PlutoTeachingTools
endmetadatashow_logsèdisabled®skip_as_script«code_folded$d70d95e9-f0e7-4225-8e06-80735278cbcacell_id$d70d95e9-f0e7-4225-8e06-80735278cbcacodemd"""
H₀ to plot: $(@bind H₀_plt_1 PlutoUI.Slider(10:100.0, default=65.0))
Show uncertainties: $(@bind plt_yerr PlutoUI.CheckBox())
Show residuals: $(@bind plt_resid PlutoUI.CheckBox())
"""metadatashow_logsèdisabled®skip_as_script«code_folded$9f4d8983-d2d7-4bc9-bb31-5460a18dac5bcell_id$9f4d8983-d2d7-4bc9-bb31-5460a18dac5bcode٭details(md"Velocity from redshift",md"By definition $1+z = \sqrt{\frac{c+v}{c-v}},$ where $c$ is the speed of light. Solving for $v$ gives $v = c - \frac{2c}{(1+z)^2+1}$.")metadatashow_logsèdisabled®skip_as_script«code_folded$22da4ae4-6aea-4d45-ab2c-6fe935783e41cell_id$22da4ae4-6aea-4d45-ab2c-6fe935783e41coderesponse_1h = missingmetadatashow_logsèdisabled®skip_as_script«code_folded$a8f89fe8-110b-452e-9ecb-53769aa02a2ecell_id$a8f89fe8-110b-452e-9ecb-53769aa02a2ecode{md"""
#### Overview
In this lab, we'll measure the Hubble constant ($H_0$), the rate of expansion of the local universe using observations of Type 1a Supernovae.
First, we'll download and read data from the Supernovae Cosmology Project and do a few transformations, so that we can plot velocity versus distance.
Second, we'll develop a statistical model to describe our data based on our understanding of the astrophysics and measurement process.
We'll compute the maximum likelihood estimate for the model parameters using conventional linear algebra techniques. This is a *very* useful and common method.
(I'd say its analogous the harmonic oscillator in physics!)
Next, we'll preview multiple other ways we could have fit the same model, but using different techniques.
You'll think about some of the results as we go and the advantages/disadvantages/limitations of each approach.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$a01c7aad-429d-447b-8943-a4e89dbac6f6cell_id$a01c7aad-429d-447b-8943-a4e89dbac6f6code.if ismissing(response_2d)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$cf6fbe84-9a51-4ddb-9def-c22893a17ed3cell_id$cf6fbe84-9a51-4ddb-9def-c22893a17ed3codeif ready_to_run_mcmc
chain_lognormal_prior = sample(model_hubble_law_lognormal_prior_given_data, NUTS(0.65), MCMCThreads(), num_mcmc_samples, num_chains)
else
chain_lognormal_prior = nothing
end;metadatashow_logsèdisabled®skip_as_script«code_folded$098afd51-2c39-42ef-91bf-cfc326827ea1cell_id$098afd51-2c39-42ef-91bf-cfc326827ea1codemd"""**Q1d:** Based on the figure above, what would you estiamte for the best-fit value of $H_0$? How did your previous estimate of the best-fit value for $H_0$ compare to the minimum in the figure above?"""metadatashow_logsèdisabled®skip_as_script«code_folded$d95d2e2d-c6ff-4cab-8b49-432ef61065bccell_id$d95d2e2d-c6ff-4cab-8b49-432ef61065bccode}let
plt = plot_hubble_diagram(df_all,H₀_bf_all, show_yerr=true)
title!(plt,"Best-fit Hubble constant for full dataset")
endmetadatashow_logsèdisabled®skip_as_script«code_folded$a1b5939b-b1ed-4d7c-907d-5bf3e7dd4479cell_id$a1b5939b-b1ed-4d7c-907d-5bf3e7dd4479codelmd"""
Now we can easily compute & plot the $\chi^2$ as a function of the value of the Hubble parameter.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$78eb6aff-3752-41d8-a4d9-0cc2b31fd56dcell_id$78eb6aff-3752-41d8-a4d9-0cc2b31fd56dcode:md"""
## Quantifying Uncertainties in Model Parameters
"""metadatashow_logsèdisabled®skip_as_script«code_folded$ac47f631-787d-4a3c-96cf-5fee86a3b2f4cell_id$ac47f631-787d-4a3c-96cf-5fee86a3b2f4code"H₀_ols_all = 1 /coef(ols_all)[1]metadatashow_logsèdisabled®skip_as_script«code_folded$42366f09-d437-4ba2-93ec-f7b8278eb09dcell_id$42366f09-d437-4ba2-93ec-f7b8278eb09dcodemd"""**Q1i:** Is the difference between the best-fit values for the two datasets larger or smaller than the formal uncertainty estimates? What does that imply for the true uncertainty in the Hubble constant?"""metadatashow_logsèdisabled®skip_as_script«code_folded$d0f003fa-f72d-4a6f-ba60-1cad3844ef6acell_id$d0f003fa-f72d-4a6f-ba60-1cad3844ef6acodemd"""
In the function `model_hubble_law_uniform_prior`, we specified the probability distribution for each random variable using the syntax `variable ~ distribution`.
For example, we specified the prior for $H_0$ and the likelihood, $\mathcal{L} = p(d_{\mathrm{obs}} | d_{\mathrm{true}}, σ_d)$.
We specify how the true distance can be computed from the velocity and Hubble constant, but we .
The input arguements for the model are the data to be analyzed, including $d_\mathrm{obs}$.
But we won't specify a value for $H_0$ or the values of $d_\mathrm{true}$.
The computer will figure out the distributions implied for both of these variables given the data we provide.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$f0220da0-6f3a-438c-b533-acb346382418cell_id$f0220da0-6f3a-438c-b533-acb346382418codedetails(md"Ready for an explanation of the syntax above?",
md"""
Before we can use Julia's CSV reader, first we need to load the [CSV.jl](https://csv.juliadata.org/stable/) package. That's done at the bottom of the notebook with the command `using CSV`. Then we can call the function `CSV.read(...)`.
The first two function arguements specify the filename and that we should read the data into a `DataFrame`.
The remaining arguements are option and allow us to fine-tune the behavior of the CSV reader.
This data file used tabs (written using the escaple code `\t`) to separate the columns rather than the standard commas, so we had to specify that with the optional `delim` parameter.
Additionally, the file included some comments/metadata in the first five rows, so we had to tell the CSV reader to skip to row 6 using the optional `skipto` parameter.
Finally, the file didn't include column heading, so we provided those manually.
""")metadatashow_logsèdisabled®skip_as_script«code_folded$8bfe0413-782b-4244-8624-a9f95e393d17cell_id$8bfe0413-782b-4244-8624-a9f95e393d17codemd"""
Packages like `GLM.jl` allow us to specify a model with a *modeling language*. The compiler figures out how to build the matrices and what calculations to peform in order to compute the MLEs. The simplest possible model formula is `d ~ v` meaning the distance is proportional to velocity (the constant to be fit is implicit). We then call the `lm` function to return the best fit model for a given model and dataset.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$3ec3e49c-04ce-4ffa-b1ca-8c58df0307ddcell_id$3ec3e49c-04ce-4ffa-b1ca-8c58df0307ddcodemd"""
In the figure above, the vertical line indicates the best-fit value for each dataset and the the horizontal bar shows the the measurement uncertainty for each dataset.
We could estimate the uncertainties by calculating the second derivative of the log likelihood.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$92238689-d9df-46c2-bd72-3fcc41a22594cell_id$92238689-d9df-46c2-bd72-3fcc41a22594codemd"""
The output table provides the best-fit value for each of the parameters in the linear model. In this case there is the the coefficient to v ($1/H_0$) as well as an intercept. That's because it fit the model $d = c_v v + c_{\mathrm{intercept}}$, rather than just $d = c_v v$. We can tell it that we don't want to include an intercept using the model formula `dₗ ~ 0 + v`.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$d5c3923b-8515-4002-9d93-6768f087aef6cell_id$d5c3923b-8515-4002-9d93-6768f087aef6code7md"""**Q2c:** Given your findings above, which of the following choices is likely to have the most significant effect on your measurement of the Hubble constant?
- Choice of supernovae to be included in analysis
- Choice of whether to estimate $H_0$ with an MLE or MAP estimator
- Choice of prior for $H_0$
"""metadatashow_logsèdisabled®skip_as_script«code_folded$4f11d0e2-7126-47f1-865a-3474298597dccell_id$4f11d0e2-7126-47f1-865a-3474298597dccodemd"""
When using a PPL like Turing, we need to specify a prior distribution for each model parameter and the likelihood function describing how the observations are generated from the model. For example, we could adopt a Normal prior for $H_0$:
$H_0 \sim \mathrm{Normal}(0,\sigma_{pr,H_0}^2)$
and a Normal likelihood for each observation
$d_{\mathrm{obs},i} \sim \mathrm{Normal}(d_{\mathrm{true},i},\sigma_{d,i}^2) \; \; ∀ \; i ∈ 1... N_{\mathrm{obs}}$
where
$d_{\mathrm{true},i} = v_i / H_0.$
It's often helpful to use vector notation, $\bf{d_{\mathrm{obs}}} \sim \mathrm{MvNormal}(\bf{d_{\mathrm{true}}},\bf{\sigma_{d}}^2).$
We implement this model using the following code
"""
metadatashow_logsèdisabled®skip_as_script«code_folded$337a560b-49f3-4fc5-926f-1689911eacf2cell_id$337a560b-49f3-4fc5-926f-1689911eacf2code٫md"""**Q2d:** Can you think of another improvement we could make to improve the quality of inferences about the expansion rate of the universe? Explain your thinking."""metadatashow_logsèdisabled®skip_as_script«code_folded$7b5b3e8a-61ac-43cf-8b85-4a8a8e4846f6cell_id$7b5b3e8a-61ac-43cf-8b85-4a8a8e4846f6codemd"""
One nice feature of using a PPL like Turing is that they make it easy to repeat analyses using slightly different assumptions. For example, we could switch to using a [LogNormal](https://en.wikipedia.org/wiki/Normal_distribution) prior distribution for $H_0$.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$d81347d3-9246-42d9-9a85-42c773ef2a1dcell_id$d81347d3-9246-42d9-9a85-42c773ef2a1dcode begin
metadata!(df_in,"source", url, style =:note)
colmetadata!(df_in, :z, "label", "redshift", style =:note)
colmetadata!(df_in, :μ, "label", "distance modulus")
colmetadata!(df_in, :σ_μ, "label", "uncertainty in distance modulus")
added_metadata = true
end;metadatashow_logsèdisabled®skip_as_script«code_folded$adf117fe-0f54-4f2a-b107-0635f56263f2cell_id$adf117fe-0f54-4f2a-b107-0635f56263f2coderesponse_2e = missingmetadatashow_logsèdisabled®skip_as_script«code_folded$2ad2563a-8a72-47e8-a36d-f6d2a49cb517cell_id$2ad2563a-8a72-47e8-a36d-f6d2a49cb517code}if chain_uniform_prior!=nothing
(H₀_posterior_mean, H₀_posterior_std) = mean_and_std(vec(chain_uniform_prior[:H₀]))
endmetadatashow_logsèdisabled®skip_as_script«code_folded$c086c05f-ad73-44c7-bf8b-a304d3533690cell_id$c086c05f-ad73-44c7-bf8b-a304d3533690code&"""
`calc_sigma_mle_linear_model(A, covar)`
Computes the uncertainty in the maximum likelihood estimate of θ for the linear model
`y = A θ`
where measurements errors of `y_obs` are normally distributed and have covariance `covar`.
Note that `y_obs` is not passed to this function, because it does not affect results.
"""
function calc_sigma_mle_linear_model(A::AbstractMatrix, covar::AbstractMatrix)
@assert size(A,1) == size(covar,1) == size(covar,2)
@assert size(A,2) >= 1
sqrt.(diag(inv(PDMat(calc_fisher_matrix_linear_model(A,covar)))))
endmetadatashow_logsèdisabled®skip_as_script«code_folded$0a4aaa8a-20da-486f-ace5-fbadcf646609cell_id$0a4aaa8a-20da-486f-ace5-fbadcf646609codeCif chain_uniform_prior!=nothing
describe(chain_uniform_prior)
endmetadatashow_logsèdisabled®skip_as_script«code_folded$eb72dda6-69f0-4538-a75a-25967653759fcell_id$eb72dda6-69f0-4538-a75a-25967653759fcode_md"""
Now we can compare the estimates of the posterior PDF under the two different models.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$2d97268e-560b-4398-982e-1cf5cd9ab7ddcell_id$2d97268e-560b-4398-982e-1cf5cd9ab7ddcode.if ismissing(response_2b)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$28b844a8-650a-47d8-b9f6-586f22d6f603cell_id$28b844a8-650a-47d8-b9f6-586f22d6f603codemd"""
Hmm... This value is significantly different that what we found above. Why?
In the calls above, we used the method of **[Ordinary Least Squares](https://en.wikipedia.org/wiki/Ordinary_least_squares)** which ignores the measurement uncertainties and weights each point equally. We can instead request it use the method of **[Weighted Least Squares](https://en.wikipedia.org/wiki/Weighted_least_squares)** and tell it to use inverse variance weighting as shown below.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$c78f0ad9-03c0-4066-a390-193e7bfa84b1cell_id$c78f0ad9-03c0-4066-a390-193e7bfa84b1codeٽmd"""**Q2e:** How could using a probabilistic programming language like Turing make it relatively easy to improve the quality of our inferences about the expansion rate of our universe?"""metadatashow_logsèdisabled®skip_as_script«code_folded$1c35257a-9b98-4af6-bde0-a945665c91cdcell_id$1c35257a-9b98-4af6-bde0-a945665c91cdcodebdetails(md"""Propagation of uncertainties""",
md"""The data file provide measurement uncertainties for the distance modulus. We'd like measurement uncertainties for the distance. Unders several assumptions (e.g., small magnitude of uncertainties, uncorrelated uncertainties), we can estimate the measurement uncertainty in the distance by applying the [Propagation of uncertainties](https://en.wikipedia.org/wiki/Propagation_of_uncertainty) method. If you haven't seen this before, then I'd suggest that you just trust me on this for now rather than worrying about the derivation of the equation used.
""")metadatashow_logsèdisabled®skip_as_script«code_folded$751eff6b-2ed2-4073-b8b5-02ad9981edb6cell_id$751eff6b-2ed2-4073-b8b5-02ad9981edb6coderesponse_1i = missingmetadatashow_logsèdisabled®skip_as_script«code_folded$5974dfc3-8416-4bae-91de-d3f0a4f0c7dbcell_id$5974dfc3-8416-4bae-91de-d3f0a4f0c7dbcode.if ismissing(response_1a)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$20f739ff-cc5b-460e-851a-141003d0c9d0cell_id$20f739ff-cc5b-460e-851a-141003d0c9d0codeSmd"""
Now the value returned should match our previous analysis very closely.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$f76f4c3c-74f1-44b0-9b3a-5895aec22e70cell_id$f76f4c3c-74f1-44b0-9b3a-5895aec22e70code!H₀_bf_all = 1/H₀inv_bf_all[1]metadatashow_logsèdisabled®skip_as_script«code_folded$a98ace81-9afd-4618-88a6-fadaf04eb364cell_id$a98ace81-9afd-4618-88a6-fadaf04eb364code.if ismissing(response_1d)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$bd1c72d5-1d0d-4e36-aba4-95c9c9e4fab1cell_id$bd1c72d5-1d0d-4e36-aba4-95c9c9e4fab1coderesponse_1e = missingmetadatashow_logsèdisabled®skip_as_script«code_folded$95bc4b33-d619-4a95-af19-5771d209336dcell_id$95bc4b33-d619-4a95-af19-5771d209336dcodek# hideall
begin
title = "Astro 416: Lab 3, Ex 1"
week = 3
topic = "Model Building I: Linear Models"
end;metadatashow_logsèdisabled®skip_as_script«code_folded$4bdd7b2d-9a67-4161-a0a0-93073d9edcf0cell_id$4bdd7b2d-9a67-4161-a0a0-93073d9edcf0codeٗmd"""
### Plot the data
It's often a good idea to make a quick plot to help validate our dataset and to understand the characteristics of the data.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$2532e9de-8fa7-44ce-8502-c2473928aa8ecell_id$2532e9de-8fa7-44ce-8502-c2473928aa8ecodeWidthOverDocs()metadatashow_logsèdisabled®skip_as_script«code_folded$4a2c109d-8e18-4561-913b-769f7eaeb6e9cell_id$4a2c109d-8e18-4561-913b-769f7eaeb6e9codemd"""
# Develop a model
In this case, it may be tempting to rush to fitting a line to the data.
But we want to use this example as a chance to think through the process of building models.
It's often useful to divide the full model for astronomical observations into two parts: a physical model and a measurement model.
### Physical model
In this case, our initial physical model is that the local universe is expanding away from us at a velocity ($v$) proportion to the distance ($d$).
The rate of proportionality is known as the Hubble constant ($H_0$).
If this model were perfectly accurate and there were no measurement errors, then we'd have
$$d = H_0^{-1} v.$$
### Measurement model
Astronomical measurements almost always have some measurement errors.
It can be useful to be explicit about when we're referring to the true value and when we're referring to the observed value.
For example, if we were measuring $v$, then we might write:
$d_{\rm obs} = d_{\rm true} + \epsilon,$
where $\epsilon$ is the measurement error.
By definition the true value of the measurement error ($\epsilon$) is almost always unknown.
Typically, astronomers design the measurements such that the expected value of the measurement is equal the true value ($E[d_{\rm obs}] = d_{\rm true}$),
and characterize the distribution of the measurement errors.
One very common model/approximation is that the measurement errors follow a [Normal (or Gaussian) distribution](https://en.wikipedia.org/wiki/Normal_distribution) with a scale parameter, $\sigma$.
We denote that with
$\epsilon \sim \mathrm{Normal}(0,\sigma^2).$
Based on the properties of the Normal distribution,
$E[\epsilon]=0$, and
$E[\epsilon^2] = \sigma^2$.
For this lab, we'll assume that's a good approximation here.
When we have many measurements, we could write the errors explicitly like
$d_{i,\rm obs} = d_{i,\rm true} + \epsilon_i = v_i/H_0 + \epsilon_i,$
where $\epsilon_i \sim \mathrm{Normal}(0,\sigma_i^2)$ for each obsevation $i\in 1...N_{\mathrm obs}$
"""metadatashow_logsèdisabled®skip_as_script«code_folded$2f3d492c-4a5e-4a56-a773-e8078e14476acell_id$2f3d492c-4a5e-4a56-a773-e8078e14476acodeBH₀inv_bf_all = calc_mle_linear_model(A_all, df_all.dₗ, Σ_all)metadatashow_logsèdisabled®skip_as_script«code_folded$c43f2686-da2d-4b67-b87e-7c22bd09257ccell_id$c43f2686-da2d-4b67-b87e-7c22bd09257ccoderesponse_1g = missingmetadatashow_logsèdisabled®skip_as_script«code_folded$b54f8319-dce6-4009-873e-bd47dc2a4619cell_id$b54f8319-dce6-4009-873e-bd47dc2a4619codeٚlet
plt = plot_hubble_diagram(df_fit,H₀_wls_fit, show_yerr=true)
title!(plt,"Fitting only SNs with v < " * string(floor(Int,max_v_plt)) * " km/s")
endmetadatashow_logsèdisabled®skip_as_script«code_folded$ddf74f60-42eb-4eab-b77b-b6e78f0e90a7cell_id$ddf74f60-42eb-4eab-b77b-b6e78f0e90a7codeEmd"""**Q1g:** How precise is our estimate of the Hubble constant?"""metadatashow_logsèdisabled®skip_as_script«code_folded$a6cf5e52-8692-4e3c-95aa-11167b118bd4cell_id$a6cf5e52-8692-4e3c-95aa-11167b118bd4code0md"""
Evaluating the posterior PDF direclty is hard (often due to the noramlization constant).
To get around the challenge, Turing allows us to generate samples from the posterior PDF.
There are many interesting algorithms for sampling from a complex posterior distribution.
Here, we'll adopt the ["NUTS" sampler](https://www.mcmchandbook.net/HandbookChapter5.pdf) which based on [Hamiltonian Monte Carlo](https://mc-stan.org/docs/2_19/reference-manual/hamiltonian-monte-carlo.html).
The algorithm is quite complex, so we won't go into the details.
But it's good to konw that this is often a good choice for models with a few to dozens of model parameters.
We'll specify that we want it to draw `num_mcmc_samples` samples and repeat the calculation `num_chains` different times (running in parallel).
"""metadatashow_logsèdisabled®skip_as_script«code_folded$0f6aa1eb-2f89-400f-b1e8-bfc07012179dcell_id$0f6aa1eb-2f89-400f-b1e8-bfc07012179dcode0md"""
## Visualize $\chi^2$ or Loss Function
"""metadatashow_logsèdisabled®skip_as_script«code_folded$0991f15c-b734-480b-8b12-6d2beda547dbcell_id$0991f15c-b734-480b-8b12-6d2beda547dbcodekbegin
# Make a new DataFrame with just data selected for new fit.
df_fit = filter(r->r.v<= max_v_plt, df_all)
# Note that the value of max_v_plt above comes from the slider below
# Fit linear model
wls_fit = lm(formula_wo_intercept, df_fit, wts=weights(df_fit.σ_dₗ.^(-2)))
# Store best-fit value of H_0 for new fit
H₀_wls_fit = 1/coef(wls_fit)[1]
end;metadatashow_logsèdisabled®skip_as_script«code_folded$5c9b0395-4d2e-4679-8d86-b9581fa0f66ccell_id$5c9b0395-4d2e-4679-8d86-b9581fa0f66ccodebmd"""
# Fitting a Model to Data
## Computing Maximum Likelihood Estimator for model parameters
"""metadatashow_logsèdisabled®skip_as_script«code_folded$340dc0ff-dd1c-44e0-907c-f357adb122d5cell_id$340dc0ff-dd1c-44e0-907c-f357adb122d5codeqmodel_hubble_law_uniform_prior_given_data = model_hubble_law_uniform_prior(df_fit.v, df_fit.dₗ, df_fit.σ_dₗ)metadatashow_logsèdisabled®skip_as_script«code_folded$0c25a4cc-9f53-4729-985b-de692278c433cell_id$0c25a4cc-9f53-4729-985b-de692278c433coderesponse_1d = missingmetadatashow_logsèdisabled®skip_as_script«code_folded$ad758d43-b135-418e-a29c-36cea3a24457cell_id$ad758d43-b135-418e-a29c-36cea3a24457code aside(tip(md"""
`MvNormal` specifies a multivariate normal distribution with one output for each observation.
The expected value for each observation is $d_\mathrm{true}$ and the standard deviation of each measurement is given by $\sigma_{d,\mathrm{obs}}$.
In principle, we could pass a covariance matrix for the second argument to describe correlated measurement uncertainties.
Instead, we've provided a vector for the second arguement, asking for a model that assume uncorrelated measurement errors."""))metadatashow_logsèdisabled®skip_as_script«code_folded$18edb530-8472-4932-946c-917d8ab9dbdecell_id$18edb530-8472-4932-946c-917d8ab9dbdecodeFif chain_lognormal_prior!= nothing
describe(chain_lognormal_prior)
endmetadatashow_logsèdisabled®skip_as_script«code_folded$2b3729fc-62b2-478d-b35c-e69e6b0efc93cell_id$2b3729fc-62b2-478d-b35c-e69e6b0efc93code"""
`plot_hubble_diagram(df, H₀; show_yerr, show_resid)`
Plots distance versus velocity.
Adds a line based on linear model of Hubble expansion with rate `H₀`.
Optionally, plots residuals to model.
"""
function plot_hubble_diagram(df::AbstractDataFrame, H₀::Real;
show_yerr::Bool=false, show_resid::Bool=false)
@assert "v" ∈ names(df)
@assert "dₗ" ∈ names(df)
@assert "σ_dₗ" ∈ names(df)
@assert size(df,1) >= 2
x_pred = range(0,stop=maximum(df.v), length=10)
y_pred = predict_linear_model(reshape(x_pred,length(x_pred),1), [1/H₀])
plt = plot(x_pred,y_pred, label="H₀ = " * string(round(H₀,digits=2)) * " km/s/Mpc")
if show_yerr
scatter!(plt,df.v,df.dₗ, yerr=df.σ_dₗ, label=:none)
else
scatter!(plt,df.v,df.dₗ, label=:none)
end
xlabel!(plt,"v (km/s)")
ylabel!(plt,"D (Mpc)")
plt
if show_resid
y_pred = predict_linear_model(reshape(df.v,size(df,1),1), [1/H₀])
plt2 = plot()
if show_yerr
scatter!(plt2,df.v, df.dₗ.-y_pred, yerr=df.σ_dₗ, label=:none)
else
scatter!(plt2,df.v, df.dₗ.-y_pred, label=:none)
end
xlabel!(plt2,"v (km/s)")
ylabel!(plt2,"Residual D (Mpc)")
return plot(plt,plt2,layout=(2,1))
else
return plt
end
endmetadatashow_logsèdisabled®skip_as_script«code_folded$92af3764-6c53-4e03-a338-81b1de202a04cell_id$92af3764-6c53-4e03-a338-81b1de202a04code.md"""
### Computing MLE via Linear Algebra
"""metadatashow_logsèdisabled®skip_as_script«code_folded$0a998eb0-4bf5-4d4e-b193-54964b54f1abcell_id$0a998eb0-4bf5-4d4e-b193-54964b54f1abcode1md"""
One disadvantage of using simple text files to share data is that they don't have a uniform way to store **metadata** about the file. For example, it would be nice to include a longer description of the of each column. So after we've read in the datafile, we can add that metadata ourselves.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$8ab41174-d7d9-42c7-8129-381c54906280cell_id$8ab41174-d7d9-42c7-8129-381c54906280code md"""
Thinking a little more broadly, sometimes the goodness-of-fit statistic that we want to minimize is more complicated. So often people refer to minimizing an **objective function** or a **loss function**.
We'll implement a $\chi^2$ loss function below.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$7f61c098-1c54-4c82-be40-f9fd154156a1cell_id$7f61c098-1c54-4c82-be40-f9fd154156a1codedetails("Choice of Priors",
md"""
The choice of prior probability distributions is often subjective.
For example, one astronomer might say that they want to assume as little as practical about the Hubble constant.
This has some philosophical appeal, but it can be tricky to put into practice.
For example, should we assume that the Hubble constant is positive?
Should we assume that there's some maximum velocity (maybe the speed of light)?
Should we assume that very large Hubble constants are less likely that small values? If so, what counts as "large"?
Another astronomer might take a different approach and want to approximation to our prior knowledge about the value of the Hubble constant based on previous studies (critically, **excluding** the data that is currently being analyzing).
This also comes with follow-up questions.
Should we approximate our prior knowledge with a simple, named distribution?
Or should we use a more complex prior to better describe the conclusions of previous research?
There are almost always some very bad choices for priors.
But often there's not a single best choice.
Rather than trying to find "the correct" prior, it's often wise to perform multiple analyses with multiple reasonable choices for priors.
You may find the the difference in outcome are so small that they don't affect your conclusions.
Or you might discover than even small differences in the choice of the priors can have a substantive impact on the conclusions.
Knowing which regime you're in is often more a valuable result than computing the result under one single choice of priors.
Here, if we assumed a broad uniform prior for $H_0$, then the maximum a posteriori estimate (MAP) would match the MLE.
Altneratively, we might adopt as [Log-Normal distribution](https://en.wikipedia.org/wiki/Normal_distribution), since it has non-zero probability for all positive values and we're able to pick a location for the peak and scale that are qualitatively similar to what we had already learned from previous studies.
""")metadatashow_logsèdisabled®skip_as_script«code_folded$ef78506f-0326-4bd0-8fb3-af6f34a06841cell_id$ef78506f-0326-4bd0-8fb3-af6f34a06841codeaside(tip(md"To type unicode characters like σ, first type `\`, then the name of the character your want (e.g., 'sigma'), then press tab. For the subscript `ₜ`, you'd type `\_t`+`TAB`."), v_offset=-180)metadatashow_logsèdisabled®skip_as_script«code_folded$316f65a1-9062-4cb0-a3fb-af54149afb4dcell_id$316f65a1-9062-4cb0-a3fb-af54149afb4dcodeSplot_hubble_diagram(df_all,H₀_plt_1; show_yerr = plt_yerr, show_resid=plt_resid )metadatashow_logsèdisabled®skip_as_script«code_folded$8ddbc921-4541-4c45-90c5-5c51ec2e391acell_id$8ddbc921-4541-4c45-90c5-5c51ec2e391acode,ols_all = lm(@formula(dₗ ~ 0 + v), df_all)metadatashow_logsèdisabled®skip_as_script«code_folded$d67f0b9e-10be-4511-93ab-edceb3aef49ccell_id$d67f0b9e-10be-4511-93ab-edceb3aef49ccode%H₀_wls_all = 1.0 / coef(wls_all)[1]metadatashow_logsèdisabled®skip_as_script«code_folded$1eb1face-a795-4cb3-9c51-f995f73e2509cell_id$1eb1face-a795-4cb3-9c51-f995f73e2509code:if chain_uniform_prior!=nothing
let
bins = range(floor(minimum(chain_uniform_prior[:H₀])),
stop=ceil(maximum(chain_uniform_prior[:H₀])), step=0.1)
plt = histogram_or_density(chain_uniform_prior, bins = bins, alpha=0.3, lw=3)
xlabel!(plt, "H₀")
title!(plt,"Posterior Samples (Uniform Prior)")
end
endmetadatashow_logsèdisabled®skip_as_script«code_folded$ccdf5418-12e4-49ee-bdea-90efc202583acell_id$ccdf5418-12e4-49ee-bdea-90efc202583acodelet
if !ismissing(response_1f)
ϵ = 1e-4
d2lossdH₀2_fit = (loss(H₀_fit+ϵ,df_fit)+loss(H₀_fit-ϵ,df_fit)-2*loss(H₀_fit,df_fit))/(ϵ^2)
σ_H₀_fit_deriv_fit = sqrt(2/(d2lossdH₀2_fit))
end
endmetadatashow_logsèdisabled®skip_as_script«code_folded$2efb7fb2-0402-4eec-b58a-a52fd530e3f5cell_id$2efb7fb2-0402-4eec-b58a-a52fd530e3f5codemd"""
## Read the Data file
"""metadatashow_logsèdisabled®skip_as_script«code_folded$3a06141a-a92d-4c31-9950-13f0bf2c023acell_id$3a06141a-a92d-4c31-9950-13f0bf2c023acode.if ismissing(response_1g)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$bf7ca6ea-4c07-46e4-baff-b205e1663505cell_id$bf7ca6ea-4c07-46e4-baff-b205e1663505codemd"""
#### $title
# $topic
"""metadatashow_logsèdisabled®skip_as_script«code_folded$2656d886-60dc-44a7-b318-100335dcf05bcell_id$2656d886-60dc-44a7-b318-100335dcf05bcode.formula_wo_intercept = @formula(dₗ ~ 0 + v);metadatashow_logsèdisabled®skip_as_script«code_folded$2c1c7d8a-eda6-4900-9be7-e4b908fed81acell_id$2c1c7d8a-eda6-4900-9be7-e4b908fed81acode4begin
num_mcmc_samples = 2_000
num_chains = 3
end;metadatashow_logsèdisabled®skip_as_script«code_folded$a1eb8ff3-ff5b-4e50-a75c-480cf2a7d861cell_id$a1eb8ff3-ff5b-4e50-a75c-480cf2a7d861code"""
`predict_linear_model(A, θ)`
Computes the predictions of a linear model with design matrix `A` and parameters `θ`.
"""
function predict_linear_model(A::AbstractMatrix, b::AbstractVector)
@assert size(A,2) == length(b)
A*b
endmetadatashow_logsèdisabled®skip_as_script«code_folded$9ef67704-5063-4f75-853b-078261e50684cell_id$9ef67704-5063-4f75-853b-078261e50684code١md"""
Now we'll call a function (implemented at the bottom of the notebook) to compute the best-fit model parameters for a linear model using linear algebra.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$095bdd00-ffa9-434d-9fc5-0978019b060bcell_id$095bdd00-ffa9-434d-9fc5-0978019b060bcodeebegin
# Scratch cell for any optional calculations to check out concerns to be described in Q1b.
endmetadatashow_logsèdisabled®skip_as_script«code_folded$3f2ef930-d367-11ef-3b5c-95b4989c2e21cell_id$3f2ef930-d367-11ef-3b5c-95b4989c2e21codeYbegin
mkpath("data")
filename = joinpath("data","SCPUnion2.1_mu_vs_z.txt")
if split(pwd(),'/')[end] == "test"
filename = "../" * filename
end
url = "http://supernova.lbl.gov/Union/figures/SCPUnion2.1_mu_vs_z.txt"
if !(filesize(filename)>0)
Downloads.download(url,filename)
@assert filesize(filename)>0
end
endmetadatashow_logsèdisabled®skip_as_script«code_folded$299e9d1f-bbb2-4fc3-8682-51c3daab0ba3cell_id$299e9d1f-bbb2-4fc3-8682-51c3daab0ba3coderesponse_1c = missingmetadatashow_logsèdisabled®skip_as_script«code_folded$ac1a64e8-7a5c-43b5-923a-633dac4fb664cell_id$ac1a64e8-7a5c-43b5-923a-633dac4fb664codemd"""
A vector containing each distance measurement,
$y_{\mathrm{obs}} =
\left(
\begin{matrix}
d_{\mathrm{obs},1} \\
d_{\mathrm{obs},1} \\
\ddots \\
d_{\mathrm{obs},N_{\mathrm{obs}}} \\
\end{matrix}
\right),$
"""metadatashow_logsèdisabled®skip_as_script«code_folded$d9058d92-b1c9-44f2-af5c-3540932dec48cell_id$d9058d92-b1c9-44f2-af5c-3540932dec48code,A_all = reshape(df_all.v, size(df_all,1), 1)metadatashow_logsèdisabled®skip_as_script«code_folded$0e42628b-95a9-45a0-b375-51f62dc0ce0dcell_id$0e42628b-95a9-45a0-b375-51f62dc0ce0dcode.if ismissing(response_1c)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$e09e30f7-838f-4571-bb9c-6d0b425807a8cell_id$e09e30f7-838f-4571-bb9c-6d0b425807a8codeX# hideall
md"""
+++
title = $title
week_num = $week
topic = $topic
+++
""" |> Base.Text;metadatashow_logsèdisabled®skip_as_script«code_folded$a2b56e1c-9735-49b0-91da-40f489823390cell_id$a2b56e1c-9735-49b0-91da-40f489823390codelaside(tip(md"""We could have created the matrix with the `diagm` function. Instead, we used `PDiagMat` to create a matrix that the computer knows will always be diagonal and positive definite. This allows it to use algorithms that are faster and/or more robust than if it couldn't assume anything about the propeties of the covariance matrix."""),v_offset=-280)metadatashow_logsèdisabled®skip_as_script«code_folded$4dcd9f29-8dd5-40ff-9ce9-50bd38f33348cell_id$4dcd9f29-8dd5-40ff-9ce9-50bd38f33348codemd"""**Q1f:** What velocity threshold did you pick? What is the results best-fit Hubble constant? How does this compare to the best-fit Hubble constant when fitting to the full data set?"""metadatashow_logsèdisabled®skip_as_script«code_folded$e4ca4972-42c8-4426-94ae-e316662513c3cell_id$e4ca4972-42c8-4426-94ae-e316662513c3codeBmd"""
### Computing MLE with package for General Linear Models
"""metadatashow_logsèdisabled®skip_as_script«code_folded$6a75dd68-5034-4b16-89a0-7388cfae0cbdcell_id$6a75dd68-5034-4b16-89a0-7388cfae0cbdcodeGmd"""
We'll package the velocities and uncertainties into matrices.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$d68dd4b4-c091-4435-a286-3de589b19009cell_id$d68dd4b4-c091-4435-a286-3de589b19009codemd"""
For a linear model, we can estimate the measurement uncertainty associated with our maximum likelihood estimators.
It's most intuitive to plot the $\Delta\chi^2$ as a function of $H_0$.
By construction, the slope of $\chi^2(H_0)$ is zero at the best-fit value.
The magnitude of the measurement uncertainty is inversely proportional to the curvature (absolute value of the second derivative) of $\chi^2(H_0)$.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$33045766-a7c5-42ff-bdf8-bde1fce9c56bcell_id$33045766-a7c5-42ff-bdf8-bde1fce9c56bcode.if ismissing(response_2a)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$9555ae09-be3d-48e8-9ee2-face2bef8ccecell_id$9555ae09-be3d-48e8-9ee2-face2bef8ccecode5md"""
## Fitting model to a subset of the dataset
"""metadatashow_logsèdisabled®skip_as_script«code_folded$466dacda-1835-42fd-b8f4-cac5c22f6756cell_id$466dacda-1835-42fd-b8f4-cac5c22f6756code)ols_try1 = lm(@formula(dₗ ~ v), df_all)metadatashow_logsèdisabled®skip_as_script«code_folded$915fc532-7d8d-4ff5-81fd-af8f6e9e64dacell_id$915fc532-7d8d-4ff5-81fd-af8f6e9e64dacode7if chain_lognormal_prior!=nothing
let
bins = range(floor(minimum(vcat(chain_uniform_prior[:H₀],chain_lognormal_prior[:H₀]))),
stop=ceil(maximum(vcat(chain_uniform_prior[:H₀],chain_lognormal_prior[:H₀]))), step=0.1)
plt1 = histogram_or_density(chain_lognormal_prior, bins=bins, alpha=0.3, lw=3)
xlabel!(plt1, "")
title!(plt1,"Posterior Samples (Uniform Prior)")
plt2 = histogram_or_density(chain_uniform_prior, bins=bins, alpha=0.3,lw=3)
xlabel!(plt2, "H₀")
title!(plt2,"Posterior Samples (LogNormal Prior)")
plot(plt1,plt2,layout=(2,1))
end
endmetadatashow_logsèdisabled®skip_as_script«code_folded$fb2545e1-d624-4669-9e07-c94a598ba8edcell_id$fb2545e1-d624-4669-9e07-c94a598ba8edcode.if ismissing(response_2c)
still_missing()
endmetadatashow_logsèdisabled®skip_as_script«code_folded$67e9aa94-2d3a-4888-97b7-c2cb019eb04acell_id$67e9aa94-2d3a-4888-97b7-c2cb019eb04acodelet
if !ismissing(response_1f)
ϵ = 1e-4
d2lossdH₀2_all = (loss(H₀_bf_all+ϵ,df_all)+loss(H₀_bf_all-ϵ,df_all)-2*loss(H₀_bf_all,df_all))/(ϵ^2)
σ_H₀_fit_deriv_all = sqrt(2/(d2lossdH₀2_all))
end
endmetadatashow_logsèdisabled®skip_as_script«code_folded$389e8239-6f62-46fb-8b13-d336516e5e2fcell_id$389e8239-6f62-46fb-8b13-d336516e5e2fcode٬let
plt = plot(LogNormal(log(100),1.0), xlims=(0,250), label=:none)
xlabel!(plt,"H_0")
ylabel!(plt,"p(H₀)")
title!(plt,"Log Normal prior (median 100, scale 1.0)")
endmetadatashow_logsèdisabled®skip_as_script«code_folded$54b51239-59b0-4767-9bc8-018d77ee5e44cell_id$54b51239-59b0-4767-9bc8-018d77ee5e44codemd"""
We can plot histograms for each of the calculations.
We can compute an estimtae of the posterior density ($p(H_0|\bf{v},\bf{d},\bf{\sigma_d})$) using the samples from each chain.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$d62cf934-a592-482b-b0b3-7fd55f9e0617cell_id$d62cf934-a592-482b-b0b3-7fd55f9e0617code~md"""
Or we can use an analytic expression for the uncertainty for parameters estimated from a general linear model,
$$E\left[\mathrm{Covar}(\theta)\right] = (A' \Sigma^{-1} A)^{-1}$$
and propagate the uncertainties for $H_0$.
We'll skip the tedious derivation, but you can see the implementation in the function `calc_sigma_mle_linear_model` at the bottom of this notebook.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$aab05339-dea2-4c7f-8b38-bc7d7971f1d9cell_id$aab05339-dea2-4c7f-8b38-bc7d7971f1d9codemd"""
## Functions used
"""metadatashow_logsèdisabled®skip_as_script«code_folded$46eb1088-0265-46a6-91ec-44dfa3bd862dcell_id$46eb1088-0265-46a6-91ec-44dfa3bd862dcodewscatter(df_all.v, df_all.dₗ, yerr=df_all.σ_dₗ, label=:none, xlabel="v (km/s)", ylabel="Luminosity Distance (Mpc)")metadatashow_logsèdisabled®skip_as_script«code_folded$25ef65bd-8b86-4f88-b0b7-32f05ce82adacell_id$25ef65bd-8b86-4f88-b0b7-32f05ce82adacoderesponse_1a = missingmetadatashow_logsèdisabled®skip_as_script«code_folded«notebook_id$48878fcc-de6d-11ef-223b-a1023bf62c62in_temp_dir¨metadata