Ț „bondsŹcell_resultsȚ eÙ$71829343-22e9-4f5f-8c54-afd4dffa826aŠqueued€logs§runningŠoutput€bodyÚEWe see that the datetime2julian function takes a single input parameter of type DateTime object (provided by the Dates module) and returns a variable of type Float64. How do we make a DateTime object? Pluto offers a convenient way to view the documentation for a function (or type, module, etc.). Click on the DateTime in the cell below and then open the Live Docs panel (probably in the bottom right of your browser).
Visualizing the Results
Now, we'd like to compare the results of the two algorithms. It will be helpful to visualize the difference as a function of the number of samples. Therefore, we'll make a function to generate a random data set with N samples and a specified true_mean for the distribution the samples are drawn from. Here true_mean is an optional, named arguement that defaults to zero.
Missing Response
The variable response_1d is still set to missing.
Under what circumstance would it be a good/poor choice to use?
Keep working on it!
The answer is not quite right.
Next, I'd like to compute the corresponding Julian date . The Dates module provides a function, datetime2julian to do that for us. Let's check how to call that function.
c. How large are the differences? Are they significant relative to the true values? Why is the difference for one quantity a larger fraction of its true value than the other?
Calling Functions
A key principle of writing code for non-trivial tasks is to organize one's code into many small functions, each of which do one thing (hopefully well). High-level languages typically come with numerous functions that allow developers to accomplish common tasks without reinventing the wheel. For example, the function sqrt(x) computes the square root of x.
Missing Response
The variable response_2h is still set to missing.
If you suceeded above, then Pluto will soon display a plot showing the absolute value of the difference between the two variance estimates below as a function of the number of observation dates in the sample. First, make a prediction for what you expect such a plot to look like.
e. What lessons does this exercise illustrate that could be important when writing similar code for your research?
Helper Code
€mime©text/htmlŹrootassigneeÀČlast_run_timestampËAÚ,YêP°persist_js_state·has_pluto_hook_features§cell_idÙ$b760fedd-41ea-4784-845f-ede0163c0d12čdepends_on_disabled_cells§runtimeÎ `”published_object_keysždepends_on_skipped_cells§erroredÂÙ$2ca37a14-1e12-41ab-b85a-b8d6e9a28ab6Šqueued€logs§runningŠoutput€bodyČ0.9615827199241083€mimeȘtext/plainŹrootassignee§s_64bitČlast_run_timestampËAÚ,+ÜÆ°persist_js_state·has_pluto_hook_features§cell_idÙ$2ca37a14-1e12-41ab-b85a-b8d6e9a28ab6čdepends_on_disabled_cells§runtimeÍE”published_object_keysždepends_on_skipped_cells§erroredÂÙ$af2bd92f-a67f-4fe3-965e-d337d57a2368Šqueued€logs§runningŠoutput€body§missing€mimeȘtext/plainŹrootassignee«response_2gČlast_run_timestampËAÚ,á_é°persist_js_state·has_pluto_hook_features§cell_idÙ$af2bd92f-a67f-4fe3-965e-d337d57a2368čdepends_on_disabled_cells§runtimeÍHh”published_object_keysždepends_on_skipped_cells§erroredÂÙ$f0c73fc1-8da9-4579-a369-a3c907fc56f4Šqueued€logs§runningŠoutput€bodyÙĐe. Consider the online 1-pass algorithm below for calculating the sample variance given below and then compare its results to the other algorithms for different data sets.
Computing Variances
Next, you will compute the variance of the above data using multiple algorithms and compare their relative merits. Algebraically, the sample mean is calculated via $m = 1/N \times \sum_{i=1}^{N} y_i$ and the sample variance can be written two ways $s^2 = 1/(N-1) \times \sum_{i=1}^N (y_i-m)^2$ or $s^2 = 1/(N-1) \times \left[ \left( \sum_{i=1}^N y_i^2 \right) - N m^2 \right] = 1/(N-1) \times \left[ \left( \sum_{i=1}^N y_i^2 \right) - \left(\sum_{i=1}^N y_i\right)^2 /N \right]$ .
In this exercise, you will consider how to calculate the sample variance accurately and efficiently. First, you'll try writing a function yourself. To get help with syntax, you can hover your mouse over the following tip boxes below. The example in the first hint box demonstrates how to write a function with a for loop and how to access elements of an array in Julia. The second hint box demonstrates using a two function calls.
datetime2julian(dt::DateTime) -> Float64Take the given DateTime and return the number of Julian calendar days since the julian epoch -4713-11-24T12:00:00 as a Float64
First, let's review the code in the cell above. The first line is a "docstring", it describes what the function below does both for developers reading the code and for users who might get the same information from Pluto's LiveDocs featuore or a website with documentation automatically extracted from the package's docstrings (using the Documenter.jl package and a GitHub Action).
The rest of the cell defines a function that takes two input parameters and returns a 1-d array of random variables. The first parameter (N) is required and must be some form of an integer. The second parameter (true_mean) could have any type and has a default value of zero. The third parameter (after the ;) is a named parameter (i.e., you must specify the name of the parameter when calling the function, instead of just using its position). Because it has a default value it is optional.
Each time the function is called, it will begin by seeding a pseudo-random number generator. This is important so that results will be reproducible when run multiple times. The function randn returns a 1-d array of standard random variables (i.e., drawn from a normal distribution with zero mean and unit variance) drawn using Julia's default pseudo-random number generator. Then the function returns the variable sample.
When you execute the code block above, julia parses the function, but does not compile or execute the code. That will only happen once the function is called. Since the last line of the cell is the end of the function, the output of the cell is the function. By ending the line with a ;, we tell Julia not to display the output. Now let's try out using this function.
Hint
"Calculate mean value of an array using a simple for loop."
function mean_demo_verbose(y::Array) # the syntax ::Array specifies that this function can only be applied if argument is an array.
n = length(y) # get the number of elements in the array y
sum = zero(first(y)) # set sum to zero. Using zero(first(y)) makes sum have the same data type as the first element of y
for i in 1:n # In Julia and Fortran, arrays start a 1, not 0 (like in C arrays and Python lists)
sum += y[i] # Short-hand for sum += sum + y[i]
end
return sum/n # return isn't necessary since functions return the last value by default
end
var_one_pass([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) â 11
var_one_pass( ) â 11
â 11
€mime©text/htmlŹrootassigneeÀČlast_run_timestampËAÚ,C°persist_js_state·has_pluto_hook_features§cell_idÙ$a83d07b9-d7b7-4274-929a-9a3474e44f08čdepends_on_disabled_cells§runtimeÎ9
<ö”published_object_keysÙ5ffb28f94-8543-11f0-120d-b33132eb5b69/7ccb0078dcd6bfa1Ù5ffb28f94-8543-11f0-120d-b33132eb5b69/b8855b077021025aždepends_on_skipped_cells§erroredÂÙ$b4d6143c-42ad-460d-8af3-a36dae1a8879Šqueued€logs§runningŠoutput€bodyÚ6Broadcasting
Look more closely at the function generate_sample above. Note the syntax .+ that tells julia the programmer wants to "broadcast" the scalar true_mean to have the same dimensions as the result of randn(...). What do you think would happen if you replaced this with true_mean+randn(N)? Try it. How does the behavior compare to what you expected?
Restore the code in generate_sample and execute the cell again, so the rest of the lab works as intended.
b. What is the advantage of julia having different syntax for arithmetic on variables with compatible dimensions from arithmetic on variables with different dimensions?
Curious to learn more?
Pro Tip
Normally, we'd use the Test module for the @test macro. Julia has a large set of modules and packages, that range from very basic functionality to complex science codes. The quality also varries widely. Several modules (like Test) are included in Julia standard library, so they're already installed for us.
However, inside Pluto, it can be helpful to instead import PlutoTest, since it displays the results particularly nicely. (It's an external package and it's still experimental, so if things break in the future, then we can revert to just using Test.
Below, I pick one based on whether we are inside a Pluto notebook session.
€mime©text/htmlŹrootassigneeÀČlast_run_timestampËAÚ,ÀÏ°persist_js_state·has_pluto_hook_features§cell_idÙ$58d4e72a-aed5-4b8f-846e-bb63e4cc54c7čdepends_on_disabled_cells§runtimeÎe[”published_object_keysždepends_on_skipped_cells§erroredÂÙ$65ef6abf-460d-49ac-80e1-1ecc1d0250ddŠqueued€logs§runningŠoutput€bodyÚùPluto: A Reactive Notebook Experience
Some of you may have experience using Jupyter notebooks. Indeed, Jupyter notebooks are a useful and commonly used for small Astronomy and Data Science projects. One big disadvantage of Jupyter notebooks is that the notebook doesn't provide a complete description of the kernel state. That's a fancy way of saying that you can run cells out of order, or change a cell and not recalculate something that depended on the results of that cell. It's suprisingly easy to confuse yourself. Indeed, the first time Astro 528 was offered, we used Jupyter notebooks for nearly all the assignments. When students encountered trouble, the most common advice they got was "Restart your notebook and step through the notebook, one cell at a time until you find where it breaks." In contrast, Pluto keeps track of all dependancies across cells. When you update a cell, it recalculates all the cells that depend on it!
Pluto can also be useful for making interactive visualizations. In the example below, we'll make a plot that depends on a variable true_mean_plt defined below. When you change the value of true_mean_plt, the plots below should automatically update itself. Try setting it to a value of 10 or 100 times larger or smaller and observed how the difference in the estimates of the sample standard deviation change.
We see that there are several different constructors to construct a DateTime object. We'll pick one below.
Missing Response
The variable response_1e is still set to missing.
Missing Response
The variable response_2e is still set to missing.
Example Function
Consider an astronomer analyzing data from a large survey or simulation. A common task is to compute the mean value ($m$ ) and sample variance ($s^2$ ) of a data set ($y_i$ ) with $i=1...N$ . The data might be of observations of some quantity or the results of performing a Monte Carlo integration. In principle, this seems very simple. In practice, floating point arithmetic can result in some suprising behavior. In this example, you will investigate some of the potential complications of performing even basic mathematical calculations.
Keep working on it!
The answer is not quite right.
Missing Response
The variable response_2c is still set to missing.
h. I've written some tests in 'test/test2.jl'. Because of Pluto's reactivity, it's tricky to run a file from inside a notebook. Instead, run julia --project test/runtests2.jl to run the code in this Pluto notebook and then the tests in 'test/test2.jl'. First, check that your functions pass my tests. If not, is it because your function has a bug? If so, fix your functions. Or is there another explanation? It may help to look at the source code for the tests to see what it means to have "passed".
Can you suggest additional tests for such functions? Feel free to add them to the tests in 'test/test2.jl' and check that your code still passes.
Missing Response
The variable response_1a is still set to missing.
Your code should pass the following tests. If it doesn't, fix your code so it does.
Missing Response
The variable response_1c is still set to missing.
Astro 528, Lab 1, Exercise 2
Introduction to Julia & Floating Point Arithmetic
€mime©text/htmlŹrootassigneeÀČlast_run_timestampËAÚ,YŐ0°persist_js_state·has_pluto_hook_features§cell_idÙ$27667e0a-8ebc-4397-8ac3-33a0f19f6987čdepends_on_disabled_cells§runtimeÎ ”published_object_keysždepends_on_skipped_cells§erroredÂÙ$48a888f3-7067-4524-b818-279e3ed2ffdcŠqueued€logs§runningŠoutput€bodyÙMissing Response
The variable response_2f is still set to missing.
The first time you execute any command in Julia (or start up a notebook), you'll notice a delay while the Julia kernel starts. Then, the first time you import a module, Julia will parse the code in the module and compile some functions. The next you import the same module it will be mucuh faster, as it won't need to reparse and recompile some of that module's code (as long as the module hasn't changed, e.g., you modify the module's code directly or due to the package being updated).
Tip
Note the cell above returns a NamedTupleTuple
Missing Response
The variable response_2d is still set to missing.
g. Don't forget that we should test your functions for accuracy. Should we expect all of these functions to return the exact same value? How can we test functions that return floating point values?
d. What considerations would affect the decision of whether to use the one-pass algorithm or the two-pass algorithm?
Calculating basic summary statistics
Now, we're going to generate a much larger sample of numbers and compute their mean and standard deviations using multiple different methods. You will compare the results. The goal is to help you to understand when floating point arithmetic is likely to be problematic, so you can anticipate potential pitfalls that might affect your own research.
Change the value of the variable num_obs defined in a cell above to smaller and larger values. How does the mangitude of the differnces depend on the number of observation dates?
Hint
The above could also be written more succinctly as
"Calculate mean value of an array using sum and length functions."
mean_demo_concise(y::Array) = sum(y)/length(y);Indeed, Julia's function Statistics.mean() that is written almost identically to this.
b. Write a function named var_two_pass take takes input similar to mean_demo_verbose and provides a two-pass algorithm to calculate the variance more accurately than the one pass algoritihm by using two loops over the $y_i$ 's.
a. Write a function named var_one_pass that takes inputs similar to mean_demo_verbose and implements a one-pass algorithm to calculate the variance, reading each value of y from the computer's main memory only once. Note that using the same element of an array repeatedly (i.e., before accessing the any other elements of the array) only counts as a single pass, since it can be reused without repeatedly copying the data from main memory.
The cell above assigns multiple variables. When writing Pluto notebooks, any cell that assigns multiple variables must be wrapped inside a begin...end block (or split into multiple cells). Note that this is different from Jupyter notebooks. The final line calls the functions mean and std (that were exported by the Statistics package) to compute the mean and standard deviation of our sample.
Using the same mean and std function as before, compute (and report) the sample mean and sample variance for each of these arrays. Compare the results by subtracting each of the results computed using Float64's and Float32's
Generate an array of simulated data using the following function:
€mime©text/htmlŹrootassigneeÀČlast_run_timestampËAÚ,Yâpm°persist_js_state·has_pluto_hook_features§cell_idÙ$c9b2756d-46cf-43a1-81a7-2aecf50fd69ečdepends_on_disabled_cells§runtimeÎ żX”published_object_keysždepends_on_skipped_cells§erroredÂÙ$66a5de37-ff4c-40f6-99fc-624ca571b881Šqueued€logs§runningŠoutput€bodyÚThe above code calls the function generate_sample, asking it to compute 10 random variables with true mean equal to the julian date for September 1, 2021. The output will be a list of floating point numbers enclosed in square brackets to denote that it's a vector, which is equivalent to a 1-dimensional array.
Look at the results above. Are the output consistent with your expectations? (If not, then try changing the inputs to generate_sample to see what happens.) Write your responce as Markdown text in the cell below and store the result as a variable named response_1a.
Full Width Mode
€mime©text/htmlŹrootassigneeÀČlast_run_timestampËAÚ,ïĄä°persist_js_state·has_pluto_hook_features§cell_idÙ$9ca08091-8906-4b2e-b965-f8dbc385623dčdepends_on_disabled_cells§runtimeÎ j B”published_object_keysždepends_on_skipped_cells§erroredÂÙ$37b12cb7-377b-48ba-b3ff-2457aa2c44a9Šqueued€logs§runningŠoutput€bodyÚFinite precision arithmetic
By default, Julia uses 64 bits of memory to store each floating point value. Often this is referred to as "double precission" (for historical reasons, although technically this is machine dependent and thus less precise) to differentiate it from "single precission" floating point values typically stored with 32 bits. To explore the effects of floating point arithmetic, let us convert the array of y values into arrays of floating point values that use fewer bits to store each value using the following code.
Tip
Julia allows unicode characters for variable and function names. This can be very useful for mathematical work. However, some programs don't display unicode characters correctly.
c. Compare the accuracy of the results using data sets of different sizes and values of the true sample mean. Under what conditions do they give results that differ by an ammount that is potentially scientifically significant?
To make Plots we'll import the Plots package. (If you're interested, you can click the eyeball to the left of the plot cells to see the plotting code.)
Writing Functions
It will be very useful to write and organize your code into many small functions. I strongly recommend you develop a habit of writing code in the form of functions. A good rule of thumb is that each function should do one specific thing. Another rule of thumb is to try to keep each function to no more than can fit on one page of paper (or one screen), even when it's complicated. If nothing else, this makes it easier for humans to debug the function. The code for most functions is considerably smaller, but sometimes a hard scientific problem demands a longer function. Often, after writing a complex function, one can refactor the code into multiple smaller functions, resulting in code that is easier to read, debug, maintain and optimize. Julia provides multiple syntaxes for writing functions, as described in the Julia manual . (I suggest stopping before the subsection on "Operators Are Functions" for now.)
It's often good to double check the return type of a function you call to make sure it's what you expected. For functions in Base Julia, this can usually be looked up in the function documentation, either in the Julia manual or using the Live Docs feature of Pluto. If you want to check a variable's type, the typeof function is quite useful.
Missing Response
The variable response_2g is still set to missing.
Using Packages
The Julia language includes many powerful features. While many of the most commonly used functions and macros are available by default (such as sqrt above), other functions are only avaliable if you import a module. For the first part of this exercise, we'll be using the Dates, Random and Statistics modules. (The Dates, Random, and Statistics are part of Julia's standard library, a set of modules that are distributed with Julia.)
To be able to access functions in a module, you execute import MyModule and then execute MyModule.fn(x) to call a function named fn with parameter x. Alternatively, using MyModule will import all the functions that the Module module has specified should be exported by default. Basically, this means you don't have to write MyModule. before every call to a function that MyModule intends for end users to call. Often, using is very convenient. For very common functions names (e.g., mean, apply) using risks creating confusion about which module is being reference. In these caeses import provides more control of exactly which functions are loaded into the current namespace.
Missing Response
The variable response_1b is still set to missing.
= 1 &&
!ismissing(var_two_pass(ones(10))) &&
var_two_pass(ones(10)) â 0 &&
var_two_pass([0,1,2,3,4,5,6,7,8,9,10]) â 11
correct()
else
keep_working()
end
end
endšmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$569df89d-5039-4a63-8396-ab595811584c§cell_idÙ$569df89d-5039-4a63-8396-ab595811584c€codeÙJdisplay_msg_if_fail(check_type_isa(:response_2c,response_2c,Markdown.MD)) šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$a387d515-82d6-4211-8934-b5f0d3b062dc§cell_idÙ$a387d515-82d6-4211-8934-b5f0d3b062dc€codeÚmd"""
h. I've written some tests in 'test/test2.jl'. Because of Pluto's reactivity, it's tricky to run a file from inside a notebook. Instead, run `julia --project test/runtests2.jl` to run the code in this Pluto notebook and then the tests in 'test/test2.jl'. First, check that your functions pass my tests. If not, is it because your function has a bug? If so, fix your functions. Or is there another explanation?
It may help to look at the source code for the tests to see what it means to have "passed".
Can you suggest additional tests for such functions? Feel free to add them to the tests in 'test/test2.jl' and check that your code still passes.
"""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$b77abd33-0214-46f4-9fde-8b38afafd224§cell_idÙ$b77abd33-0214-46f4-9fde-8b38afafd224€codeÙJdisplay_msg_if_fail(check_type_isa(:response_1a,response_1a,Markdown.MD)) šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$72474fca-4bc7-471e-9116-c48023f147dd§cell_idÙ$72474fca-4bc7-471e-9116-c48023f147dd€codeÙ^md"""
Your code should pass the following tests. If it doesn't, fix your code so it does.
"""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$da5572db-5df8-4753-876f-e1b3a186f8a8§cell_idÙ$da5572db-5df8-4753-876f-e1b3a186f8a8€code·m_32bit = mean(y_32bit)šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$a200eb3c-7041-47e4-89d2-d077dccc18c2§cell_idÙ$a200eb3c-7041-47e4-89d2-d077dccc18c2€codeÙJdisplay_msg_if_fail(check_type_isa(:response_1c,response_1c,Markdown.MD)) šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$86a442a6-fb6e-45c7-9ab9-83aee71b028c§cell_idÙ$86a442a6-fb6e-45c7-9ab9-83aee71b028c€codeÙSfunction var_two_pass(y::Array)
# INSERT CODE for var_two_pass
return missing
endšmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$27667e0a-8ebc-4397-8ac3-33a0f19f6987§cell_idÙ$27667e0a-8ebc-4397-8ac3-33a0f19f6987€codeÙ^md"""
# Astro 528, Lab 1, Exercise 2
## Introduction to Julia & Floating Point Arithmetic
"""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$48a888f3-7067-4524-b818-279e3ed2ffdc§cell_idÙ$48a888f3-7067-4524-b818-279e3ed2ffdc€codeÙJdisplay_msg_if_fail(check_type_isa(:response_2f,response_2f,Markdown.MD)) šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$a802faec-6745-4f9d-820a-bb8c1aa25fc1§cell_idÙ$a802faec-6745-4f9d-820a-bb8c1aa25fc1€codeœtrue_mean_plt = jd_sept1_2021šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$60dc2204-db4f-4038-8158-d0694dd720ba§cell_idÙ$60dc2204-db4f-4038-8158-d0694dd720ba€codeÚómd"""
The first time you execute any command in Julia (or start up a notebook), you'll notice a delay while the Julia kernel starts. Then, the first time you import a module, Julia will parse the code in the module and compile some functions. The next you import the same module it will be mucuh faster, as it won't need to reparse and recompile some of that module's code (as long as the module hasn't changed, e.g., you modify the module's code directly or due to the package being updated).
"""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$a4a9b516-5473-478e-a390-9e4f715310eb§cell_idÙ$a4a9b516-5473-478e-a390-9e4f715310eb€codeÚ_tip(md"Note the cell above returns a [`NamedTuple`](https://docs.julialang.org/en/v1/base/base/#Core.NamedTuple) that contains two Float64's. Naming the two elements of the [`Tuple`](https://docs.julialang.org/en/v1/manual/types/#Tuple-Types) can be useful for preventing silly mistakes when order is the only way to distinguish the two numbers.
")šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$cd09046b-753e-4843-a4e4-3be0b0c7fb97§cell_idÙ$cd09046b-753e-4843-a4e4-3be0b0c7fb97€codeÚif ready_to_plot && !ismissing(compare_var_calcs(3)) && !ismissing(compare_var_calcs_online(3))
local N_list = [2^n for n in 2:20]
local plt = plot()
local true_mean = jd_sept1_2021
local label_1 = "1-pass minus 2-pass"
local label_online ="online minus 2-pass"
scatter!(plt,N_list, compare_var_calcs.(N_list,true_mean), x_scale=:log10, color=:red,label=label_1, legend=:topleft)
scatter!(plt,N_list, compare_var_calcs_online.(N_list,true_mean), x_scale=:log10,color=:green,label=label_online)
xlabel!(plt,"Number of observation dates")
ylabel!(plt,"|Î Estimated Sample Var| (daysÂČ)")
title!(plt,"Sample mean = $true_mean")
endšmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$48d96e9d-b34b-4899-a976-a92602156981§cell_idÙ$48d96e9d-b34b-4899-a976-a92602156981€codeÙJdisplay_msg_if_fail(check_type_isa(:response_2d,response_2d,Markdown.MD)) šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$339639fc-77d8-4e88-85f3-59c7821cd01f§cell_idÙ$339639fc-77d8-4e88-85f3-59c7821cd01f€codeÙÔmd"""
g. Don't forget that we should test your functions for accuracy. Should we expect all of these functions to return the exact same value? How can we test functions that return floating point values?
"""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$3cb23c2a-2611-47e6-9ee5-4d1d4f0a84dc§cell_idÙ$3cb23c2a-2611-47e6-9ee5-4d1d4f0a84dc€codeÙSbegin
y = generate_sample(num_obs, jd_sept1_2021)
(ÎŒ = mean(y), Ï = std(y))
endšmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$66d9bc94-3e61-41e8-a81d-88e307d97653§cell_idÙ$66d9bc94-3e61-41e8-a81d-88e307d97653€codeÙymd"d. What considerations would affect the decision of whether to use the one-pass algorithm or the two-pass algorithm?"šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$61503eb8-8a70-4ff7-b5bb-0a73c501d6c7§cell_idÙ$61503eb8-8a70-4ff7-b5bb-0a73c501d6c7€codeÚmd"""
### Calculating basic summary statistics
Now, we're going to generate a much larger sample of numbers and compute their
mean and standard deviations using multiple different methods. You will compare the results. The goal is to help you to understand when floating point arithmetic is likely to be problematic, so you can anticipate potential pitfalls that might affect your own research.
"""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$731f047f-0f26-4ab7-8810-398659642b0c§cell_idÙ$731f047f-0f26-4ab7-8810-398659642b0c€codeÙżmd"""
Change the value of the variable `num_obs` defined in a cell above to smaller and larger values.
How does the mangitude of the differnces depend on the number of observation dates?"""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$0b768a24-e97a-47e8-925f-b9e75601ceae§cell_idÙ$0b768a24-e97a-47e8-925f-b9e75601ceae€codeÚ hint(md"""
The above could also be written more succinctly as
```julia
"Calculate mean value of an array using sum and length functions."
mean_demo_concise(y::Array) = sum(y)/length(y);
```
Indeed, Julia's function `Statistics.mean()` that is written almost identically to this.
""") šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$8b3572b8-e571-43b7-ab45-68eabecace69§cell_idÙ$8b3572b8-e571-43b7-ab45-68eabecace69€codeÙôfunction compare_var_calcs(N::Integer, true_mean::Real = 0.0)
@assert N > 2
@assert !isnan(true_mean)
@assert !isinf(true_mean)
input_data = generate_sample(N,true_mean)
Îvar = abs(var_one_pass(input_data) - var_two_pass(input_data))
end;šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$4ab6efe2-271c-4574-898e-ce0817fc5033§cell_idÙ$4ab6efe2-271c-4574-898e-ce0817fc5033€codeÙ3function var_one_pass(y::Array)
return missing
endšmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$58d0e74a-d4f6-4aab-97aa-18d305e888e1§cell_idÙ$58d0e74a-d4f6-4aab-97aa-18d305e888e1€codeÙïmd"""b. Write a function named `var_two_pass` take takes input similar to `mean_demo_verbose` and provides a two-pass algorithm to calculate the variance more accurately than the one pass algoritihm by using two loops over the $y_i$'s."""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$c337d59a-cb1c-4542-ae16-830bf8a2afc5§cell_idÙ$c337d59a-cb1c-4542-ae16-830bf8a2afc5€codeÚÉmd"""a. Write a function named `var_one_pass` that takes inputs similar to `mean_demo_verbose` and implements a **one-pass algorithm** to calculate the variance, reading each value of y from the computer's main memory only once. Note that using the same element of an array repeatedly (i.e., before accessing the any other elements of the array) only counts as a single pass, since it can be reused without repeatedly copying the data from main memory."""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$32831892-6da7-4f85-80ff-73c60a638382§cell_idÙ$32831892-6da7-4f85-80ff-73c60a638382€codeÚ md"""The cell above assigns multiple variables. When writing Pluto notebooks, any cell that assigns multiple variables must be wrapped inside a `begin`...`end` block (or split into multiple cells). Note that this is different from Jupyter notebooks. The final line calls the functions `mean` and `std` (that were exported by the Statistics package) to compute the mean and standard deviation of our sample.
"""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$a1699fca-90ec-418b-a675-3982dd4c11ff§cell_idÙ$a1699fca-90ec-418b-a675-3982dd4c11ff€code©DateTime;šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$76cd6fbe-9be4-4e29-a3e9-4fac87d1a0c8§cell_idÙ$76cd6fbe-9be4-4e29-a3e9-4fac87d1a0c8€codeÙ}md"Once you've completed the questions above and made your prediction, **click this box**: $(@bind ready_to_plot CheckBox())"šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$192e6360-5eba-4a4d-b203-363caba8af64§cell_idÙ$192e6360-5eba-4a4d-b203-363caba8af64€codeÙ"@test var_one_pass(ones(10)) â 0šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$84f0ce6c-9ae0-4016-91e2-d436d4385366§cell_idÙ$84f0ce6c-9ae0-4016-91e2-d436d4385366€codeÙćmd"Using the same mean and std function as before, compute (and report) the sample mean and sample variance for each of these arrays. Compare the results by subtracting each of the results computed using Float64's and Float32's"šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$119e51ef-ed7b-4f9b-b7dd-67ae70bf934a§cell_idÙ$119e51ef-ed7b-4f9b-b7dd-67ae70bf934a€codeÚuif ready_to_plot && !ismissing(compare_var_calcs(3))
local N_list = [2^n for n in 2:20]
local true_mean = jd_sept1_2021
local plt = scatter(N_list, compare_var_calcs.(N_list,true_mean), x_scale=:log10, legend=:none)
xlabel!(plt,"Number of observation dates")
ylabel!(plt,"|Î Estimated Sample Var| (daysÂČ)")
title!(plt,"Sample mean = $true_mean (Static Plot)")
endšmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$10dca248-2004-4719-9e30-eb3025da0513§cell_idÙ$10dca248-2004-4719-9e30-eb3025da0513€code”response_2c = missingšmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$8c2cddbf-3a02-4969-952e-4d76ca23f95b§cell_idÙ$8c2cddbf-3a02-4969-952e-4d76ca23f95b€code”response_1b = missingšmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$c9b2756d-46cf-43a1-81a7-2aecf50fd69e§cell_idÙ$c9b2756d-46cf-43a1-81a7-2aecf50fd69e€codeÙOmd"""
### Generate an array of simulated data using the following function:
"""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$66a5de37-ff4c-40f6-99fc-624ca571b881§cell_idÙ$66a5de37-ff4c-40f6-99fc-624ca571b881€codeÚRmd"""The above code calls the function generate_sample, asking it to compute 10
random variables with true mean equal to the julian date for September 1, 2021. The output will be a list of floating point numbers enclosed in square brackets to denote that it's a vector, which is equivalent to a 1-dimensional array.
Look at the results above. Are the output consistent with your expectations? (If not, then try changing the inputs to `generate_sample` to see what happens.) Write your responce as Markdown text in the cell below and store the result as a variable named `response_1a`.
"""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$0c9c629c-cb43-4fa7-8322-0b47718f2a9a§cell_idÙ$0c9c629c-cb43-4fa7-8322-0b47718f2a9a€code«using Plotsšmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$9ca08091-8906-4b2e-b965-f8dbc385623d§cell_idÙ$9ca08091-8906-4b2e-b965-f8dbc385623d€codeŻWidthOverDocs()šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$37b12cb7-377b-48ba-b3ff-2457aa2c44a9§cell_idÙ$37b12cb7-377b-48ba-b3ff-2457aa2c44a9€codeÚ*md"""
## Finite precision arithmetic
By default, Julia uses 64 bits of memory to store each floating point value. Often this is referred to as "double precission" (for historical reasons, although technically this is machine dependent and thus less precise) to differentiate it from "single precission" floating point values typically stored with 32 bits. To explore the effects of floating point arithmetic, let us convert the array of y values into arrays of floating point values that use fewer bits to store each value using the following code."""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$a6cfeba5-d7c7-4493-aeae-379569932ef0§cell_idÙ$a6cfeba5-d7c7-4493-aeae-379569932ef0€codeÙœtip(md"Julia allows unicode characters for variable and function names. This can be very useful for mathematical work. However, some programs don't display unicode characters correctly.")šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$1907ec31-4b3f-4db6-a42f-fffb3b722d7e§cell_idÙ$1907ec31-4b3f-4db6-a42f-fffb3b722d7e€codeżusing Dates, Random, Statisticsšmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$64dcecf3-1561-411a-8759-a4ccb219e303§cell_idÙ$64dcecf3-1561-411a-8759-a4ccb219e303€code”response_2f = missingšmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$d470aea9-b69e-4f02-bbaf-4d61cb5244b9§cell_idÙ$d470aea9-b69e-4f02-bbaf-4d61cb5244b9€codeÙémd"c. Compare the accuracy of the results using data sets of different sizes and values of the true sample mean. Under what conditions do they give results that differ by an ammount that is potentially scientifically significant?"šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$43f717de-f5c8-43b5-9229-254b9cb89ca7§cell_idÙ$43f717de-f5c8-43b5-9229-254b9cb89ca7€codeÙmd"To make Plots we'll import the `Plots` package. (If you're interested, you can click the eyeball to the left of the plot cells to see the plotting code.)"šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$98861118-c1c1-48d3-bb7b-1dc8f8e52604§cell_idÙ$98861118-c1c1-48d3-bb7b-1dc8f8e52604€codeÚúmd"""
## Writing Functions
It will be very useful to write and organize your code into many small functions. I strongly recommend you develop a habit of writing code in the form of functions. A good rule of thumb is that each function should do one specific thing. Another rule of thumb is to try to keep each function to no more than can fit on one page of paper (or one screen), even when it's complicated. If nothing else, this makes it easier for humans to debug the function. The code for most functions is considerably smaller, but sometimes a hard scientific problem demands a longer function. Often, after writing a complex function, one can refactor the code into multiple smaller functions, resulting in code that is easier to read, debug, maintain and optimize. Julia provides multiple syntaxes for writing functions, [as described in the Julia manual](https://docs.julialang.org/en/v1/manual/functions/index.html). (I suggest stopping before the subsection on "Operators Are Functions" for now.)
"""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$a11ffb3d-310f-4e28-b8f2-724aab7006a0§cell_idÙ$a11ffb3d-310f-4e28-b8f2-724aab7006a0€code”response_2h = missingšmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$9f4c6a89-754d-4dee-813f-b2d902439ea1§cell_idÙ$9f4c6a89-754d-4dee-813f-b2d902439ea1€code”typeof(jd_sept1_2021)šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$8b757575-9b7e-4154-8222-d024cb62f08f§cell_idÙ$8b757575-9b7e-4154-8222-d024cb62f08f€codeÙ+jd_sept1_2021 = datetime2julian(sept1_2021)šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$29ab9da1-2130-4ea1-aa4d-af08f8011bd0§cell_idÙ$29ab9da1-2130-4ea1-aa4d-af08f8011bd0€codeÚmd"""
It's often good to double check the return type of a function you call to make sure it's what you expected. For functions in Base Julia, this can usually be looked up in the function documentation, either in the [Julia manual](https://docs.julialang.org/en/v1/) or using the Live Docs feature of Pluto. If you want to check a variable's type, the `typeof` function is quite useful."""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$691632b7-3435-4c0d-aff4-3bdc87d77e7b§cell_idÙ$691632b7-3435-4c0d-aff4-3bdc87d77e7b€code”y_32bit = Float32.(y)šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$f362498a-fe8e-440d-afab-c817545b3144§cell_idÙ$f362498a-fe8e-440d-afab-c817545b3144€codeÙJdisplay_msg_if_fail(check_type_isa(:response_2g,response_2g,Markdown.MD)) šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$5bcf3076-f31b-47e8-8297-8cf406ff71ab§cell_idÙ$5bcf3076-f31b-47e8-8297-8cf406ff71ab€code±m_64bit = mean(y)šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$9ad35a1b-3cc1-475a-afd6-8f340a86cdd0§cell_idÙ$9ad35a1b-3cc1-475a-afd6-8f340a86cdd0€codeÚmd"""
## Using Packages
The Julia language includes many powerful features. While many of the most commonly used functions and macros are available by default (such as `sqrt` above), other functions are only avaliable if you *import* a module. For the first part of this exercise, we'll be using the Dates, Random and Statistics modules. (The Dates, Random, and Statistics are part of Julia's standard library, a set of modules that are distributed with Julia.)
To be able to access functions in a module, you execute `import MyModule` and then execute `MyModule.fn(x)` to call a function named `fn` with parameter `x`.
Alternatively, `using MyModule` will import all the functions that the Module module has specified should be *exported* by default. Basically, this means you don't have to write `MyModule.` before every call to a function that `MyModule` intends for end users to call. Often, `using` is very convenient. For very common functions names (e.g., `mean`, `apply`) `using` risks creating confusion about which module is being reference. In these caeses `import` provides more control of exactly which functions are loaded into the current namespace.
"""šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedĂÙ$fe4601cb-0cc3-4ac1-ae18-aa5fc7c35bb5§cell_idÙ$fe4601cb-0cc3-4ac1-ae18-aa5fc7c35bb5€code”response_1c = missingšmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$5f84a3cc-dab6-4ad0-9644-7ea803f43475§cell_idÙ$5f84a3cc-dab6-4ad0-9644-7ea803f43475€codeÙJdisplay_msg_if_fail(check_type_isa(:response_1b,response_1b,Markdown.MD)) šmetadata©show_logsĂšdisabledÂźskip_as_script«code_foldedÂÙ$ba1bac88-5b7d-4b21-991e-948fb00fc2bb§cell_idÙ$ba1bac88-5b7d-4b21-991e-948fb00fc2bb€codenum_obs = 100šmetadata©show_logsĂšdisabledÂźskip_as_script«code_folded«notebook_idÙ$ffb28f94-8543-11f0-120d-b33132eb5b69«in_temp_diršmetadata 