6 Control Structures Homework
Allen Aguas
7/6/24
Instructions
There are four exercises below. Three exercises are worth 10 points each, and are labelled (Either R, Python or Julia). For these exercises, you can implement a solution in the language of your choice. I’ve included code chunk templates for all three languages, but you only need to fill in the chunks for one language.
One exercise is labelled (R, Python and Julia). This exercise is worth 20 points, and you must provide a solution using all three languages. This exercise also includes code chunk templates for all three languages.
Some of the exercises have questions relating to the results. Write your answers in the narrative part of the text, formatted as italicized text, bold text or
quoted text.
Submit both Rmd and typeset results using file name
formats as in previous exercises.
Do not use libraries in R that are not part of the base R
installation (with one exception as stated below); this includes popular
libraries for data tables like data.table or tidyverse
libraries like dplyr and tidyr. You may use
any Python or Julia libraries introduced in lecture, except for graphs.
You are free to use any Python or Julia graphics library, as long as I
can install it on my machine if necessary.
Reuse
We will be using your NormPDF function for two of these
exercises. Define that function here.
R
NormPDF <- function(x, mu = 0, sigma = 1) {
return((1 / sqrt(2 * pi * sigma^2)) * exp(-((x - mu)^2) / (2 * sigma^2)))
}Python
import math
def NormPdf(x, mu=0, sigma=1):
coefficient = 1 / (sigma * math.sqrt(2 * math.pi))
exponent = math.exp(- ((x - mu) ** 2) / (2 * sigma ** 2))
return coefficient * exponentJulia
function NormPDF(x, mu=0.0, sigma=1.0)
return (1 / sqrt(2 * π * sigma^2)) * exp(-((x - mu)^2) / (2 * sigma^2))
end## NormPDF (generic function with 3 methods)Exercise 1 (Either R, Python or Julia)
Part a.
Write a function MyMoments to compute mean, standard
deviation, skewness and kurtosis from a single vector of numeric values.
You can use library mean functions, but should use one (and
only one) for loop to compute the rest. Note that computationally
efficient implementations of moments take advantage of \((Y_i-\bar{Y})^4 = (Y_i-\bar{Y}) \times
(Y_i-\bar{Y})^3\), etc.
See https://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm for formula for skewness and kurtosis. This reference gives several definitions for both skewness and kurtosis, you only need to implement one formula for each. Note that for computing skewness and kurtosis, standard deviation is computed using \(N\) as a divisor, not \(N-1\).
Your function should return a list or a tuple with Mean,
SD, Skewness and Kurtosis.
Your function should check for missing values, and accumulate sums
only if values are non-missing. You might need to pass arguments to
specify how to handle missing values if you use library functions to
compute a mean (na.rm in R, nanmean in
numpy or skipmissing in Julia)
R
Python
import numpy as np
def MyMoments(x, na_rm=True):
if na_rm:
x = x[~np.isnan(x)]
n = len(x)
mean_x = np.mean(x)
sum_sq_diff = 0
sum_cub_diff = 0
sum_quart_diff = 0
for i in range(n):
diff = x[i] - mean_x
sum_sq_diff += diff**2
sum_cub_diff += diff**3
sum_quart_diff += diff**4
sd_x = np.sqrt(sum_sq_diff / n)
skewness_x = (sum_cub_diff / n) / (sd_x**3)
kurtosis_x = (sum_quart_diff / n) / (sd_x**4) - 3
return mean_x, sd_x, skewness_x, kurtosis_xJulia
Part b.
Test your function by computing moments for MPG from the
file mpg.csv. This file will be available on D2L under Week
5 Lecture Materials. Note that there are missing values in this column,
so will be testing your functions ability to process a vector with
missing values.
R
Python
import pandas as pd
from scipy.stats import skew, kurtosis
# Load the data
mpg_data = pd.read_csv(r"C:\Users\Allen\OneDrive - Dakota State University\Summer 24\Statistical Programming 600\Week 6\HW\mpg.csv")
# Extract the 'MPG' column
mpg = mpg_data['MPG'].to_numpy()
# Compute moments using MyMoments
mean, sd, skewness, kurtosis_val = MyMoments(mpg)
# Print the results
print("MyMoments Results:")## MyMoments Results:print(f"Mean: {mean}")## Mean: 33.977000000000004print(f"Standard Deviation: {sd}")## Standard Deviation: 2.7808220007760287print(f"Skewness: {skewness}")## Skewness: 0.38610542371158807print(f"Kurtosis: {kurtosis_val}")## Kurtosis: -0.5258174866779437Julia
If you wish, compare your function results with the
skewness and kurtosis in the R
moments package (you may need to use
na.rm=TRUE), or the skew and
kurtosis functions from scipy.stats (you may
need to call .to_list() and use
nan_policy='omit'), or the StatsBase functions
skewness and kurtosis in Julia (You may need
to use collect(skipmissing(...)). Note that some functions
may return excess kurtosis (Pearson’s kurtosis minus 3).
R
Python
# Compare with scipy.stats functions
skewness_scipy = skew(mpg, nan_policy='omit')
kurtosis_scipy = kurtosis(mpg, nan_policy='omit')
print("\nscipy.stats Results:")##
## scipy.stats Results:print(f"Skewness: {skewness_scipy}")## Skewness: 0.38610542371158785print(f"Kurtosis: {kurtosis_scipy}")## Kurtosis: -0.5258174866779455Julia
Exercise 2 (Either R, Python or Julia)
Consider Newton’s method to find a minimum or maximum value attained by a function over an interval. Given a function \(f\), we wish to find
\[ \max_{x \in [a,b]} f(x) \]
Start with an initial guess, \(x_0\), then generate a sequence of guesses using the formula
\[ x_{k+1} = x_{k} - \frac{f'(x_{k})}{f''(x_{k})} \]
where \(f'\) and \(f''\) are first and second derivatives. We won’t be finding derivatives analytically, instead, we will be using numerical approximations (central finite differences), given by
\[ \begin{aligned} f' & \approx \frac{f(x+\frac{h}{2}) - f(x-\frac{h}{2})}{h} \\ f'' & \approx \frac{f(x+h) - 2f(x)+f(x-h)}{h^2} \end{aligned} \] where \(h\) is some arbitrary small value.
We will work with the normal pdf, \(f (x ; \mu, \sigma^2) = \frac{1}{\sigma \sqrt{2 \pi}^{}} e^{- \frac{(x - \mu)^2}{2 \sigma^2}}\). Let \(\mu = m_{1936}\) be the mean Calories per Serving from 1936, and let \(\sigma = s_{1936}\) be the corresponding standard deviation. We will wish to find the \(x_*\) that maximizes
\[ \max pdf (x ; m_{1936}, s_{1936}^2) \]
Let the initial guess be \(x_0 = 180\) and let \(h = 0.1\). Calculate 10 successive \(x_k\), saving each value in a vector. Print the final \(x_k\). Why does this value maximize the likelihood function?
R
Python
import numpy as np
import matplotlib.pyplot as plt
# Constants
mu = 1936
sigma = 1936
x0 = 180
h = 0.1
tol = 1e-6
max_iter = 100
# Normal PDF function
def normal_pdf(x, mu, sigma):
return (1 / (np.sqrt(2 * np.pi * sigma**2))) * np.exp(-((x - mu)**2) / (2 * sigma**2))
# First derivative using central finite difference
def first_derivative(f, x, h, mu, sigma):
return (f(x + h, mu, sigma) - f(x - h, mu, sigma)) / (2 * h)
# Second derivative using central finite difference
def second_derivative(f, x, h, mu, sigma):
return (f(x + h, mu, sigma) - 2 * f(x, mu, sigma) + f(x - h, mu, sigma)) / (h**2)
# Newton's method for maximization
xk = x0
iterations = [xk]
for i in range(max_iter):
f_prime = first_derivative(normal_pdf, xk, h, mu, sigma)
f_double_prime = second_derivative(normal_pdf, xk, h, mu, sigma)
if f_double_prime == 0:
break
xk_new = xk - (f_prime / f_double_prime)
if abs(xk_new - xk) < tol:
break
xk = xk_new
iterations.append(xk)Julia
Part b.
Plot the sequence of \(x\) versus iteration number (\(k\)) as the independent variable. Add a horizontal line corresponding to \(m_{1936}=268.1\). How many iterations are required until \(|x_{k+1} - x_{k}| < 10^{-6}\)?
It would be a number of 5 iterations, for the difference between consecutive X values to be less than 10^-6,indicating convergence.
R
Python
# Plotting
plt.figure(figsize=(10, 6))
plt.plot(iterations, marker='o')
plt.axhline(y=268.1, color='r', linestyle='--', label='y = 268.1')
plt.xlabel('Iteration number (k)')
plt.ylabel('x value')
plt.title('x value vs Iteration number')
plt.legend()
plt.grid(True)
plt.show()
print(f"Final xk value: {xk}")## Final xk value: 30622.631850966536print(f"Number of iterations: {len(iterations)}")## Number of iterations: 101Julia
Exercise 3 (Either R, Python or Julia)
Consider the Trapezoidal Rule for integration. From “Analysis by Its History” (https://books.google.com/books/about/Analysis_by_Its_History.html?id=E2IhMXPZMNIC)
On the interval \(\left[ x_i, x_{i+1}\right]\) the function \(f(x)\) is replaced by a straight line passing through \(\left(x_i,f(x_i)\right)\) and \(\left(x_{i+1},f(x_{i+1})\right)\). The integral between \(x_i\) and \(x_{i+1}\) is then approximated by the trapezoidal area \(h \cdot \left(f(x_i)+f(x_{i+1})\right)/2\) and we obtain
\[ \int _{a} ^{b} f(x) dx = F(x) \approx \sum _{i=1} ^{N-1} \frac{h}{2} \left(f(x_i)+f(x_{i+1})\right) \]
We will calculate the integral for the normal pdf
\[ \int _{-3} ^{3} f (x ; \mu, \sigma^2) dx = \int _{-3} ^{3} \frac{1}{\sigma \sqrt{2 \pi}^{}} e^{- \frac{(x - \mu)^2}{2 \sigma^2}} dx \]
with \(\mu=0\) and \(sigma=1\), using your NormPdf
function. We will do this by creating a sequence of approximations, each
more precise than the preceding approximation.
Part a.
Calculate a first approximation of step size \(h_0=1\), using the sequence of \(x_i = \left\{-3.0,-2.0,-1.0,0.0,1.0,2.0,3.0\right\}\). Let this approximation be \(F_0\). Print the first approximation.
R
Python
import numpy as np
# Normal PDF function
def normal_pdf(x, mu, sigma):
return (1 / (np.sqrt(2 * np.pi * sigma**2))) * np.exp(-((x - mu)**2) / (2 * sigma**2))
# Trapezoidal rule function
def trapezoidal_rule(f, a, b, n, mu, sigma):
x = np.linspace(a, b, n + 1)
h = (b - a) / n
integral = 0.5 * f(x[0], mu, sigma) + 0.5 * f(x[-1], mu, sigma)
integral += np.sum(f(x[1:-1], mu, sigma))
integral *= h
return integral
# Parameters
mu = 0
sigma = 1
a = -3
b = 3
# Part a: First approximation
h0 = 1
x0 = np.arange(a, b + h0, h0)
F0 = trapezoidal_rule(normal_pdf, a, b, len(x0) - 1, mu, sigma)
print(f"F0: {F0}")## F0: 0.9952975108780335Julia
Part b.
Continue to calculate a series of approximations \(F_0, F_1, F_2, \dots\) such that \(F_{k+1}\) improves on \(F_k\) by increasing \(N\). Do this by decreasing the step size by
2, \(h_{k+1} = h_{k}/2\). Thus, the
sequence used to calculate \(F_1\) will
be of the form \(x_i = \left\{-3.0, -2.5,
-2.0, -1.5, -1.0, \dots, 1.5, 2.0, 2.5, 3.0 \right\}\). If you
use Python, you may need to use numpy.arange to produce a
range of float values.
Calculate the first 10 approximations in the series and print the final approximation.
R
Python
# Part b: Successive approximations
approximations = [F0]
for k in range(1, 11):
n = 2**k * 3 # Number of intervals
Fk = trapezoidal_rule(normal_pdf, a, b, n, mu, sigma)
approximations.append(Fk)
print(f"F{k}: {Fk}")## F1: 0.9952975108780335
## F2: 0.9967599783627767
## F3: 0.997162572936611
## F4: 0.997265634200434
## F5: 0.9972915513637273
## F6: 0.9972980401595793
## F7: 0.997299662952827
## F8: 0.9973000686882851
## F9: 0.9973001701244715
## F10: 0.9973001954836631Julia
Part c.
Plot the successive approximations \(F_i\) against iteration number (you will
need to define an array to store each approximation). Add a horizontal
line for the expected value
(pnorm(3, lower.tail = TRUE)-pnorm(-3, lower.tail = TRUE)).
Set y-axis limits for this plot to be \([0.92,1]\) to best view the progression of
approximations.
It is common practice to terminate a sequence of approximations when the difference between successive approximations is less than some small value. What is the difference between your final two approximations (It should be less than \(10^{-6}\))?
R
Python
import matplotlib.pyplot as plt
# Part c: Plotting
plt.figure(figsize=(10, 6))
plt.plot(approximations, marker='o', label='Approximations $F_i$')
plt.axhline(y=0.9973, color='r', linestyle='--', label='Expected value')
plt.xlabel('Iteration number (i)')
plt.ylabel('Approximation $F_i$')
plt.title('Successive Approximations of the Integral')
plt.ylim([0.925, 1.25])## (0.925, 1.25)plt.yticks(np.arange(0.925, 1.275, 0.025)) # Set y-axis ticks from 0.925 to 1.25 with increments of 0.025## ([<matplotlib.axis.YTick object at 0x00000140562222E0>, <matplotlib.axis.YTick object at 0x000001405621CC40>, <matplotlib.axis.YTick object at 0x000001405621C160>, <matplotlib.axis.YTick object at 0x00000140562756D0>, <matplotlib.axis.YTick object at 0x0000014056275E20>, <matplotlib.axis.YTick object at 0x000001405627BC70>, <matplotlib.axis.YTick object at 0x0000014056281760>, <matplotlib.axis.YTick object at 0x0000014056288250>, <matplotlib.axis.YTick object at 0x0000014056288D00>, <matplotlib.axis.YTick object at 0x000001405628D7F0>, <matplotlib.axis.YTick object at 0x00000140562932E0>, <matplotlib.axis.YTick object at 0x0000014056293D90>, <matplotlib.axis.YTick object at 0x0000014056299880>, <matplotlib.axis.YTick object at 0x000001405629E370>], [Text(0, 0.925, '0.925'), Text(0, 0.9500000000000001, '0.950'), Text(0, 0.9750000000000001, '0.975'), Text(0, 1.0, '1.000'), Text(0, 1.0250000000000001, '1.025'), Text(0, 1.0500000000000003, '1.050'), Text(0, 1.0750000000000002, '1.075'), Text(0, 1.1, '1.100'), Text(0, 1.1250000000000002, '1.125'), Text(0, 1.1500000000000004, '1.150'), Text(0, 1.1750000000000003, '1.175'), Text(0, 1.2000000000000002, '1.200'), Text(0, 1.2250000000000003, '1.225'), Text(0, 1.2500000000000004, '1.250')])plt.legend()
plt.grid(True)
plt.show()
# Difference between the final two approximations
difference = abs(approximations[-1] - approximations[-2])
print(f"Difference between the final two approximations: {difference}")## Difference between the final two approximations: 2.5359191635487832e-08Julia
Exercise 4. (R, Python and Julia)
You will be simulating a series of Binomial experiments. A simple
example of a Binomial experiment is the number of heads in a series of
coin flips. Your simulation will require two nested for
loops.
Part a.
First, define a vector of length 11, populated with 0s. This will
represent a count of the number of times heads appear in a series of 10
flips. Call this vector Successes.
R
Successes <- rep(0, 11)
print(Successes)## [1] 0 0 0 0 0 0 0 0 0 0 0Python
Successes = [0] * 11
print(Successes)## [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]Julia
Successes = fill(0, 11)## 11-element Vector{Int64}:
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0println(Successes)## [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]Part b.
Next, write two loops. The outer loop will run from 1 to 1000. This
will allow use to repeat the simple experiment 1000 times. The inner
loop will perform the experiment, and will run from 1 to 10. Remember to
adjust the indices for Python. You might use i as counter
for the inner loop, and j as a counter for the inner loop,
but we won’t need those variables.
Inside the outer loop, but outside the inner loop, define an
accumulator variable Count, setting this value to 0. Inside
the inner loop, perform a random coin toss by drawing from a uniform
random sample (see this week’s lecture notes for code to simulate a coin
toss). If the random value is greater than 0.5, count this as a success,
and increment Count by 1. Remember to set
Count to 0 before the start of the inner loop.
After the inner loop executes, Count will have a value
between 0 (no heads) and 10 (all heads). You can now use
Count as an index into Successes. After each
iteration of the inner loop, increment Successes[Count] by
1. (You may need to use Count+1 as an index; remember
Count can be 0).
We will be using Successes to sum the number of times
each possible value of Count occurs in the 1000 repetitions
of the simple experiment. Increment Successes outside the
inner loop, but inside the outer loop.
R
set.seed(123) # For reproducibility
Successes <- rep(0, 11)
for (i in 1:1000) {
Count <- 0
for (j in 1:10) {
if (runif(1) > 0.5) {
Count <- Count + 1
}
}
Successes[Count + 1] <- Successes[Count + 1] + 1
}
print(Successes)## [1] 0 12 46 115 224 236 201 125 32 9 0Python
import random
random.seed(123) # For reproducibility
Successes = [0] * 11
for _ in range(1000):
Count = 0
for _ in range(10):
if random.random() > 0.5:
Count += 1
Successes[Count] += 1
print(Successes)## [1, 12, 48, 106, 204, 231, 224, 111, 55, 8, 0]Julia
using Random
Random.seed!(123) # For reproducibility## TaskLocalRNG()Successes = fill(0, 11)## 11-element Vector{Int64}:
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
for i in 1:1000
Count = 0
for j in 1:10
if rand() > 0.5
Count += 1
end
end
Successes[Count + 1] += 1
end
println(Successes)## [1, 7, 37, 114, 194, 231, 222, 133, 49, 8, 4]Part c.
Divide each element of Successes by 1000.
Successes will now represent the proportion of success out
of \(n=10\) attempts. Plot
Successes against a vector of values from 0 to 10 as
points. Compare this with your BinomPMF function, where
x takes values from 0 to 10, n = 10 and
pi = 0.5. Add the values produced by BinomPMF
to the plot of Successes as a line. How well did the
simulation compare to the theoretical values from
BinomPMF?
The points and the line are closely aligned; in other words the simulation and theoretical values from the ‘BinomPMF’ accurately represents the theoretical distribution.
R
Successes <- Successes / 1000
x <- 0:10
plot(x, Successes, type="b", col="blue", xlab="Number of Heads", ylab="Proportion", main="Simulation vs Binomial PMF")
# Binomial PMF
n <- 10
p <- 0.5
BinomPMF <- dbinom(x, n, p)
lines(x, BinomPMF, type="b", col="red")
legend("topright", legend=c("Simulation", "Binomial PMF"), col=c("blue", "red"), lty=1)
Python
import matplotlib.pyplot as plt
from scipy.stats import binom
Successes = [s / 1000 for s in Successes]
x = list(range(11))
plt.plot(x, Successes, 'bo-', label='Simulation')## [<matplotlib.lines.Line2D object at 0x00000140568A1D60>]# Binomial PMF
n = 10
p = 0.5
BinomPMF = [binom.pmf(k, n, p) for k in x]
plt.plot(x, BinomPMF, 'ro-', label='Binomial PMF')## [<matplotlib.lines.Line2D object at 0x00000140568C23D0>]plt.xlabel('Number of Heads')## Text(0.5, 0, 'Number of Heads')plt.ylabel('Proportion')## Text(0, 0.5, 'Proportion')plt.title('Simulation vs Binomial PMF')## Text(0.5, 1.0, 'Simulation vs Binomial PMF')plt.legend()## <matplotlib.legend.Legend object at 0x00000140568CDFD0>plt.show()
Julia
using Plots
using Distributions
Successes = Successes ./ 1000## 11-element Vector{Float64}:
## 0.001
## 0.007
## 0.037
## 0.114
## 0.194
## 0.231
## 0.222
## 0.133
## 0.049
## 0.008
## 0.004x = 0:10## 0:10
plot(x, Successes, label="Simulation", seriestype=:scatter, color=:blue, xlabel="Number of Heads", ylabel="Proportion", title="Simulation vs Binomial PMF")
# Binomial PMF
n = 10## 10p = 0.5## 0.5BinomPMF = pdf.(Binomial(n, p), x)## 11-element Vector{Float64}:
## 0.0009765625
## 0.00976562500000001
## 0.04394531249999999
## 0.1171875000000004
## 0.20507812500000033
## 0.24609375000000022
## 0.20507812500000033
## 0.1171875000000004
## 0.04394531249999999
## 0.00976562500000001
## 0.0009765625plot!(x, BinomPMF, label="Binomial PMF", seriestype=:line, color=:red)