My main website is https://brandonrozek.com. Though I like the customization that Wordpress allows, it does come with additional overhead.
+Therefore, I decided to use the Tilda space offered by the Computer Science department to host content pertaining to my Academic Life.
+As such, my research and other class documents will appear here.
+For the curious, this website was built with Pico CMS using the BitsAndPieces theme.
+Probabilities must be between zero and one, i.e., $0≤P(A)≤1$ for any event A.
+Probabilities add to one, i.e., $\sum{P(X_i)} = 1$
+The complement of an event, $A^c$, denotes that the event did not happen. Since probabilities must add to one, $P(A^c) = 1 - P(A)$
+If A and B are two events, the probability that A or B happens (this is an inclusive or) is the probability of the union of the events: +$$ +P(A \cup B) = P(A) + P(B) - P(A\cap B) +$$ +where $\cup$ represents union ("or") and $\cap$ represents intersection ("and"). If a set of events $Ai$ are mutually exclusive (only one event may happen), then +$$ +P(\cup{i=1}^n{Ai}) = \sum{i=1}^n{P(A_i)} +$$
+The odds for event A, denoted $\mathcal{O}(A)$ is defined as $\mathcal{O}(A) = P(A)/P(A^c)$
+This is the probability for divided by probability against the event
+From odds, we can also compute back probabilities +$$ +\frac{P(A)}{P(A^c)} = \mathcal{O}(A) +$$
+$$ +\frac{P(A)}{1-P(A)} = \mathcal{O}(A) +$$
+$$ +\frac{1 -P(A)}{P(A)} = \frac{1}{\mathcal{O}(A)} +$$
+$$ +\frac{1}{P(A)} - 1 = \frac{1}{\mathcal{O}(A)} +$$
+$$ +\frac{1}{P(A)} = \frac{1}{\mathcal{O}(A)} + 1 +$$
+$$ +\frac{1}{P(A)} = \frac{1 + \mathcal{O}(A)}{\mathcal{O}(A)} +$$
+$$ +P(A) = \frac{\mathcal{O}(A)}{1 + \mathcal{O}(A)} +$$
+The expected value of a random variable X is a weighted average of values X can take, with weights given by the probabilities of those values. +$$ +E(X) = \sum_{i=1}^n{x_i * P(X=x_i)} +$$
+Classical -- Outcomes that are equally likely have equal probabilities
+Frequentist -- In an infinite sequence of events, what is the relative frequency
+Bayesian -- Personal perspective (your own measure of uncertainty)
+In betting, one must make sure that all the rules of probability are followed. That the events are "coherent", otherwise one might construct a series of bets where you're guaranteed to lose money. This is referred to as a Dutch book.
+$$ +P(A|B) = \frac{P(A\cup B)}{P(B)} +$$
+Where $A|B$ denotes "A given B"
+Example from lecture:
+Suppose there are 30 students, 9 of which are female. From the 30 students, 12 are computer science majors. 4 of those 12 computer science majors are female +$$ +P(Female) = \frac{9}{30} = \frac{3}{10} +$$
+$$ +P(CS) = \frac{12}{30} = \frac{2}{5} +$$
+$$ +P(F\cap CS) = \frac{4}{30} = \frac{2}{15} +$$
+$$ +P(F|CS) = \frac{P(F \cap CS)}{P(CS)} = \frac{2/15}{2/5} = \frac{1}{3} +$$
+An intuitive way to think about a conditional probability is that we're looking at a subsegment of the original population, and asking a probability question within that segment +$$ +P(F|CS^c) = \frac{P(F\cap CS^c)}{PS(CS^c)} = \frac{5/30}{18/30} = \frac{5}{18} +$$ +The concept of independence is when one event does not depend on another. +$$ +P(A|B) = P(A) +$$ +It doesn't matter that B occurred.
+If two events are independent then the following is true +$$ +P(A\cap B) = P(A)P(B) +$$ +This can be derived from the conditional probability equation.
+Suppose we don't know what $P(A|B)$ is but we do know what $P(B|A)$ is. We can then rewrite $P(A|B)$ in terms of $P(B|A)$ +$$ +P(A|B) = \frac{P(B|A)P(A)}{P(B|A)P(A) + P(B|A^c)P(A^c)} +$$ +Let's look at an example of an early test for HIV antibodies known as the ELISA test. +$$ +P(+ | HIV) = 0.977 +$$
+$$ +P(- | NO_HIV) = 0.926 +$$
+As you can see over 90% of the time, this test was accurate.
+The probability of someone in North America having this disease was $P(HIV) = .0026$
+Now let's consider the following problem: the probability of having the disease given that they tested positive $P(HIV | +)$ +$$ +P(HIV|+) = \frac{P(+|HIV)P(HIV)}{P(+|HIV)P(HIV) + P(+|NO_HIV){P(NO_HIV)}} +$$
+$$ +P(HIV|+) = \frac{(.977)(.0026)}{(.977)(.0026) + (1-.977)(1-.0026)} +$$
+$$ +P(HIV|+) = 0.033 +$$
+This example looked at Bayes Theorem for the two event case. We can generalize it to n events through the following formula +$$ +P(A|B) = \frac{P(B|A_1){(A1)}}{\sum{i=1}^{n}{P(B|A_i)}P(A_i)} +$$
+~ means 'is distributed as'
+We'll be first studying the Bernoulli Distribution. This is when your event has two outcomes, which is commonly referred to as a success outcome and a failure outcome. The probability of success is $p$ which means the probability of failure is $(1-p)$ +$$ +X \sim B(p) +$$
+$$ +P(X = 1) = p +$$
+$$ +P(X = 0) = 1-p +$$
+The probability of a random variable $X$ taking some value $x$ given $p$ is +$$ +f(X = x | p) = f(x|p) = p^x(1-p)^{1 - x}I +$$ +Where $I$ is the Heavenside function
+Recall the expected value +$$ +E(X) = \sum_{x_i}{x_iP(X=x_i)} = (1)p + (0)(1-p) = p +$$ +We can also define the variance of Bernoulli +$$ +Var(X) = p(1-p) +$$
+The binomial distribution is the sum of n independent Bernoulli trials +$$ +X \sim Bin(n, p) +$$
+$$ +P(X=x|p) = f(x|p) = {n \choose x} p^x (1-p)^{n-x} +$$
+$n\choose x$ is the combinatoric term which is defined as +$$ +{n \choose x} = \frac{n!}{x! (n - x)!} +$$
+$$ +E(X) = np +$$
+$$ +Var(X) = np(1-p) +$$
+Let's say X is uniformally distributed +$$ +X \sim U[0,1] +$$
+$$ +f(x) = \left{ +\begin{array}{lr} +1 & : x \in [0,1]\ +0 & : otherwise +\end{array} +\right. +$$
+$$ +P(0 < x < \frac{1}{2}) = \int_0^\frac{1}{2}{f(x)dx} = \int_0^\frac{1}{2}{dx} = \frac{1}{2} +$$
+$$ +P(0 \leq x \leq \frac{1}{2}) = \int_0^\frac{1}{2}{f(x)dx} = \int_0^\frac{1}{2}{dx} = \frac{1}{2} +$$
+$$ +P(x = \frac{1}{2}) = 0 +$$
+$$ +\int_{-\infty}^\infty{f(x)dx} = 1 +$$
+$$ +f(x) \ge 0 +$$
+$$ +E(X) = \int_{-\infty}^\infty{xf(x)dx} +$$
+$$ +E(g(X)) = \int{g(x)f(x)dx} +$$
+$$ +E(aX) = aE(X) +$$
+$$ +E(X + Y) = E(X) + E(Y) +$$
+If X & Y are independent +$$ +E(XY) = E(X)E(Y) +$$
+$$ +X \sim Exp(\lambda) +$$
+Where $\lambda$ is the average unit between observations +$$ +f(x|\lambda) = \lambda e^{-\lambda x} +$$
+$$ +E(X) = \frac{1}{\lambda} +$$
+$$ +Var(X) = \frac{1}{\lambda^2} +$$
+$$ +X \sim [\theta_1, \theta_2] +$$
+$$ +f(x|\theta_1,\theta_2) = \frac{1}{\theta_2 - \theta1}I{\theta_1 \le x \le \theta_2} +$$
+$$ +X \sim N(\mu, \sigma^2) +$$
+$$ +f(x|\mu,\sigma^2) = \frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{1}{2\sigma^2}(x-\mu)^2} +$$
+$$ +E(X) = \mu +$$
+$$ +Var(X) = \sigma^2 +$$
+Variance is the squared distance from the mean +$$ +Var(X) = \int_{-\infty}^\infty {(x - \mu)^2f(x)dx} +$$
+The geometric distribution is the number of trails needed to get the first success, i.e, the number of Bernoulli events until a success is observed. +$$ +X \sim Geo(p) +$$
+$$ +P(X = x|p) = p(1-p)^{x-1} +$$
+$$ +E(X) = \frac{1}{p} +$$
+Multinomial is like a binomial when there are more than two possible outcomes.
+$$ +f(x_1,...,x_k|p_1,...,p_k) = \frac{n!}{x_1! ... x_k!}p_1^{x_1}...p_k^{x_k} +$$
+The Poisson distribution is used for counts. The parameter $\lambda > 0$ is the rate at which we expect to observe the thing we are counting. +$$ +X \sim Pois(\lambda) +$$
+$$ +P(X=x|\lambda) = \frac{\lambda^xe^{-\lambda}}{x!} +$$
+$$ +E(X) = \lambda +$$
+$$ +Var(X) = \lambda +$$
+If $X_1, X_2, ..., X_n$ are independent and identically distributed Exponentials,waiting time between success events, then the total waiting time for all $n$ events to occur will follow a gamma distribution with shape parameter $\alpha = n$ and rate parameter $\beta = \lambda$ +$$ +Y \sim Gamma(\alpha, \beta) +$$
+$$ +f(y|\alpha,\beta) = \frac{\beta^n}{\Gamma(\alpha)}y^{n-1}e^{-\beta y}I_{y\ge0}(y) +$$
+$$ +E(Y) = \frac{\alpha}{\beta} +$$
+$$ +Var(Y) = \frac{\alpha}{\beta^2} +$$
+Where $\Gamma(x)$ is the gamma function. The exponential distribution is a special case of the gamma distribution with $\alpha = 1$. As $\alpha$ increases, the gamma distribution more closely resembles the normal distribution.
+The beta distribution is used for random variables which take on values between 0 and 1. For this reason, the beta distribution is commonly used to model probabilities. +$$ +X \sim Beta(\alpha, \beta) +$$
+$$ +f(x|\alpha,\beta) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)}x^{n -1}(1 - x)^{\beta - 1}I_{{0 < x < 1}} +$$
+$$ +E(X) = \frac{\alpha}{\alpha + \beta} +$$
+$$ +Var(X) = \frac{\alpha\beta}{(\alpha + \beta)^2(\alpha+\beta+1)} +$$
+The standard uniform distribution is a special case of the beta distribution with $\alpha = \beta = 1$
+$$ +f(\theta|y) = \frac{f(y|\theta)f(\theta)}{\int{f(y|\theta)f(\theta)d\theta}} +$$
+Under the frequentest paradigm, you view the data as a random sample from some larger, potentially hypothetical population. We can then make probability statements i.e, long-run frequency statements based on this larger population.
+Let's suppose we flip a coin 100 times and we get 44 heads and 56 tails. +$$ +n = 100 +$$ +We can view these 100 flips as a random sample from a much larger infinite hypothetical population of flips from this coin.
+Let's say that each flip follows a Bournelli distribution with some probability p. In this case $p$ is unknown, but we're assuming it's fixed because we have a particular physical coin.
+We can ask what's our best estimate of the probability of getting a head, or an estimate of $p$. We can also ask about how confident we are about that estimate.
+Let's start by applying the Central Limit Theorem. The Central Limit Theorem states that the sum of 100 flips will follow approximately a Normal distribution with mean 100p and variance 100p(1-p) +$$ +\sum^n_{i=1}{x_i} \sim N(np, np(1-p)) +$$
+$$ +\sum_{i = 1}^{100}{x_i} \sim N(100p, 100p(1-p)) +$$
+By the properties of a Normal distribution, 95% of the time we'll get a result within 1.96 standard deviations of the mean. Our estimate is $100\hat{p}$ and our error is 1.96 times the standard deviation. +$$ +n\hat{p} \pm 1.96\sqrt{n\hat{p}(1-\hat{p})} +$$
+$$ +100\hat{p} \pm 1.96\sqrt{100\hat{p}(1-\hat{p})} +$$
+This is referred to as a Confidence Interval. Confidence Intervals are commonly abbreviated as CI. In our example $\hat{p} = \frac{44}{n} = \frac{44}{100}$. Therefore, the 95% Confidence Interval in the true number of heads after flipping a coin 100 times is: +$$ +100(.44) \pm 1.96\sqrt{100(.44)(1-.44)} +$$
+$$ +44 \pm 1.96\sqrt{44(.56)} +$$
+$$ +44\pm 1.96\sqrt{24.64} +$$
+$$ +(34.27, 53.73) +$$
+We can divide this by 100 to get the 95% Confidence Interval for $p$ +$$ +(0.34, 0.53) +$$ +Let's step back and ask, what does it mean when I say we're 95% confident?
+Under the frequentest paradigm, what this means is we have to think back to our infinite hypothetical sequence of events. So if we were to repeat this trial an infinite number of times, or an arbitrary large number of times. Each time we create a confidence interval based on the data we observe, than on average 95% of the intervals we make will contain the true value of p.
+On the other hand, we might want to know something about this particular interval. Does this interval contain the true p? What's the probability that this interval contains a true p? Well, we don't know for this particular interval. But under the frequentest paradigm, we're assuming that there is a fixed answer for p. Either p is in that interval or it's not in that interval. The probability that p is in that interval is either 0 or 1.
+Consider a hospital where 400 patients are admitted over a month for heart attacks, and a month later 72 of them have died and 328 of them have survived.
+We can ask, what's our estimate of the mortality rate?
+Under the frequentest paradigm, we must first establish our reference population. What do we think our reference population is here? One possibility is we could think about heart attack patients in the region.
+Another reference population we can think about is heart attack patients that are admitted to this hospital, but over a longer period of time.
+Both of these might be reasonable attempts, but in this case our actual data are not random sample from either of those populations. We could sort of pretend they are and move on, or we could also try to think harder about what a random sample situation might be. We can think about all the people in the region who might possibly have a heart attack and might possibly get admitted to this hospital.
+It's a bit of an odd hypothetical situation, and so there are some philosophical issues with the setup of this whole problem with the frequentest paradigm, In any case, let's forge forward and think about how we might do some estimation.
+Moving on, we can say each patient comes from a Bernoulli distribution with an unknown parameter $\theta$. +$$ +Y_i \sim B(\theta) +$$
+$$ +P(Y_i = 1) = \theta +$$
+In this case, let's call the "success" a mortality.
+The probability density function for the entire set of data we can write in vector form. Probability of all the Y's take some value little y given a value of theta. +$$ +P(Y = y | \theta) = P(Y_1 = y_1, Y_2, = y_2,\dots, Y_n=y_n|\theta) +$$ +Since we're viewing these as independent events, then the probability of each of these individual ones we can write in product notation. +$$ +P(Y = y | \theta) = P(Y_1 = y_1|\theta)\dots P(Y_n = y_n | \theta) +$$
+$$ +P(Y = y | \theta) = \prod_{i = 1}^n{P(Y_i =yi | \theta)} = \prod{i = 1}^n{(\theta^{y_i}(1-\theta)^{1-y_i})} +$$
+This is the probability of observing the actual data that we collected, conditioned on a value of the parameter $\theta$. We can now think about this expression as a function of theta. This is a concept of a likelihood. +$$ +L(\theta|y) = \prod_{i = 1}^n{(\theta^{y_i}(1-\theta)^{1-y_i})} +$$ +It looks like the same function, but here it is a function of y given $\theta$. And now we're thinking of it as a function of $\theta$ given y.
+This is not a probability distribution anymore, but it is still a function for $\theta$.
+One we to estimate $\theta$ is that we choose the $\theta$ that gives us the largest value of the likelihood. It makes the data the most likely to occur for the particular data we observed.
+This is referred to as the maximum likelihood estimate (MLE),
+We're trying to find the $\theta$ that maximizes the likelihood.
+In practice, it's usually easier to maximize the natural logarithm of the likelihood, commonly referred to as the log likelihood. +$$ +\mathcal{L}(\theta) = \log{L(\theta|y)} +$$ +Since the logarithm is a monotonic function, if we maximize the logarithm of the function, we also maximize the original function. +$$ +\mathcal{L(\theta)} = \log{\prod_{i = 1}^n{(\theta^{y_i}(1-\theta)^{1-y_i})}} +$$
+$$ +\mathcal{L}(\theta) = \sum_{i = 1}^n{\log{(\theta^{y_i}(1-\theta)^{1-y_i})}} +$$
+$$ +\mathcal{L}(\theta) = \sum_{i = 1}^n{(\log{(\theta^{y_i}}) + \log{(1-\theta)^{1-y_i}})} +$$
+$$ +\mathcal{L}(\theta) = \sum_{i = 1}^n{(y_i\log{\theta} + (1 - y_i)\log{(1-\theta)})} +$$
+$$ +\mathcal{L}(\theta) = \log{\theta}\sum_{i = 1}^n{yi} + \log{(1-\theta)}\sum{i = 1}^n{(1-y_i)} +$$
+How do we find the theta that maximizes this function? Recall from calculus that we can maximize a function by taking the derivative and setting it equal to 0. +$$ +\mathcal{L}^\prime(\theta) = \frac{1}{\theta}\sum{y_i} - \frac{1}{1-\theta}\sum{(1 - y_i)} +$$ +Now we need to set it equal to zero and solve for $\hat{\theta}$ +$$ +\frac{\sum{y_i}}{\hat{\theta}} = \frac{\sum{(1-y_i)}}{1-\hat{\theta}} +$$
+$$ +\hat{\theta}\sum{(1-y_i)} = (1-\hat{\theta})\sum{y_i} +$$
+$$ +\hat{\theta}\sum{(1-y_i)} + \hat{\theta}\sum{y_i} = \sum{y_i} +$$
+$$ +\hat{\theta}(\sum^n{(1 - y_i + y_i)}) = \sum{y_i} +$$
+$$ +\hat{\theta} = \frac{1}{n}\sum{y_i} = \hat{p} = \frac{72}{400} = 0.18 +$$
+Maximum likelihood estimators have many desirable mathematical properties. They're unbiased, they're consistent, and they're invariant.
+In general, under certain regularity conditions, we can say that the MLE is approximately normally distributed with mean at the true value of $\theta$ and the variance $\frac{1}{I(\hat{\theta})}$ where $I(\hat{\theta})$ is the Fisher information at the value of $\hat{\theta}$. The Fisher information is a measure of how much information about $\theta$ is in each data point. +$$ +\hat{\theta} \sim N(\theta, \frac{1}{I(\hat{\theta})}) +$$ +For a Bernoulli random variable, the Fisher information turns out to be +$$ +I(\theta) = \frac{1}{\theta(1-\theta)} +$$ +So the information is larger, when theta is near zero or near one, and it's the smallest when theta is near one half.
+This makes sense, because if you're flipping a coin, and you're getting a mix of heads and tails, that tells you a little bit less than if you're getting nearly all heads or nearly all tails.
+Let's say $X_i$ are distributed so +$$ +Xi \sim Exp(\lambda) +$$ +Let's say the data is independent and identically distributed, therefore making the overall density function +$$ +f(x|\lambda) = \prod{i = 1}^n{\lambda e^{-\lambda x_i}} +$$
+$$ +f(x|\lambda) = \lambda^ne^{-\lambda \sum{x_i}} +$$
+Now the likelihood function is +$$ +L(\lambda|x) = \lambda^n e^{-\lambda \sum{x_i}} +$$
+$$ +\mathcal{L}(\lambda) = n\ln{\lambda} - \lambda\sum{x_i} +$$
+Taking the derivative +$$ +\mathcal{L}^\prime(\lambda) = \frac{n}{\lambda} - \sum{x_i} +$$ +Setting this equal to zero +$$ +\frac{n}{\hat{\lambda}} =\sum{x_i} +$$
+$$ +\hat{\lambda} = \frac{n}{\sum{x_i}} = \frac{1}{\bar{x}} +$$
+$$ +X_i \sim U[0, \theta] +$$
+$$ +f(x|\theta) = \prod{i = 1}^n{\frac{1}{\theta}I{0 \le x_i \le \theta}} +$$
+Combining all the indicator functions, for this to be a 1, each of these has to be true. These are going to be true if all the observations are bigger than 0, as in the minimum of the x is bigger than or equal to 0. The maximum of the x's is also less than or equal to $\theta$. +$$ +L(\theta|x) = \theta^{-1} I_{0\le min(x_i) \le max(x_i) \le \theta} +$$
+$$ +L^\prime(\theta) = -n\theta^{-(n + 1)}I_{0 \le min(x_i) \le max(x_i)\le \theta} +$$
+So now we can ask, can we set this equal to zero and solve for $\theta$? Well it turns out, this is not equal to zero for any $\theta$ positive value. We need $\theta$ to be strictly larger than zero.
+However, we can also note that for $\theta$ positive, this will always be negative. The derivative is negative, that says this is a decreasing function. So this funciton will be maximized when we pick $\theta$ as small as possible. What's the smallest possible value of $\theta$ we can pick? Well we need in particular for $\theta$ to be larger than all of the $X_i$ . And so, the maximum likelihood estimate is the maximum of $X_i$ +$$ +\hat{\theta} = max(x_i) +$$
+Because 0 * 1 = 0, the product of indicator functions can be combined into a single indicator function with a modified condition.
+Example: $I{x < 5} * I{x \ge 0} = I_{0 \le x < 5}$
+Example: $\prod{i = 1}^n{I{xi < 2}} = I{xi < 2 for all i} = I{max(x_i) < 2}$
+R has some nice functions that one can use for analysis
+mean(z) gives the mean of some row vector $z$
var(z) reports the variance of some row vector
sqrt(var(z)) gives the standard deviation of some row vector
seq(from=0.1, to = 0.9, by = 0.1) creates a vector that starts from $0.1$ and goes to $0.9$ incrementing by $0.1$
seq(from=0.1, to = 0.9, by = 0.1)
+[1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
+seq(1, 10)
+[1] 1 2 3 4 5 6 7 8 9 10names(x) gives the names of all the columns in the dataset.
names(trees)
+[1] "Girth" "Height" "Volume"hist(x) provides a histogram based on a vector
The more general plot function tries to guess at which type of plot to make. Feeding it two numerical vectors will make a scatter plot.
The R function pairs takes in a data frame and tries to make all possible Pairwise scatterplots for the dataset.
The summary command gives the five/six number summary (minimum, first quartile, median, mean, third quartile, maximum)
Going back to the hospital example
+## Likelihood function
+likelihood = function(n, y, theta) {
+ return(theta^y * (1 - theta)^(n - y))
+}
+theta = seq(from = 0.01, to = 0.99, by = 0.01)
+plot(theta, likelihood(400, 72, theta))You can also do this with log likelihoods. This is typically more numerically stable to compute
+loglike = function(n, y, theta) {
+ return(y * log(theta) + (n - y) * log(1 - theta))
+}
+plot(theta, loglike(400, 72, theta))Having these plotted as points makes it difficult to see, let's plot it as lines
+plot(theta, loglike(400, 72, theta), type = "l")The cumulative distribution function (CDF) exists for every distribution. We define it as $F(x) = P(X \le x)$ for random variable $X$.
+If $X$ is discrete-valued, then the CDF is computed with summation $F(x) = \sum_{t = -\infty}^x {f(t)}$. where $f(t) = P(X = t)$ is the probability mass function (PMF) which we've already seen.
+If $X$ is continuous, the CDF is computed with an integral $F(x) = \int_{-\infty}^x{f(t)dt}$
+The CDF is convenient for calculating probabilities of intervals. Let $a$ and $b$ be any real numbers with $a < b$. Then the probability that $X$ falls between $a$ and $b$ is equal to $P(a < X < b) = P(X \le b) - P(X \le a) = F(b) - F(a)$
+The CDF takes a value for a random variable and returns a probability. Suppose instead we start with a number between $0$ and $1$, which we call $p$, and we wish to find a value $x$ so that $P(X \le x) = p$. The value $x$ which satisfies this equation is called the $p$ quantile. (or $100p$ percentile) of the distribution of $X$.
+Each of the distributions introduced in Lesson 3 have convenient functions in R which allow you to evaluate the PDF/PMF, CDF, and quantile functions, as well as generate random samples from the distribution. To illustrate, Table I list these functions for the normal distribution
+| Function | +What it does | +
|---|---|
dnorm(x, mean, sd) |
+Evaluate the PDF at $x$ (mean = $\mu$ and sd = $\sqrt{\sigma^2}$) | +
pnorm(q, mean, sd) |
+Evaluate the CDF at $q$ | +
qnorm(p, mean, sd) |
+Evaluate the quantile function at $p$ | +
rnorm(n, mean, sd) |
+Generate $n$ pseudo-random samples from the normal distribution | +
These four functions exist for each distribution where d... function evaluates the density/mass, p... evaluates the CDF, q... evaluates the quantile, and r... generates a sample. Table 2 lists the d... functions for some of the most popular distributions. The d can be replaced with p, q, or r for any of the distributions, depending on what you want to calculate.
For details enter ?dnorm to view R's documentation page for the Normal distribution. As usual , replace the norm with any distribution to read the documentation for that distribution.
| Distribution | +Function | +Parameters | +
|---|---|---|
| $Binomial(n,p)$ | +dbinom(x, size, prob) |
+size = $n$, prob = $p$ | +
| $Poisson(\lambda)$ | +dpois(x, lambda) |
+lambda = $\lambda$ | +
| $Exp(\lambda)$ | +dexp(x, rate) |
+rate = $\lambda$ | +
| $Gamma(\alpha, \beta)$ | +dgamma(x, shape, rate) |
+shape = $\alpha$, rate = $\beta$ | +
| $Uniform(a, b)$ | +dunif(x, min, max) |
+min = $a$, max = $b$ | +
| $Beta(\alpha, \beta)$ | +dbeta(x, shape1, shape2) |
+shape1 = $\alpha$, shape2 = $\beta$ | +
| $N(\mu, \sigma^2)$ | +dnorm(x, mean, sd) |
+mean = $\mu$, sd = $\sqrt{\sigma^2}$ | +
| $t_v$ | +dt(x, df) |
+df = $v$ | +
Suppose your brother has a coin which you know to be loaded so that it comes up heads 70% of the time. He then comes to you with some coin, you're not sure which one and he wants to make a bet with you. Betting money that it's going to come up heads.
+You're not sure if it's the loaded coin or if it's just a fair one. So he gives you a chance to flip it 5 times to check it out.
+You flip it five times and get 2 heads and 3 tails. Which coin do you think it is and how sure are you about that?
+We'll start by defining the unknown parameter $\theta$, this is either that the coin is fair or it's a loaded coin. +$$ +\theta = {fair ,loaded} +$$
+$$ +X \sim Bin(5, ?) +$$
+$$ +f(x|\theta) = \begin{cases} +{5 \choose x}(\frac{1}{2})^5 & \theta = fair \ +{5 \choose x} (.7)^x (.3)^{5 - x} & \theta = loaded\ +\end{cases} +$$
+We can also rewrite $f(x|\theta)$ with indicator functions +$$ +f(x|\theta) = {5\choose x}(.5)^5I{{\theta= fair}} + {5 \choose x}(.7)^x(.3)^{5 - x}I{{\theta = loaded}} +$$ +In this case, we observed that $x = 2$ +$$ +f(\theta | x = 2) = \begin{cases} +0.3125 & \theta = fair \ +0.1323 & \theta = loaded +\end{cases} +$$ +MLE $\hat{\theta} = fair$
+That's a good point estimate, but then how do we answer the question, how sure are you?
+This is not a question that's easily answered in the frequentest paradigm. Another question is that we might like to know what is the probability that theta equals fair, give, we observe two heads. +$$ +P(\theta = fair|x = 2) = ? +$$ +In the frequentest paradigm, the coin is a physical quantity. It's a fixed coin, and therefore it has a fixed probability of coining up heads. It is either the fair coin, or it's the loaded coin. +$$ +P(\theta = fair) = {0,1} +$$
+An advantage of the Bayesian approach is that it allows you to easily incorporate prior information, when you know something in advance of the looking at the data. This is difficult to do under the Frequentest paradigm.
+In this case, we're talking about your brother. You probably know him pretty well. So suppose you think that before you've looked at the coin, there's a 60% probability that this is the loaded coin.
+This case, we put this into our prior. Our prior is that the probability the coin is loaded is 0.6. We can update our prior with the data to get our posterior beliefs, and we can do this using the bayes theorem.
+Prior: $P(loaded) = 0.6$ +$$ +f(\theta|x) = \frac{f(x|\theta)f(\theta)}{\sum_\theta{f(x|\theta)f(\theta)}} +$$
+$$ +f(\theta|x) = \frac{{5\choose x} [(\frac{1}{2})^5(.4)I{{\theta = fair}} + (.7)^x (.3)^{5-x}(.6)I{{\theta = loaded}} ] } +{{5\choose x} [(\frac{1}{2})^5(.4) + (.7)^x (.3)^{5-x}(0.6) ] } +$$
+$$ +f(\theta|x=2)= \frac{0.0125I{{\theta=fair}} + 0.0079I{{\theta=loaded}} }{0.0125+0.0079} +$$
+$$ +f(\theta|x=2) = 0.612I{{\theta=fair}} + 0.388I{{\theta = loaded}} +$$
+As you can see in the calculation here, we have the likelihood times the prior in the numerator, and in the denominator, we have a normalizing constant, so that when we divide by this, we'll get answer that add up to one. These numbers match exactly in this case, because it's a very simple problem. But this is a concept that goes on, what's in the denominator here is always a normalizing constant. +$$ +P(\theta = loaded | x = 2) = 0.388 +$$ +This here updates our beliefs after seeing some data about what the probability might be.
+We can also examine what would happen under different choices of prior. +$$ +P(\theta = loaded) = \frac{1}{2} \implies P(\theta = loaded | x = 2) = 0.297 +$$
+$$ +P(\theta = loaded) = 0.9 \implies P(\theta = loaded | x = 2) = 0.792 +$$
+In this case, the Bayesian approach is inherently subjective. It represents your own personal perspective, and this is an important part of the paradigm. If you have a different perspective, you will get different answers, and that's okay. It's all done in a mathematically vigorous framework, and it's all mathematically consistent and coherent.
+And in the end, we get results that are interpretable
+$$ +f(\theta | y) = \frac{f(y | \theta)f(\theta)}{f(y)} = \frac{f(y|\theta)f(\theta)}{\int{f(y|\theta)f(\theta)d\theta}} = \frac{likelihood prior}{normalization} \propto likelihood prior +$$
+In practice, sometimes this integral can be a pain to compute. And so, we may work with looking at saying this is proportional to the likelihood times the prior. And if we can figure out what this looks like and just put the appropriate normalizing constant on at the end, we don't necessarily have to compute this integral.
+So for example, suppose we're looking at a coin and it has unknown probability $\theta$ of coming up heads. Suppose we express ignorance about the value of $\theta$ by assigning it a uniform distribution. +$$ +\theta \sim U[0, 1] +$$
+$$ +f(\theta) = I_{{0 \le \theta\le 1}} +$$
+$$ +f(\theta | y = 1) = \frac{\theta^1(1-\theta)^0I{{0 \le \theta\le1}}}{\int{-\infty}^\infty{\theta^1(1-\theta)^0I_{{0\le \theta \le 1}}}} +$$
+$$ +f(\theta | y = 1) = \frac{\theta I_{{0\le\theta\le1}}}{\int0^1{\theta d\theta}} = 2\theta I{{0\le\theta\le1}} +$$
+Now if we didn't want to take the integral we could've done this approach +$$ +f(\theta | y) \propto f(y|\theta)f(\theta) \propto \theta I_{{0\le\theta\le1}} +$$ +Which then we need to find the constant such that it's a proper PMF. In this case, it's $2$.
+Since it's a proper PMF, we can perform interval probabilities as well. This is called Posterior interval estimates. +$$ +P(0.025 < \theta < 0.975) = \int_{0.025}^{0.975}{2\theta d \theta} = (0.975)^2 - (0.025)^2 = 0.95 +$$
+$$ +P(\theta > 0.05) = 1 - (0.05)^2 = 0.9975 +$$
+These are the sort of intervals we would get from the prior and asking about what their posterior probability is.
+In other cases, we may want to ask, what is the posterior interval of interest? What's an interval that contains 95% of posterior probability in some meaningful way? This would be equivalent then to a frequentest confidence interval. We can do this in several different ways, 2 main ways that we make Bayesian Posterior intervals or credible intervals are equal-tailed intervals and highest posterior density intervals.
+In the case of an equal-tailed interval, we put the equal amount of probability in each tail. So to make a 95% interval we'll put 0.025 in each tail.
+To be able to do this, we're going to have to figure out what the quantiles are. So we're going to need some value, $q$, so that +$$ +P(\theta < q | Y = 1) = \int_0^9{2\theta d\theta} = q^2 +$$
+$$ +P(\sqrt{0.025} < \theta < \sqrt{0.975}) = P(0.158 < \theta < 0.987) = 0.95 +$$
+This is an equal tailed interval in that the probability that $\theta$ is less than 0.18 is the same as the probability that $\theta$ is greater than 0.987. We can say that under the posterior, there's a 95% probability that $\theta$ is in this interval.
+Here we want to ask where in the density function is it highest? Theoretically this will be the shortest possible interval that contains the given probability, in this case a 95% probability. +$$ +P(\theta > \sqrt{0.05} | Y = 1) = P(\theta > 0.224 | Y = 1) = 0.95 +$$ +This is the shortest possible interval, that under the posterior has a probability 0.95. it's $\theta$ going from 0.224 up to 1.
+The posterior distribution describes our understanding of our uncertainty combinbing our prior beliefs and the data. It does this with a probability density function, so at the end of teh day, we can make intervals and talk about probabilities of data being in the interval.
+This is different from the frequentest approach, where we get confidence intervals. But we can't say a whole lot about the actual parameter relative to the confidence interval. We can only make long run frequency statements about hypothetical intervals.
+In this case, we can legitimately say that the posterior probability that $\theta$ is bigger than 0.05 is $0.9975$. We can also say that we believe there's a 95% probability that $\theta$ is in between 0.158 and 0.987.
+Bayesians represent uncertainty with probabilities, so that the coin itself is a physical quantity. It may have a particular value for $\theta$.
+It may be fixed, but because we don't know what that value is, we represent our uncertainty about that value with a distribution. And at the end of the day, we can represent our uncertainty, collect it with the data, and get a posterior distribution and make intuitive statements.
+ +Frequentest confidence intervals have the interpretation that "If you were to repeat many times the process of collecting data and computing a 95% confidence interval, then on average about 95% of those intervals would contain the true parameter value; however, once you observe data and compute an interval the true value is either in the interval or it is not, but you can't tell which."
+Bayesian credible intervals have the interpretation that "Your posterior probability that the parameter is in a 95% credible interval is 95%."
+How do we choose a prior?
+Our prior needs to represent our personal perspective, beliefs, and our uncertainties.
+Theoretically, we're defining a cumulative distribution function for the parameter +$$ +\begin{cases} +P(\theta \le c) & c \in \mathbb{R} +\end{cases} +$$ +This is true for an infinite number of possible sets. This isn't practical to do, and it would be very difficult to do coherently so that all the probabilities were consistent.
+In practice, we work with a convenient family that's sufficiently flexible such that a member of a family represents our beliefs.
+Generally if one has enough data, the information in the data will overwhealm the information in the prior. And so, the prior is not particularly important in terms of what you get for the posterior. Any reasonable choice of prior will lead to approximately the same posterior. However, there are some things that can go wrong.
+Suppose we chose a prior that says the probability of $P(\theta = \frac{1}{2}) = 1$
+And thus, the probability of $\theta$ equaling any other value is $0$. If we do this, our data won't make a difference since we only put a probability of $1$ at a single point. +$$ +f(\theta|y) \propto f(y|\theta)f(\theta) = f(\theta) = \delta(\theta) +$$
+In the basic context, events with prior probability of zero have a posterior probability of zero. Events with a prior probability of one, have a posterior probability of one.
+Thus a good Bayesian will not assign probability of zero or one to any event that has already occurred or already known not to occur.
+A useful concept in terms of choosing priors is that of the calibration of predictive intervals.
+If we make an interval where we're saying we predict 95% of new data points will occur in this interval. It would be good if in reality 95% of new data points actually did fall in that interval.
+How do we calibrate to reality? This is actually a frequentest concept but this is important for practical statistical purposes that our results reflect reality.
+We can compute a predictive interval, this is an interval such that 95% of new observations are expected to fall into it. It's an interval for the data rather than an interval for $\theta$ +$$ +f(y) = \int{f(y|\theta)f(\theta)d\theta} = \int{f(y, \theta)d\theta} +$$ +Where $f(y,\theta)$ is the joint density of Y and $\theta$.
+This is the prior predictive before any data is observed.
+Side Note: From this you can say that $f(y, \theta) = f(y|\theta)f(\theta)$
+Suppose we're going to flip a coin ten times and count the number of heads we see. We're thinking about this in advance of actually doing it, so we're interested in the predictive distribution. How many heads do we predict we're going to see? +$$ +X = \sum_{i = 1}^{10}{Y_i} +$$ +Where $Y_i$ is each individual coin flip.
+If we think that all possible coins or all possible probabilities are equally likely, then we can put a prior for $\theta$ that's flat over the interval from 0 to 1. +$$ +f(\theta) = I_{{0 \le \theta \le 1}} +$$
+$$ +f(x) = \int{f(x|\theta)f(\theta)d\theta} = \int_0^1{\frac{10!}{x!(10-x)!}\theta^x(1-\theta)^{10 -x}(1)d\theta} +$$
+Note that because we're interested in $X$ at the end, it's important that we distinguish between a binomial density and a Bernoulli density. Here we just care about the total count rather than the exact ordering which would be Bernoulli's
+For most of the analyses we're doing, where we're interested in $\theta$ rather than x, the binomial and the Bernoulli are interchangeable because the part in here that depends on $\theta$ is the same.
+To solve this integral let us recall some facts +$$ +n! = \Gamma(n + 1) +$$
+$$ +Z \sim Beta(\alpha, \beta) +$$
+$$ +f(z) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)}z^{\alpha - 1}(1-z)^{\beta - 1} +$$
+Let us rewrite $f(x)$ +$$ +f(x) = \int_0^1{\frac{\Gamma(11)}{\Gamma(x + 1)\Gamma(11 - x)}\theta^{(x + 1)-1}(1-\theta)^{(11-x)-1}d\theta} +$$
+$$ +f(x) = \frac{\Gamma(11)}{\Gamma(12)}\int_0^1{\frac{\Gamma(12)}{\Gamma(x + 1)\Gamma(11 - x)}\theta^{(x + 1)-1}(1-\theta)^{(11-x)-1}d\theta} +$$
+The integral above is a beta density, all integrals of valid beta densities equals to one. +$$ +f(x) = \frac{\Gamma(11)}{\Gamma(12)} = \frac{10!}{11!} = \frac{1}{11} +$$ +For $x \in {0, 1, 2, \dots, 10}$
+Thus we see that if we start with a uniform prior, we then end up with a discrete uniform predictive density for $X$. If all possible probabilities are equally likely, then all possible $X$ outcomes are equally likely.
+What about after we've observed data? What's our posterior predictive distribution?
+Going from the previous example, let us observe after one flip that we got a head. We want to ask, what's our predictive distribution for the second flip, given we saw a head on the first flip. +$$ +f(y_2|y_1) = \int{f(y_2|\theta,y_1)f(\theta|y_1)}d\theta +$$ +We're going to assume that $Y_2$ is independent of $Y_1$. Therefore, +$$ +f(y_2 |y_1) = \int{f(y_2|\theta)f(\theta|y_1)d\theta} +$$ +Suppose we're thinking of a uniform distribution for $\theta$ and we observe the first flip is a heads. What do we predict for the second flip?
+This is no longer going to be a uniform distribution like it was before, because we have some data. We're going to think it's more likely that we're going to get a second head. We think this because since we observed a head $\theta$ is now likely to be at least $\frac{1}{2}$ possibly larger. +$$ +f(y_2 | Y_1 = 1) = \int_0^1{\theta^{y_2}(1-\theta)^{1-y_2}2\theta d\theta} +$$
+$$ +f(y_2|Y_1 = 1) = \int_0^1{2\theta^{y_2 + 1}(1-\theta)^{1-y_2}d\theta} +$$
+We could work this out in a more general form, but in this case, $Y_2$ has to take the value $0$ or $1$. The next flip is either going to be heads or tails so it's easier to just plop in a particular example. +$$ +P(Y_2 = 1|Y_1 = 1) = \int_0^1{2\theta^2d\theta} = \frac{2}{3} +$$
+$$ +P(Y_2 = 0 | Y_1 = 1) = 1 - P(Y_2 = 1 | Y_1 = 1) = 1 - \frac{2}{3} = \frac{1}{3} +$$
+We can see here that the posterior is a combination of the information in the prior and the information in the data. In this case, our prior is like having two data points, one head and one tail.
+Saying we have a uniform prior for $\theta$ is equivalent in an information sense as saying we have observed one head and one tail.
+So then when we observe one head, it's like we now have seen two heads and one tail. So our predictive distribution for the second flip says if we have two heads and one tail, then we have a $\frac{2}{3}$ probability of getting another head and a $\frac{1}{3}$ probability of getting another tail.
+Likelihood of y given theta is +$$ +f(y|\theta) = \theta^{\sum{y_i}}(1-\theta)^{n - \sum{y_i}} +$$
+Our prior for theta is just a uniform distribution +$$ +f(\theta) = I_{{0 \le \theta \le 1}} +$$ +Thus our posterior for $\theta$ is +$$ +f(\theta | y) = \frac{f(y|\theta)f(\theta)}{\int{f(y|\theta)f(\theta)d\theta}} = \frac{\theta^{\sum{y_i}}(1-\theta)^{n - \sum{yi}} I{{0 \le \theta \le 1}}}{\int_0^1{\theta^{\sum{y_i}}(1-\theta)^{n - \sum{yi}} I{{0 \le \theta \le 1}}d\theta}} +$$ +Recalling the form of the beta distribution we can rewrite our posterior as +$$ +f(\theta | y) = \frac{\theta^{\sum{y_i}}(1-\theta)^{n - \sum{yi}} I{{0 \le \theta \le 1}}}{\frac{\Gamma(\sum{y_i} + 1)\Gamma(n - \sum{y_i} + 1)}{\Gamma(n + 2)}\int_0^1{\frac{\Gamma(n + 2)}{\Gamma(\sum{y_i} + 1)\Gamma(n - \sum{y_i} + 1)}\theta^{\sum{y_i}}(1-\theta)^{n - \sum{y_i}}d\theta}} +$$ +Since the beta density integrates to $1$, we can simplify this as +$$ +f(\theta | y) = \frac{\Gamma(n + 2)}{\Gamma(\sum{y_i}+ 1)\Gamma(n - \sum{y_i}+ 1)}\theta^{\sum{y_i}}(1-\theta)^{n-\sum{yi}}I{{0 \le \theta \le 1}} +$$ +From here we can see that the posterior follows a beta distribution +$$ +\theta | y \sim Beta(\sum{y_i} + 1, n - \sum{y_i} + 1) +$$
+The uniform distribution is $Beta(1, 1)$
+Any beta distribution is conjugate for the Bernoulli distribution. Any beta prior will give a beta posterior. +$$ +f(\theta) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha - 1}(1-\theta)^{\beta -1}I_{{0 \le \theta \le 1}} +$$
+$$ +f(\theta | y) \propto f(y|\theta)f(\theta) = \theta^{\sum{y_i}}(1-\theta)^{n - \sum{yi}}\frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha - 1}(1 - \theta)^{\beta - 1}I{{0 \le \theta \le 1}} +$$
+$$ +f(y|\theta)f(\theta) \propto \theta^{\alpha + \sum{y_i}-1}(1-\theta)^{\beta + n - \sum{y_i} - 1} +$$
+Thus we see that this is a beta distribution +$$ +\theta | y \sim Beta(\alpha + \sum{y_i}, \beta + n - \sum{y_i}) +$$ +When $\alpha$ and $\beta$ is one like in the uniform distribution, then we get the same result as earlier.
+This whole process where we choose a particular form of prior that works with a likelihood is called using a conjugate family.
+A family of distributions is referred to as conjugate if when you use a member of that family as a prior, you get another member of that family as your posterior.
+The beta distribution is conjugate for the Bernoulli distribution. It's also conjugate for the binomial distribution. The only difference in the binomial likelihood is that there is a combinatoric term. Since that does not depend on $\theta$, we get the same posterior.
+We often use conjugate priors because they make life much more simpler, sticking to conjugate families allows us to get closed form solutions easily.
+If the family is flexible enough, then you can find a member of that family that closely represents your beliefs.
+Returning to the beta posterior model it is clear how both the prior and the data contribute to the posterior.
+We can say that the effect size of the prior is $\alpha + \beta$
+Recall that the expected value or mean of a beta distribution is $\frac{\alpha}{\alpha + \beta}$
+Therefore we can derive the posterior mean as +$$ +posterior_{mean} = \frac{\alpha + \sum{y_i}}{\alpha + \sum{y_i}+\beta + n - \sum{y_i}}= \frac{\alpha+\sum{yi}}{\alpha + \beta + n} +$$ +We can further decompose this as +$$ +posterior{mean} = \frac{\alpha + \beta}{\alpha + \beta + n}\frac{\alpha}{\alpha + \beta} + \frac{n}{\alpha + \beta + n}\frac{\sum{y_i}}{n} +$$ +We can describe this as the (prior weight prior mean) + (data weight data mean)
+The posterior mean is a weighted average of the prior mean and the data mean.
+This effective sample size gives you an idea of how much data you would need to make sure that your prior doesn't have much influence on your posterior.
+If $\alpha + \beta$ is small compared to $n$ then the posterior will largely just be driven by the data. If $\alpha + \beta$ is large relative to $n$ then the posterior will be largely driven by the prior.
+We can make a 95% credible interval using our posterior distribution for $\theta$. We can find an interval that actually has 95% probability of containing $\theta$.
+Using Bayesian Statistics we can chain together dong a sequential update every time we get new data. We can get a new posterior, and we just use our previous posterior as a prior to do another update using Baye's theorem.
+Suppose we're giving two students a multiple-choice exam with 40 questions, where each question has four choices. We don't know how much the students have studied for this exam, but we think that they'll do better than just guessing randomly
+1) What are the parameters of interest?
+The parameters of interests are $\theta_1 = true$ the probability that the first student will answer a question correctly, $\theta_2 = true$ the probability that the second student will answer a question correctly.
+2) What is our likelihood?
+The likelihood is $Binomial(40, \theta)$, if we assume that each question is independent and that the probability a student gets each question right is the same for all questions for that student.
+3) What prior should we use?
+The conjugate prior is a beta prior. We can plot the density with dbeta
theta = seq(from = 0, to = 1, by = 0.1)
+# Uniform
+plot(theta, dbeta(theta, 1, 1), type = 'l')
+# Prior mean 2/3
+plot(theta, dbeta(theta, 4, 2), type = 'l')
+# Prior mean 2/3 but higher effect size (more concentrated at mean)
+plot(theta, dbeta(theta, 8, 4), type = 'l')4 ) What are the prior probabilities $P(\theta > 0.25)$? $P(\theta > 0.5)$? $P(\theta > 0.8)$?
+1 - pbeta(0.25, 8, 4)
+#[1] 0.998117
+1 - pbeta(0.5, 8, 4)
+#[1] 0.8867188
+1 - pbeta(0.8, 8, 4)
+#[1] 0.161133925) Suppose the first student gets 33 questions right. What is the posterior distribution for $\theta_1$? $P(\theta > 0.25)$? $P(\theta > 0.5)$? $P(\theta > 0.8)$? What is the 95% posterior credible interval for $\theta_1$? +$$ +Posterior \sim Beta(8 + 33, 4 + 40 - 33) = Beta(41, 11) +$$ +With a posterior mean of $\frac{41}{41+11} = \frac{41}{52}$
+We can plot the posterior distribution with the prior
+plot(theta, dbeta(theta, 41, 11), type = 'l')
+lines(theta, dbeta(theta, 8 ,4), lty = 2) #Dashed line for priorPosterior probabilities
+1 - pbeta(0.25, 41, 11)
+#[1] 1
+1 - pbeta(0.5, 41, 11)
+#[1] 0.9999926
+1 - pbeta(0.8, 41, 11)
+#[1] 0.4444044Equal tailed 95% credible interval
+qbeta(0.025, 41, 11)
+#[1] 0.6688426
+qbeta(0.975, 41, 11)
+#[1] 0.887109495% confidence that $\theta_1$ is between 0.67 and 0.89
+6) Suppose the second student gets 24 questions right. What is the posterior distribution for $\theta_2$? $P(\theta > 0.25)$? $P(\theta > 0.5)$? $P(\theta > 0.8)$? What is the 95% posterior credible interval for $\theta_2$ +$$ +Posterior \sim Beta(8 + 24, 4 + 40 - 24) = Beta(32, 20) +$$ +With a posterior mean of $\frac{32}{32+20} = \frac{32}{52}$
+We can plot the posterior distribution with the prior
+plot(theta, dbeta(theta, 32, 20), type = 'l')
+lines(theta, dbeta(theta, 8 ,4), lty = 2) #Dashed line for priorPosterior probabilities
+1 - pbeta(0.25, 32, 20)
+#[1] 1
+1 - pbeta(0.5, 32, 20)
+#[1] 0.9540427
+1 - pbeta(0.8, 32, 20)
+#[1] 0.00124819Equal tailed 95% credible interval
+qbeta(0.025, 32, 20)
+#[1] 0.4808022
+qbeta(0.975, 32, 20)
+#[1] 0.741556495% confidence that $\theta_1$ is between 0.48 and 0.74
+7) What is the posterior probability that $\theta_1 > \theta_2$? i.e., that the first student has a better chance of getting a question right than the second student?
+Estimate by simulation: draw 1,000 samples from each and see how often we observe $\theta_1 > \theta_2$
+theta1 = rbeta(100000, 41, 11)
+theta2 = rbeta(100000, 32, 20)
+mean(theta1 > theta2)
+#[1] 0.975In mass produced chocolate chip cookies, they make a large amount of dough. They mix in a large number of chips, mix it up really well and then chunk out individual cookies. In this process, the number of chips per cookie approximately follow a Poisson distribution.
+If we were to assume that chips have no volume, then this would be exactly a Poisson process and follow exactly a Poisson distribution. In practice, however, chips aren't that big so they follow approximately a Poisson distribution for the number of chips per cookie. +$$ +Y_i \sim Poisson(\lambda) +$$
+$$ +f(y|\lambda) = \frac{\lambda^{\sum{yi}}e^{-n\lambda}}{\prod{i = 1}^n{y_i!}} +$$
+This is for $\lambda > 0$
+What type of prior should we put on $\lambda$? It would be convenient if we could put a conjugate prior. What distribution looks like lambda raised to a power and e raised to a negative power?
+For this, we're going to use a Gamma prior. +$$ +\lambda \sim \Gamma(\alpha, \beta) +$$
+$$ +f(\lambda) = \frac{\beta^\alpha}{\Gamma(\alpha)}\lambda^{\alpha - 1}e^{-\beta\lambda} +$$
+$$ +f(\lambda | y) \propto f(y|\lambda)f(\lambda) \propto \lambda^{\sum{y_i}}e^{-n\lambda}\lambda^{\alpha - 1}e^{-\beta \lambda} +$$
+$$ +f(\lambda | y) \propto \lambda^{\alpha + \sum{y_i} - 1}e^{-(\beta + n)\lambda} +$$
+Thus we can see that the posterior is a Gamma Distribution +$$ +\lambda|y \sim \Gamma(\alpha + \sum{y_i}, \beta + n) +$$ +The mean of Gamma under this parameterization is $\frac{\alpha}{\beta}$
+The posterior mean is going to be +$$ +posterior_{mean} = \frac{\alpha + \sum{y_i}}{\beta + n} = \frac{\beta}{\beta + n}\frac{\alpha}{\beta} + \frac{n}{\beta + n}\frac{\sum{y_i}}{n} +$$ +As you can see here the posterior mean of the Gamma distribution is also the weighted average of the prior mean and the data mean.
+Let us present two strategies on how to choose our hyper parameters $\alpha$ and $\beta$
+After this, we need some other piece of knowledge to pin point both parameters. Here are some options.
+This would create a distribution with a mean of 1 and a huge standard deviation across the whole space. Hence the posterior will be largely driven by the data and very little by the prior.
+Suppose you're waiting for a bus that you think comes on average once every 10 minutes, but you're not sure exactly how often it comes. +$$ +Y \sim Exp(\lambda) +$$ +Your waiting time has a prior expectation of $\frac{1}{\lambda}$
+It turns out the gamma distribution is conjugate for an exponential likelihood. We need to specify a prior, or a particular gamma in this case. If we think that the buses come on average every ten minutes, that's a rate of one over ten. +$$ +prior_{mean} = \frac{1}{10} +$$ +Thus, we'll want to specify a gamma distribution so that the first parameter divded by the second parameter is $\frac{1}{10}$
+We can now think about our variability. Perhaps you specify +$$ +\Gamma(100, 1000) +$$ +This will indeed have a prior mean of $\frac{1}{10}$ and it'll have a standard deviation of $\frac{1}{100}$. If you want to have a rough estimate of our mean plus or minus two standard deviations then we have the following +$$ +0.1 \pm 0.02 +$$ +Suppose that we wait for 12 minutes and a bus arrives. Now you want to update your posterior for $\lambda$ about how often this bus will arrive. +$$ +f(\lambda | y) \propto f(y|\lambda)f(\lambda) +$$
+$$ +f(\lambda | y) \propto \lambda e^{-\lambda y}\lambda^{\alpha - 1}e^{-\beta \lambda} +$$
+$$ +f(\lambda | y) \propto \lambda^{(\alpha + 1) - 1}e^{-(\beta + y)\lambda} +$$
+$$ +\lambda | y \sim \Gamma(\alpha + 1, \beta + y) +$$
+Plugging in our particular prior gives us a posterior for $\lambda$ which is +$$ +\lambda | y \sim \Gamma(101, 1012) +$$ +Thus our posterior mean is going to be $\frac{101}{1012} Which is equal to 0.0998.
+This one observation doesn't contain a lot of data under this likelihood. When the bus comes and it takes 12 minutes instead of 10, it barely shifts our posterior mean up. One data point doesn't have a big impact here.
+Let's suppose the standard deviation or variance is known and we're only interested in learning about the mean. This is the situation that often arises in monitoring industrial production processes. +$$ +X_i \sim N(\mu, \sigma^2) +$$ +It turns out that the Normal distribution is conjugate for itself when looking for the mean parameter
+Prior +$$ +\mu \sim N(m_0,S_0^2) +$$
+$$ +f(\mu |x ) \propto f(x|\mu)f(\mu) +$$
+$$ +\mu | x \sim N(\frac{n\bar{x}/\sigma_0^2 + m_0/s_0^2}{n/\sigma_0^2 + 1/s_0^2}, \frac{1}{n/\sigma_0^2 + 1/s_0^2}) +$$
+Let's look at the posterior mean +$$ +posterior_{mean} = \frac{n/\sigma_0^2}{n/\sigma_0^2 + 1/s_0^2}\bar{x} + \frac{1/s_0^2}{n/\sigma_0^2 + 1/s_0^2} +$$
+$$ +posterior_{mean} = \frac{n}{n + \sigma_0^2/s_0^2}\bar{x} + \frac{\sigma_0^2/s_0^2}{n + \sigma_0^2/s_0^2} +$$
+Thus we see, that the posterior mean is a weighted average of the prior mean and the data mean. And indeed that the effective sample size for this prior is the ratio of the variance for the data to the variance in the prior.
+This makes sense, because the larger the variance of the prior, the less information that's in it.
+The marginal distribution for Y is +$$ +N(m_0, s_0^2 + \sigma^2) +$$
+$$ +X_i | \mu, \sigma^2 \sim N(\mu, \sigma^2) +$$
+A prior from $\mu$ conditional on the value for $\sigma^2$ +$$ +\mu | \sigma^2 \sim N(m, \frac{\sigma^2}{w}) +$$ +$w$ is going to be the ratio of $\sigma^2$ and some variance for the Normal distribution. This is the effective sample size of the prior.
+Finally, the last step is to specify a prior for $\sigma^2$. The conjugate prior here is an inverse gamma distribution with parameters $\alpha$ and $\beta$. +$$ +\sigma^2 \sim \Gamma^{-1}(\alpha, \beta) +$$ +After many calculations... we get the posterior distribution +$$ +\sigma^2 | x \sim \Gamma^{-1}(\alpha + \frac{n}{2}, \beta + \frac{1}{2}\sum_{i = 1}^n{(x-\bar{x}^2 + \frac{nw}{2(n+2)}(\bar{x} - m)^2)}) +$$
+$$ +\mu | \sigma^2,x \sim N(\frac{n\bar{x}+wm}{n+w}, \frac{\sigma^2}{n + w}) +$$
+Where the posterior mean can be written as the weighted average of the prior mean and the data mean. +$$ +\frac{n\bar{x}+wm}{n+w} = \frac{w}{n + w}m + \frac{n}{n + w}\bar{x} +$$ +In some cases, we really only care about $\mu$. We want some inference on $\mu$ and we may want it such that it doesn't depend on $\sigma^2$. We can marginalize that $\sigma^2$ integrating it out. The posterior for $\mu$ marginally follows a $t$ distribution. +$$ +\mu | x \sim t +$$ +Similarly the posterior predictive distribution also is a $t$ distribution.
+Finally, note that we can extend this in various directions, this can be extended to the multivariate normal case that requires matrix vector notations and can be extended in a hierarchial fashion if we want to specify priors for $m$, $w$ and $\beta$
+We've seen examples of choosing priors that contain a significant amount of information. We've also seen some examples of choosing priors where we're attempting to not put too much information in to keep them vague.
+Another approach is referred to as objective Bayesian statistics or inference where we explicitly try to minimize the amount of information that goes into the prior.
+This is an attempt to have the data have maximum influence on the posterior
+Let's go back to coin flipping +$$ +Y_i \sim B(\theta) +$$ +How do we minimize our prior information in $\theta$? One obvious intuitive approach is to say that all values of $\theta$ are equally likely. So we could have a piror for $\theta$ which follows a uniform distribution on the interval $[0, 1]$
+Saying all values of $\theta$ are equally likely seems like it would have no information in it.
+Recall however, that a $Uniform(0, 1)$ is the same as $Beta(1, 1)$
+The effective sample size of a beta prior is the sum of its two parameters. So in this case, it has an effective sample size of 2. This is equivalent to data, with one head and one tail already in it.
+So this is not a completely non informative prior.
+We could think about a prior that has less information. For example $Beta(\frac{1}{2}, \frac{1}{2})$, this would have half as much information with an effective sample size of one.
+We can take this even further. Think about something like $Beta(0.001, 0.001)$ This would have much less information, with the effective sample size fairly close to zero. In this case, the data would determine the posterior and there would be very little influence from the prior.
+Can we go even further? We can think of the limiting case. Let's think of $Beta(0,0)$, what would that look like? +$$ +f(\theta) \propto \theta^{-1}(1-\theta)^{-1} +$$ +This is not a proper density. If you integrate this over $(0,1)$, you'll get an infinite integral, so it's not a true density in the sense of it not integrating to 1.
+There's no way to normalize it, since it has an infinite integral. This is what we refer to as an improper prior.
+It's improper in the sense that it doesn't have a proper density. But it's not necessarily imporper in the sense that we can't use it. If we collect data, we use this prior and as long as we observe one head and one tail, or at least one success and one failure. Then we can get a posterior +$$ +f(\theta|y) \propto \theta^{y-1}(1-\theta)^{n-y-1} \sim Beta(y, n-y) +$$ +With a posterior mean of $\frac{y}{n} =\hat{\theta}$, which you should recognize as the maximum likelihood estimate. So by using this improper prior, we get a posterior which gives us point estimates exactly the same as the frequentest approach.
+But in this case, we can also think of having a full posterior. From this, we can make interval statements, probability statements, and we can actually find an interval and say that there's 95% probability that $\theta$ is in this interval. Which is not something you can do under the frequentest approach even though we may get the same exact interval.
+Improper priors are okay as long as the posterior itself is proper. There may be some mathematical things that need to be checked and you may need to have certain restrictions on the data. In this case, we needed to make sure that we observed at least one head and one tail to get a proper posterior.
+But as long as the posterior is proper, we can go forwards and do Bayesian inference even with an improper prior.
+The second point is that for many problems there does exist a prior, typically an improper prior that will lead to the same point estimates as you would get under the frequentest paradigm. So we can get very similar results, results that are fully dependent on the data, under the Bayesian approach.
+But in this case, we can also have continue to have a posterior and make posterior interval estimates and talk about posterior probabilities of the parameter.
+Another example is thinking about the normal case. +$$ +Y_i \stackrel{iid}\sim N(\mu, \sigma^2) +$$ +Let's start off by assuming that $\sigma^2$ is known and we'll just focus on the mean $\mu$.
+We can think about a vague prior like before and say that +$$ +\mu \sim N(0, 1000000^2) +$$ +This would just spread things out across the real line. That would be a fairly non informative prior covering a lot of possibilities. We can then think about taking the limit, what happens if we let the variance go to $\infty$. In that case, we're basically spreading out this distribution across the entire real number line. We can say that the density is just a constant across the whole real line. +$$ +f(\mu) \propto 1 +$$ +This is an improper prior because if you integrate the real line you get an infinite answer. However, if we go ahead and find the posterior +$$ +f(\mu|y) \propto f(y|\mu)f(\mu) \propto exp(-\frac{1}{2\sigma^2}\sum{(y_i - \mu)^2})(1) +$$
+$$ +f(\mu | y) \propto exp(-\frac{1}{2\sigma^2/n}(\mu - \bar{y})^2) +$$
+$$ +\mu | y \sim N(\bar{y}, \frac{\sigma^2}{n}) +$$
+This should look just like the maximum likelihood estimate.
+In the case that $\sigma^2$ is unknown, the standard non informative prior is +$$ +f(\sigma^2) \propto \frac{1}{\sigma^2} +$$
+$$ +\sigma^2 \sim \Gamma^{-1}(0,0) +$$
+This is an improper prior and it's uniform on the log scale of $\sigma^2$.
+In this case, we'll end up with a posterior for $\sigma^2$ +$$ +\sigma^2|y \sim \Gamma^{-1}(\frac{n-1}{2}, \frac{1}{2}\sum{(y_i - \bar{y})^2}) +$$ +This should also look reminiscent of quantities we get as a frequentest. For example, the samples standard deviation
+Choosing a uniform prior depends upon the particular parameterization.
+Suppose I used a prior which is uniform on the log scale for $\sigma^2$ +$$ +f(\sigma^2) \propto \frac{1}{\sigma^2} +$$ +Suppose somebody else decides, that they just want to put a uniform prior on $\sigma^2$ itself. +$$ +f(\sigma^2) \propto 1 +$$ +These are both uniform on certain scales or certain parameterizations, but they are different priors. So when we compute the posteriors, they will be different as well.
+The key thing is that uniform priors are not invariant with respect to transformation. Depending on how you parameterize the problem, you can get different answers by using a uniform prior
+One attempt to round this out is to use Jeffrey's Prior
+Jeffrey's Prior is defined as the following +$$ +f(\theta) \propto \sqrt{\mathcal{I(\theta)}} +$$ +Where $\mathcal{I}(\theta)$ is the fisher information of $\theta$. In most cases, this will be an improper prior.
+For the example of Normal Data +$$ +Y_i \sim N(\mu, \sigma^2) +$$
+$$ +f(\mu) \propto 1 +$$
+$$ +f(\sigma^2) \propto \frac{1}{\sigma^2} +$$
+Where $\mu$ is uniform and $\sigma^2$ is uniform on the log scale.
+This prior will then be transformation invariant. We will end up putting the same information into the prior even if we use a different parameterization for the Normal.
+$$ +Y_i \sim B(\theta) +$$
+$$ +f(\theta) \propto \theta^{-\frac{1}{2}}(1-\theta)^{-\frac{1}{2}} \sim Beta(\frac{1}{2},\frac{1}{2}) +$$
+This is a rare example of where the Jeffreys prior turns out to be a proper prior.
+You'll note that this prior actually does have some information in it. It's equivalent to an effective sample size of one data point. However, this information will be the same, not depending on the parameterization we use.
+In this case, we have $\theta$ as a probability, but another alternative which is sometimes used is when we model things on a logistics scale.
+By using the Jeffreys prior, we'll maintain the exact same information.
+Other possible approaches to objective Bayesian inference includes priors such as reference priors and maximum entropy priors.
+A related concept to this is called empirical Bayesian analysis. The idea in empirical Baye's is that you use the data to help inform your prior; such as by using the mean of the data to set the mean of the prior distribution. This approach often leads to reasonable point estimates in your posterior. However, it's sort of cheating since you're using your data twice and as a result may lead to improper uncertainty estimates.
+The Fisher information (for one parameter) is defined as +$$ +\mathcal{I}(\theta) = E[(\frac{d}{d\theta}log{(f(X|\theta))})^2] +$$ +Where the expectation is taken with respect to $X$ which has PDF $f(x|\theta)$. This quantity is useful in obtaining estimators for $\theta$ with good properties, such as low variance. It is also the basis for Jeffreys prior.
+Example: Let $X | \theta \sim N(\theta, 1)$. Then we have +$$ +f(x|\theta) = \frac{1}{\sqrt{2\pi}}exp[-\frac{1}{2}(x-\theta)^2] +$$
+$$ +\log{(f(x|\theta))} = -\frac{1}{2}\log{(2\pi)}-\frac{1}{2}(x-\theta)^2 +$$
+$$ +(\frac{d}{d\theta}log{(f(x|\theta))})^2 = (x-\theta)^2 +$$
+and so $\mathcal{I}(\theta) = E[(X - \theta)^2] = Var(X) = 1$
+Recall that linear regression is a model for predicting a response or dependent variable ($Y$, also called an output) from one or more covariates or independent variables ($X$, also called explanatory variables, inputs, or features). For a given value of a single $x$, the expected value of $y$ is +$$ +E[y] = \beta_0 + \beta_1x +$$ +or we could say that $Y \sim N(\beta_0 + \beta_1x, \sigma^2)$. For data $(x_1, y_1), \dots , (x_n, y_n)$, the fitted values for the coefficients, $\hat{\beta_0}$ and $\hat{\beta1}$ are those that minimize the sum of squared errors $\sum{i = 1}^n{(y_i - \hat{y_i})^2}$, where the predicted values for the response are $\hat{y} = \hat{\beta_0} + \hat{\beta_1}x$. We can get these values from R. These fitted coefficients give the least-squares line for the data.
+This model extends to multiple covariates, with one $\beta_j$ for each $k$ covariates +$$ +E[y_i] = \beta_0 + \beta1x{i1} + \dots + \betakx{ik} +$$ +Optionally, we can represent the multivariate case using vector-matrix notation.
+In the Bayesian framework, we treat the $\beta$ parameters as unknown, put a prior on them, and then find the posterior. We might treat $\sigma^2$ as fixed and known, or we might treat it as an unknown and also put a prior on it. Because the underlying assumption of a regression model is that the errors are independent and identically normally distributed with mean $0$ and variance $\sigma^2$, this defines a normal likelihood.
+Sometimes we may know the value of the error variance $\sigma^2$. This simplifies calculations. The conjugate prior for the $\beta$'s is a normal prior. In practice, people typically use a non-informative prior, i.e., the limit as the variance of the normal prior goes to infinity, which has the same mean as the standard least-squares estimates. If we are only estimating $\beta$ and treating $\sigma^2$ as known, then the posterior for $\beta$ is a (multivariate) normal distribution. If we just have a single covariate, then the posterior for the slope is +$$ +\beta1 | y \sim N(\frac{\sum{i = 1}^n{(x_i-\bar{x})(yi - \bar{y})}}{\sum{i=1}^n{(xi-\bar{x})^2}}, \frac{\sigma^2}{\sum{i=1}^n{(x_i - \bar{x})^2}}) +$$ +If we have multiple covariates, then using a matrix-vector notation, the posterior for the vector of coefficients is +$$ +\beta | y \sim N((X^tX)^{-1}X^ty, (X^tX)^{-1}\sigma^2) +$$ +where $X$ denotes the design matrix and $X^t$ is the transpose of $X$. The intercept is typically included in $X$ as a column of $1$'s. Using an improper prior requires us to have at least as many data points as we have parameters to ensure that the posterior is proper.
+If we treat both $\beta$ and $\sigma^2$ as unknown, the standard prior is the non-informative Jeffreys prior, $f(\beta, \sigma^2) \propto \frac{1}{\sigma^2}$. Again, the posterior mean for $\beta$ will be the same as the standard least-squares estimates. The posterior for $\beta$ conditional on $\sigma^2$ is the same normal distributions as when $\sigma^2$ is known, but the marginal posterior distribution for $\beta$, with $\sigma^2$ integrated out is a $t$ distribution, analogous to the $t$ tests for significance in standard linear regression. The posterior $t$ distribution has mean $(X^tX)^{-1}X^ty$ and scale matrix (related to the variance matrix) $s^2(X^tX)^{-1}$, where $s^2 = \sum_{i = 1}^n{(y_i - \hat{yi})^2/(n - k - 1)}$. The posterior distribution for $\sigma^2$ is an inverse gamma distribution +$$ +\sigma^2 | y \sim \Gamma^{-1}(\frac{n - k - 1}{2}, \frac{n - k - 1}{2}s^2) +$$ +In the simple linear regression case (single variable), the marginal posterior for $\beta$ is a $t$ distribution with mean $\frac{\sum{i = 1}^n{(x_i-\bar{x})(yi - \bar{y})}}{\sum{i=1}^n{(xi-\bar{x})^2}}$ and scale $ \frac{s^2}{\sum{i=1}^n{(x_i - \bar{x})^2}}$. If we are trying to predict a new observation at a specified input $x^$, that predicted value has a marginal posterior predictive distribution that is a $t$ distribution, with mean $\hat{y} = \hat{\beta_0} + \hat{\beta_1}x^$ and scale $se_r\sqrt{1 + \frac{1}{n} + \frac{(x^* - \bar{x})^2}{(n - 1)s_x^2}}$. $se_r$ is the residual standard error of the regression, which can be found easily in R. $s_x^2$ is the sample variance of $x$. Recall that the predictive distribution for a new observation has more variability than the posterior distribution for $\hat{y}$ , because individual observations are more variable than the mean.
+We'll be looking at the Challenger dataset. It contains 23 past launches where it has the temperature at the day of launch and the O-ring damage index
+http://www.randomservices.org/random/data/Challenger2.txt +Read in the data
+oring=read.table("http://www.randomservices.org/random/data/Challenger2.txt",
+ header=T)
+# Note that attaching this masks T which is originally TRUE
+attach(oring)Now we'll see the plot
+plot(T, I)
Fit a linear model
+oring.lm=lm(I~T)
+summary(oring.lm)Output of the summary
+Call:
+lm(formula = I ~ T)
+
+Residuals:
+ Min 1Q Median 3Q Max
+-2.3025 -1.4507 -0.4928 0.7397 5.5337
+
+Coefficients:
+ Estimate Std. Error t value Pr(>|t|)
+(Intercept) 18.36508 4.43859 4.138 0.000468
+T -0.24337 0.06349 -3.833 0.000968
+
+(Intercept) ***
+T ***
+---
+Signif. codes:
+0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+
+Residual standard error: 2.102 on 21 degrees of freedom
+Multiple R-squared: 0.4116, Adjusted R-squared: 0.3836
+F-statistic: 14.69 on 1 and 21 DF, p-value: 0.0009677Add the fitted line into the scatterplot
+lines(T,fitted(oring.lm)) 
Create a 95% posterior interval for the slope
+-0.24337 - 0.06349*qt(.975,21)
+# [1] -0.3754047
+-0.24337 + 0.06349*qt(.975,21)
+# [1] -0.1113353Note: These are the same as the frequentest confidence intervals
+If the challenger launch was at 31 degrees Fahrenheit, how much O-Ring damage would we predict?
+coef(oring.lm)[1] + coef(oring.lm)[2]*31
+# [1] 10.82052 Let's make our posterior prediction interval
+predict(oring.lm,data.frame(T=31),interval="predict")Output of predict
fit lwr upr
+1 10.82052 4.048269 17.59276We can calculate the lower bound through the following formula
+10.82052-2.102*qt(.975,21)*sqrt(1+1/23+((31-mean(T))^2/22/var(T)))What's the posterior probability that the damage index is greater than zero?
+1-pt((0-10.82052)/(2.102*sqrt(1+1/23+((31-mean(T))^2/22/var(T)))),21)We're looking at Galton's seminal data predicting the height of children from the height of the parents
+http://www.randomservices.org/random/data/Galton.txt +Read in the data
+heights=read.table("http://www.randomservices.org/random/data/Galton.txt",
+ header=T)
+attach(heights)What are the columns in the dataset?
+names(heights)
+# [1] "Family" "Father" "Mother" "Gender" "Height" "Kids" Let's look at the relationship between the different variables
+pairs(heights)
First let's start by creating a linear model taking all of the columns into account
+summary(lm(Height~Father+Mother+Gender+Kids))Output of summary
Call:
+lm(formula = Height ~ Father + Mother + Gender + Kids)
+
+Residuals:
+ Min 1Q Median 3Q Max
+-9.4748 -1.4500 0.0889 1.4716 9.1656
+
+Coefficients:
+ Estimate Std. Error t value Pr(>|t|)
+(Intercept) 16.18771 2.79387 5.794 9.52e-09
+Father 0.39831 0.02957 13.472 < 2e-16
+Mother 0.32096 0.03126 10.269 < 2e-16
+GenderM 5.20995 0.14422 36.125 < 2e-16
+Kids -0.04382 0.02718 -1.612 0.107
+
+(Intercept) ***
+Father ***
+Mother ***
+GenderM ***
+Kids
+---
+Signif. codes:
+0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+
+Residual standard error: 2.152 on 893 degrees of freedom
+Multiple R-squared: 0.6407, Adjusted R-squared: 0.6391
+F-statistic: 398.1 on 4 and 893 DF, p-value: < 2.2e-16As you can see here, the Kids column is not significant. Let's look at a model with it removed.
summary(lm(Height~Father+Mother+Gender))Output of summary
Call:
+lm(formula = Height ~ Father + Mother + Gender)
+
+Residuals:
+ Min 1Q Median 3Q Max
+-9.523 -1.440 0.117 1.473 9.114
+
+Coefficients:
+ Estimate Std. Error t value Pr(>|t|)
+(Intercept) 15.34476 2.74696 5.586 3.08e-08
+Father 0.40598 0.02921 13.900 < 2e-16
+Mother 0.32150 0.03128 10.277 < 2e-16
+GenderM 5.22595 0.14401 36.289 < 2e-16
+
+(Intercept) ***
+Father ***
+Mother ***
+GenderM ***
+---
+Signif. codes:
+0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+
+Residual standard error: 2.154 on 894 degrees of freedom
+Multiple R-squared: 0.6397, Adjusted R-squared: 0.6385
+F-statistic: 529 on 3 and 894 DF, p-value: < 2.2e-16This model looks good, let's go ahead and save it to a variable
+heights.lm=lm(Height~Father+Mother+Gender)From this we can tell that for each extra inch of height in a father is correlated with an extra 0.4 inches extra to the height of a child.
+We can also tell that each extra inch of height in a mother is correlated with an extra 0.3 inches extra to the height of the child.
+A male child is on average 5.2 inches taller than a female child.
+Let's create a 95% posterior interval for the difference in height by gender
+5.226 - 0.144*qt(.975,894)
+# [1] 4.943383
+5.226 + 0.144*qt(.975,894)
+# [1] 5.508617Let's make a posterior prediction interval for a male and female with a father whose 68 inches and a mother whose 64 inches.
+predict(heights.lm,data.frame(Father=68,Mother=64,Gender="M"),
+ interval="predict")
+
+# fit lwr upr
+# 1 68.75291 64.51971 72.9861predict(heights.lm,data.frame(Father=68,Mother=64,Gender="F"),
+ interval="predict")
+
+# fit lwr upr
+# 1 63.52695 59.29329 67.76062This concludes this course, if you want to go further with Bayesian statistics, the next topics to explore would be hierarchal modeling and fitting of non conjugate models with Markov Chain Monte Carlo or MCMC.
+In the Winter of 2017, I took a course on Bayesian Statistics on Coursera offered by Dr. Herbert Lee.
+Below are the notes for each of the four weeks.
+ + + + +The ultimate standard for strengthening scientific evidence is replication of finding and conducting studies with independent
+Replication is particularly important in studies that can impact broad policy or regulatory decisions
+Some studies cannot be replicated
+Reproducible Research: Make analytic data and code available so that others may reproduce findings
+Reproducibility bridges the gap between replication which is awesome and doing nothing.
+New technologies increasing data collection throughput; data are more complex and extremely high dimensional
+Existing databases can be merged into new "megadatabases"
+Computing power is greatly increased, allowing more sophisticated analyses
+For every field "X" there is a field "Computational X"
+Measured Data -> Analytic Data -> Computational Results -> Figures/Tables/Numeric Summaries -> Articles -> Text
+Data/Metadata used to develop test should be made publically available
+The computer code and fully specified computational procedures used for development of the candidate omics-based test should be made sustainably available
+"Ideally, the computer code that is released will encompass all of the steps of computational analysis, including all data preprocessing steps. All aspects of the analysis needs to be transparently reported" -- IOM Report
+Authors:
+Readers:
+Authors:
+Readers:
+An article is a stream of text and code
+Analysis code is divided into text and code "chunks"
+Each code chunk loads data and computes results
+Presentation code formats results (tables, figures, etc.)
+Article text explains what is going on
+Literate programs can be weaved to produce human-readable documents and tagled to produce machine-readable documents
+Literate programming is a general concept that requires
+Knitr is an R package that brings a variety of documentation languages such as Latex, Markdown, and HTML
+Reproducible research is important as a minimum standard, particularly for studies that are difficult to replicate
+Infrastructure is needed for creating and distributing reproducible document, beyond what is currently available
+There is a growing number of tools for creating reproducible documents
+Golden Rule of Reproducibility: Script Everything
+"Ask yourselves, what problem have you solved, ever, that was worth solving, where you knew all of the given information in advance? Where you didn't have a surplus of information and have to filter it out, or you had insufficient information and have to go find some?" -- Dan Myer
+Defining a question is the kind of most powerful dimension reduction tool you can ever employ.
+Start with a general question
+Can I automatically detect emails that are SPAM or not?
+Make it concrete
+Can I use quantitative characteristics of emails to classify them as SPAM?
+The data set may depend on your goal
+Sometimes you can find data free on the web
+Other times you may need to buy the data
+Be sure to respect the terms of use
+If the data don't exist, you may need to generate it yourself.
+Try to obtain the raw data
+Be sure to reference the source
+Polite emails go a long way
+If you load the data from an Internet source, record the URL and time accessed
+Raw data often needs to be processed
+If it is pre-processed, make sure you understand how
+Understand the source of the data (census, sample, convenience sample, etc)
+May need reformatting, subsampling -- record these steps
+Determine if the data are good enough -- If not, quit or change data
+Look at summaries of the data
+Check for missing data
+-> Why is there missing data?
+Look for outliers
+Create exploratory plots
+Perform exploratory analyses such as clustering
+If it's hard to see your plots since it's all bunched up, consider taking the log base 10 of an axis
+plot(log10(trainSpan$capitalAve + 1) ~ trainSpam$type)
Should be informed by the results of your exploratory analysis
+Exact methods depend on the question of interest
+Transformations/processing should be accounted for when necessary
+Measures of uncertainty should be reported.
+Use the appropriate language
+Gives an explanation
+Interpret Coefficients
+Interpret measures of uncertainty
+Challenge all steps:
+Challenge measures of uncertainty
+Challenge choices of terms to include in models
+Think of potential alternative analyses
+Lead with the question
+Summarize the analyses into the story
+Don't include every analysis, include it
+Lead with the question
+ Can I use quantitative characteristics of the emails to classify them as SPAM?
+Describe the approach
+ Collected data from UCI -> created training/test sets
+ Explored Relationships
+ Choose logistic model on training set by cross validation
+ Applied to test, 78% test set accuracy
+Interpret results
+ Number of dollar signs seem reasonable, e.g. "Make more money with Viagra $ $ $ $"
+Challenge Results
+ 78% isn't that great
+ Could use more variables
+ Why use logistic regression?
+Data
+Figures
+R Code
+Text
+Should be stored in the analysis folder
+If accessed from the web, include URL, description, and date accessed in README
+Processed data should be named so it is easy to see which script generated the data
+The processing script -- processed data mapping should occur in the README
+Processed data should be tidy
+Figures made during the course of your analysis, not necessarily part of your final report
+They do not need to be "pretty"
+Usually a small subset of the original figures
+Axes/Colors set to make the figure clear
+Possibly multiple panels
+May be less commented (but comments help you!)
+May be multiple versions
+May include analyses that are later discarded
+Clearly commented
+Small comments liberally - what, when, why, how
+Include processing details
+Only analyses that appear in the final write-up
+R Markdown files can be used to generate reproducible reports
+Text and R code are integrated
+Very easy to create in RStudio
+Not necessary if you use R Markdown
+Should contain step-by-step instructions for analysis
+It should contain a title, introduction (motivation), methods (statistics you used), results (including measures of uncertainty), and conclusions (including potential problems)
+It should tell a story
+It should not include every analysis you performed
+References should be included for statistical methods
+Markdown is a text-to-HTML conversion tool for web writers. Markdown allows you to write using an easy-to-read, easy-to-write plain text format, then convert it structurally to valid XHTML/HTML
+*This text will appear italicized!*
This text will appear italicized!
+**This text will appear bold!**
This text will appear bold
+## This is a secondary heading
###This is a tertiary heading
Unordered Lists
+- first item in list
- second item in list
Ordered lists
+1. first item in list
2. second item in list
3. third item in list
Create links
+[Download R](http://www.r-project.org/)
Advanced linking
+I spent so much time reading [R bloggers][1] and [Simply Statistics][2]!
[1]: http://www.r-bloggers.com/ "R bloggers"
[2]: http://simplystatistics.org/ "Simply Statistics"
I spent so much time reading R bloggers and Simply Statistics!
+Newlines require a double space after the end of a line
+Created by John Gruber and Aaron Swartz. It is a simplified version of "markup" languages. It allows one to focus on writing as opposed to formatting. Markdown provides a simple, minimal, and intuitive way of formatting elements.
+You can easily convert Markdown to valid HTML (and other formats) using existing tools.
+R Markdown is the integration of R code with markdown. It allows one to create documents containing "live" R code. R code is evaluated as part of the processing of the markdown and its results are inserted into the Markdown document. R Markdown is a core tool in literate statistical programming
+R Markdown can be converted to standard markdown using knitr package in R. Markdown can then be converted to HTML using the markdown package in R. This workflow can be easily managed using R Studio. One can create powerpoint like slides using the slidify package.
One of the ways to resolve this is to simply put the data and code together in the same document so that people can execute the code in the right order, and the data are read at the right times. You can have a single document that integrates the data analysis with all the textual representations.
+knitr supports a variety of documentation languageslibrary(knitr)
+setwd(<working directory>)
+knit2html('document.Rmd')
+browseURL('document.html')knitr will fill a new document with filler text; delete it
+Code chunks begin with ```{r} and ends with ```
All R code goes in between these markers
+Code chunks can have names, which is useful when we start making graphics
+```{r firstchunk}
## R code goes here
```
You can reference variable in RMarkdown through the following
+`The current time is `r time`. My favorite random number is `r rand`Example for suppressing all code chunks
+```{r setoptions, echo=FALSE}
+opts_chunk$set(echo=False, results = "hide")
+```cache=TRUE option can be set on a chunk-by-chunk basis to store results of computationPeople are busy, especially managers and leaders. Results of data analyses are sometimes presented in oral form, but often the first cut is presented via email.
+It is often useful therefore, to breakdown the results of an analysis into different levels of granularity/detail
+Things done by hand need to precisely documented (this is harder than it sounds!)
+If something needs to be done as part of your analysis / investigation, try to teach your computer to do it (even if you only need to do it once)
+In order to give your computer instructions, you need to write down exactly what you mean to do and how it should be done. Teaching a computer almost guarantees reproducibility
+For example, by, hand you can
+ 1. Go to the UCI Machine Learning Repository at http://archive.ics.uci.edu/mil/
+ 2. Download the Bike Sharing DatasetOr you can teach your computer to do it using R
+download.file("http://archive.ics.uci.edu/ml/machine-learning-databases/00275/Bike-Sharing-Dataset.zip", "ProjectData/Bike-Sharing-Dataset.zip")Notice here that:
+It helps you slow things down by adding changes into small chunks. (Don't just do one massive commit). It allows one to track / tag snapshots so that one can revert back to older versions of the project. Software like Github / Bitbucket / SourceForge make it easy to publish results.
+If you work on a complex project involving many tools / datasets, the software and computing environment can be critical for reproducing your analysis.
+Computer Architecture: CPU (Intel, AMD, ARM), CPU Architecture, GPUs
+Operating System: Windows, Mac OS, Linux / Unix
+Software Toolchain: Compilers, interpreters, command shell, programming language (C, Perl, Python, etc.), database backends, data analysis software
+Supporting software / infrastructure: Libraries, R packages, dependencies
+External dependencies: Websites, data repositories, remote databases, software repositories
+Version Numbers: Ideally, for everything (if available)
+This function in R helps report a bunch of information relating to the software environment
+sessionInfo()Avoid saving data analysis output (tables, figures, summaries, processed data, etc.), except perhaps temporarily for efficiency purposes.
+If a stray output file cannot be easily connected with the means by which it was created, then it is not reproducible
+Save the data + code that generated the output, rather than the output itself.
+Intermediate files are okay as long as there is clear documentation of how they were created.
+Random number generators generate pseudo-random numbers based on an initial seed (usually a number or set of numbers)
+ In R, you can use the set.seed() function to set the seed and to specify the random number generator to use
Setting the seed allows for the stream of random numbers to be exactly reproducible
+Whenever you generate random numbers for a non-trivial purpose, always set the seed.
+Replication
+Reproducibility
+What we get:
+What we do NOT get
+An analysis can be reproducible and still be wrong
+We want to know 'can we trust this analysis
+Does requiring reproducibility deter bad analysis?
+The premise of reproducible research is that with data/code available, people can check each other and the whole system is self-correcting
+cacher Package for RThe value of this is so other people can get the analysis or clone the analysis and look at subsets of the code. Or maybe more specifically data objects. People who want to run your code may not necessarily have the resources that you have. Because of that, they may not want to run the entire Markov chain Monte Carlo simulation that you did to get the posterior distribution or the histogram that you got at the end.
+But the idea is that you peel the onion a little bit rather than just go straight to the core.
+cacher as an Authorcachepackage function creates a cacher package storing
+cacher packagecacher as a ReaderA journal article says
+ "... the code and data for this analysis can be found in the cacher package 092dcc7dda4b93e42f23e038a60e1d44dbec7b3f"
+library(cacher)
+clonecache(id = "092dcc7dda4b93e42f23e038a60e1d44dbec7b3f")
+clonecache(id = "092d") ## Same as above
+# Created cache directory `.cache`
+showfiles()
+# [1] "top20.R"
+sourcefile("top20.R")graphcode() gives a node graph representing the code
runcode function executes code in the source filecheckcode function evaluates all expressions from scratch (no lazy-loading)checkobjects function to check for possible corruption of data perhaps during transit.You can inspect data objects with loadcache. This loads in pointers to each of the data objects into the workspace. Once you access the object, it will transfer it from the cache.
cacher Summarycacher package can be used by authors to create cache packages from data analyses for distributioncacher package to inspect others' data analyses by examing cached computationscacher is mindful of readers' resources and efficiently loads only those data objects that are needed.Particulate Matter -- PM
+When doing air pollution studies you're looking at particulate matter pollution. The dust is not just one monolithic piece of dirt or soot but it's actually composed of many different chemical constituents.
+Metals inert things like salts and other kinds of components so there's a possibility that a subset of those constituents are really harmful elements.
+PM is composed of many different chemical constituents and it's important to understand that the Environmental Protection Agency (EPA) monitors the chemical constituents of particulate matter and has been doing so since 1999 or 2000 on a national basis.
+The National Morbidity, Mortality, and Air Pollution Study (NMMAPS) was a national study of the short-term health effects of ambient air pollution
+Focused primarily on particulate matter ($PM_{10}$) and Ozone ($O_3$)
+Health outcomes include mortality from all causes and hospitalizations for cardiovascular and respiratory diseases
+Key publications
+ +Funded by the Heath Effects Institute
+
+http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1665439
+http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2137127
+One of the most important things about those three points to the right is those are called high leverage points. So the regression line can be very senstive to high leverage points. Removing those three points from the dataset brings the regression line's slope down a little bit. Which then produces a line that is no longer statistical significant (p-value about 0.31)
+In the Winter of 2017, I took a coursera course on Reproducible Research taught by Dr. Peng from John Hopkins University
+Below are my notes for each of the four weeks.
+ + + + +In the Fall of 2017, I took the course STAT 381 with Dr. Debra Hydorn. Below I included the interesting labs we worked on in the class.
+Please note that these reports were not formatted for this site. So equations and images may not show up.
+ + + + +Brandon Rozek
+The Central Limit Theorem tells us that if the sample size is large, then the distribution of sample means approach the Normal Distribution. For distributions that are more skewed, a larger sample size is needed, since that lowers the impact of extreme values on the sample mean.
+Skewness can be determined by the following formula +$$ +Sk = E((\frac{X - \mu}{\sigma})^3) = \frac{E((X - \mu)^3)}{\sigma^3} +$$ +Uniform distributions have a skewness of zero. Poisson distributions however have a skewness of $\lambda^{-\frac{1}{2}}$.
+In this lab, we are interested in the sample size needed to obtain a distribution of sample means that is approximately normal.
+In this lab, we will test for normality using the Shapiro-Wilk test. The null hypothesis of this test is that the data is normally distributed. The alternative hypothesis is that the data is not normally distributed. This test is known to favor the alternative hypothesis for a large number of sample means. To circumvent this, we will test normality starting with a small sample size $n$ and steadily increase it until we obtain a distribution of sample means that has a p-value greater than 0.05 in the Shapiro-Wilk test.
+This tells us that with a false positive rate of 5%, there is no evidence to suggest that the distribution of sample means don't follow the normal distribution.
+We will use this test to look at the distribution of sample means of both the Uniform and Poisson distribution in this lab.
+The Uniform distribution has a mean of $0.5$ and a standard deviation of $\frac{1}{\sqrt{12n}}$ and the Poisson distribution has a mean of $\lambda$ and a standard deviation of $\sqrt{\frac{\lambda}{n}}$.
+For the first part of the lab, we will sample means from a Uniform distribution and a Poisson distribution of $\lambda = 1$ both with a sample size $n = 5$.
+Doing so shows us how the Uniform distribution is roughly symmetric while the Poisson distribution is highly skewed. This begs the question: what sample size $(n)$ do I need for the Poisson distribution to be approximately normal?
+The maximum number of mean observations that the Shapiro-Wilk test allows is 5000 observations. Therefore, we will obtain n observations separately from both the Uniform or Poisson distribution and calculate the mean from it. Repeating that process 5000 times.
The mean can be calculated from the following way +$$ +Mean = \frac{\sum x_i}{n} +$$ +Where $x_i$ is the observation obtained from the Uniform or Poisson distribution
+Having a sample size of a certain amount doesn't always guarantee that it will fail to reject the Shapiro-Wilk test. Therefore, it is useful to run the test multiple times so that we can create a 95th percentile of sample sizes that fails to reject the Shapiro-Wilk test.
+The issue with this is that lower lambda values result in higher skewness's. Which is described by the skewness formula. If a distribution has a high degree of skewness, then it will take a larger sample size n to make the sample mean distribution approximately normal.
+Finding large values of n result in longer computational time. Therefore, the code takes this into account by starting at a larger value of n and/or incrementing by a larger value of n each iteration. Incrementing by a larger value of n decreases the precision, though that is the compromise I'm willing to take in order to achieve faster results.
+Finding just the first value of $n$ that generates the sample means that fails to reject the Shapiro-Wilk test doesn't tell us much in terms of the sample size needed for the distribution of sample means to be approximately normal. Instead, it is better to run this process many times, finding the values of n that satisfy this condition multiple times. That way we can look at the distribution of sample sizes required and return back the 95th percentile.
+Returning the 95th percentile tells us that 95% of the time, it was the sample size returned or lower that first failed to reject the Shapiro-Wilk test. One must be careful, because it can be wrongly interpreted as the sample size needed to fail to reject the Shapiro-Wilk test 95% of the time. Using that logic requires additional mathematics outside the scope of this paper. Returning the 95th percentile of the first sample size that failed to reject the Shapiro-Wilk test, however, will give us a good enough estimate for a sample size needed.
+Once a value for n is determined, we sample the means of the particular distribution (Uniform/Poisson) and create histograms and Q-Q plots for each of the parameters we're interested in. We're looking to verify that the histogram looks symmetric and that the points on the Q-Q Plot fit closely to the Q-Q Line with one end of the scattering of points on the opposite side of the line as the other.
Sampling the mean of the uniform distribution with $n = 5$ results in a mean of $\bar{x} = 0.498$ and standard deviation of $sd = 0.1288$. The histogram and Q-Q Plot can be seen in Figure I and Figure II respectively.
+$\bar{x}$ isn't far from the theoretical 0.5 and the standard deviation is also close to +$$ +\frac{1}{\sqrt{12(5)}} \approx 0.129 +$$ +Looking at the histogram and Q-Q plot, it suggests that data is approximately normal. Therefore we can conclude that a sample size of 5 is sufficient for the sample mean distribution coming from the normal distribution to be approximately normal.
+Sampling the mean of the Poisson distribution with $n = 5$ and $\lambda = 1$ results in a mean of $\bar{x} = 0.9918$ and a standard deviation of $sd = 0.443$. The histogram and Q-Q Plot can be seen in Figures III and IV respectively.
+$\bar{x}$ is not too far from the theoretical $\lambda = 1$, the standard deviation is a bit farther from the theoretical +$$ +\sqrt{\frac{\lambda}{n}} = \sqrt{\frac{1}{5}} = 0.447 +$$ +Looking at the Figures, however, shows us that the data does not appear normal. Therefore, we cannot conclude that 5 is a big enough sample size for the Poisson Distribution of $\lambda = 1$ to be approximately normal.
+Running the algorithm I described, I produced the following table
+| $\lambda$ | +Skewness | +Sample Size Needed | +Shapiro-Wilk P-Value | +Average of Sample Means | +Standard Deviation of Sample Means | +Theoretical Standard Deviation of Sample Means | +
|---|---|---|---|---|---|---|
| 0.1 | +3.16228 | +2710 | +0.05778 | +0.099 | +0.0060 | +0.0061 | +
| 0.5 | +1.41421 | +802 | +0.16840 | +0.499 | +0.0250 | +0.0249 | +
| 1 | +1.00000 | +215 | +0.06479 | +1.000 | +0.0675 | +0.0682 | +
| 5 | +0.44721 | +53 | +0.12550 | +4.997 | +0.3060 | +0.3071 | +
| 10 | +0.31622 | +31 | +0.14120 | +9.999 | +0.5617 | +0.5679 | +
| 50 | +0.14142 | +10 | +0.48440 | +50.03 | +2.2461 | +2.2361 | +
| 100 | +0.10000 | +6 | +0.47230 | +100.0027 | +4.1245 | +4.0824 | +
The skewness was derived from the formula in the first section while the sample size was obtained by looking at the .95 blue quantile line in Figures XVIII-XIV. The rest of the columns are obtained from the output of the R Code function show_results.
Looking at the histograms and Q-Q Plots produced by the algorithm, the distribution of sample means distributions are all roughly symmetric. The sample means are also tightly clustered around the Q-Q line, showing that the normal distribution is a good fit. This allows us to be confident that using these values of n as the sample size would result in the distribution of sample means of Uniform or Poisson (with a given lambda) to be approximately normal.
All the values of the average sampling means are within 0.001 of the theoretical average of sample means. The standard deviation of sample means slightly increase as the value of $\lambda$ increases, but it still is quite low.
+The table in the results section clearly show that as the skewness increases, so does the sample size needed to make the distribution of sample means approximately normal. This shows the central limit theorem in action in that no matter the skewness, if you obtain a large enough sample, the distribution of sample means will be approximately normal.
+These conclusions pave the way for more interesting applications such as hypothesis testing and confidence intervals.
+
rm(list=ls())
+library(ggplot2)
+
+sample_mean_uniform = function(n) {
+ xbarsunif = numeric(5000)
+ for (i in 1:5000) {
+ sumunif = 0
+ for (j in 1:n) {
+ sumunif = sumunif + runif(1, 0, 1)
+ }
+ xbarsunif[i] = sumunif / n
+ }
+ xbarsunif
+}
+
+sample_mean_poisson = function(n, lambda) {
+ xbarspois = numeric(5000)
+ for (i in 1:5000) {
+ sumpois = 0
+ for (j in 1:n) {
+ sumpois = sumpois + rpois(1, lambda)
+ }
+ xbarspois[i] = sumpois / n
+ }
+ xbarspois
+}
+
+poisson_n_to_approx_normal = function(lambda) {
+ print(paste("Looking at Lambda =", lambda))
+ ns = c()
+
+ # Speed up computation of lower lambda values by starting at a different sample size
+ # and/or lowering the precision by increasing the delta sample size
+ # and/or lowering the number of sample sizes we obtain from the shapiro loop
+ increaseBy = 1;
+ iter = 3;
+ startingValue = 2
+ if (lambda == 0.1) {
+ startingValue = 2000;
+ iter = 3;
+ increaseBy = 50;
+ } else if (lambda == 0.5) {
+ startingValue = 200;
+ iter = 5;
+ increaseBy = 10;
+ } else if (lambda == 1) {
+ startingValue = 150;
+ iter = 25;
+ } else if (lambda == 5) {
+ startingValue = 20;
+ iter = 50;
+ startingValue = 10;
+ } else if (lambda == 10) {
+ iter = 100;
+ } else {
+ iter = 500;
+ }
+
+ progressIter = 1
+ for (i in 1:iter) {
+
+ # Include a progress indicator for personal sanity
+ if (i / iter > .1 * progressIter) {
+ print(paste("Progress", i / iter * 100, "% complete"))
+ progressIter = progressIter + 1
+ }
+
+ n = startingValue
+
+ dist = sample_mean_poisson(n, lambda)
+ p.value = shapiro.test(dist)$p.value
+ while (p.value < 0.05) {
+ n = n + increaseBy
+ dist = sample_mean_poisson(n, lambda)
+ p.value = shapiro.test(dist)$p.value
+
+ # More sanity checks
+ if (n %% 10 == 0) {
+ print(paste("N =", n, " p.value =", p.value))
+ }
+ }
+ ns = c(ns, n)
+ }
+
+ print(ggplot(data.frame(ns), aes(x = ns)) +
+ geom_histogram(fill = 'white', color = 'black', bins = 10) +
+ geom_vline(xintercept = ceiling(quantile(ns, .95)), col = '#0000AA') +
+ ggtitle(paste("Histogram of N needed for Poisson distribution of lambda =", lambda)) +
+ xlab("N") +
+ ylab("Count") +
+ theme_bw())
+
+ ceiling(quantile(ns, .95)) #95% of the time, this value of n will give you a sampling distribution that is approximately normal
+}
+
+uniform_n_to_approx_normal = function() {
+ ns = c()
+ progressIter = 1
+
+ for (i in 1:500) {
+
+ # Include a progress indicator for personal sanity
+ if (i / 500 > .1 * progressIter) {
+ print(paste("Progress", i / 5, "% complete"))
+ progressIter = progressIter + 1
+ }
+
+ n = 2
+ dist = sample_mean_uniform(n)
+ p.value = shapiro.test(dist)$p.value
+ while (p.value < 0.05) {
+ n = n + 1
+ dist = sample_mean_uniform(n)
+ p.value = shapiro.test(dist)$p.value
+
+ if (n %% 10 == 0) {
+ print(paste("N =", n, " p.value =", p.value))
+ }
+
+ }
+
+ ns = c(ns, n)
+ }
+
+ print(ggplot(data.frame(ns), aes(x = ns)) +
+ geom_histogram(fill = 'white', color = 'black', bins = 10) +
+ geom_vline(xintercept = ceiling(quantile(ns, .95)), col = '#0000AA') +
+ ggtitle("Histogram of N needed for Uniform Distribution") +
+ xlab("N") +
+ ylab("Count") +
+ theme_bw())
+
+ ceiling(quantile(ns, .95)) #95% of the time, this value of n will give you a sampling distribution that is approximately normal
+}
+
+show_results = function(dist) {
+ print(paste("The mean of the sample mean distribution is:", mean(dist)))
+ print(paste("The standard deviation of the sample mean distribution is:", sd(dist)))
+ print(shapiro.test(dist))
+ print(ggplot(data.frame(dist), aes(x = dist)) +
+ geom_histogram(fill = 'white', color = 'black', bins = 20) +
+ ggtitle("Histogram of Sample Means") +
+ xlab("Mean") +
+ ylab("Count") +
+ theme_bw())
+ qqnorm(dist, pch = 1, col = '#001155', main = "QQ Plot", xlab = "Sample Data", ylab = "Theoretical Data")
+ qqline(dist, col="#AA0000", lty=2)
+}
+
+## PART I
+uniform_mean_dist = sample_mean_uniform(n = 5)
+poisson_mean_dist = sample_mean_poisson(n = 5, lambda = 1)
+
+show_results(uniform_mean_dist)
+show_results(poisson_mean_dist)
+
+## PART II
+
+print("Starting Algorithm to Find Sample Size Needed for the Poisson Distribution of Lambda = 0.1");
+n.01 = poisson_n_to_approx_normal(0.1)
+show_results(sample_mean_poisson(n.01, 0.1))
+print("Starting Algorithm to Find Sample Size Needed for the Poisson Distribution of Lambda = 0.5");
+n.05 = poisson_n_to_approx_normal(0.5)
+show_results(sample_mean_poisson(n.05, 0.5))
+print("Starting Algorithm to Find Sample Size Needed for the Poisson Distribution of Lambda = 1");
+n.1 = poisson_n_to_approx_normal(1)
+show_results(sample_mean_poisson(n.1, 1))
+print("Starting Algorithm to Find Sample Size Needed for the Poisson Distribution of Lambda = 5");
+n.5 = poisson_n_to_approx_normal(5)
+show_results(sample_mean_poisson(n.5, 5))
+print("Starting Algorithm to Find Sample Size Needed for the Poisson Distribution of Lambda = 10");
+n.10 = poisson_n_to_approx_normal(10)
+show_results(sample_mean_poisson(n.10, 10))
+print("Starting Algorithm to Find Sample Size Needed for the Poisson Distribution of Lambda = 50");
+n.50 = poisson_n_to_approx_normal(50)
+show_results(sample_mean_poisson(n.50, 50))
+print("Starting Algorithm to Find Sample Size Needed for the Poisson Distribution of Lambda = 100");
+n.100 = poisson_n_to_approx_normal(100)
+show_results(sample_mean_poisson(n.100, 100))
+
+print("Starting Algorithm to Find Sample Size Needed for the Uniform Distribution")
+n.uniform = uniform_n_to_approx_normal()
+show_results(sample_mean_uniform(n.uniform))Written by Brandon Rozek
+Confidence intervals expands the concept of a point estimation by giving a margin of error such that one can be confident that a certain percentage of the time the true parameter falls within that interval.
+In this lab, we will look at confidence intervals for a mean. This lab focuses on a certain method of confidence intervals that depends on the distribution of sample means being Normal. We will show how the violation of this assumption impacts the probability that the true parameter falls within the interval.
+The observed level of confidence tells us the proportion of times the true mean falls within a confidence interval. To show how the violation of the Normality assumption affects this, we will sample from both a Normal distribution, T distribution, and exponential distribution with different sample sizes.
+The normal and T distributions are sampled with a mean of 5 and a standard deviation of 2. The exponential deviation is sampled with a lambda of 2 or mean of 0.5.
+From the samples, we obtain the mean and the upper/lower bounds of the confidence interval. This is performed 100,000 times. That way we obtain a distribution of these statistics.
+We know that a confidence interval is valid, if the lower bound is no more than the true mean and the upper bound is no less than the true mean. From this definition, we can compute a proportion of observed confidence from the simulations
+From the distributions of statistics, we can create visualizations to support the understanding of confidence intervals.
+The first one is a scatterplot of lower bounds vs upper bounds. This plot demonstrates the valid confidence intervals in blue and the invalid ones in red. It demonstrates how confidence intervals that are invalid are not located inside the box.
+The second visualization is a histogram of all the sample means collected. The sample means that didn't belong to a valid confidence interval are shaded in red. This graphic helps demonstrate the type I errors on both sides of the distribution.
+In this lab, we're interested in seeing how our observed level of confidence differs from our theoretical level of confidence (95%) when different sample sizes and distributions are applied.
+We can see from the table section in the Appendix that sampling from a Normal or t distribution does not adversely affect our observed level of confidence. The observed level of confidence varies slightly from the theoretical level of confidence of 0.95.
+When sampling from the exponential distribution, however, the observed level of confidence highly depends upon the sample size.
+Looking at Table III, we can see that for a sample size of 10, the observed level of confidence is at a meager 90%. This is 5% off from our theoretical level of confidence. This shows how the normality assumption is vital to the precision of our estimate.
+This comes from the fact that using this type of confidence interval on a mean from a non-normal distribution requires a large sample size for the central limit theorem to take affect.
+The central limit theorem states that if the sample size is "large", the distribution of sample means approach the normal distribution. You can see how in Figure XVIII, the distribution of sample means is skewed, though as the sample size increases, the distribution of sample means become more symmetric (Figure XIX).
+From this, we can conclude that violating the underlying assumption of normality decreases the observed level of confidence. We can mitigate the decrease of the observed level of confidence when sampling means from a non-normal distribution by having a larger sample size. This is due to the central limit theorem.
+| Sample Size | +Proportion of Means Within CI | +
|---|---|
| 10 | +0.94849 | +
| 20 | +0.94913 | +
| 50 | +0.95045 | +
| 100 | +0.94955 | +
| Sample Size | +Proportion of Means Within CI | +
|---|---|
| 10 | +0.94966 | +
| 20 | +0.94983 | +
| 50 | +0.94932 | +
| 100 | +0.94999 | +
| Sample Size | +Proportion of Means Within CI | +
|---|---|
| 10 | +0.89934 | +
| 20 | +0.91829 | +
| 50 | +0.93505 | +
| 100 | +0.94172 | +
rm(list=ls())
+library(ggplot2)
+library(functional) # For function currying
+
+proportion_in_CI = function(n, mu, dist) {
+
+ # Preallocate vectors
+ lower_bound = numeric(100000)
+ upper_bound = numeric(100000)
+ means = numeric(100000)
+
+ number_within_CI = 0
+
+ ME = 1.96 * 2 / sqrt(n) ## Normal Margin of Error
+
+ for (i in 1:100000) {
+ x = numeric(n)
+
+ # Sample from distribution
+ if (dist == "Normal" | dist == "t") {
+ x = rnorm(n,mu,2)
+ } else if (dist == "Exponential") {
+ x = rexp(n, 1 / mu)
+ }
+
+ ## Correct ME if non-normal
+ if (dist != "Normal") {
+ ME = qt(0.975,n-1)*sd(x)/sqrt(n)
+ }
+
+ ## Store statistics
+ means[i] = mean(x)
+ lower_bound[i] = mean(x) - ME
+ upper_bound[i] = mean(x) + ME
+
+ # Is Confidence Interval Valid?
+ if (lower_bound[i] < mu & upper_bound[i] > mu) {
+ number_within_CI = number_within_CI + 1
+ }
+ }
+
+ # Prepare for plotting
+ lbub = data.frame(lower_bound, upper_bound, means)
+ lbub$col = ifelse(lbub$lower_bound < mu & lbub$upper_bound > mu, 'Within CI', 'Outside CI')
+ print(ggplot(lbub, aes(x = lower_bound, y = upper_bound, col = col)) +
+ geom_point(pch = 1) +
+ geom_hline(yintercept = mu, col = '#000055') +
+ geom_vline(xintercept = mu, col = '#000055') +
+ ggtitle(paste("Plot of Lower Bounds vs Upper Bounds with Sample Size of ", n)) +
+ xlab("Lower Bound") +
+ ylab("Upper Bounds") +
+ theme_bw()
+ )
+ print(ggplot(lbub, aes(x = means, fill = col)) +
+ geom_histogram(color = 'black') +
+ ggtitle(paste("Histogram of Sample Means with Sample Size of ", n)) +
+ xlab("Sample Mean") +
+ ylab("Count") +
+ theme_bw()
+ )
+
+ # Return proportion within CI
+ number_within_CI / 100000
+}
+
+sample_sizes = c(10, 20, 50, 100)
+
+### PART I
+proportion_in_CI_Normal = Curry(proportion_in_CI, dist = "Normal", mu = 5)
+p_norm = sapply(sample_sizes, proportion_in_CI_Normal)
+sapply(p_norm, function(x) {
+ cat("The observed proportion of intervals containing mu is", x, "\n")
+ invisible(x)
+})
+
+### PART II
+proportion_in_CI_T = Curry(proportion_in_CI, dist = "t", mu = 5)
+p_t = sapply(sample_sizes, proportion_in_CI_T)
+sapply(p_t, function(x) {
+ cat("The observed proportion of intervals containing mu is", x, "\n")
+ invisible(x)
+})
+
+### PART III
+proportion_in_CI_Exp = Curry(proportion_in_CI, dist = "Exponential", mu = 0.5)
+p_exp = sapply(sample_sizes, proportion_in_CI_Exp)
+sapply(p_exp, function(x) {
+ cat("The observed proportion of intervals containing mu is", x, "\n")
+ invisible(x)
+})The generation of random numbers have a variety of applications including but not limited to the modeling of stochastic processes. It is important, therefore, to be able to generate any random number following a given distribution. One of the many ways to achieve this is to transform numbers sampled from a random uniform distribution.
+In this lab, we will compare the effectiveness between the linear congruence method (LCM), runif, and rexp to generate random numbers. The programming language R will be used and different plots and summary statistics are drawn upon to analyze the effectiveness of the methods.
The linear congruential method (LCM) is an algorithm that produces a sequence of pseudo-random numbers using the following recursive function
+$$
+X_{n + 1} = (aX_n + c) \mod m
+$$
+The R code we use has a c value of 0 which is a special case known as the multiplicative congruential generator (MCG).
There are several conditions for using a MCG. First, the seed value must be co-prime to m. In other words, the greatest common denominator between the two must be 1. A function was written in R to check this. Secondly, m must be a prime number or a power of a prime number.
An element of periodicity comes into play when dealing with pseudo-random number generators. Once the generator has produced a certain number of terms, it is shown to loop back to the first term of the sequence. It is advantageous, therefore, to keep the period high.
+The period in a MCG is at most m - 1. In this lab, m has a value of $2^{31}$ and only 100 numbers were sampled from the LCM. This allows us to avoid the issue of periodicity entirely in our analysis.
The default pseudo-random number generator (pseudo RNG) in R is the Mersenne-Twister with the default seed value of the current time and the process id of the current R instance. Mersenne-Twister is part of a class of pseudo RNG called the generalized feedback register. This class of pseudo RNGs provide a very long period (VLP) or in other words, a long cycle before the values start repeating. VLPs are often useful when executing large simulations in R.
+Since this method is already coded in the base R-package, we will leave the implementation for the curious to research.
The uniform distribution function $Unif(\theta_1, \theta_2)$ has a probability mass function (PMF) of the following. +$$ +f(x) = \frac{1}{\theta_2 - \theta_1} +$$ +Where x is valid between [$\theta_1$, $\theta_2$].
+In our case, we are producing numbers from [0,1]. We can therefore reduce the probability mass function to the following +$$ +f(x) = 1 +$$
+The expected value can be defined as +$$ +E(X) = \int_a^b xf(x)dx +$$ +For the uniform distribution we're using that becomes +$$ +E(X) = \int_0^1 xdx = \frac{1}{2}[x^2]_0^1 = \frac{1}{2} +$$ +We should, therefore, expect the mean of the generated random number sequence to roughly equal $\frac{1}{2}$.
+To find the median of the distribution, we need to find the point at which the cumulative density function (CDF) is equal to .5. +$$ +\int_0^{Median(X)} f(x)dx = \frac{1}{2} +$$ +$$ +\int_0^{Median(X)} dx = \frac{1}{2} +$$
+$$ +[X]_0^{Median(X)} = \frac{1}{2} +$$
+$$ +Median(X) = \frac{1}{2} +$$
+This shows us that the median of the distribution should be roughly equal to $\frac{1}{2}$ as well.
+The histogram of a uniform distribution of data should look rectangular in shape. This means that the counts of all the sub-intervals should be about the same.
+Another way to test whether or not the distribution is a good fit is to use what is called a Quantile-Quantile plot (Q-Q Plot). This is a probability plot that compares the quantiles from the theoretical distribution, this is the uniform distribution, to that of the sample.
+For a precise Q-Q Plot, we need a lot of quantiles to compare. In this lab, we will compare 100 different quantiles. The quantiles from the theoretical distribution can be easily derived from the fact that all value ranges are equally probable. The theoretical quantiles in this case is 0.01, 0.02, ..., 0.98, 0.99, 1. The sample quantiles are obtianed by receiving 100 observations from runif or the LCM.
In the Q-Q plot, a good fit is shown when the points hug closely to the Q-Q line. It is also important to have symmetry in the points. This means that the points are ending on opposite sides of the Q-Q line.
+For sake of exploration, we will use 5 different seed values for the linear congruential method (while making sure that the seed values are co-prime). We will then use the results of these and compare LCM to how the standard runif method generates random numbers.
The exponential distribution is a type of continuous distribution that is defined by the following PMF +$$ +f(x) = \lambda e^{-\lambda t} +$$ +We can find the CDF by taking the integral of the PMF. +$$ +F(x) = \int f(x)dx = \lambda \int e^{-\lambda x} dx = \lambda * (\frac{-1}{\lambda}) e^{-\lambda x} + C = -e^{-\lambda x} + C +$$ +One of the conditions for the cumulative density function is that as $x \to \infty$, $F(x) \to 1$. +$$ +1 = \lim{x \to \infty} (-e^{-\lambda x} + C) = \lim{x \to \infty} (-e^{-\lambda x}) + lim{x \to \infty} ( C) = \lim{x \to \infty}(-e^{\lambda x}) + C +$$ +This shows that $C = 1$, since $lim_{x \to \infty} (-e^{-\lambda x})$ is equal to 0.
+Substituting $C$ gives us the following. +$$ +F(x) = 1 - e^{-\lambda x} +$$
+We can find the expected value using the formula described in the previous Expected Value section. Reflected in the bounds of integration is the domain of the exponential distribution $[0, \infty)$. +$$ +E(X) = \lim_{a \to \infty}\int_0^a x \lambda e^{-\lambda x} dx +$$ +The integral above features a product of two functions. So as an aside, we will derive a formula so we can take the integral above.
+The total derivative is defined as +$$ +d(uv) = udv + vdu +$$
+After taking the integral of both sides, we can solve for a formula that gives the following +$$ +\int d(uv) = \int udv + \int vdu +$$
+$$ +\int udv = uv - \int vdu +$$
+The formula above is called Integration by Parts. We will make use of this formula by defining $u = \lambda x$ and $dv = e^{-\lambda x} dx$
+This implies that $du = \lambda dx$ and $v= -\frac{1}{\lambda}e^{-\lambda x}$
+$$
+E(X) = -\lim_{a \to \infty}[\lambda x \frac{1}{\lambda}e^{-\lambda x}]0^a - \lim{b \to \infty}\int_0^b -\frac{1}{\lambda}e^{-\lambda x}\lambda dx
+$$
$$ +E(X) = -\lim_{a \to \infty} [xe^{-\lambda x}]0^a - \lim{b \to \infty}\int_0^b -e^{-\lambda x}dx +$$
+$$ +E(X) = -\lim{a \to \infty}(ae^{-\lambda a}) - \lim{b \to \infty}[\frac{1}{\lambda}e^{-\lambda x}]_0^b +$$
+$$ +E(X) = 0 - \frac{1}{\lambda}[\lim_{b \to \infty}(e^{-\lambda b}) - e^{-0\lambda}] +$$
+$$ +E(X) = -\frac{1}{\lambda}(-1) = \frac{1}{\lambda} +$$
+For the purposes of this lab, we will make the rate ($\lambda$) equal to 3. Which means we should expect our mean to be roughly equal to $\frac{1}{3}$.
+The median can be found by setting the CDF equal to $\frac{1}{2}$ +$$ +1- e^{-\lambda Median(X)} = \frac{1}{2} +$$
+$$ +e^{-\lambda Median(X)} = \frac{1}{2} +$$
+$$ +-\lambda Median(X) = \ln(\frac{1}{2}) +$$
+$$ +Median(X) = \frac{\ln(2)}{\lambda} +$$
+This shows that we should expect our median to be around $\frac{\ln(2)}{3} \approx 0.231$.
+Once we have a uniform distribution of values, we can then use these values to map to an exponential distribution. This is done through a technique called inverse transform sampling, the technique works as follows:
+Using these steps, we'll derive the exponential transform below. +$$ +F(X) = 1 - e^{-\lambda x} = u +$$ +Then proceeding to solve for X, we obtain the following. +$$ +e^{-\lambda X} = 1 - u +$$
+$$ +-\lambda X = \ln(1 - u) +$$
+$$ +X = \frac{-\ln(1 - u)}{\lambda} +$$
+With this formula, we can now find values for X which is a random variable that follows an exponential distribution given random uniform data $u$.
+In this lab, we are looking at the exponential distribution with a rate of 3. Therefore, the probability mass function, the cumulative distribution function, and the exponential transform can be redefined as these respectively. +$$ +f(x) = 3e^{-3x} +$$
+$$ +F(x) = 1 - e^{-3x} +$$
+$$ +X = \frac{-\ln(1 - u)}{3} +$$
+The R code makes use of the concepts above. The purpose of this code is to output the summary statistics, histograms, and Q-Q plots are used to compare the different methods.
+First the LCM is executed four times with three different seed values. The LCM is used to create a uniform distribution of data that is then compared to the standard runif function in the R language.
Then, transformations of a LCM, runif, are executed and compared with the standard rexp to create an exponential distribution of data with $\lambda = 3$.
For a small sample of 100 values, the data seems evenly distributed for all seeds and methods used. The peaks and troughs are more pronounced in the LCM methods suggesting that the runif command creates more evenly distributed data than the LCM.
Deviations from the mean and median are lowest for the LCM rather than the standard runif command with the exception of the LCM with the seed of 93463.
The Q-Q plots for all of the methods follow the Q-Q Line tightly and appear symmetric.
+A bin size of 20 is used to make the histograms easily comparable to each other. One interesting thing to note is that in Figure XI, there is an observation located far out from the rest of the data. For the purpose of a histogram, which is to show us the shape of the distribution, this one far observation skews what we can see. Therefore the next figure, Figure XII has that single outlier removed and shows the histogram of the rest of the data.
+All of the histograms appear exponential in shape. Looking at the Q-Q plots, the LCM transformation and the rexp distribution of data fit closely to the QQ-Line. All of the Q-Q Plots have the points getting further away from the line as you get higher up in the percentiles. The runif transformation quantiles diverge from the line after about the 50th percentile.
Deviations from the mean and median were about the same for both transformed techniques (runif and LCM). The rexp command did better when it came to the deviation from the mean, but performed worse when it came to deviation from the median.
The linear congruential methods performed better when it came to simulating the distribution than the standard R functions. It more accurately captured the median for both the standard uniform data, and the exponential data. Against the runif command, it also performed better at simulating the mean.
This can especially be seen when comparing the two transform techniques. The transformed LCM distribution of data followed the Q-Q line more tightly than the transformed runif.
I speculate that this is due to the seed value used. The Mersenne-Twist method performs better when the seed value doesn't have many zeros in it. Since the seed value is determined by the system time and process id, we don't know for sure if it's a proper input for the Mersenne-Twist. For the LCM, seeds were properly tested to make sure that it was co-prime with one of its parameters. This condition allowed proper seeds to work well with the LCM.
+Further research can be done on standardizing the seed values used across all the different pseudo random number generators and looking at the 6 other pseudo RNG that comes built into R. Changing the seed and random number generator can be achieved through the set.seed function.
| Method | +Mean ($\bar{x}$) | +Median ($\tilde{x}$) | +$\mu - \bar{x}$ | +$m - \tilde{x}$ | +
|---|---|---|---|---|
| LCM with seed 55555 | +0.505 | +0.521 | +-0.005 | +-0.021 | +
| LCM with seed 93463 | +0.452 | +0.402 | +0.048 | +0.098 | +
| LCM with seed 29345 | +0.520 | +0.502 | +-0.020 | +-0.002 | +
| LCM with seed 68237 | +0.489 | +0.517 | +0.011 | +-0.017 | +
| R Random Uniform | +0.480 | +0.471 | +0.020 | +0.029 | +
| Method | +Mean | +Median | +$\mu - \bar{x}$ | +$m - \tilde{x}$ | +
|---|---|---|---|---|
| LCM Transform | +0.380 | +0.246 | +-0.047 | +-0.015 | +
| RUnif Transform | +0.283 | +0.218 | +0.050 | +0.013 | +
| R Exponential | +0.363 | +0.278 | +-0.030 | +-0.047 | +
library(ggplot2)
+linear_congruential = function(seed) {
+ a = 69069
+ c = 0
+ m = 2^31
+ x = numeric(100)
+ x[1] = seed
+ for (i in 2:100) {
+ x[i] = (a * x[i-1] + c) %% m
+ }
+ xnew = x / (max(x) + 1)
+ xnew
+}
+
+gcd = function(x,y) {
+ r = x %% y;
+ return(ifelse(r, gcd(y, r), y))
+}
+check_seed = function(seed) {
+ if (gcd(seed, 2^31) == 1) {
+ print(paste("The seed value of", seed, "is co-prime"))
+ } else {
+ print(paste("The seed value of", seed, "is NOT co-prime"))
+ }
+}
+
+check_data = function(data, distribution, title) {
+ print(paste("Currently looking at", title))
+ print(summary(data))
+ print(ggplot(data.frame(data), aes(x = data)) +
+ geom_histogram(fill = 'white', color = 'black', bins = 10) +
+ xlab("Date") +
+ ylab("Frequency") +
+ ggtitle(paste("Histogram of", title)) +
+ theme_bw())
+ uqs = (1:100) / 100
+ if (distribution == "Uniform") {
+ qqplot(data, uqs, pch = 1, col = '#001155', main = paste("QQ Plot of", title), xlab = "Sample Data", ylab = "Theoretical Data")
+ qqline(uqs, distribution = qunif, col="red", lty=2)
+ } else if (distribution == "Exponential") {
+ eqs = -log(1-uqs) / 3
+ qqplot(data, eqs, pch = 1, col = '#001155', main = paste("QQ Plot of", title), xlab = "Sample Data", ylab = "Theoretical Data")
+ qqline(eqs, distribution=function(p) { qexp(p, rate = 3)}, col="#AA0000", lty=2)
+ }
+}
+
+seed1 = 55555
+seed2 = 93463
+seed3 = 29345
+seed4 = 68237
+
+check_seed(seed1)
+lin_cong = linear_congruential(seed1)
+check_data(lin_cong, "Uniform", paste("LCM Random Data with seed", seed1))
+
+check_seed(seed2)
+lin_cong2 = linear_congruential(seed2)
+check_data(lin_cong2, "Uniform", paste("LCM Random Data with seed", seed2))
+
+check_seed(seed3)
+lin_cong3 = linear_congruential(seed3)
+check_data(lin_cong3, "Uniform", paste("LCM Random Data with seed", seed3))
+
+check_seed(seed4)
+lin_cong4 = linear_congruential(seed4)
+check_data(lin_cong4, "Uniform", paste("LCM Random Data with seed", seed4))
+
+# Using the standard built-in R function
+unif = runif(100, 0, 1)
+check_data(unif, "Uniform", "R Random Uniform Data")
+
+# Building exponential from linear congruence method
+exp1 = -log(1 - lin_cong) / 3
+check_data(exp1, "Exponential", "Exponential Data from LCM Random")
+
+# Building exponential from runif
+exp2 = -log(1 - unif) / 3
+check_data(exp2, "Exponential", "Exponential Data from R Random Uniform")
+
+# Building exponential from rexp
+exp3 = rexp(100, rate = 3)
+check_data(exp3, "Exponential", "R Generated Exponential Data")This lab features a random walk simulation that answers the following question: Given the radius of a circular room, how many random steps does it take to reach to reach the wall from the center? To explore this, we will look at the distribution of steps required and describe it in terms of a discrete probability distribution. Later on, we will look at the radius parameter and find relationships between it and the different summary statistics that we compute (mean, median, range, variance, skewness).
+The paper will start off with the methods used in this lab. The programming language R was used and a simulation was created to simulate a random walk according to the description above. It will also discuss the methods used to find the discrete probability distribution that best fits the data and how the relationship between the radii and the summary statistics were found using regression.
+After discussing the methods, the paper goes into the results. The discrete probability distribution of the data is revealed and the equations of the different summary statistics based on a regression of the radius is shared.
+The simulation written performs the following steps given the size of a room:
+The number of steps required is then recorded into a vector and the simulation is performed 999 more times. This gives us the distribution of steps required to reach the wall given the room size.
+Given the distribution vector pertaining to the room size, we then find the following information:
+For every room size selected, a total of three trials of 1000 simulations each is conducted.
+Since the simulation counts the steps required before reaching the walls of the room, it describes the process of a geometric distribution. To confirm the suspicion, the function fitdistr inside the package MASS was used to find the parameter (probability) of the geometric distribution to best fit the simulation data. In this case the simulation data used is of room size 100.
The probability density histogram of the data is then shown with the overlay of the geometric distribution that best fit the data.
+Another more robust way to visually see the goodness of fit of the distribution is to use what is called a Quantile-Quantile Plot (Q-Q Plot), Using the quantiles from the theoretical distribution, we compare it to the quantiles of the simulated data and check to see that it follows the Q-Q line closely.
+From the simulations performed, vectors have been saved containing the summary statistics for the different room sizes. Vectors are then created from these vectors as columns inside a data frame. Having this data frame allows us to use ggplot to graph the data using scatter plots.
ggplot was used as opposed to the conventional base package since it allows us to add layers to the graphs, which was used to show the geometric distribution on top of the density histogram in the previous section. ggplot is also aesthetically more pleasing as it decreases the margins and makes smart choices on the shadings of the axes on the graph.
After plotting the data, regressions are then taken with respect to different polynomial degrees of radii and analyzed for best fit using the adjusted r squared value as a means for comparison.
+The adjusted r squared value is a good comparator since it tells us a proportion of how much the variation in the data is accounted for by the model.
+The shape of the distribution matches the geometric distribution as shown in Figure IX & X. Figure IX shows how the probability density histogram nicely fits the overlay of the geometric distribution with the fitted parameters. The description of the problem closely matches the description of the geometric probability distribution and both distributions are skewed right. These observations increase our confidence that the distribution is a good fit. This can further be verified with the Q-Q Plot.
+The Q-Q Plot in Figure X shows how the theoretical quantiles compare with the quantiles from the distribution. Since the observed sample quantiles start off in one side of the Q-Q line and ends up on the other side at the end, we can say that the data is not skewed with respect to the theoretical distribution. This along with the fact that the sample quantiles closely fit the Q-Q line supports the proposition that the geometric distribution is a good fit for the steps required in the random walk.
+Regressions were then performed for each summary statistic and the models that had the lowest degree polynomial with respect to radii and with relatively small standard error were chosen.
+For the mean, quadratic regression was performed and it produces the following model +$$ +\mu(Steps) = -31.409 + 7.800(radii) + 0.837(radii^2) +$$ +With this model, the adjusted $r^2$ value is 99.938%, which tells us that 99.938% of the variation in means is accounted for by the quadratic regression above. It's fit to the observed values can be seen in Figure XI.
+For the median, quadratic regression was performed and produces the following model +$$ +\widetilde{Steps} = -19.525 + 5.096(radii) + 0.687(radii^2) +$$ +The model above has an adjusted $r^2$ value of 99.957% which tells us that 99.957% of the variation in medians is accounted for by the regression above. It's fit to the observed values can be seen in Figure XIII.
+For the range, the quadratic regression performed creates the following model +$$ +Range(Steps) = 63.873 - 9.073(radii) + 5.572(radii^2) +$$ +The model has an adjusted $r^2$ value of 99.278% which tells us that 99.278% of the variation in the ranges of steps is accounted for by the regression above. It's fit to the observed values can be seen in Figure XV.
+For the variance, the cubic regression performed creates the following model +$$ +\sigma^2(Steps) = -88471.866 + 32256.076(radii) - 2418.835(radii^2) + 61.650(radii^3) +$$ +The model has an adjusted $r^2$ value of 99.779% which tells us that 99.779% of the variation in the variations of steps is accounted for by the regression above. It's fit to the observed values can be seen in Figure XVII
+Looking at the scatter plot of radii vs skewness of steps shown in Figure XVIII, there appears to be no relationship between the radii and skewness of steps. What the scatter plot suggests that the skewness is uniformly distributed.
+In summation, the distribution of steps required to reach the end of the wall follow a geometric distribution. This was backed up using the probability density histogram and the Q-Q Plot.
+Most of the summary statistics follow a quadratic regression while the variation follows a cubic regression. Skewness are uniformly distributed across the different simulations.
+In other words, the mean, median, and range follow a function that is a second degree polynomial based on the size of the room ($r^2$) and the variation follows a function that is a third degree polynomial based on the size of the room $(r^3)$
+| Trial | +Mean | +Median | +Variance | +Range | +Skewness | +
|---|---|---|---|---|---|
| 1 | +5.827 | +5.000 | +13.471 | +36.000 | +2.446 | +
| 2 | +5.761 | +5.000 | +10.635 | +22.000 | +1.780 | +
| 3 | +5.805 | +5.000 | +11.430 | +21.000 | +1.882 | +
| Trial | +Mean | +Median | +Variance | +Range | +Skewness | +
|---|---|---|---|---|---|
| 1 | +11.566 | +9.000 | +59.055 | +60.000 | +1.854 | +
| 2 | +11.940 | +10.000 | +74.487 | +68.000 | +2.467 | +
| 3 | +11.242 | +9.000 | +53.417 | +48.000 | +1.922 | +
| Trial | +Mean | +Median | +Variance | +Range | +Skewness | +
|---|---|---|---|---|---|
| 1 | +19.306 | +15.000 | +171.688 | +90.000 | +1.721 | +
| 2 | +18.127 | +15.000 | +143.284 | +74.000 | +1.712 | +
| 3 | +18.993 | +15.000 | +180.892 | +111.000 | +2.089 | +
| Trial | +Mean | +Median | +Variance | +Range | +Skewness | +
|---|---|---|---|---|---|
| 1 | +29.272 | +23.000 | +405.817 | +150.000 | +1.737 | +
| 2 | +28.382 | +22.000 | +402.801 | +146.000 | +1.862 | +
| 3 | +27.891 | +22.000 | +390.221 | +211.000 | +2.308 | +
| Trial | +Mean | +Median | +Variance | +Range | +Skewness | +
|---|---|---|---|---|---|
| 1 | +110.234 | +87.000 | +6070.976 | +547.000 | +1.780 | +
| 2 | +109.594 | +88.000 | +5928.037 | +542.000 | +1.864 | +
| 3 | +108.398 | +86.000 | +6469.467 | +621.000 | +2.191 | +
| Trial | +Mean | +Median | +Variance | +Range | +Skewness | +
|---|---|---|---|---|---|
| 1 | +623.260 | +492.000 | +197161.8 | +3262 | +1.965 | +
| 2 | +637.890 | +511.000 | +212826.6 | +3394 | +2.009 | +
| 3 | +653.661 | +528.5 | +211126.2 | +3694 | +1.841 | +
| Trial | +Mean | +Median | +Variance | +Range | +Skewness | +
|---|---|---|---|---|---|
| 1 | +2507.935 | +1939.0 | +3231980 | +11471 | +1.726 | +
| 2 | +24675.615 | +1992.0 | +3248443 | +14150 | +2.060 | +
| 3 | +2488.795 | +2001.5 | +3038915 | +14689 | +1.843 | +
| Trial | +Mean | +Median | +Variance | +Range | +Skewness | +
|---|---|---|---|---|---|
| 1 | +8983.406 | +7257.0 | +38271895 | +52549 | +2.027 | +
| 2 | +9372.917 | +7537.0 | +41658620 | +51810 | +1.867 | +
| 3 | +8974.549 | +7294.5 | +41868765 | +60326 | +2.420 | +
rm(list=ls())
+library(moments)
+library(ggplot2)
+
+simulateRoom = function(r) {
+ n = 1000
+ results = numeric(n)
+ x = 0
+ y = 0
+ dist = 0
+ num = 0
+ for (i in 1:n) {
+ while(dist < r) {
+ deg = runif(1, min=0, max=360)
+ newx = cos(deg)
+ newy = sin(deg)
+ x = x + newx
+ y = y + newy
+ num = num+ 1
+ dist=(x^2 + y^2)^0.5
+ }
+ results[i] = num
+ x= 0
+ y = 0
+ dist = 0
+ num = 0
+ i = i + 1
+ }
+ return(results)
+}
+
+repetitions = 3
+radii = c(2, 3, 4, 5, 10, 25, 50, 100)
+
+median_stepcount = c()
+mean_stepcount = c()
+variance_stepcount = c()
+range_stepcount = c()
+skewness_stepcount = c()
+for (r in radii) {
+ room_data = data.frame()
+ print(paste("Currently simulating r =", r))
+ for (i in 1:repetitions) {
+ simulation = simulateRoom(r)
+
+ ## Summary Statistics
+ variance = var(simulation)
+ skewness_ = skewness(simulation)
+ range_data = range(simulation)
+
+ room_data_temp = data.frame(i, mean(simulation), median(simulation), variance, range_data[2] - range_data[1], skewness_)
+ names(room_data_temp) = c("Trial", "Mean", "Median", "Variance", "Range", "Skewness")
+ room_data = rbind(room_data, room_data_temp)
+
+ median_stepcount = c(median_stepcount, median(simulation))
+ mean_stepcount = c(mean_stepcount, mean(simulation))
+ variance_stepcount = c(variance_stepcount, variance)
+ range_stepcount = c(range_stepcount, range_data[2] - range_data[1])
+ skewness_stepcount = c(skewness_stepcount, skewness_)
+
+ simulation = data.frame(simulation)
+ print(ggplot(simulation, aes(x = simulation)) +
+ geom_histogram(color = 'black', fill = 'white') +
+ xlab("Steps") +
+ ylab("Count") +
+ ggtitle(paste("Distribution of r =", r)) +
+ theme_bw())
+ }
+ print(room_data)
+}
+
+## Show the probability distribution
+library(MASS)
+distribution = fitdistr(simulation$simulation, "geometric")
+probability_simulation = data.frame(simulation,
+ rgeom(1:1000, distribution$estimate['prob']))
+names(probability_simulation) = c("simulation", "geometric")
+ggplot(data = probability_simulation, aes(x = simulation)) +
+ geom_histogram(aes(y = ..density..), color = 'black', fill = 'white') +
+ geom_density(aes(x = geometric), fill = '#0000AA', alpha = .5) +
+ xlab("Steps") +
+ ylab("Probability Density") +
+ ggtitle("Probability Density Histogram with Geometric Overlay") +
+ theme_bw()
+
+## Calculate QQ Plot of Geometric Distribution and Sample
+yquantiles = quantile(simulation$simulation, c(0.25, 0.75))
+xquantiles = qgeom(prob = distribution$estimate['prob'], c(0.25, 0.75))
+slope = diff(yquantiles) / diff(xquantiles)
+intercept = yquantiles[1] - slope * xquantiles[1]
+ggplot() + aes(sample=simulation$simulation) +
+ stat_qq(distribution=qgeom, dparams = distribution$estimate['prob']) +
+ geom_abline(intercept=intercept, slope=slope) +
+ ggtitle("QQ Plot of Geometric Distribution") +
+ xlab("Theoretical Quantiles") +
+ ylab("Sample Quantiles") +
+ theme_bw()
+
+## Create a dataframe
+room_simulation_data = data.frame(as.vector(sapply(radii, function(x) { rep(x, repetitions)})), mean_stepcount, median_stepcount, range_stepcount, skewness_stepcount, variance_stepcount)
+names(room_simulation_data) = c("radii", "mean_stepcount", "median_stepcount", "range_stepcount", "skewness_stepcount", "variance_stepcount")
+######### ANALYSIS OF MEANS
+ggplot(room_simulation_data, aes(x = radii, y = mean_stepcount)) +
+ geom_point(shape = 21, size = 2) +
+ xlab("Radii") +
+ ylab("Mean Steps") +
+ ggtitle("Radii vs Mean Steps") +
+ theme_bw()
+
+## Quadratic Regression of Means
+quadratic_mean_model = lm(mean_stepcount ~ poly(radii, 2, raw = T), data = room_simulation_data)
+quadratic_mean_model_summary = summary(quadratic_mean_model)
+ggplot(room_simulation_data, aes(x = radii, y = mean_stepcount)) +
+ geom_point(shape = 21, size = 2) +
+ geom_smooth(method = "lm", formula = y ~ poly(x, 2, raw = T), se = T) +
+ xlab("Radii") +
+ ylab("Mean Steps") +
+ ggtitle("Radii vs Mean Steps") +
+ theme_bw()
+
+## Analysis of Quadratic Regression of Means
+quadratic_mean_coefficients = as.vector(quadratic_mean_model$coefficients)
+print(paste("The quadratic regression equation is: Mean Steps =", quadratic_mean_coefficients[1], "+", quadratic_mean_coefficients[2], "radii", quadratic_mean_coefficients[3], "radii^2"))
+print(paste(quadratic_mean_model_summary$adj.r.squared * 100, "% of the variation of the means is accounted for by the Quadratic Model", sep = ""))
+
+######### ANALYSIS OF MEDIANS
+ggplot(room_simulation_data, aes(x = radii, y = median_stepcount)) +
+ geom_point(shape = 21, size = 2) +
+ xlab("Radii") +
+ ylab("Median Steps") +
+ ggtitle("Radii vs Median Steps") +
+ theme_bw()
+
+## Quadratic Regression of Medians
+quadratic_median_model = lm(median_stepcount ~ poly(radii, 2, raw = T), data = room_simulation_data)
+quadratic_median_model_summary = summary(quadratic_median_model)
+ggplot(room_simulation_data, aes(x = radii, y = median_stepcount)) +
+ geom_point(shape = 21, size = 2) +
+ geom_smooth(method = "lm", formula = y ~ poly(x, 2, raw = T), se = T) +
+ xlab("Radii") +
+ ylab("Median Steps") +
+ ggtitle("Radii vs Median Steps") +
+ theme_bw()
+
+## Analysis of Quadratic Regression of Medians
+quadratic_median_coefficients = as.vector(quadratic_median_model$coefficients)
+print(paste("The Quadratic Regression Equation is: Median Steps =", quadratic_median_coefficients[1], "+", quadratic_median_coefficients[2], "radii", quadratic_median_coefficients[3], "radii^2"))
+print(paste(quadratic_median_model_summary$adj.r.squared * 100, "% of the variation of the medians is accounted for by the Quadratic Model", sep = ""))
+
+######### ANALYSIS OF RANGES
+ggplot(room_simulation_data, aes(x = radii, y = range_stepcount)) +
+ geom_point(shape = 21, size = 2) +
+ xlab("Radii") +
+ ylab("Range of Steps") +
+ ggtitle("Radii vs Range Steps") +
+ theme_bw()
+
+## Quadratic Regression of Ranges
+quadratic_range_model = lm(range_stepcount ~ poly(radii, 2, raw = T), data = room_simulation_data)
+quadratic_range_model_summary = summary(quadratic_range_model)
+ggplot(room_simulation_data, aes(x = radii, y = range_stepcount)) +
+ geom_point(shape = 21, size = 2) +
+ geom_smooth(method = "lm", formula = y ~ poly(x, 2, raw = T), se = T) +
+ xlab("Radii") +
+ ylab("Range of Steps") +
+ ggtitle("Radii vs Range of Steps") +
+ theme_bw()
+
+## Analysis of Quadratic Regression of Ranges
+quadratic_range_coefficients = as.vector(quadratic_range_model$coefficients)
+print(paste("The Quadratic Regression Equation is: Range of Steps =", quadratic_range_coefficients[1], "+", quadratic_range_coefficients[2], "radii", quadratic_range_coefficients[3], "radii^2"))
+print(paste(quadratic_range_model_summary$adj.r.squared * 100, "% of the variation of the ranges is accounted for by the Quadratic Model", sep = ""))
+
+######### ANALYSIS OF SKEWNESS
+ggplot(room_simulation_data, aes(x = radii, y = skewness_stepcount)) +
+ geom_point(shape = 21, size = 2) +
+ xlab("Radii") +
+ ylab("Skewness") +
+ ggtitle("Radii vs Skewness Steps") +
+ theme_bw()
+
+######### ANALYSIS OF VARIANCES
+ggplot(room_simulation_data, aes(x = radii, y = variance_stepcount)) +
+ geom_point(shape = 21, size = 2) +
+ xlab("Radii") +
+ ylab("Variance of Steps") +
+ ggtitle("Radii vs Variance Steps") +
+ theme_bw()
+
+## Cubic Regression of Variance
+cubic_variance_model = lm(variance_stepcount ~ poly(radii, 3, raw = T), data = room_simulation_data)
+cubic_variance_model_summary = summary(cubic_variance_model)
+ggplot(room_simulation_data, aes(x = radii, y = variance_stepcount)) +
+ geom_point(shape = 21, size = 2) +
+ geom_smooth(method = "lm", formula = y ~ poly(x, 3, raw = T), se = T) +
+ xlab("Radii") +
+ ylab("Variance of Steps") +
+ ggtitle("Radii vs Variance of Steps") +
+ theme_bw()
+
+## Analysis of Cubic Regression
+cubic_variance_coefficients = as.vector(cubic_variance_model$coefficients)
+print(paste("The Cubic Regression Equation is: Variance of Steps =", cubic_variance_coefficients[1], "+", cubic_variance_coefficients[2], "radii", cubic_variance_coefficients[3], "radii^2", cubic_variance_coefficients[4], "radii^3"))
+print(paste(cubic_variance_model_summary$adj.r.squared * 100, "% of the variation of the variations is accounted for by the Cubic Model", sep = ""))Below are the courses that I have documents up online for.
+Probability and Statistical Inference I
+Occassionally during breaks, I like to take some coursera courses. I will post my notes from the lecture here as well.
+ + +My main website is https://brandonrozek.com. Though I like the customization that Wordpress allows, it does come with additional overhead.
+Therefore, I decided to use the Tilda space offered by the Computer Science department to host content pertaining to my Academic Life.
+As such, my research and other class documents will appear here.
+For the curious, this website was built with Pico CMS using the BitsAndPieces theme.
+Here is how you can create a class called "Employee" with a non-default constructor (a constructor that takes parameters) and a getter and setter
+public class Employee {
+ // Our private variables
+ private String name;
+ private double salary;
+ // Non-default constructor
+ public Employee(String name, double salary) {
+ this.name = name;
+ this.salarly = salary;
+ }
+ // This is a getter
+ public string getName() {
+ return name;
+ }
+ public double setSalarly(double salary) {
+ this.salary = salary;
+ }
+}For loops are constructed in the following way
+for (initialization; condition to stop; increment to get closer to condition to stop)
//Assume an array with variable name array is declared before
+for (int i = 0; i < array.length; i++) {
+ // This code will loop through every entry in the array
+}Note that you don't always have to start from zero, you can start anywhere from the array.
+This is an example of how you can take in an array in a method
+public static boolean isEven(int[] array) { // <-- Note the int[] array
+ for (int i = 0; i < array.length; i++) { // Iterate through the entire array
+ // If you find an odd number, return false
+ if (array[i] % 2 == 1) {
+ return false;
+ }
+ }
+ // If you didn't find any odd numbers, return true
+ return true;
+}Let's say that you have the following file
+4
+chicken
+3
+saladAnd you want to make it so that you take the number, and print the word after it a certain number of times. For this example we would want to see the following
+chicken chicken chicken chicken
+salad salad saladHere is the code to write it
+public static void printStrings() {
+ FileInputStream file = new FileInputStream("stuff.txt"); // File contents are in stuff.txt
+ Scanner scnr = new Scanner(file); // Scanner takes in a file to read
+ while (scnr.hasNext()) { // While there are still stuff in the file to read
+ int number = scnr.nextInt(); // Grab the number
+ String word = scnr.next(); // Grab the word after the number
+ // Print the word number times
+ for (int i = 0; i < number; i++) {
+ System.out.print(word);
+ }
+ // Put a new line here
+ System.out.println();
+ }
+
+}Look at handout and carefully trace recursion problems
+Declare a 2D int array with 4 rows and 7 columns
+int[][] dataVals = new int[4][7];A 2D array with 4 rows and 7 columns has 7 * 4 = 28 entries.
+If you want to sum up all the numbers in a 2 dimensional array, do the following
+// Assume numbers is declared beforehand
+int sum = 0;
+for (int i = 0; i < numbers.length; i++) { // Loop through every row in the 2D array
+ for (int j = 0; j < numbers[i].length; j++) { // Loop through every column in a row
+ // This code now looks at one entry in a 2D array
+ sum += numbers[i][j];
+ }
+}Write a Java method to sum values of an array
+// Assume variable numbers is an array of 10 integers
+int sum = 0;
+for (int i = 0; i < numbers.length; i++) {
+ sum += numbers[i];
+}Write a Java method to test if an array contains a specific value
+// Assume variable numbers is an array of 10 integers
+int numWanted = 4;
+boolean found = false;
+for (int i = 0; i < numbers.length; i++) {
+ if (numbers[i] == numWanted) {
+ System.out.println("The number " + numWanted + " was found!");
+ found = true;
+ }
+}
+if (!found) {
+ System.out.println("The number " + numWanted + " was not found!");
+}Write a Java method to find the index of an array element
+// Assume variable numbers is an array of 10 integers
+int numWanted = 4;
+boolean found = false;
+for (int i = 0; i < numbers.length; i++) {
+ if (numbers[i] == numWanted) {
+ System.out.println("The index of number " + numWanted + " is " + i);
+ }
+}How many lines will the following loop print when it is run?
+int i = 0;
+while (i <= 6) {
+ System.out.println("i is now " + (i));
+ i++;
+}i is now 0
+i is now 1
+i is now 2
+i is now 3
+i is now 4
+i is now 5
+i is now 6Last class we talked about reading from files, we can also write to files.
+First you must import all of the necessary libraries
+// To read
+import java.util.Scanner;
+import java.io.FileOutputStream;
+// To Write
+import java.io.FileReader;
+import java.io.PrintWriter;
+// For Exception Handling
+import java.io.IOException;Then in your main, declare a FileOutputStream and PrintWriter
FileOutputStream file;
+PrintWriter print;Create a try block to open a file for writing
+try {
+ // If the file doesn't exist, it'll create it
+ file = new FileOutputStream("output.txt");
+ print = new PrintWriter(file);
+} catch (IOException except) {
+ // Prints out the error message
+ System.out.println("File error " + except.getMessage());
+} Adding a finally block allows the program to clean up before it closes
+ try {
+ file = new FileOutputStream("output.txt");
+ print = new PrintWriter(file);
+ } catch (IOException except) {
+ System.out.println("File error " + except.getMessage());
+ } finally { // It starts here!
+ delete file;
+ delete print;
+ file.close();
+ return;
+ }Then you can write the to file!
+// Do you notice the following methods?
+ print.println("Your number is");
+ print.print("My name is..\n");
+ print.printf("%s %d", "Hello ", 5);
+ print.flush(); //Clears the output stream
+ file.close(); //Closes the file
+Extra Note: Disk fragmentation is a way of cleaning up memory that isn't being used by any of the code on your computer.
+You need to import all the necessary libraries first
+import javax.swing.*;
+import java.awt.FlowLayout;
+import java.awt.event.ActionListener;Your class file needs to extend JFrame that way it can use a bunch of the already existent code written for Swing applications
public class firstGUi extends JFrame {
+ //....Java Swing makes use of what is called Swing Components. These are basic building blocks of GUI items you can add to the screen. For example, a checkbox, a radio button, text field. These are all swing components.
+I wrote a blog post back in the summer which is an overview of them. You can check it out here: https://brandonrozek.com/2017/06/java-swing-components/
+Inside your firstGUI class, declare some Swing components you would like to use in the program
public class firstGUI extends JFrame {
+ JButton button1;
+ JTextArea area;
+ JTextField text;
+ // ....
+You need to create a constructor for this class that way you can initiate all of the swing component values.
+// ...
+JButton button1;
+JTextArea area;
+JTextField text;
+
+// Begin Constructor
+firstGUI() {
+ // Define the components
+ JLabel name = new JLabel("Enter in your name:");
+ text = new JTextField("Jennifer", 20); // 20 characters long, default value: Jennifer
+ area = new JTextArea(10, 10); //Width and Height is 10 characters big
+ JScrollPane sp = new JScrollPane(area); //Adds a scroll bar for the text area
+ button1 = new JButton("Press Me");
+
+ // Set the Layout
+ // FlowLayout organizes each of the components top to bottom, left to right
+ setLayout(new FlowLayout());
+
+ // Add the components to the screen
+ add(name);
+ add(text);
+ add(sp); // never add the textarea when surrounded by a ScrollPane
+ add(button1);
+}Finally, you need to create the Main method which will initiate it all
+public static void main(String[] args) {
+ firstGUI myFrame = new firstGUI();
+
+ // End the program when the 'x' button (not keyboard) is pressed
+ myFrame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
+
+ myFrame.setTitle("My title"); // Titles the window
+ myFrame.pack(); // Packs it all into the frame
+ myFrame.setVisible(true); // Makes the frame appear on the screen
+}You need to change your class header to the following
+public class firstGUI extends JFrame implements ActionListener {
+ // ...Then in your class, add the following method
+@Override public void actionPerformed(ActionEvent event) {
+ // Do stuff as a result of an event here
+ area.append("You Pressed the Button");
+}To make it actually activate as a result of an event. You need to attach it to a swing component.
+For example, I want the code in actionPerformed to activate in result of a button press.
Add the following line to your code in the constructor.
+//...
+button1 = new JButton("Press Me");
+button1.addActionListener(this); // New Code
+//....How do you know which button was pressed in the actionPerformed method?
You can use event.getSource() to find out.
Example:
+@Override public void actionPerformed(ActionEvent event) {
+ if (event.getSource() == button1) { // Replace button1 with appropriate variable name
+ // Do things as a result of a specific button being pressed
+ }
+}To use Swing, do the following steps
+ActionListener to the components you wish to monitorImport the necessary libraries
+import java.awt.*;
+import java.awt.event.*;
+import javax.swing.*;This time we'll be using the null layout. This will require you to set the x, and y positions of each of the elements on the screen
+First we'll make a class called GUILayout which will extend JPanel and contain our layout and components
import java.awt.*;
+import java.awt.event.*;
+import javax.swing.*;
+
+public class GUILayout extends JPanel {
+ private JButton one, two;
+ private JTextField text;
+ private JTextArea area;
+ private int count; // Used to keep track of button pushes
+ private int number = 0;
+
+ GUILayout() {
+ // Set the layout to the null layout
+ setLayout(null);
+ // Tells the computer that you wish to have the size of your program to be 500x500
+ setPreferredSize(new Dimension(500, 500));
+
+ // Create the buttons
+ one = new JButton("Change color");
+ two = new JButton("Count Pushes");
+
+ // Set the x, y, width, height inside parenthesis
+ // This places the buttons on the specified locations on the screen
+ one.setBounds(10, 10, 100, 24);
+ two.setBounds(120, 10, 100, 24);
+
+ // Sets the color of the text of button 1 to blue
+ one.setForeground(Color.BLUE);
+
+ // Add the components to the screen
+ add(one);
+ add(two);
+
+ text = new JTextField("Today is Wednesday, a photographer is here.");
+ text.setBounds(10, 40 ,480, 24);
+ add(text);
+
+ area = new JTextArea(20, 20);
+ // Wrap the text area into the scroll pane
+ // This adds the scroll bars (vertical/horizontal)
+ JScrollPane sp = new JScrollPane(area);
+ sp.setBounds(10, 75, 490, 400);
+ add(sp);
+
+ // Bind the Listener class to the button one
+ one.addActionListener(new Listener());
+
+ // Bind it to two and text
+ two.addActionListener(new Listener());
+ text.addActionListener(new Listener());
+
+ }
+
+ // Create the class for the listener
+ private class Listener implements ActionListener {
+ // Define what you want to occur after an event
+ public void actionPerformed(ActionEvent event) {
+ if (event.getSource() == one) {
+
+ } else if (event.getSource() == two) {
+ count++;
+ // \n is needed to prevent scrolling
+ area.append("You pressed the button " + count + " times!");
+ } else if (event.getSource() == text) {
+ // Grab the number inside the text area and set it to the variable number
+ number = Integer.paseInt(text.getText());
+ number += 10;
+ // add "" to number so that it converts number to a String
+ text.setText(number + "");
+ // Convert the number to a String
+ String temp = Integer.toString(number);
+ area.setText(temp);
+ }
+
+ }
+ }
+}The main code is really short, you just need to extend JFrame and some short code seen last class to set it up.
public class mainTester extends JFrame {
+ public static void main(String[] args) {
+ JFrame frame = new JFrame("Sample GUI"); // Sets the title to "Sample GUI"
+ frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
+ // Instantiate the GUILayout file you made before
+ frame.getContentPane().add(new GUILayout());
+ }
+}You can't put swing and graphics together, therefore you need to make a seperate JPanels for swing and graphics.
+Add necessary libraries
+import java.awt.*;
+import java.awt.Graphics;
+import java.awt.event.*;
+import javax.swing.*;public class drawingWindow extends JPanel {
+ JTextField field;
+ JButton draw;
+ DrawingPanel drawingPanel;
+
+ public drawingWindow() {
+ // Each new component would be vertically stacked upon each other
+ setLayout(new BoxLayout(this, BoxLayout.Y_AXIS));
+
+ JPanel swingSTuff = new JPanel();
+
+ // Add things to the screen
+ draw = new JButton("Draw");
+ field = new JTextField();
+ swingStuff.add(field);
+ swingStuff.add(draw)
+
+ // Add the drawing panel onto the screen
+ drawingPanel = new DrawingPanel(200, 400);
+ add(drawingPanel);
+
+ // Activate the listener if the button was pressed
+ draw.addActionListener(new Listener());
+ }
+
+ // Add the listener to respond to events
+ private class listener implements ActionListener {
+ public void actionPerformed(ActionEvent event) {
+ if (event.getSource() == draw) {
+ drawingPanel.setFlag(1);
+ // Repaints the screen so the oval can appear
+ drawingPanel.repaint();
+ }
+ }
+ }
+
+ // Create the draw panel so we can add it to the screen
+ private class DrawingPanel extends JPanel {
+ private int width, height;
+ DrawingPanel(int width, int height) {
+ this.width = width;
+ this.height = height;
+ setPreferredSize(new Dimension(width, height));
+ }
+ public void setFlag(int flag) {
+ this.flag = flag;
+ }
+ public void paintComponent(Graphics g) {
+ super.paintComponent(g);
+
+ // Every time the flag is set, draw an oval at a random location and color
+ if (flag == 1) {
+ Random rand = new Random();
+ int x = rand.nextInt(width);
+ int y = rand.nextInt(height);
+ g.setColor(Color.RED);
+ g.fillOval(x, y, 20, 30);
+ }
+ }
+ }
+}There are a myriad of different methods you can use.
+// Assume width, height, y, x, etc are defined above
+public void paintComponent(Graphics g) {
+ //....
+ g.dispose(); // Flushes the graphics buffer
+}You have the traditional fill and draw methods. Fill creates the shape shaded in with a color. Draw creates an outline of the shape.
+
+// ...
+g.fillRect(x ,y, width, height);
+ g.drawRect(x, y, width, height);
+ g.fillOval(x, y, width, height);
+ g.drawOval(x, y, width, height);
+ //g.drawPoly(parematers...);
+ //g.fillPoly(parameters...);
+ g.drawArc(x, y, width, height, startingAngle, sweepingAngle);
+ g.fillArc(x, y, width, height, startingAngle, sweepingAngle);You can also create complex shapes like a polygon. When adding points, you need to make sure you add them Clockwise or Counterclockwise (but NOT both)
+ Polygon tri = new Polygon();
+ tri.addPoint(150, 10);
+ tri.addPoint(175, 100);
+ tri.addPoint(125, 100);
+ // Add points clockwise or counterclockwise (NOT BOTH)When doing recursion, make sure not to use loops.
+Recursion is when a function makes a function call to itself.
+It consists of two parts:
+You can have one or more base cases or caller cases.
+To teach recursion, we'll start with a problem that should be written iteratively (with loops) but we'll show how to do it with recursion.
+void CountDown(int number) {
+ while (number > 0) {
+ System.out.println(number + " ");
+ number = number - 1;
+ }
+ System.out.println("Blast Off")
+}Base Case: It stops when the number equals 0
+// Base Case
+if (number == 0) {
+ System.out.println("Blast Off");
+ return;
+}Smaller Caller Case: It decreases number by 1
+// Smaller Caller Case
+System.out.print(number + " ");
+countDownRecursive(number - 1);Put it together...
+void countDownRecursive(int number) {
+ if (number == 0) {
+ System.out.println("Blast Off");
+ } else {
+ System.out.print(number + " ");
+ countDownRecursive(number - 1);
+ }
+}Prints 10 9 8 7 6 5 4 3 2 1 Blast Off
Every time you make a recursive call, it keeps track of all the local variables at the current function call and pushes it onto the stack.
+That means if you make 10 recursive calls, you'll have 10 slots added onto the stack. If you want to return back to the beginning, you would need to return 10 times.
+Whether you do the recursive step first or some other action, it completely changes the output. Look at the following example that's similar to countDownRecursive.
void countUpRecursive(int number) {
+ if (number == 0) {
+ System.out.println("Blast Off");
+ } else {
+ // Notice the swapping of the next two lines
+ countUpRecursive(number - 1);
+ System.out.print(number + " ");
+ }
+}This would print out Blast Off 1 2 3 4 5 6 7 8 9 10
This would be our iterative solution
+int sum(int number) {
+ int sum = 0;
+ while (number > 0) {
+ sum += number;
+ number--;
+ }
+ return sum;
+}How does this loop stop?
+ Same as before. Think about it, if the number you pass in is zero, what should be the result of sum? Zero. Since adding up from 0 to 0 is just 0.
+if (number == 0) {
+ return 0;
+}How does this loop change each time through?
+ You want to update your sum, so return the sum of the next call to the current sum.
+return (number + sum(number - 1));Putting it all together
+int sumRecursive(int number) {
+ if (number == 0) {
+ return 0;
+ } else {
+ return number + sumRecursive(number - 1);
+ }
+}How to do it iteratively.
+void linearSearch(int[] array, int number) {
+ int i = 0;
+ while (i < array.length && number != array[i]) {
+ i++;
+ }
+ if (i == array.length) {
+ System.out.println("Not Found");
+ } else {
+ System.out.println("Found");
+ }
+}Base Case: There are two base cases, when it reaches the end of the array and when it finds the number
+if (array.length == i) {
+ System.out.println("Not Found");
+} else (array[i] == number) {
+ System.out.println(number + " found at index " + i);
+}Smaller Caller Case: Check the next element
+linearSearchRecursion(number, array, i + 1);Putting it all together...
+void linearSearchRecursion(int[] array, int number) {
+ if (array.length == i) {
+ System.out.println("Not Found");
+ } else if (array[i] == number) {
+ System.out.println(number + " found at index " + index);
+ } else {
+ linearSearchRecursion(number, array, i + 1);
+ }
+}This is a much more efficient search than the linear search we have been doing. The only condition is that your array must be sorted beforehand.
+A regular linear search O(n) -- Check at most all of the elements of the array.
Binary Search O(log(n)) -- Checks at most ceil(log(n)) elements of an array.
To demonstrate this with real numbers, let's take an array of 100 elements
+Crazy right?
+Iterative approach
+void binarySearch(int number, int[] array) {
+ int lower = 0;
+ int upper = array.length - 1;
+ int mid = (lower + upper) / 2
+ while (lower <= upper && array[mid] != number) {
+ mid = (lower + upper) / 2;
+ if (array[mid] < number) {
+ lower = mid + 1;
+ } else if () {
+ upper = mid -1;
+ }
+ }
+ if (lower > upper) {
+ System.out.println("Not Found");
+ } else {
+ System.out.println(number + " found at index " + mid);
+ }
+}Recursive Approach
+Base Case: There are two cases, you either found it, or you made it to the end of the array without finding it
+if (lower > upper) {
+ System.out.println("Not Found");
+} else if (array[mid] == number) {
+ System.out.println(number + " found at index " + mid);
+}Smaller Caller Case: Change the lower or upper bounds to cut the search in half
+if (array[mid] < number) {
+ lower = mid + 1;
+ binarySearchRecursive(number, array, lower, upper);
+} else if (array[mid] > number) {
+ upper = mid - 1;
+ binarySearchRecursive(number, array, lower, upper);
+}Putting it together....
+binarySearch(int number, int[] array, int lower, int upper) {
+ if (lower > upper) {
+ System.out.println("Not Found");
+ } else if (array[mid] == number) {
+ System.out.println(number + " found at index " + mid);
+ }
+ else if (array[mid] < number) {
+ lower = mid + 1;
+ binarySearchRecursive(number, array, lower, upper);
+ } else if (array[mid] > number) {
+ upper = mid - 1;
+ binarySearchRecursive(number, array, lower, upper);
+ }
+}Continuing from last week's topic, writing a function called compare in your custom class.
For example, for your project:
+public bool compare(ResearchProject right, int type) {
+ if (type == 0) { // Is last lastname greater than right lastname alphabeticallly
+ if (lastname.compareTo(right.lastname) > 0) {
+ return true;
+ } else {
+ return false;
+ }
+ } else if (type == 1) { // Implement a different type of comparison
+
+ }
+}You can then use the sorting algorithms discussed previously to sort the ResearchProjects using the compare method.
First start by importing the required libraries.
+import java.io.FileInputStream;
+import java.io.IOException; // For error handling
+import java.io.FileNotFoundException; // For error handling
+import java.util.Scanner;Then inside your main method which throws an IOException (see example "Reading the filename from the keyboard")
+FileStream file;
+Scanner in;
+try { // Try something
+ file = new FileInputStream("test.txt");
+ in = new Scanner(file);
+} catch (FileNotFoundException e) { // Catch IF it fails
+ System.out.println("File not found.");
+ in = new Scanner(System.in); // Read from the keyboard instead
+}
+// Read the file now....
+String name = in.nextLine();Before we had linked the Scanner to the keyboard, now we're doing it to a file.
+In the code above you see what's called a try-catch block. That means that if something was to happen in the execution of the code. Instead of just crashing as we have been seeing normally. We can deal with the error. In the example above you saw that the program reads from the keyboard instead of the file
+public static void main(String[] args) throws IOException {
+ FileInputStream file;
+ Scanner in, scnr;
+ // Connect one scanner to the keyboard to get the filename
+ scnr = new Scanner(System.in);
+
+ System.out.print("Enter file name: ");
+ String filename = in.nextLine();
+
+ // Repeat code from previous example
+ FileStream file;
+ Scanner in;
+ try {
+ file = new FileInputStream(filename); // Only this line changed
+ in = new Scanner(file);
+ } catch (FileNotFoundException e) {
+ System.out.println("File not found.");
+ in = new Scanner(System.in);
+ }
+ String name = in.nextLine();
+ // ....
+}The main throws the IOException since we don't deal with it in any of the catch blocks.
+Let's say we have a file with the following in it
+3
+Paul
+Peter
+Mary3 is to indicate the number of names in th file.
+Let's write code in order to read these names into an array!
+import java.io.FileInputStream;
+import java.io.IOException;
+import java.io.FileNotFoundException;
+import java.util.Scanner;
+public static void main(String[] args) throws IOException {
+ FileInputStream file;
+ Scanner in, scnr;
+ scnr = new Scanner(System.in);
+
+ System.out.print("Enter file name: ");
+ String filename = in.nextLine();
+
+ FileStream file;
+ Scanner in;
+ try {
+ file = new FileInputStream(filename); // Only this line changed
+ in = new Scanner(file);
+ } catch (FileNotFoundException e) {
+ System.out.println("File not found.");
+ in = new Scanner(System.in);
+ }
+
+ // For the size of the array, get the first number in the file
+ int size = in.nextInt();
+ String[] nameArray = new String[size];
+ for (int i = 0; i < size; i++) {
+ namesArray[i] = in.nextLine();
+ }
+
+ // Now all the names in the file is in the namesArray
+}Global variables is where you don't declare a variable inside a method. This is generally not a recommended practice. It is recommended to declare variables inside methods so that it is easier to reuse code.
+public class mainDriver {
+ public static int globalVariable = 5; // I am a global variable
+ public static void main(String[] args) {
+ int localVariable = 4; // I am a local variable
+ }
+}You can format strings in java by using the printf method in System.out.
Format strings work by using placeholders starting with % to specify where to place the value of a variable.
General format of command
+//System.out.printf(FormatString, variable1, variable2, ....)
+String s = "The number is";
+int x = 5;
+System.out.printf("%s %d", s, x); // Prints "The number is 5"If you want to print out money, you can do it through the following
+float tax = 0.45698;
+System.out.printf("The tax is %.2f"); //prints "The tax is 0.46"Due to how computers store non-integers, math can be non-precise after some mathematical operations. It is therefore advantageous to do operations on integers to the level of precision you want.
+For example, instead of working in dollars, work in pennies since it's the smallest division we care about.
+Standard arrays are static, meaning they have no ability to grow. Instead of doing the operations we described last class, in order to add something to an array. We can use something called an ArrayList instead.
The methods in the ArrayList library are useful abstractions to make the life of a programmer easier.
ArrayLists can also hold only one type. The type cannot be a primitive like a standard array. Instead it must be a class representing the desired type.
+int -> Integer
+double -> Double
+char -> Character
+float -> String
+To declare and initialize an ArrayList
import java.util.ArrayList;
+public class ArrayListTest {
+ public static void main(String[] args) {
+ ArrayList<Integer> numbers = new ArrayList<Integer>();
+ ArrayList<String> names = new ArrayList<String>();
+ }
+}ArrayLists has a variety of different methods that can be used to access/manipulate it
If you want to add the numbers 1 through 10 into the ArrayList of numbers
for (int i = 1; i < 11; i++) {
+ numbers.add(i);
+}To print out the entire contents of the ArrayList
for (int i = 0; i < numbers.size(); i++) {
+ System.out.println(numbers.get(i));
+}How can you replace a value?
+int n = 5; // Make this the number you wish to replace
+int i = 0;
+while (i < numbers.size() && numbers.get(i) != n) {
+ i++;
+}
+if (i == numbers.size()) {
+ System.out.println(n + " not found.");
+} else {
+ int r = 10; // Make this the value you want to replace n with
+ numbers.set(i, r);
+}The remove method removes an item given an index. This shifts all the elements above the removed index down.
numbers.remove(3); // Removes the element in the third indexThe add method can also take an index. It pushes all the elements at the index specified up.
numbers.add(3, 5); // Adds the number 5 to the third indexYou can clone an array using the clone method
ArrayList<Integer> numbers2 = new ArrayList<Integer>();
+numbers2.clone(numbers); // Clones numbers into numbers2public class Numbers {
+ private ArrayList<Integer> used;
+ private ArrayList<Integer> unused;
+ numbers () {
+ // Debug: System.out.println("1. Constructor Entry Point");
+ used = new ArrayList<Integer>();
+ unused = new ArrayList<Integer>();
+ // Debug: System.out.println("2, Constructor Size of ArrayLists" + used.size() + " " + unused.size())
+ }
+ // Adds the numbers 1-5 into the used ArrayList
+ public void fillUsedArrayList() {
+ for (int i = 0; i < 5; i++) {
+ used.add(i + 1);
+ }
+ }
+ // Move an item from the unused ArrayList to the used ArrayList
+ public int moveIt(int index) {
+ int temp = used.get(index);
+ unused.add(temp);
+ // Debug: System.out.println("Adding to used array:" + (i + 1));
+ used.remove(index);
+ return temp;
+ }
+ // The method below is created for debugging purposes
+ public void printBothArrayLists() {
+ // Print the used arraylist
+ System.out.println("Used ArrayList");
+ for (int i = 0; i < used.size(); i++) {
+ System.out.println(used.get(i));
+ }
+
+ // Print the unused arraylist
+ System.out.println("Unused ArrayList");
+ for (int i = 0; i < unused.size(); i ++) {
+ System.out.println(unused.get(i));
+ }
+ }
+}Recall that you can compile the code above but you cannot run it. To run code, you must have a main method.
+public class NumberTester {
+ public static void main(String[] args) {
+ Numbers list;
+ list = new Numbers();
+ list.fillUsedArrayList();
+ list.printBothArrayLists();
+ }
+}An Array is a static structure of contiguous memory of a single type.
+An ArrayList is a dynamic structure of contiguous memory of a single type
+To get the size of an array you use .length
To get the size of an ArrayList you use .size()
array is not a reserved word, it's a concept. Arrays are able to hold multiple values under one name of the same type.
For instance, you can have an array of integers.
+Properties of an array
+You can declare arrays by saying the type '[]' and the name of the array
+int[] numbers;
+double[] gpas;
+float[] grades;Before you can use the array, you must new it
numbers = new int[10];Where 10 is the size of the array.
+You can combine both the declaration and initialization
+double[] points = new double[7];You can access individual elements of the array by using its index. Indexes start from zero
+points[0] = 5.4; // The first element in the point array is 5.4The .length property of an array gives the size of the array
You can print out each element in the array using a for loop
+for (int i = 0; i < numbers.length; i++) {
+ System.out.println(numbers[i]);
+}You can ask a user to input a value to each element in the array
+for (int i = 0; i < points.length; i++) {
+ System.out.print("Enter a number: ");
+ points[i] = scnr.nextInt();
+}You can use a while loop to search the array
+int i = 0;
+int number = 5;
+// While the index is within the array size and the number isn't found
+while (i != number.length && number != numbers[i]) {
+ i++
+}
+if (i == numbers.length) {
+ System.out.println(number + " was not found.")
+} else {
+ System.out.println(number + " was found at index " + i)
+}If you don't include the i != number.length you will obtain an IndexOutOfBounds error.
The example above is called a Linear Search.
+Linear searches work on an unsorted and sorted arrays.
+You can pass an array into a method
+public static void exampleMethod(int[] sample) {
+ // Do something
+}
+public static void main(String[] args) {
+ int[] s = new int[30];
+ exampleMethod(s);
+}For-loops can run 0 or more times. If you want something to execute at least once. Use a do-while loop.
+do {
+ // Code
+} while (condition);For example, to search at least once and asking whether the user wants to search again
+// Assume linearSearch and array are defined
+char answer;
+Scanner input = new Scanner(System.in);
+do {
+ linearSearch(array, input);
+ System.out.print("Do you want to search again? (Y/N) ");
+ input.nextLine();
+ answer = input.next().charAt(0);
+} while( answer != 'N');You can create any type of loop just by using a while loop.
+You can find the max of an array using the following method
+double max = arrayName[0];
+for (int i = 1; i < arrayName.length; i++) {
+ if (max < arrayName[i]) {
+ max = arrayName[i];
+ }
+}
+System.out.println("The max is " + max);You can sum the array using the following method
+double sum = 0;
+for (int i = 0; i < arrayName.length; i++) {
+ sum += arrayName[i];
+}
+System.out.println("The sum is " + sum);You can think of a two dimensional array as a grid that is organized by rows then columns.
+To declare a two dimensional array do the following
+int[][] grid = new int[5][5]; // Declares a 2D array of 5 rows and 5 columnsYou can have as many dimensions as you want. For graphics a 3-dimensional array would render 3D points.
+It doesn't have to be inherently visual though, you can have the n-dimensional array look at the interaction between n different variables. For example, the relationships to different questions in a survey.
+Strings are essentially a char array with extra methods attached. We can imitate an array of strings with a 2D char array.
+char[][] helloWorld = new char[5][5];
+hello[0][0] = 'h';
+hello[0][1] = 'e';
+hello[0][2] = 'l';
+hello[0][3] = 'l';
+hello[0][4] = 'o';
+hello[1][0] = 'w';
+hello[1][1] = 'o';
+hello[1][2] = 'r';
+hello[1][3] = 'l';
+hello[1][4] = 'd';
+To access the elements in a 2D array, you need to use a nested for-loop.
+Here is how you print out the hello world example above
+for (int row = 0; row < helloWorld.length; row++) {
+ for (int col = 0; col < helloWorld[row].length; col++) {
+ System.out.print(helloWorld[row][col]);
+ }
+ System.out.print(" ")
+}The code above prints out "hello world"
+You can write a get like method in the following way
+public static void get(int[][] array, int row, int col) {
+ return array[row][col];
+}Arrays in Java are pass by reference not pass by value. Meaning that if you change the array within the method then it will change outside the method.
+public static void insert(int[][] array, int row, int col, int numToInsert) {
+ array[row][col] = numToInsert;
+}To make it not be able to change the array inside the method, use the keyword const inside the method header. To code below will throw a compiler error.
public static void insert(const int[][] array, int row, int col, int numToInsert) {
+ array[row][col] = numToInsert; // This line will throw an errror
+}You can not only do a two dimensional array of primitive types, but you can also do two dimensional arrays of objects/classes.
+animalLocation[][] map;
+map = new animalLocation[5][4];Since we are dealing with classes, you cannot use this array right away. The code above creates the space in memory to store the objects. To have the animalLocation objects in the array, you must new each instance of the object.
for (int i = 0; i < map.length; i++) {
+ for (int j = 0; j < map[i].length; j++) {
+ map[i][j] = new animalLocation();
+ }
+}These instructions are to sort in descending order, to sort in ascending order just negate the condition.
+This sort is a series of iterations, for each iteration you go to the bottom of the array. Then compare if the value is greater than the element before it. If it is, then you
+Once you go through an iteration, the last thing being compared is the greatest value of the entire array. That means you don't have to check it every time anymore.
+Keep going through all the iterations until n, where n is the size of the array, iterations have been completed.
+If you try to swap variables by saying
+y = x;
+x = y;then you'll end up overwriting y's value with x and both variable would have x's value.
+If you want to actually swap variables, you must create a temporary variable that saves y's value so that it can be properly assigned to x.
+int temp;
+temp = y;
+y = x;
+x = temp;// Each iteration
+for (int j = 0; j < array.length - 1; j++) {
+ // Each element in the list
+ for (int i = 0; i < array.length - 1; i++) {
+ // Is the element greater than the one after it?
+ if (array[i] > array[i + 1]) {
+ // Swap
+ int temp = array[i + 1];
+ array[i + 1] = array[i];
+ array[i] = temp;
+ }
+ }
+}This here compares each of the values each time. If you remember before, I said that you don't have to compare the topmost each time.
+To change this, just change the second loop condition
+// Each iteration
+for (int j = 0; j < array.length - 1; j++) {
+ // Each element in the list
+ for (int i = 0; i < array.length - 1 - i; i++) { // Note this line
+ // Is the element greater than the one after it?
+ if (array[i] > array[i + 1]) {
+ // Swap
+ int temp = array[i + 1];
+ array[i + 1] = array[i];
+ array[i] = temp;
+ }
+ }
+}In Java, you can compare numbers, strings, and even your own customized objects. To compare your own customize object, you must write a method called compare in your class.
// Each iteration
+for (int j = 0; j < array.length - 1; j++) {
+ // Each element in the list
+ for (int i = 0; i < array.length - 1 - i; i++) {
+ // Is the element greater than the one after it?
+ if (array[i].compare(array[i + 1])) { // Note this line
+ // Swap
+ int temp = array[i + 1];
+ array[i + 1] = array[i];
+ array[i] = temp;
+ }
+ }
+}When you pass a primitive type into a method (int, char, double, float, etc), it makes a copy of the value of the variable and brings it into the method
+When you pass an array into a method (int[], char[], double[], etc[]), it passes in the reference of the variable into the method. In other words, you give the actual array into the method and allows the method to change it.
+If you change the primitive in a method, it doesn't actually change the value of the variable.
+If you pass in an array and change it in the method, it has been permanently changed outside the method as well.
+Use the finalkeyword in the method header
public static void printAll(final int[] array) {
+ for (int i = 0; i < array.length; i++) {
+ System.out.println("Number " + (i + 1) + " is " + array[i])
+ }
+}String[] names = {"Jennifer", "Noodle", "Fluffy", "Rosie", "Cinnamon", "Brianne", "Oliver"}Everything between the {} is the initial values in the names array in the order that it is written.
Recall that arrays are of a fixed size. The names array above has 7 elements.
Do the following steps:
+// (1)
+String[] temp = new String[7];
+// (2)
+temp.clone(names);
+// (3)
+names = new String[20]; // Now it can hold up to 20 names
+// (4)
+names.clone(temp);
+// (5)
+names[7] = "New name!";public class car { // begin car
+ private int year;
+ private String make;
+ public int getYear(void) {
+ return year;
+ }
+ public String getMake() {
+ return make;
+ }
+ public void setYear(int y) {
+ if (y > 1890) {
+ year = y;
+ } else {
+ System.out.println(y + " is not a valid year.");
+ }
+ }
+ public void setMake(String m) {
+ make = m;
+ }
+}Local variables are only existent within the curly braces that it is defined in.
+Boolean expressions return a boolean
+1 < 4; // 1 is less than 4: TRUE
+3 > 5; // 3 is greater than 5: FALSE
+5 == 5; // 5 is equal to 5: TRUE
+5 != 5; // 5 is not equal to 5: FALSE
+1 >= 1; // 1 is greater than or equal to 1: TRUE
+5 <= 1; // 5 is less than or equal to 1: FALSEIf statements only occur if the boolean expression is true, otherwise the else block is executed.
if (true) {
+ System.out.println("I am always printed");
+} else {
+ System.out.println("I am never printed");
+}You can only have one else per if. If you have an if you don't necessarily need an else
If you have a local variable and the class variable sharing the same name, then the local variable is always used first.
+public class car { // begin car
+ private int year;
+ public void setYear(int year) {
+ year = year;
+ }
+}This is a redundant statement, it makes the argument that is passed in equal to itself.
+To avoid this situation, use the keyword this to access the class variable
public class car {
+ private int year;
+ public void setYear(int year) {
+ this.year = year;
+ }
+}The code above runs as expected.
+Rewriting our class with this
public class car { // begin car
+ private int year;
+ private String make;
+ public int getYear(void) {
+ return year;
+ }
+ public String getMake() {
+ return make;
+ }
+ public void setYear(int year) {
+ if (y > 1890) {
+ this.year = year;
+ } else {
+ System.out.println(y + " is not a valid year.");
+ }
+ }
+ public void setMake(String make) {
+ this.make = make;
+ }
+}When the code hits a return statement, it stops executing the rest of the code in the method. Also throws an Unreachable Code Error.
public int add(int x, int y) {
+ return x + y;
+ System.out.println("x + y = " + x + y);
+}
+add();
+System.out.println("Hello");Here the code above will not compile, though assuming the error doesn't exist then it would only print out "Hello"
+You cannot have a private or protected constructors.
+Constructors are used to initialize your objects.
+You want to have the class variables to the left of the assignment statement.
+public class car {
+ private int year;
+ private String make;
+ car() {
+ year = 1890;
+ make = "Ford";
+ }
+ car(int year, String make) {
+ this.year = year;
+ this.make = make;
+ }
+}Testers are useful to check that the class is implemented correctly. Both the tester and the class have to be in the same folder/directory.
+public class carTester {
+ public static void main(String[] args) {
+ Car myCar; // Declaration
+ myCar = new Car(); // Initilization
+ Car yourCar = new Car(2009, "Hyundai"); // Declaration + Initialization
+ }
+}public class Car {
+ private String name;
+ private int odometer;
+ public void setOdometer(int od) {
+ odometer = od;
+ }
+ public void setName(String n) {
+ this.name = n;
+ }
+ public void changeOilRequest(String name, int od) {
+ if (name == this.name) {
+ int difference = od - this.odometer;
+ if (difference > = 3000) {
+ // You can call other methods in the class
+ setOdo(od); // Equivalent to "this.setOdo(od);"
+ this.odometer = od;
+ System.out.println("Ready for oil change.");
+ } else {
+ System.out.println(name + " not ready for oil change.")
+ }
+ } // end if
+ } // end changeOil request
+} // end classTo call public methods outside the class use the variable name to do so.
+public class CarTester {
+ public static void main(String[] args) {
+ Car myCar = new Car();
+ myCar.setName("Honda")
+ myCar.changeOilRequest("Honda", 3400);
+ }
+}The ceil method rounds up while the floor method runs down.
Math.ceil(3.2); // 4
+Math.ceil(4.1); // 4Looking at the following example code
+int i;
+for (i = 0; i < 3; i++) { //begin for
+ System.out.println("i = " + i); //body
+} //end for
+System.out.println("After loop, i = " + i);i = 0 is the initializing statement
i < 3 is the conditional, that is when the loop ends
i++ is the increment/decrement
i++ is synonymous with i = i + 1
The initialization statement only occurs once at the beginning of the loop.
+Let us go through this for loop example
+i = 0i < 3? Yes execute the body
+"i = 0"i ++, i is now 1i < 3? Yes, 1 is less than 3. Execute body
+"i = 1"i++ i is now 2i < 3? Yes 2 is less than 3. Execute body
+"i = 2"i++, i is now 3i < 3? No 3 is not less than 3
+Exit loop. Print "After loop, i = 3"
You can also do the declaration in the initialization statement
+for (int i = 0; i < 3; i++) {
+ System.out.println("i = " + i);
+}This now runs like above without the "After loop, i = 3" print. You cannot access the variable i outside the for loop since in this example, it belongs to the for loop's scope.
With the AND operator && both the left and right side needs to be true for the expression to be true.
true && true // true
+true && false // false
+false && true // false
+false && false // falseWith the OR operator || either the left or right side needs to be true for the expression to be true.
true || true // true
+true || false // true
+false || true // true
+false || false // falseExample: Print out the number n if it is between 10 and 20 (inclusive)
if (n >= 10 && n <= 20) {
+ System.out.println(n);
+}Example: Print out the age if it is not of young adult age. Young adult range is from 18 to 39 (inclusive)
if (!(age >= 18 && age <= 39)) {
+ System.out.println(age);
+}Or you can use De Morgan's Law (for the curious)
+if (age < 18 || age > 39) {
+ System.out.println(age);
+}You can use the loop to count backwards
+for (int i = 10; i > -1; i--) {
+ System.out.println(i);
+}This prints the following
+10
+9
+8
+7
+6
+5
+4
+3
+2
+0You can make rows and columns of asterisks
+for (int j = 0; j < someNumber; j++) { // Corresponds to rows
+ for (int i = 0; i < someNumber2; i++) { // Corresponds to columns
+ System.out.print("*");
+ }
+ System.out.println(""); // Goes to the next row
+}If someNumber equals someNumber2, then we have the same amount of rows as columns.
Let someNumber equal to 2 and someNumber2 equal to 2
Output:
+**
+**You can make a right triangle of Tilda with the following code
+for (int i = 1; i <= num; i++) { // Corresponds to the row
+ for (int j = 0; j < i; j++) { // Corresponds to the column and stops at the current row number
+ System.out.print("~");
+ }
+ System.out.println(""); // Moves to next row
+}Adding the keyword final in front of a variable declaration makes the variable constant. Meaning you cannot later change it in the code.
final int MAX = 10;By industry norm, the style of the variable is to have it in all caps.
+You CANNOT do the following
+final int MAX = 10;
+MAX = 15;Example:
+import java.util.Math;
+public class MathTest {
+ public static void main(String[] args) {
+ double answer = Math.ceil(5.4);
+ System.out.println(Math.ceil(4.5));
+ }
+}You can only cast a variable if you are casting it to a type that is larger than the one it was previously. The expression Polack used is that you cannot fit into super skinny jeans, but you can fit into bigger pants.
+double dnum;
+float fnum;
+int inum;
+dnum = (double)fnum * (double)inum;Strings are initialized in Java with double quotes while chars are initialized with single quotes
char initial = 'j';
+String name = "Jennifer";Characters can be read in as an integer.
+java.util.Random;newimport java.util.Random;
+public class RandTest {
+ public static void main(String[] args) {
+ Random rand;
+ rand = new Random();
+ int number = rand.nextInt(100); // Random generates number between 0-99
+ }
+}How do you generate numbers in a different range? [50, 150]
+rand.nextInt(100); // 0-99
+rand.nextInt(101) // 0 - 100
+rand.nextInt(101) + 50 //50-150In more general terms
+rand.nextInt(end - start + 1) + startVariable -- Storage of information
+The type cannot change in a variable.
+Examples of types include
+Declaration: int num;
Initialization: num = 5;
Declaration + Initialization: int num = 5;
You cannot declare a variable multiple times.
+Undefined variables are when you do not declare a variable before attempting to use it.
+You need to cast if you are attempting to lose data or store a larger memory type into a smaller one.
+double -> float -> int (casting required)
+double gpa = 3.2;
+int num1 = 10 * (int)gpa // 30The basic number operations are
+Examples
+0 % 2 // 0
+1 % 2 // 1
+2 % 2 // 0
+
+3 % 2 // 1
+4 % 2 // 0
+5 % 2 // 1
+
+3 % 5 // 3
+7 % 5 // 2You can test if something is even using modulus %
+// Assuming i was initiliazed to a value earlier
+if (i % 2 == 0) {
+ System.out.println("i is even");
+} else {
+ System.out.println("i is odd");
+}Here is sample code using a Scanner as input
+import java.util.Scanner;
+public class ScannerExample {
+ public static void main(String[] args) {
+ Scanner in;
+ in = new Scanner(System.in);
+
+ // Grab numerical values
+ int num = in.nextInt();
+ float gpa = in.nextFloat();
+ double weight = in.nextDouble();
+
+ // Grab a single character
+ in.nextLine()
+ char initial = in.next().charAt(0);
+
+ // To get the entire line of a string
+ in.nextLine();
+ String name = in.nextLine();
+ }
+}You need to use in.nextLine() to grab the carriage return that is left after grabbing a numeric value.
Classes are a new type that you can have multiple things of.
+These classes are blueprints that are made up of primitives or more basic types.
+First create a Pet.java file (Name of the class must match the name of the file)
+public class Pet {
+ private String name;
+ private int years;
+}You can then use the Pet class in your main program. The terminology here is that you can create instances or objects of the class.
+In PetTester.java
+public class PetTester {
+ public static void main(String[] args) {
+ Pet myPet;
+ myPet = new Pet();
+ }
+}Both Pet.java and PetTester.java must be in the same directory/folder
+Since the variables are private we cannot access them in the main program. To work around this, we can write what is called a mutator method.
+public class Pet {
+ private String name;
+ private int years;
+
+ // Mutators
+ public void setName(String n) {
+ name = n;
+ }
+ public void setYears(int y) {
+ if (y >= 0) {
+ years = y;
+ } else {
+ System.out.println("No one is less than 0 years old.")
+ }
+ }
+}Now let's use these new methods
+public class PetTester {
+ public static void main(String[] args) {
+ Pet mypet;
+ myPet = new Pet();
+ myPet.setName("Fred");
+ myPet.setAge(20);
+ }
+}We need a method that will allow us to access the data type. Let's add accessors to our pet class.
+public class Pet {
+ private String name;
+ private int years;
+
+ // Mutators
+ public void setName(String n) {
+ name = n;
+ }
+ public void setYears(int y) {
+ if (y >= 0) {
+ years = y;
+ } else {
+ System.out.println("No one is less than 0 years old.")
+ }
+ }
+
+ // Accessors
+ public String getName() {
+ return name;
+ }
+ public int getYears() {
+ return years;
+ }
+}Now let's get some information from the pet object, such as the age.
+public class PetTester {
+ public static void main(String[] args) {
+ Pet mypet;
+ myPet = new Pet();
+ myPet.setName("Fred");
+ myPet.setYears(20);
+
+ int year = myPet.getYears();
+ }
+}Constructors lets us initialize variables in the class without having to use mutator methods.
+public class Pet {
+ private String name;
+ private int years;
+
+ // Default Constructor
+ public Pet() {
+ name = "";
+ years = 0;
+ }
+ // Non-Default Constructor
+ public Pet(int y, String n) {
+ name = n;
+ years = y;
+ }
+
+ // Mutator Methods
+ public void setName(String n) {
+ name = n;
+ }
+ public void setYears(int y) {
+ if (y >= 0) {
+ years = y;
+ } else {
+ System.out.println("No one is less than 0 years old.")
+ }
+ }
+
+ // Accessor Methods
+ public String getName() {
+ return name;
+ }
+ public int getYears() {
+ return years;
+ }
+}Now let us see this in action.
+public class PetTester {
+ public static void main(String[] args) {
+ Pet yourPet = new Pet(10, "Fluffy");
+ }
+}You can have as many constructors as you want, but they must be different.
+Example:
+public class Pet {
+ ...
+ pet() {
+ name = "";
+ years = 0;
+ }
+ pet(int y, String n) {
+ name = n;
+ years = y;
+ }
+ pet(String n) {
+ years = 1;
+ name = n;
+ }
+ ...
+}The base class, super class, or parent class is the initial class that we are working with. Let's say that you want to extend the class, or add additional functionality. The class that inherits from the parent class is called the child class, subclass or derived class.
+public class Truck extends Car {
+ // Truck Appropriate Fields
+ // Necessary methods for truck
+}This code adds all the methods from Car into the Truck class. You can then add methods that is specific to a Truck into the Truck class.
+A child class has all parent fields and access to all parent methods!
+Recall the words public and private
The public modifier makes the field/method accessible by any class
The private modifier makes the field/method only accessible within the method itself
The protected modifier makes the field/method accessible within the same class or any subclasses.
+You can override a parent class method by declaring a method in the child class with the same...
+but this method would have different behavior!
+In this class we will talk about three types of control structures
+Sequential is what is most familiar to us. Write the lines from top to bottom and it executes it in that order
+Selection depends on the question of if.
If it is raining, wear boots
+if (raining) {
+ wearingBoots = true;
+}If you want something to happen also when it is not true, consider an if-else statement
If the light is off, turn it on.
+Otherwise, turn it on
+if (lightIsOn) {
+ lightIsOn = false;
+} else {
+ lightIsOn = true;
+}Sometimes you can have multiple branches depending on a condition. Let us take a stop light as an example
+if (light == "red") {
+ car.stop()
+} else if (light == "yellow") {
+ car.slow()
+} else {
+ car.go()
+}There is a specific method in the String class when it comes to checking for string equality
boolean equals(String s)Let us look at an example
+String word = "hello";
+boolean ans = word.equals("hello"); // Returns true
+boolean ans2 = word.equals("Hello"); // Returns falseWhile some check is true, repeat the work.
+While the cookies aren't baked, keep baking
+while (expression) {
+ // Loop body executes if expression is true
+}
+// Statements execute after expression is false// Ask user about first song
+while (user says play next song) {
+ // play next song
+ // ask user about next song
+}You can have loops inside loops
+int outer = 1;
+while (outer < 4) {
+ int inner = 1;
+ while (inner < 4) {
+ System.out.println(outer + ":" + inner);
+ inner++;
+ }
+ outer++;
+}This code does the following
+1:1
+1:2
+1:3
+2:1
+2:2
+2:3
+3:1
+3:2
+3:3Never write the entire program at once! This makes it incredibly hard to debug. Instead break it into small parts.
+Write one part -> debug until it works
+Write second part -> debug until it works
+This way you know which part of your code failed, instead of having everything fail.
+You can get input from a file instead of from the terminal
+FileInputStream fileIn = new FileInputStream("myFile.txt");
+// Our familiar Scanner
+Scanner scnr = new Scanner(fileIn);
+// We can use our usual Scanner methods
+String line = scnr.nextLine();
+fileIn.close(); // Remember to close the file when you're finished with it!To understand some of the Scanner methods we need to be aware of the "newline" character. This character is equivalent to the Enter button on the keyboard.
scnr.nextLine() This get's all the characters in the buffer up to the newline character.
scnr.next() Grabs the characters in the next "token". Tokens are usually separated by any whitespace type character (spaces, enters, tabs, etc.)
Prints information to a file instead of to the screen
+FileOutputStream fileOut = new FileOutputStream("myOutfile.txt");
+PrintWriter out = new PrintWriter(fileOut);
+out.println("Print this as the first line.");
+out.flush(); // Pushes the file changes to the file
+fileOut.close(); // If you forget this then it won't remember your changesArrays are containers of fixed size. It contains a fixed number of values of the same type. (Ex: 10 integers, 2 strings, 5 booleans)
+Declaration
+int[] array; // This declares an integer arrayInitialization
+array = new int[7]; // This states that this array can hold up to 7 integersStoring a value in an array
+int[] array = new int[7];
+array[0] = 5; // Stores 5 into the first slotNow let us attempt to retrieve the value
+int temp = array[0];
+System.out.println(temp); // Prints "5"Let's say we have the following array
+int[] numbers = {3, 5, 2, 7, 9};Let's print out each of the values in the array
+for (int i = 0; i < numbers.length; i++) {
+ System.out.print("value in " + i " is " + numbers[i]);
+}int highest = numbers[0];
+for (int i = 0; i < numbers.length; i++) {
+ if (numbers[i] > highest) {
+ highest = numbers[x];
+ }
+}Chapter 1 -- Code Style, API
+Chapter 2 -- Variables & Assignments, strings
+Chapter 3 -- input & output
+Chapter 4 -- branches (if, if/else, switch)
+Chapter 5 -- loops (while, for), scope
+Chapter 6 -- File Reading and Writing
+What is the output of this code?
+String preferredLanguage = "Spanish";
+if (preferredLanguage.equals("Chinese")) {
+ System.out.println("Ni hao!");
+}
+if (preferredLanguage.equals("Spanish")) {
+ System.out.println("Hola!");
+}
+if (preferredLanguage.equals("French")) {
+ System.out.println("Bonjour!");
+}
+if (preferredLanguage.equals("German")) {
+ System.out.println("Gutentag!")
+} else {
+ System.out.println("Hello!")
+}The output is
+Hola!
+Hello!This is because each of the if statements are independent from each other. Whether or not the if statement gets check is not affected by the if statements around it.
+Since the preferred language equals Spanish it outputs Hola! But since the language is also not German it prints out Hello! as well.
Square brackets notation is used to access elements, array slots can be used as variables
+int[] array = new int[7]; // Creates an integer array of size 7
+array[0] = 5;You can swap x and y in the following way with a temporary variable
int x = 6;
+int y = 1;
+
+int temp = x;
+
+x = y;
+y = temp;// Creates a 2D array of two rows and three columns
+int[][] a = new int[2][3]You can access an element of this 2D array using the conventional square bracket notation
+a[0][0] = 5;Inside the parenthesis of the if statement must be a boolean expression. This is an expression that evaluates to either true or false. We can do more complex boolean expressions through logical operators.
NOT !a this is true when a is false
AND a && b this is true when both operands are true
OR a || b this is true when either a is true OR b is true
Let's look at the most simplest case. Not.
+| a | +!a | +
|---|---|
| true | +false | +
| false | +true | +
| a | +b | +a && b | +
|---|---|---|
| true | +true | +true | +
| true | +false | +false | +
| false | +true | +false | +
| false | +false | +false | +
Notice here that a && b is only true when both a and b are true.
| a | +b | +a \ | +\ | +b | +
|---|---|---|---|---|
| true | +true | +true | +||
| true | +false | +true | +||
| false | +true | +true | +||
| false | +false | +false | +
Notice here that a || b is only false when both a and b are false.
| + | + |
|---|---|
| Parenthesis | +() |
+
| Logical Not | +! |
+
| Arithmetic Operators | +* / % + - |
+
| Relational Operators | +< <= > >= |
+
| Equality and Inequality operators | +== != |
+
| Logical AND | +&& |
+
| Logical OR | +|| |
+
| a | +b | +!b | +a && !b | +
|---|---|---|---|
| true | +true | +false | +false | +
| true | +false | +true | +true | +
| false | +true | +false | +false | +
| false | +false | +true | +false | +
| a | +b | +!a | +!a \ | +\ | +b | +
|---|---|---|---|---|---|
| true | +true | +false | +true | +||
| true | +false | +false | +false | +||
| false | +true | +true | +true | +||
| false | +false | +true | +true | +
| a | +b | +c | +b && c | +a \ | +\ | +(b && c) | +!(a \ | +\ | +b && c) | +
|---|---|---|---|---|---|---|---|---|---|
| true | +true | +true | +true | +true | +false | +||||
| true | +true | +false | +false | +true | +false | +||||
| true | +false | +true | +false | +true | +false | +||||
| false | +true | +true | +true | +true | +false | +||||
| true | +true | +false | +false | +true | +false | +||||
| true | +false | +true | +false | +true | +false | +||||
| false | +true | +true | +true | +true | +false | +||||
| false | +false | +false | +false | +false | +true | +
| a | +b | +c | +!a | +b && c | +!a \ | +\ | +b && c | +
|---|---|---|---|---|---|---|---|
| true | +true | +true | +false | +true | +true | +||
| true | +true | +false | +false | +false | +false | +||
| true | +false | +true | +false | +false | +false | +||
| false | +true | +true | +true | +true | +true | +||
| true | +false | +false | +false | +false | +false | +||
| false | +true | +false | +true | +false | +true | +||
| false | +false | +true | +true | +false | +true | +||
| false | +false | +false | +true | +false | +true | +
The following statements are equivalent
+!(a && b) is equivalent to !a || !b
Notice how when you distribute the ! you have to flip the operand as well. && becomes ||
Same is true for the following example
+!(a || b) is equivalent to !a && !b
!(a || b && c) is equivalent to !a && (!b || !c)
In an && (AND) statement, if the left side is false, there is no need to evaluate the right side. Since it's going to be false anyways!!
false && true; // FALSE no matter what the right side isIn an || (OR) statement, if the left side is `true, there is no need to evaluate the right side. Since it's going to be true by default!!
true || false; // TRUE no matter what the right side isJava takes this shortcut by default for efficiency reasons
+Another way to perform multiway branching. Comparing a variable and constant values (String, int, char)
Switch statements cannot be used with boolean, double, or floats
switch (variable) {
+ case value1:
+ // Do something
+ break;
+ case value2:
+ // Do something else
+ break;
+ //...
+ default:
+ // If all else fails do this
+ break;
+}case is a reserved word that means "when our variable in consideration is equal to..."
If you forget the break keyword, then the program will keep doing the work of all the statements until it hits a break keyword.
switch (birthday) {
+ case 1:
+ birthstone = "garnet";
+ break;
+ case 2:
+ birthstone = "amethyst";
+ break;
+ // ....
+ default:
+ System.out.println("Not valid");
+ break;
+}Comparing strings are based on the ASCII value of characters
+Sorting strings will result in strings being in alphabetical or reverse alphabetical order. The values of the strings are compared character by character from the left with each ASCII value.
+To compare strings use the compareTo() method. Here is the format of the call
str1.compareTo(str2)This returns a negative number when str1 is less than str2
This returns 0 when str1 is equal to str1
This returns a positive number when str1 is greater than str2
String a = "apple";
+String b = "banana";
+
+int x = a.compareTo(b); // x = -1
+
+int y = b.compareTo(a); // y = 1With a ternary operator, you can shorten statements where a value is determined by an if statement
+String output = "";
+if (movieRating > 4) {
+ output = "Fan favorite";
+} else {
+ output = "Alright";
+}Is equivalent to
+String output = "";
+output = (movieRating > 4)? "Fan favorite": "Alright";double shipping;
+if (isPrimeMember) {
+ shipping = 0;
+} else {
+ shipping = 3.99;
+}Is equivalent to
+double shipping = (isPrimeMember)? 0: 3.99;You can use multi-line comments or single line comments to note about a piece of code. This helps you so that you don't forget what was happening in the code
+/* Multi Line Comment
+I am multiple lines woo!
+*/
+
+System.out.println("Hello"); // I am an inline commentThis is a standardized method in Java to describe your functions and classes
+/**
+ * This is a Javadoc comment. A description of your class should appear here
+ * @author Brandon
+ * @version 2018
+ */
+public class Example{
+ /** Another Javadoc comment.
+ * A description of your program should go here
+ */
+ public static void main(String[] args) {
+
+ }
+}Convention in Java is for all of your functions/method names to be lowercase.
+Class names should be title cased, each word is capitalized ex: IntroExample
The Java API is publicly accessible and is extremely useful for finding methods in existing classes.
+The API documentation for the String class is located here: https://docs.oracle.com/javase/8/docs/api/index.html?java/lang/String.html
Think about variables as buckets that hold information. Once the bucket is created, only one type of item can go in the bucket.
+sand bucket1;We can say that bucket1 is of type sand, that means the only thing that can go in the bucket is sand.
int bucket1;
+double bucket2;From the two lines above, we have declared the variable.
+Variables store state, they are a name for a location in memory.
+Always remember to initialize your variables. Otherwise there's nothing in the bucket!
+bucket1 = 5;You can combine both the declaration and initialization
+int count = 15;Remember when dealing with variables to stay true with the type, don't mix a bucket of water with a bucket of sand.
+We can update count to contain a true value
count = 55;count no longer has the value of 15 in it. There's no record of it! It has been overwritten with the value 55
There are 8 primitive types in Java
+byte through double are all numeric types
+boolean can only be equal to true or false
boolean student = true;Stores a single character from the Unicode set
+There are 16 bits per character which adds up to 65,536 characters
+ASCII is the US subset of the characters. You can look this up online when needing to deal with ASCII values
+char firstLetter = 'A';The different numeric types determine the precision of your number. Since numbers are not represented the same in the computer as they are in real life, there are some approximations.
+The default type you can use your code is int for integers and double for numbers with a decimal point
There are certain types of operations you can perform on numeric type
+| Symbol | +Meaning | +Example | +Value | +
|---|---|---|---|
| + | +addition | +43 + 8 | +51 | +
| - | +subtraction | +43.0-8.0 | +35.0 | +
| * | +multiplication | +43 * 8 | +344 | +
| / | +division | +43.0 / 8.0 | +5.375 | +
| % | +remainder / mod | +43 % 8 | +3 | +
| - | +unary minus | +-43 | +-43 | +
There are two types of in/decrementers postfix and prefix
+Postfix:
+int x = 0;
+int y = 7;
+x++; // Shortcut for x = x + 1
+y--; // Shortcut for y = y - 1Prefix
+int x = 0, y = 7, z;
+z = y * x++; // Equivalent to (y * x) + 1 = 7 * 0
+z = y * ++x; // Equivalent to y * (x + 1) = 7 * 1There are two types of data conversion, implicit and explicit
+The compiler can perform implicit data conversion automatically.
+Performing an explicit data conversion requires additional work on the programmer's part
+A conversion is implicit if you do not lose any information in it
+double price = 6.99;
+int sale = 3;
+double total = price - sale;A cast is an explicit data conversion. This is requested by a programmer, this can lead to loss of information
+int nextChar = 'b';
+Character.isAlphabetic( (char) nextChar); // Let's you print the actual letter 'b' instead of the number corresponding to it
+
+float price = 6.99;
+int cost = (int) price; // cost is now 6You can print the values of variables using System.out.println and System.out.print
The difference is that System.out.println adds a new line at the end. Meaning the next print out will be on the next line.
int cost = 5;
+double sale = .30;
+System.out.print(cost);
+System.out.print(sale);
+// Prints out '5.30`
+System.out.println(cost);
+System.out.println(sale);
+// Prints out '5'
+// Prints out '0.30'To add a space between two variables in a print, add " " to the expression in between the two variables
System.out.println("The total cost is " + 5 " dollars and" + " " + 93 + " cents");
+// The total cost is 5 dollars and 94 centsYou can get import from the user, we can do this using the Scanner class.
First import it at the top of your file
+import java.util.Scanner;All you can do with Scanner is outlined in the Java API at this link https://docs.oracle.com/javase/8/docs/api/index.html?java/util/Scanner.html
Create a Scanner object
+Scanner input = new Scanner(System.in);
+System.out.print("Please enter an integer: ");
+price = input.nextInt(); // The integer that the user inputs is now stored in price
+System.out.println("Your input: " + price);In Java, your code must live in a class.
+public class NameOfClass {
+ public static void main(String[] args) {
+ // All program code
+ }
+}It is important that NameOfClass is named meaningfully for the code. It is convention to use CamelCase when using classes. (Capitalize your class names!)
All methods have a method signature, it is unique to it. For main it is the words public static void and the argument String[] args.
public means that any other piece of code can reference it.
void means that the method returns nothing
main is the name of the method. It is important to have main since that tells the Java Interpreter where to start in your program.
String[] args is the command line arguments inputted into the program. For this part of the class, we don't need to worry about it.
If you noticed String is a class, it is not a primitive type. This is denoted in Java by having it capitalized.
There is an order of operations in programming as well. It goes like this:
+And from there you read from left to right.
+These are variables that can never be changed
+final int MINUTES_PER_HOUR = 60The keyword final indicates to the Java compiler that it is a constant variable.
By convention, constants are in all caps with underscores being separated between the words
+There are some arithmetic expressions that we want to be able to do and we cannot achieve that simply with the standard operations
+| Method | +Description | +
|---|---|
| Math.sqrt(x) | +square root | +
| Math.abs(x) | +absolute value | +
| Math.pow(a, b) | +exponentiation $a^b$ | +
| Math.max(a, b) | +returns the maximum of a or b | +
| Math.min(a, b) | +returns the minimum of a or b | +
| Math.round(x) | +rounds to the nearest integer | +
public class MoreVariables
+ public static void main(String[] args) {
+ // Decrate a variable
+ int x;
+
+ // Initialize ia variable
+ x = 5;
+
+ // Area of a square
+ int squareArea = x * x;
+ System.out.println("Area of a square: " + squareArea);
+ double newSquare = Math.pow(x, 2);
+ System.out.println("Area of square: " + newSquare);
+
+ // Area of Circle
+ final double PI = 3.14159;
+ double radius = 3;
+ double circleArea = radius * radius * PI;
+ System.out.println("Area of circle: " + circleArea);
+
+}These are concatenated chars
+'d' + 'o' + 'g' // equivalent to "dog""straw" + "berry" // strawberryStrings are denoted by double quotes "" rather than a string which is denoted by single quotes ''
String is not a primitive type, it is a class. Hence, why it is capitalized in Java.
+The java.lang.String is automatically imported in Java.
To declare and initialize a String
+String name = "Henry";In memory it appears as
+| H | +'e' | +'n' | +'r' | +'y' | +
|---|---|---|---|---|
| + | + | + | + | + |
int length()boolean equals(String another)boolean startsWith(String prefix)boolean endsWith(String suffix)String substring(int start, int end)int indexOf(int ch)String toLowerCase()char first = name.charAt(0);Remember in Java, that it starts counting from zero! If you try to access a letter that doesn't exist, it will produce an IndexOutOfBounds error.
There are two types of errors, compile-type errors and run-time errors. Later we will talk about debugging skills such as making "breakpoints" in your code so you can analyze the different variable values.
+Compile time errors are generated due to syntax errors. Forgot a semicolon? Missing a brace?
+These are logic errors. Not derived from syntax errors. An example of one that was discussed earlier is the IndexOutOfBounds error.
Let's talk about input right now. Let's say you have the following scenario
+Scanner input = new Scanner(System.in);
+System.out.println("Enter pet's age: ");
+int age = input.nextInt();
+System.out.println("Enter pet's name: ");
+String name = input.nextLine();
+System.out.println("Enter pet's breed: ");
+String breed = input.next();Then when we start to run the program...
+Enter pet's age:
+14
+Enter pet's name:
+Enter pet's breed:
+LabradoodleWhy did it skip pet's name? Let's run through the process again
+Enter pet's age:
+14 [ENTER]
+Enter pet's name:
+Enter pet's breed:
+LabradoodleHere the [ENTER] key gets saved into name.
+To resolve this, just use an input.nextLine() to throw away that [ENTER]
Scanner input = new Scanner(System.in);
+System.out.println("Enter pet's age: ");
+int age = input.nextInt();
+System.out.println("Enter pet's name: ");
+input.nextLine();
+String name = input.nextLine();
+System.out.println("Enter pet's breed: ");
+String breed = input.next();One of the ways you can do a random number generator is through this method:
+Import a class called random
+import java.util.Random;Then you need to create a Random object
Random rand = new Random();After this you can call the nextInt() method to get a random number between 0 and $2^{32}$
int randInt = rand.nextInt();If you don't want a random number between 0 and $2^{32}$ but instead to another maximum value, then you can call the nextInt method inserting the max integer as a parameter.
Random Integer from 0-10 (not including 10)
+int randInt2 = rand.nextInt(10);We have already encountered System.out.println and System.out.print but let us go over the differences again.
System.out.println() prints the contents inside the parenthesis and appends a newline character afterwards so that the next output is on a new line
System.out.print() prints the contents inside the parenthesis and does not output a newline character
If you want more control on how your output is displayed, it is recommended that you use System.out.printf to format your output
First, you need to specify your type using the % instruction
+Example:
+int sum = 50;
+System.out.printf("Total = %d", sum);This outputs
+Total = 50Notice here that there is no concatenation required like the previous two methods, instead you insert the variables as parameters
+Let us deconstruct the % instruction
+% .
+The first underline is the + - 0 space (sometimes we want to pad the money with zeros)
+The second underline is the width of the text
+The third underline is the number of decimal places
+The the final underline is the specifier f for decimal and d for integer
Example
+double amount = 0.5;
+System.out.printf("Total Due: %0.2f")This outputs
+Total Due: 0.50Methods are small blocks of statements that make it easier to solve a problem. It usually focuses on solving a small part of the overall problem.
+Usually in methods you provide some sort of input and get some output out of it.
+Consists of a method name, input and output and the block of statements.
+Usually this is succinctly written using JavaDocs which is what you see in the JavaAPI
+A method call is the execution of the method. The statements defined in the method is what will execute.
+Recall from method definition the parts of the method definition. Now look at the following method
+String[] split(String s)The output here is String[]
The method name is split
The input is String s
Let us look at the following example:
+The program should have a list of grocery prices. It should be able to calculate the total cost of groceries. The store gives a student discount of 5%. The program should calculate this discount and update the total, it should calculate and add the 2.5% tax.
+public static int timesTwo(int num) {
+ int two = num * 2;
+ return two;
+}It first starts off by declaring the visibility public
The return type if int
The method name is timesTwo
The input parameter is int num
Between the curly braces is the body of the method
+int a = 5;
+int b = 3;
+
+int ans = multiply(a, b)The method call is multiply(a, b) and the result is stored in the variable ans
With unit testing you are able to test small parts as you go. With unit testing, you can test many examples and border cases (extreme inputs) to see how your code reacts to different input.
+A variable declared inside a method is only accessible inside that method.
+A good rule of thumb is that if a variable is declared within curly braces {} then it does not exist outside that curly brace.
+Sometimes we need to have methods with the same name but different input parameters.
+public static int multiply(int a, int b) {
+ return a * b;
+}This method only works with integers, what if we want to multiply two doubles?
+We can overload the method by declaring another method with the same name.
+public static double multiply(double a, double b) {
+ return a * b;
+}Java has a special way that you can document your methods such that it will create documentation for you if you follow the convention.
+The Java API actually uses this technique to produce its own documentation.
+To create this, indicate a method with special comments that begin with /** and ends with */
It contains block tags that describe input and output parameters
+@param and @return
/**
+ * @param y an integer to sum
+ * @param x an integer to sum
+ * @return the sum of x and y
+ */
+public int multiply(int x, int y) {
+ return x + y;
+}We only want to create one user input scanner per program, we also only want one file input scanner per program.
+If a method needs a scanner, you can pass the one you already created in as an input parameter.
+Primitive types (int, char, double, etc.) are passed by value. Modifications made inside a method cannot be seen outside the method.
Arrays on the other hand, is pass by reference. Changes made to an array inside the method can be seen outside the method.
+public static void main(String[] args) {
+ int[] nums = {1, 3, 5, 7, 9};
+
+ timesTwo(nums);
+}
+public static void timesTwo(int[] arr) {
+ for (int i = 0; i < arr.length; i++) {
+ arr[i] *= 2;
+ }
+}At the end of the timesTwo method call, the variable nums would have {2, 6, 10, 14, 18}
When we declare an array, Java automatically fills every slot of the array with the type in memory. So if you know that you need exactly 8 slots, then you only ask for 8.
+This is when we don't know how many slots we need. Therefore, we ask for more than we think we'll need. That way we don't go out of bounds.
+If we do this, then we don't know how many elements we have already inserted into the array. Since the length is the number of slots.
+So we can create another variable, which will keep track of the index in where we can add the next element.
+We use oversized arrays when the size of the array is unknown or may change.
+Objects are known for having characteristics. A car has on average 4 wheels, 2-4 doors, a steering wheel.
+Objects can perform actions. A car can drive, hold cargo, and honk.
+Object-Oriented Programming
+public class Classname {
+ // Fields
+ // Constructors
+ // Methods
+}Fields are instance variables, they store values, help define state, and exist in memory for the lifetime of the object.
+public class Car {
+ private double price;
+ private double gas;
+}We can build an object through a constructor. It is a special kind of method, this method requires that you not have a return type and that you name it the same as the class itself.
+Constructors help set default field values for the different properties of our class.
+public class Car {
+ // Begin Constructor
+ public Car(double cost) {
+ this.price = cost;
+ this.gas = 0;
+ }
+ // End Constructor
+ private double price;
+ private double gas;
+}Note: The this keyword refers to the object's fields. This helps keep it separate from other variables you can create in the method and the input parameters you receive.
We like to classify methods into two types, accessors and mutators.
+Getter methods return a copy of an instance field. It does not change the state of the object.
+public double getPrice() {
+ return this.price;
+}This type of method modifies an instance field. It does not return anything and changes the state of the object.
+public void setPrice(double cost) {
+ this.price = cost;
+}public class Car {
+ // Instance Variables
+ private int mpg;
+ private double price;
+
+ // Constructors
+ public Car() {
+ this.price = 0;
+ this.mpg = 0;
+ }
+ public Car(double cost, int mpg) {
+ this.price = cost;
+ this.mpg = mpg;
+ }
+
+ // Accessors
+ public double getPrice() {
+ return this.price''
+ }
+ public int getMpg() {
+ return this.mpg;
+ }
+
+ // Mutators
+ public void setPrice(double cost) {
+ this.price = cost;
+ }
+ public void setMpg(int mpg) {
+ this.mpg = mpg;
+ }
+}Just like how we used the Scanner class, we can also use our new Car class.
public class TestCar {
+ public static void main(String[] args) {
+ // Declare an object reference
+ Car c;
+
+ // Initialize the object
+ c = new Car();
+
+ // Update the fields of the object
+ c.setPrice(3000);
+ c.setMpg(22);
+
+ // Print object information
+ System.out.println("Price is " + c.getPrice() )
+ }
+}These represent a fixed set of constants and include all possible values within them.
+Let's look at coins. On a daily basis in the US, we use the following coins:
+Other examples include the days of the week, clothes sizes, etc.
+Let's define an enum type
public enum Coin { PENNY, NICKEL, DIME, QUARTER}Now declare and initialize a variable
+Coin myCoin = Coin.PENNYArrays require you to say upfront, how many slots you need. ArrayLists are more flexible since you can change the length of the array during Runtime.
+Arrays can store objects and primitives such as int, char, boolean, etc.
ArrayLists can only store objects.
+ArrayList<objectType> list = new ArrayList<objectType>();Array
+int length = array.length;ArrayList
+int length = array.size();This is a special loop in where you tell it to go through all the elements of the array, without specifying an index.
+for (String b : buildings) {
+ System.out.print(b);
+}Let us suppose that we have the following array
+int[] b = {11, 12, 15, 16, 21}for (int i = 0; i < b.length; i++) {
+ b[i] += 2;
+}for (int i = 0; i < b.length; i++) {
+ System.out.println(b[i]);
+}int sum = 0;
+for (int i = 0; i < b.length; i++) {
+ sum += b[i];
+}System.out.println(b[b.length - 1]);System.out.println(b[b.length / 2]);For the summer session, I didn't write in depth lecture notes. Instead I wrote short complementary material on my old website to help out with the labs.
+ + + + + +I am the lab aide one of her sections on Computer Programming and Problem Solving (CPSC 220).
+My role in this class is to help students during lab time and to create lecture notes based on the lectures given in class.
+ +I've been the lab aide for two of Polack's classes so far: Introduction to Computer Science (CPSC 110) and Computer Programming and Problem Solving (CPSC 220).
+My role involves helping students debug during the lab time. I also create lecture notes for students to look over while working on projects.
+ +Fall 2017/CPSC 220 Lecture Notes
+I was the lab aide for Zacharski's DATA 101 class. My role in his class mostly involved helping students debug during lab time and grading student's demonstrations of what they've done in a lab.
+I was invited to give a talk to the Data Science Group at the University of Mary Washington called the Data Mavens. It was a 25 minute presentation on the bootstrap resampling technique with code on how to get started.
+ +At the VCU RamHacks hackathon, Clare Arrington, Harrison Crosse, and I demoed our product Similyrics. It's a web application that takes your favorite song, grabs the lyrics, and finds a song from a database that closely matches to the song you have chosen lyric wise.
+ +To identify clusters of observations we need to know how close individuals are to each other or how far apart they are.
+Two individuals are 'close' when their dissimilarity of distance is small and their similarity large.
+Special attention will be paid to proximity measures suitable for data consisting of repeated measures of the same variable, for example taken at different time points.
+Measures are generally scaled to be in the interval $[0, 1]$, although occasionally they are expressed as percentages in the range $0-100\%$
+Similarity value of unity indicates that both observations have identical values for all variables
+Similarity value of zero indicates that the two individuals differ maximally for all variables.
+An extensive list of similarity measures for binary data exist, the reason for such is that a large number of possible measures has to do with the apparent uncertainty as to how to deal with the count of zero-zero matches
+In some cases, zero-zero matches are equivalent to one-one matches and therefore should be included in the calculated similarity measure
+Example: Gender, where there is no preference as to which of the two categories should be coded as zero or one
+In other cases the inclusion or otherwise of the matches is more problematic
+Example: When the zero category corresponds to the genuine absence of some property, such as wings in a study of insects
+The question that then needs to be asked is do the co-absences contain useful information about the similarity of the two objects?
+Attributing a high degree of similarity to a pair of individuals simply because they both lack a large number of attributes may not be sensible in many situations
+The following table below will help when it comes to interpreting the measures
+
Measure that ignore the co-absence (lack of both objects having a zero) are Jaccard's Coefficient (S2), Sneath and Sokal (S4)
+When co-absences are considered informative, the simple matching coefficient (S1) is usually employed.
+Measures S3 and S5 are further examples of symmetric coefficients that treat positive matches (a) and negative matches (d) in the same way.
+
+Categorical data where the variables have more than two levels (for example, eye color) could be dealt with in a similar way to binary data, with each level of a variable being regarded as a single binary variable.
+This is not an attractive approach, however, simply because of the large number of ‘negative’ matches which will inevitably be involved.
+A superior method is to allocate a score of zero or one to each variable depending on whether the two observations are the same on that variable. These scores are then averaged over all p variables to give the required similarity coefficient as +$$ +s{ij} = \frac{1}{p}\sum{k = 1}^p{s_{ik}} +$$
+A metric on a set $X$ is a distance function +$$ +d : X \times X \to [0, \infty) +$$ +where $[0, \infty)$ is the set of non-negative real numbers and for all $x, y, z \in X$, the following conditions are satisfied
+Conditions 1 and 2 define a positive-definite function
+All distance measures are formulated so as to allow for differential weighting of the quantitative variables $w_k$ denotes the nonnegative weights of $p$ variables
+
+Proposed dissimilarity measures can be broadly divided into distance measures and correlation-type measures.
+The Minkowski distance is a metric in normed vector space which can be considered as a generalization of both the Euclidean distance and the Manhattan distance +$$ +D(X, Y) = (\sum_{i = 1}^n{w_i^p|x_i - y_i|^p})^{\frac{1}{p}} +$$ +This is a metric for $p > 1$
+This is the case in the Minkowski distance when $p = 1$ +$$ +d(X, Y) = \sum_{i = 1}^n{w_i|x_i - y_i|} +$$ +Manhattan distance depends on the rotation of the coordinate system, but does not depend on its reflection about a coordinate axis or its translation +$$ +d(x, y) = d(-x, -y) +$$
+$$ +d(x, y) = d(x + a, y + a) +$$
+Shortest paths are not unique in this metric
+This is the case in the Minkowski distance when $p = 2$. The Euclidean distance between points X and Y is the length of the line segment connection them. +$$ +d(X, Y) = \sqrt{\sum_{i = 1}^n{w_i^2(x_i - y_i)^2}} +$$ +There is a unique path in which it has the shortest distance. This distance metric is also translation and rotation invariant
+The standard Euclidean distance can be squared in order to place progressively greater weight on objects that are farther apart. In this case, the equation becomes +$$ +d(X, Y) = \sum_{i = 1}^n{w_i^2(x_i - y_i)^2} +$$ +Squared Euclidean Distance is not a metric as it does not satisfy the triangle inequality, however, it is frequently used in optimization problems in which distances only have to be compared.
+The Chebyshev distance is where the distance between two vectors is the greatest of their differences along any coordinate dimension.
+It is also known as chessboard distance, since in the game of chess the minimum number of moves needed by a king to go from one square on a chessboard to another equals the Chebyshev distance +$$ +d(X, Y) = \lim{p \to \infty}{(\sum{i = 1}^n{|x_i - y_i|^p})}^\frac{1}{p} +$$
+$$ += max_i(|x_i - y_i|) +$$
+Chebyshev distance is translation invariant
+The Canberra distance (D4) is a weighted version of the $L1$ Manhattan distance. This measure is very sensitive to small changes close to $x{ik} = x_{jk} = 0$.
+It is often regarded as a generalization of the dissimilarity measure for binary data. In this context the measure can be divided by the number of variables, $p$, to ensure a dissimilarity coefficient in the interval $[0, 1]$
+It can then be shown that this measure for binary variables is just one minus the matching coefficient.
+It has often been suggested that the correlation between two observations can be used to quantify the similarity between them.
+Since for correlation coefficients we have that $-1 \le \phi{ij} \le 1$ with the value ‘1’ reflecting the strongest possible positive relationship and the value ‘-1’ the strongest possible negative relationship, these coefficients can be transformed into dissimilarities, $d{ij}$, within the interval $[0, 1]$
+The use of correlation coefficients in this context is far more contentious than its noncontroversial role in assessing the linear relationship between two variables based on $n$ observations.
+When correlations between two individuals are used to quantify their similarity the rows of the data matrix are standardized, not its columns.
+Disadvantages
+When variables are measured on different scales the notion of a difference between variable values and consequently that of a mean variable value or a variance is meaningless.
+In addition, the correlation coefficient is unable to measure the difference in size between two observations.
+Advantages
+However, the use of a correlation coefficient can be justified for situations where all of the variables have been measured on the same scale and precise values taken are important only to the extent that they provide information about the subject's relative profile
+Example: In classifying animals or plants, the absolute size of the organisms or their parts are often less important than their shapes. In such studies the investigator requires a dissimilarity coefficient that takes the value zero if and only if two individuals' profiles are multiples of each other. The angular separation dissimilarity measure has this property.
+Further considerations
+The Pearson correlation is sensitive to outliers. This has prompted a number of suggestions for modifying correlation coefficients when used as similarity measures; for example, robust versions of correlation coefficients such as jackknife correlation or altogether more general association coefficients such as mutual information distance measure
+Mahalanobis distance is a measure of distance between a point P and a distribution D. It is a multi-dimensional generalization of the idea of measuring how many standard deviations away P is from the mean of D
+Mahalanobis distance is unitless and scale-invariant and takes into account the correlations of the data set +$$ +D(\vec{x}) = \sqrt{(\vec{x} - \vec{\mu})^T S^{-1}(\vec{x}-\vec{\mu})} +$$ +Where $\mu$ is a set of mean observations and $S$ is the covariance matrix
+If the covariance matrix is diagonal then the resulting distance measure is called a normalized Euclidean distance. +$$ +d(\vec{x}, \vec{y}) = \sqrt{\sum_{i = 1}^N{\frac{(x_i - y_i)^2}{s^2_i}}} +$$ +Where $s_i$ is the standard deviation of the $x_i$ and $y_i$ over the sample set
+This metric describes whether or not two observations are equivalent +$$ +\rho(x, y) = \begin{cases} +1 & x \not= y \ +0 & x = y +\end{cases} +$$
+There are a number of approaches to constructing proximities for mixed-mode data, that is, data in which some variables are continuous and some categorical.
+Most general-purpose statistical software implement a number of measurs for converting two-mode data matrix into a one-mode dissimilarity matrix.
+R has cluster, clusterSim, or proxy
We'll be looking at data that consists of repeated measures of the same outcome variable but under different conditions.
+The simplest and perhaps most commonly used approach to exploiting the reference variable is in the construction of a reduced set of relevant summaries per object which are then used as the basis for defining object similarity.
+Here we will look at some approaches for choosing summary measures and resulting proximity measures for the most frequently encountered reference vectors (e.g. time, experimental condition, and underlying factor)
+Structured data arise when the variables can be assumed to follow a known factor model. Under confirmatory factor analysis model each variable or item can be allocated to one of a set of underlying factors or concepts. The factors cannot be observed directly but are 'indicated' by a number of items that are all measured on the same scale.
+Note that the summary approach, while typically used with continuous variables, is not limited to variables on an interval scale. The same principles can be applied to dealing with categorical data. The difference is that summary measures now need to capture relevant aspects of the distribution of categorical variables over repeated measures.
+Rows of $X$ which represent ordered lists of elements, that is all the variables provide a categorical outcome and these variables can be aligned in one dimension, are more generally referred to as sequences. Sequence analysis is an area of research that centers on problems of events and actions in their temporal context and includes the measurements of similarities between sequences.
+The most popular measure of dissimilarity between two sequences is the Levenhstein distance and counts the minimum number of operations needed to transform one sequence of categories into another, where an operation is an insertion, a deletion, or a substitution of a single category. Each operation may be assigned a penalty weight (a typical choice would be to give double the penalty to a substitution as opposed to an insertion or deletion. The measure is sometimes called the 'edit distance' due to its application in spell checkers.
+Optimal matching algorithms (OMAs) need to be employed to find the minimum number of operations required to match one sequence to another. One such algorithm for aligning sequences is the Needleman-Wunsch algorithm, which is commonly used in bioinformatics to align proteins.
+The Jary similarity measure is a related measure of similarity between sequences of categories often used to delete duplicates in the area of record linkage.
+In clustering applications, it becomes necessary to consider how to measure the proximity between groups of individuals.
+For deriving inter-group proximities from a matrix of inter-individual proximities, there are a variety of possibilities
+One method for constructing inter-group dissimilarity measures for continuous data is to simply substitute group means (also known as the centroid) for the variable values in the formulae for inter-individual measures
+More appropriate, however, might be measures which incorporate in one way or another, knowledge of within-group variation. One possibility is to use Mahallanobis's generalized distance.
+Mahalanobis distance is a measure of distance between a point P and a distribution D. It is a multi-dimensional generalization of the idea of measuring how many standard deviations away P is from the mean of D
+Mahalanobis distance is unitless and scale-invariant and takes into account the correlations of the data set +$$ +D(\vec{x}) = \sqrt{(\vec{x} - \vec{\mu})^T S^{-1}(\vec{x}-\vec{\mu})} +$$ +Where $\mu$ is a set of mean observations and $S$ is the covariance matrix
+If the covariance matrix is diagonal then the resulting distance measure is called a normalized Euclidean distance. +$$ +d(\vec{x}, \vec{y}) = \sqrt{\sum_{i = 1}^N{\frac{(x_i - y_i)^2}{s^2_i}}} +$$ +Where $s_i$ is the standard deviation of the $x_i$ and $y_i$ over the sample set
+Thus, the Mahalanobis distance incraeses with increasing distances between the two group centers and with decreasing within-group variation.
+By also employing within-group correlations, the Mahalanobis distance takes account the possibly non-spherical shapes of the groups.
+The use of Mahalanobis implies that the investigator is willing to assume that the covariance matrices are at least approximately the same in the two groups. When this is not so, this measure is an inappropriate inter-group measure. Other alternatives exist such as the one proposed by Anderson and Bahadur
+Another alternative is the normal information radius suggested by Jardine and Sibson
+Approaches for measuring inter-group dissimilarities between groups of individuals for which categorical variables have been observed have been considered by a number of authors. Balakrishnan and Sanghvi (1968), for example, proposed a dissimilarity index of the form
+where $p{Akl}$ and $p{Bkl}$ are the proportions of the lth category of the kth variable in group A and B respectively, , ck + 1 is the number of categories for the kth variable and p is the number of variables.
Kurczynski (1969) suggested adapting the generalized Mahalanobis distance, with categorical variables replacing quantitative variables. In its most general form, this measure for inter-group distance is given by
+where contains sample proportions in group A and
is defined in a similar manner, and
is the m × m common sample covariance matrix, where
.
To weight a variable means to give it greater or lesser importance than other variables in determining the proximity between two objects.
+The question is 'How should the weights be chosen?' Before we discuss this question, it is important to realize that the selection of variables for inclusion into the study already presents a form of weighting, since the variables not included are effectively being given the weight $0$.
+The weights chosen for the variables reflect the importance that the investigator assigns the variables for the classification task.
+There are several approaches to this
+The second approach assumed the importance of a variable to be inversely proportional to the total variability of that variable. The total variability of a variable comprises variation both within and between groups which may exist within the set of individuals. The aim of cluster analysis is typically to identify such groups. Hence it can be argued that the importance of a variable should not be reduced because of between-group variation (on the contrary, one might wish to assign more importance to a variable that shows larger between-group variation.)
+Gnanadesikan et al. (1995) assessed the ability of squared distance functions based on data-determined weights, both those described above and others, to recover groups in eight simulated and real continuous data sets in a subsequent cluster analysis. Their main findings were:
+In many clustering applications, the variables describing the objects to be clustered will not be measured in the same units. A number of variability measures have been used for this purpose
+The second is shown to outperform auto-scaling in many clustering applications. As pointed out in the previous section, standardization of variables to unit variance can be viewed as a special case of weighting.
+Firstly, the nature of the data should strongly influence the choice of the proximity measure.
+Next, the choice of measure should depend on the scale of the data. Similarity coefficients should be used when the data is binary. For continuous data, distance of correlation-type dissimilarity measure should be used according to whether 'size' or 'shape' of the objects is of interest.
+Finally, the clustering method to be used might have some implications for the choice of the coefficient. For example, making a choice between several proximity coefficients with similar properties which are also known to be monotonically related can be avoided by employing a cluster method that depends only on the ranking of the proximities, not their absolute values.
+This technique validates the consistency within clusters of data. It provides a succinct graphical representation of how well each object lies in its cluster.
+The silhouette ranges from -1 to 1 where a high value indicates that the object is consistent within its own cluster and poorly matched to neighboring clustesr.
+A low or negative silhouette value can mean that the current clustering configuration has too many or too few clusters.
+For each datum $i$, let $a(i)$ be the average distance of $i$ with all other data within the same cluster.
+$a(i)$ can be interpreted as how well $i$ is assigned to its cluster. (lower values mean better agreement)
+We can then define the average dissimilarity of point $i$ to a cluster $c$ as the average distance from $i$ to all points in $c$.
+Let $b(i)$ be the lowest average distance of $i$ to all other points in any other cluster in which i is not already a member.
+The cluster with this lowest average dissimilarity is said to be the neighboring cluster of $i$. From here we can define a silhouette: +$$ +s(i) = \frac{b(i) - a(i)}{max{a(i), b(i)}} +$$ +The average $s(i)$ over all data of a cluster is a measure of how tightly grouped all the data in the cluster are. A silhouette plot may be used to visualize the agreement between each of the data and its cluster.
+
Recall that $a(i)$ is a measure of how dissimilar $i$ is to its own cluster, a smaller value means that it's in agreement to its cluster. For $s(i)$ to be close to 1, we require $a(i) << b(i)$ .
+If $s(i)$ is close to negative one, then by the same logic we can see that $i$ would be more appropriate if it was clustered in its neighboring cluster.
+$s(i)$ near zero means that the datum is on the border of two natural clusters.
+This can also be used in helping to determine the number of clusters in a dataset. The ideal number of cluster is one that produces the highest silhouette value.
+Also a good indication that one has too many clusters is if there are clusters with the majority of observations being under the mean silhouette value.
+ +
In centroid-based clustering, clusters are represented by some central vector which may or may not be a member of the dataset. In practice, the number of clusters is fixed to $k$ and the goal is to solve some sort of optimization problem.
+The similarity of two clusters is defined as the similarity of their centroids.
+This problem is computationally difficult so there are efficient heuristic algorithms that are commonly employed. These usually converge quickly to a local optimum.
+This aims to partition $n$ observations into $k$ clusters in which each observation belongs to the cluster with the nearest mean which serves as the centroid of the cluster.
+This technique results in partitioning the data space into Voronoi cells.
+Given a set of observations $x$, k-means clustering aims to partition the $n$ observations into $k$ sets $S$ so as to minimize the within-cluster sum of squares (i.e. variance). More formally, the objective is to find +$$ +argmins{\sum{i = 1}^k{\sum_{x \in S_i}{||x-\mui||^2}}}= argmin{s}{\sum_{i = 1}^k{|S_i|Var(S_i)}} +$$ +where $\mu_i$ is the mean of points in $S_i$. This is equivalent to minimizing the pairwise squared deviations of points in the same cluster +$$ +argmins{\sum{i = 1}^k{\frac{1}{2|Si|}\sum{x, y \in S_i}{||x-y||^2}}} +$$
+Given an initial set of $k$ means, the algorithm proceeds by alternating between two steps.
+Assignment step: Assign each observation to the cluster whose mean has the least squared euclidean distance.
+Update Step: Calculate the new means to be the centroids of the observations in the new clusters
+The algorithm is known to have converged when assignments no longer change. There is no guarantee that the optimum is found using this algorithm.
+The result depends on the initial clusters. It is common to run this multiple times with different starting conditions.
+Using a different distance function other than the squared Euclidean distance may stop the algorithm from converging.
+Commonly used initialization methods are Forgy and Random Partition.
+Forgy Method: This method randomly chooses $k$ observations from the data set and uses these are the initial means
+This method is known to spread the initial means out
+Random Partition Method: This method first randomly assigns a cluster to each observation and then proceeds to the update step.
+This method is known to place most of the means close to the center of the dataset.
+A Voronoi diagram is a partitioning of a plan into regions based on distance to points in a specific subset of the plane.
+The set of points (often called seeds, sites, or generators) is specified beforehand, and for each seed there is a corresponding region consisting of all points closer to that seed than any other.
+Different metrics may be used and often result in different Voronoi diagrams
+Euclidean
+
Manhattan
+
K-means++ is an algorithm for choosing the initial values or seeds for the k-means clustering algorithm. This was proposed as a way of avoiding the sometimes poor clustering found by a standard k-means algorithm.
+The intuition behind this approach involves spreading out the $k$ initial cluster centers. The first cluster center is chosen uniformly at random from the data points that are being clustered, after which each subsequent cluster center is chosen from the remaining data points with probability proportional to its squared distance from the point's closest existing cluster center.
+The exact algorithm is as follows
+A medoid can be defined as the object of a cluster whose average dissimilarity to all the objects in the cluster is minimal.
+The K-medoids algorithm is related to k-means and the medoidshift algorithm. Both the k-means and k-medoids algorithms are partition and both attempt to minimize the distance between points in the cluster to it's center. In contrast to k-means, it chooses data points as centers and uses the Manhattan Norm to define the distance between data points instead of the Euclidean.
+This method is known to be more robust to noise and outliers compared to k-means since it minimizes the sum of pairwise dissimilarities instead of the sum of squared Euclidean distances.
+There are several algorithms that have been created as an optimization to an exhaustive search. In this section, we'll discuss PAM and Voronoi iteration method.
+This is a variant of the PAM algorithm that relies on the sampling approach to handle large datasets. The cost of a particular cluster configuration is the mean cost of all the dissimilarities.
+Both PAM and CLARA are defined in the cluster package in R.
clara(x, k, metric = "euclidean", stand = FALSE, samples = 5,
+ sampsize = min(n, 40 + 2 * k), trace = 0, medoids.x = TRUE,
+ keep.data = medoids.x, rngR = FALSE)pam(x, k, metric = "euclidean", stand = FALSE)This is a variation of k-means clustering where instead of calculating the mean for each cluster to determine its centroid we are going to calculate the median instead.
+This has the effect of minimizing error over all the clusters with respect to the Manhattan norm as opposed to the Euclidean squared norm which is minimized in K-means
+Given an initial set of $k$ medians, the algorithm proceeds by alternating between two steps.
+Assignment step: Assign each observation to the cluster whose median has the leas Manhattan distance.
+Update Step: Calculate the new medians to be the centroids of the observations in the new clusters
+The algorithm is known to have converged when assignments no longer change. There is no guarantee that the optimum is found using this algorithm.
+The result depends on the initial clusters. It is common to run this multiple times with different starting conditions.
+In density-based clustering, clusters are defined as areas of higher density than the remainder of the data sets. Objects in more sparse areas are considered to be outliers or border points. This helps discover clusters of arbitrary shape.
+Given a set of points in space, it groups together points that are closely packed together while marking points that lie alone in low-density regions as outliers.
+Non core points can be part of a cluster, but they form its "edge", since they cannot be used to reach more points.
+Reachability is not a symmetric relation since, by definition, no point may be reachable from a non-core point, regardless of distance.
+Two points $p$ and $q$ are density-connected if there is a point $o$ such that both $p$ and $q$ are reachable from $o$. Density-connectedness is symmetric.
+A cluster then satisfies two properties:
+$k$: $k$ must be larger than $(D + 1)$ where $D$ is the number of dimensions in the dataset. Normally $k$ is chosen to be twice the number of dimensions.
+$\epsilon$: Ideally the $k^{th}$ nearest neighbors are at roughly the same distance. Plot the sorted distance of every point to it's $k^{th}$ nearest neighbor
+Example of Run Through
+https://www.cse.buffalo.edu/~jing/cse601/fa12/materials/clustering_density.pdf
+I made an attempt to find out in what situations people use different distance measures. Looking around in the Internet usually produces the results "It depends on the problem" or "I typically just always use Euclidean"
+Which as you might imagine, isn't a terribly useful answer. Since it doesn't give me any examples of which types of problems different distances solve.
+Therefore, let's think about it in a different way. What properties do different distance measures have that make them desirable?
+Variables can be on different scales. The Mahalanobis formula has a built in variance-covariance matrix which allows you to rescale your variables to make distances of different variables more comparable.
+https://stats.stackexchange.com/questions/50949/why-use-the-mahalanobis-distance#50956
+In higher dimensions, the points essentially become uniformly distant from one another. This is a problem observed in most distance metrics but it's more obvious with the Euclidean one.
+ +Hopefully in this course, we'll discover more properties as to why it makes sense to use different distance measures since it can have a impact on how our clusters are formed.
+Principal component analysis is a statistical procedure that performs an orthogonal transformation to convert a set of variables into a set of linearly uncorrelated variables called principle components.
+Number of distinct principle components equals $min(# Variables, # Observations - 1)$
+The transformation is defined in such a way that the first principle component has the largest possible variance explained in the data.
+Each succeeding component has the highest possible variance under the constraint of having to be orthogonal to the preceding components.
+PCA is sensitive to the relative scaling of the original variables.
+Results are discussed in terms of component scores which is the transformed variables and loadings which is the weight by which each original variable should be multiplied to get the component score.
+pcal = function(data) {
+ centered_data = scale(data)
+ covariance = cov(centered_data)
+ eigen_stuff = eigen(covariance)
+ sorted_indices = sort(eigen_stuff$values,
+ index.return = T,
+ decreasing = T)$ix
+ loadings = eigen_stuff$values[sorted_indices]
+ components = eigen_stuff$vectors[sorted_indices,]
+ combined_list = list(loadings, components)
+ names(combined_list) = c("Loadings", "Components")
+ return(combined_list)
+}An additional use case for Manhatten distance is when dealing with binary vectors. This approach, otherwise known as the Hamming distance, is the number of bits that are different between two binary vectors.
+Ordinal variables can be treated as if they were on a interval scale.
+First, replace the ordinal variable value by its rank ($r_{if}$) Then map the range of each variable onto the interval $[0, 1]$ by replacing the $fi$ where f is the variable and i is the object by +$$ +z{if} = \frac{r_{if} - 1}{M_f - 1} +$$ +Where $M_f$ is the maximum rank.
+Freshman = $0$ Sophmore = $\frac{1}{3}$ Junior = $\frac{2}{3}$ Senior = $1$
+$d(freshman, senior) = 1$
+$d(junior, senior) = \frac{1}{3}$
+This is the assessment of the suitability of clustering. Cluster Tendency determines whether the data has any inherent grouping structure.
+This is a hard task since there are so many different definitions of clusters (portioning, hierarchical, density, graph, etc.) Even after fixing a cluster type, this is still hard in defining an appropriate null model for a data set.
+One way we can go about measuring cluster tendency is to compare the data against random data. On average, random data should not contain clusters.
+There are some clusterability assessment methods such as Spatial histogram, distance distribution and Hopkins statistic.
+Let $X$ be the set of $n$ data points in $d$ dimensional space. Consider a random sample (without replacement) of $m << n$ data points. Also generate a set $Y$ of $m$ uniformly randomly distributed data points.
+Now define two distance measures $u_i$ to be the distance of $y_i \in Y$ from its nearest neighbor in X and $w_i$ to be the distance of $x_i \in X$ from its nearest neighbor in X
+We can then define Hopkins statistic as +$$ +H = \frac{\sum_{i = 1}^m{ui^d}}{\sum{i = 1}^m{ui^d} + \sum{i =1}^m{w_i^d}} +$$
+With this definition, uniform random data should tend to have values near 0.5, and clustered data should tend to have values nearer to 1.
+However, data containing a single Gaussian will also score close to one. As this statistic measures deviation from a uniform distribution. Making this statistic less useful in application as real data is usually not remotely uniform.
+For this method, I'm not too sure how this works, but here are some key points I found.
+Divide each dimension in equal width bins, and count how many points lie in each of the bins and obtain the empirical joint probability mass function.
+Do the same for the randomly sampled data
+Finally compute how much they differ using the Kullback-Leibler (KL) divergence value. If it differs greatly than we can say that the data is clusterable.
+There are a number of ways to maximize the variance of a principal component. To create an unique solution we should impose a constraint. Let us say that the sum of the square of the coefficients must equal 1. In vector notation this is the same as +$$ +a_i^Ta_i = 1 +$$ +Every future principal component is said to be orthogonal to all the principal components previous to it. +$$ +a_j^Tai = 0, i < j +$$ +The total variance of the $q$ principal components will equal the total variance of the original variables +$$ +\sum{i = 1}^q {\lambda_i} = trace(S) +$$ +Where $S$ is the sample covariance matrix.
+The proportion of accounted variation in each principle component is +$$ +P_j = \frac{\lambdaj}{trace(S)} +$$ +From this, we can generalize to the first $m$ principal components where $m < q$ and find the proportion $P^{(m)}$ of variation accounted for +$$ +P^{(m)} = \frac{\sum{i = 1}^m{\lambda_i}}{trace(S)} +$$ +You can think of the first principal component as the line of best fit that minimizes the residuals orthogonal to it.
+As a reminder to the last lecture, PCA is not scale-invariant. Therefore, transformations done to the dataset before PCA and after PCA often lead to different results and possibly conclusions.
+Additionally, if there are large differences between the variances of the original variables, then those whose variances are largest will tend to dominate the early components.
+Therefore, principal components should only be extracted from the sample covariance matrix when all of the original variables have roughly the same scale.
+But it is rare in practice to have a scenario when all of the variables are of the same scale. Therefore, principal components are typically extracted from the correlation matrix $R$
+Choosing to work with the correlation matrix rather than the covariance matrix treats the variables as all equally important when performing PCA.
+Let $R$ be the correlation matrix +$$ +R = \begin{pmatrix} +1 & r \ +r & 1 +\end{pmatrix} +$$ +Let us find the eigenvectors and eigenvalues of the correlation matrix +$$ +det(R - \lambda I) = 0 +$$
+$$ +(1-\lambda)^2 - r^2 = 0 +$$
+$$ +\lambda_1 = 1 + r, \lambda_2 = 1 - r +$$
+Let us remember to check the condition "sum of the principal components equals the trace of the correlation matrix": +$$ +\lambda_1 + \lambda_2 = 1+r + (1 - r) = 2 = trace(R) +$$
+Looking back at the characteristic equation +$$ +Ra_1 = \lambda a1 +$$ +We can get the following two formulas +$$ +a{11} + ra{12} = (1+r)a{11} \tag{1} +$$
+$$ +ra{11} + a{12} = (1 + r)a_{12} \tag{2} +$$
+Now let us find out what $a{11}$ and $a{12}$ equal. First let us solve for $a{11}$ using equation $(1)$ +$$ +ra{12} = (1+r)a{11} - a{11} +$$
+$$ +ra{12} = a{11}(1 + r - 1) +$$
+$$ +ra{12} = ra{11} +$$
+$$ +a{12} = a{11} +$$
+Where $r$ does not equal $0$.
+Now we must apply the condition of sum squares +$$ +a_1^Ta_1 = 1 +$$
+$$ +a{11}^2 + a{12}^2 = 1 +$$
+Recall that $a{12} = a{11}$ +$$ +2a_{11}^2 = 1 +$$
+$$ +a_{11}^2 = \frac{1}{2} +$$
+$$ +a_{11} =\pm \frac{1}{\sqrt{2}} +$$
+For sake of choosing a value, let us take the principal root and say $a_{11} = \frac{1}{\sqrt{2}}$
+Recall the fact that each subsequent eigenvector is orthogonal to the first. This means +$$ +a{11}a{21} + a{12}a{22} = 0 +$$ +Substituting the values for $a{11}$ and $a{12}$ calculated in the previous section +$$ +\frac{1}{\sqrt{2}}a{21} + \frac{1}{\sqrt{2}}a{22} = 0 +$$
+$$ +a{21} + a{22} = 0 +$$
+$$ +a{21} = -a{22} +$$
+Since this eigenvector also needs to satisfy the first condition, we get the following values +$$ +a{21} = \frac{1}{\sqrt{2}} , a{22} = \frac{-1}{\sqrt{2}} +$$
+From this, we can say that the first principal components are given by +$$ +y_1 = \frac{1}{\sqrt{2}}(x_1 + x_2), y_2 = \frac{1}{\sqrt{2}}(x_1-x_2) +$$ +With the variance of the first principal component being given by $(1+r)$ and the second by $(1-r)$
+Due to this, as $r$ increases, so does the variance explained in the first principal component. This in turn, lowers the variance explained in the second principal component.
+Principal Component Analysis is typically used in dimensionality reduction efforts. Therefore, there are several strategies for picking the right number of principal components to keep. Here are a few:
+Hierarchical algorithms combine observations to form clusters based on their distance.
+Hierarchal Clustering techniques can be subdivided depending on the method of going about it.
+First there are two different methods in forming the clusters Agglomerative and Divisive
+Agglomerative is when you combine the n individuals into groups through each iteration
+Divisive is when you are separating one giant group into finer groupings with each iteration.
+Hierarchical methods are an irrevocable algorithm, once it joins or separates a grouping, it cannot be undone. As Kaufman and Rousseeuw (1990) colorfully comment: "A hierarchical method suffers from the defect that it can never repair what was done in previous steps".
+It is the job of the statistician to decide when to stop the agglomerative or decisive algorithm, since having one giant cluster containing all observations or having each observation be a cluster isn't particularly useful.
+At different distances, different clusters are formed and are more readily represented using a dendrogram. These algorithms do not provide a unique solution but rather provide an extensive hierarchy of clusters that merge or divide at different distances.
+Apart from the method of forming clusters, the user also needs to decide on a linkage criterion to use. Meaning, how do you want to optimize your clusters.
+Do you want to group based on the nearest points in each cluster? Nearest Neighbor Clustering
+Or do you want to based on the farthest observations in each cluster? Farthest neighbor clustering.
+
This method is not very robust towards outliers, which will either show up as additional clusters or even cause other clusters to merge depending on the clustering method.
+As we go through this section, we will go into detail about the different linkage criterion and other parameters of this model.
+First let us consider the single linkage (nearest neighbor) approach. The clusters can be found through the following algorithm
+Single linkage can operate directly on a proximity matrix, the actual data is not required.
+A wonderful visual representation can be found in Everitt Section 4.2
+This is another criterion measure that requires both the data and proximity matrix. These are the following steps of the algorithm. Requires Euclidean distance measure to preserve geometric correctness
+This is like Single Linkage, except now we're taking the farthest distance. The algorithm can be adjusted to the following
+In this criterion, we are no longer summarizing each cluster before taking distances, but instead comparing each observation in the cluster to the outside point and taking the average
+This approach is similar to the UPGMA approach, except now we're taking the median instead of the mean
+This one I didn't look too far into but here's the description: With Ward's linkage method, the distance between two clusters is the sum of squared deviations from points to centroids. The objective of Ward's linkage is to minimize the within-cluster sum of squares.
+According to the following two stack overflow posts: https://stats.stackexchange.com/questions/195446/choosing-the-right-linkage-method-for-hierarchical-clustering and https://stats.stackexchange.com/questions/195456/how-to-select-a-clustering-method-how-to-validate-a-cluster-solution-to-warran/195481#195481
+These are the following ways you can justify a linkage type.
+Cluster metaphor. "I preferred this method because it constitutes clusters such (or such a way) which meets with my concept of a cluster in my particular project."
+Data/method assumptions. "I preferred this method because my data nature or format predispose to it."
+Internal validity. "I preferred this method because it gave me most clear-cut, tight-and-isolated clusters."
+External validity. "I preferred this method because it gave me clusters which differ by their background or clusters which match with the true ones I know."
+Cross-validity. "I preferred this method because it is giving me very similar clusters on equivalent samples of the data or extrapolates well onto such samples."
+Interpretation. "I preferred this method because it gave me clusters which, explained, are most persuasive that there is meaning in the world."
+Let us explore the idea of cluster metaphors now.
+Single Linkage or Nearest Neighbor is a spectrum or chain.
+Since single linkage joins clusters by the shortest link between them, the technique cannot discern poorly separated clusters. On the other hand, single linkage is one of the few clustering methods that can delineate nonelipsodial clusters.
+Complete Linkage or Farthest Neighbor is a circle.
+Between-Group Average linkage (UPGMA) is a united *class
+Centroid method (UPGMC) is proximity of platforms (commonly used in politics)
+A dendrogram is a tree diagram frequently used to illustrate the arrangement of the clusters produced by hierarchical clustering. It shows how different clusters are formed at different distance groupings.
+Divisive methods work in the opposite direction of agglomerative methods. They take one large cluster and successively splits it.
+This is computationally demanding if all $2^{k - 1} - 1$ possible divisions into two sub-clusters of a cluster of $k$ objects are considered at each stage.
+While less common than Agglomerative methods, divisive methods have the advantage that most users are interested in the main structure in their data, and this is revealed from the outset of a divisive method.
+For data consisting of $p$ binary variables, there is a computationally efficient method known as monothetic divisive methods available.
+Monothetic divisions divide clusters according to the presence or absence of each of the $p$ variables, so that at each stage, clusters contain members with certain attributes that are either all present or all absent.
+The term 'monothetic' refers to the use of a single variable on which to base the split on. Polythetic methods use all the variables at each stage.
+The choice of the variable on which a split is made depends on optimizing a criterion reflecting either cluster homogeneity or association with other variables.
+This tends to minimize the number of splits that have to be made.
+An example of an homogeneity criterion is the information content $C$
+This is defined with $p$ variables and $n$ objections where $fk$ is the number of individuals with the $k$ attribute +$$ +C = pn\log{n}-\sum{k = 1}^p{(f_k\log{f_k} - (n-f_k)\log{(n-f_k)})} +$$
+Recall that another way to split is based on the association with other variables. The attribute used at each step can be chosen according to its overall association with all attributes remaining at each step.
+This is sometimes termed association analysis.
+| + | V1 | +V2 | +
|---|---|---|
| + | 1 | +0 | +
| 1 | +a | +b | +
| 0 | +c | +d | +
$$ +|ad-bc| \tag{4.6} +$$
+$$ +(ad-bc)^2 \tag{4.7} +$$
+$$ +\frac{(ad-bc)^2n}{(a+b)(a+c)(b+d)(c+d)} \tag{4.8} +$$
+$$ +\sqrt{\frac{(ad-bc)^2n}{(a+b)(a+c)(b+d)(c+d)}} \tag{4.9} +$$
+$$ +\frac{(ad-bc)^2}{(a+b)(a+c)(b+d)(c+d)} \tag{4.10} +$$
+$(4.6)$ and $(4.7)$ have the advantage that there is no danger of computational problems if any marginal totals are near zero.
+The last three, $(4.8)$, $(4.9)$, $(4.10)$, are all related to the $\chi^2$ statistic, its square root, and the Pearson correlation coefficient respectively.
+Appealing features of monothetic divisive methods are the easy classification of new members and the including of cases with missing values.
+A further advantage of monothetic divisive methods is that it is obvious which variables produce the split at any stage of the process.
+A disadvantage with these methods is that the possession of a particular attribute which is either rare or rarely found in combination with others may take an individual down a different path.
+Recall in the previous section that we spoke about Monothetic and Polythetic methods. Monothetic methods only looks at a single variable at a time while Polythetic looks at multiple variables simultaneously. In this section, we will speak more about polythetic divisive methods.
+Polythetic methods operate via a distance matrix.
+This procedure avoids considering all possible splits by
+This process has been developed into a program named DIANA (DIvisive ANAlysis Clustering) which is implemented in R.
This somewhat resembles a way a political party might split due to inner conflicts.
+Firstly, the most discontented member leaves the party and starts a new one, and then some others follow him until a kind of equilibrium is attained.
+There are two common hierarchical methods used for large data sets BIRCH and CURE. Both of these algorithms employ a pre-clustering phase in where dense regions are summarized, the summaries being then clustered using a hierarchical method based on centroids.
CURE starts with a random sample of points and represents clusters by a smaller number of points that capture the shape of the clusterCURE has been shown to be able to cope with arbitrary-shaped clusters and in that respect may be superior to BIRCH, although it does require judgment as to the number of clusters and also a parameter which favors either more or less compact clusters.
In order to decide where clusters should be combined (for agglomerative), or where a cluster should be split (for divisive), a measure of dissimilarity between sets of observations is required.
+In most methods of hierarchical clustering this is achieved by a use of an appropriate
+Clustering using Representatives is a Hierarchical clustering technique in which you can represent a cluster using a set of well-scattered representative points.
+This algorithm has a parameter $\alpha$ which defines the factor of the points in which to shrink towards the centroid.
+CURE is known to be robust to outliers and able to identify clusters that have a non-spherical shape and size variance.
+The clusters with the closest pair of representatives are the clusters that are merged at each step of CURE's algorithm.
+This algorithm cannot be directly applied to large datasets due to high runtime complexity. Several enhancements were added to address this requirement
+Youtube Video: https://www.youtube.com/watch?v=JrOJspZ1CUw
+Steps
+TSNE allows us to reduce the dimensionality of a dataset to two which allows us to visualize the data.
+It is able to do this since many real-world datasets have a low intrinsic dimensionality embedded within the high-dimensional space.
+Since the technique needs to conserve the structure of the data, two corresponding mapped points must be close to each other distance wise as well. Let $|x_i - x_j|$ be the Euclidean distance between two data points, and $|y_i - yj|$ he distance between the map points. This conditional similarity between two data points is: +$$ +p{j|i} = \frac{exp(-|x_i-x_j|^2 / (2\sigmai^2))}{\sum{k \ne i}{exp(-|x_i-x_k|^2/(2\sigma_i^2))}} +$$ +Where we are considering the Gaussian distribution surrounding the distance between $x_j$ from $x_i$ with a given variance $\sigma_i^2$. The variance is different for every point; it is chosen such that points in dense areas are given a smaller variance than points in sparse areas.
+Now the similarity matrix for mapped points are +$$ +q_{ij} = \frac{f(|x_i - xj|)}{\sum{k \ne i}{f(|x_i - x_k)}} +$$ +Where $f(z) = \frac{1}{1 + z^2}$
+This has the same idea as the conditional similarity between two data points, except this is based on the Cauchy distribution.
+TSNE works at minimizing the Kullback-Leiber divergence between the two distributions $p{ij}$ and $q{ij}$ +$$ +KL(P || Q) = \sum{i,j}{p{i,j} \log{\frac{p{ij}}{q{ij}}}} +$$ +To minimize this score, gradient descent is typically performed +$$ +\frac{\partial KL(P||Q)}{\partial y_i} = 4\sumj{(p{ij} - q_{ij})} +$$
+There are multiple approaches to validating your cluster models
+Internal evaluation measures suffer form the problem that they represent functions that are objectives for many clustering algorithms. So of course the result of the clustering algorithm will be such that the objective would be minimized.
+External evaluation suffers from the fact that if we had labels to begin with then we wouldn't need to cluster. Practical applications of clustering occur usually when we don't have labels. On the other hand, possible labeling can reflect one possible partitioning of the data set. There could exist different, perhaps even better clustering.
+We like to see a few qualities in cluster models
+Let us focus on the second and third bullet point for now. Internal evaluation measures are best suited to get some insight into situations where one algorithm performs better than another, this does not imply that one algorithm produces more valid results than another.
+$$ +DB = \frac{1}{n}\sum{i=1}^n{max{j\ne i}{(\frac{\sigma_i + \sigma_j}{d(c_i,c_j)})}} +$$
+Where $n$ is the number of clusters, $c$ indicates a centroid, and $\sigma$ represents the deviation from the centroid.
+Better clustering algorithms are indicated by smaller DB values.
+$$ +D= \frac{min{1\le i < j \le n}{d(i,j)}}{max{1\le k \le n}{d^\prime(k)}} +$$
+The Dunn index aims to identify dense and well-separated clusters. This is defined as the ratio between the minimal inter-cluster distance to maximal intra-cluster distance.
+High Dunn Index values are more desirable.
+In terms of robustness we can measure uncertainty in each of the individual clusters. This can be examined using a bootstrapping approach by Suzuki and Shimodaira (2006). The probability or "p-value" is the proportion of bootstrapped samples that contain the cluster. Larger p-values in this case indicate more support for the cluster.
+This is available in R via Pvclust
One approach to assess the effects of perturbations of the data is by randomly dividing the data into two subsets and performing an analysis on each subset separately. This method was proposed by McIntyre and Blashfield in 1980; their method involves the following steps
+Using internal evaluation metrics, you can see the impact of each point by doing a "leave one out" analysis. Here you evaluate the dataset minus one point for each of the points. If a positive difference is found, the point is regarded as a facilitator, whereas if it is negative then it is considered an inhibitor. once an influential inhibitor is found, the suggestion is to normally omit it from the clustering.
+clValid contains a variety of internal validation measures.
Paper: https://cran.r-project.org/web/packages/clValid/vignettes/clValid.pdf
+Lecture 1: Measures of Similarity
+Lecture 2.1: Distance Measures Reasoning
+Lecture 2.2: Principle Component Analysis Pt. 1
+Lecture 3: Discussion of Dataset
+Lecture 4: Principal Component Analysis Pt. 2
+Lecture 4.2: Revisiting Measures
+ +Lecture 5: Introduction to Connectivity Based Models
+Lecture 6: Agglomerative Methods
+Lecture 7: Divisive Methods Part 1: Monothetic
+Lecture 8: Divisive Methods Part 2: Polythetic
+ +Lecture 9.2: Cluster Validation Part I
+Lecture 10.1: Silhouette Coefficient
+Lecture 10.2: Centroid-Based Clustering
+ + + + + +Garson Textbook Chapter 3
+A Tutorial on Principal Component Analysis
+No Readings
+An Introduction to Applied Multivariate Analysis with R by Brian Evveritt and Torsten Horthorn.
+Sections 3.0-3.9 Everitt
+Section 4.1 Everitt
+Section 4.2 Everitt
+Applied Multivariate Statistical Analysis Johnson Section 12.3
+Linkage Methods for Cluster Observations
+Section 4.3 Everitt
+Introduction to the TSNE Algorithm
+Section 9.5 Everitt
+Cluster Analysis 5th Edition
+Cluster Analysis: 2014 Edition (Statistical Associates Blue Book Series 24)
+In an ideal world, the topics below I estimated being a certain time period for learning them. Of course you have more experience when it comes to how long it actually takes to learn these topics, so I'll leave this mostly to your discretion.
+Distance, Dimensionality Reduction, and Tendency -- 3 Weeks
+Validating Cluster Models -- 1 Week
+Connectivity Models -- 2 Weeks
+Cluster Evaluation -- 1 Week
+Centroid Models -- 3 Weeks
+Density Models -- 3 Weeks
+Analysis of Model Appropriateness -- 1 Week
+The schedule above accounts for 14 weeks, so there is a week that is free as a buffer.
+Creating this document got me really excited for this independent study. Feel free to give me feedback :)
+Cluster Analysis is the art of finding inherent structures in data to form groups of similar observations. This has a myriad of applications from recommendation engines to social network analysis.
+This is an independent study, meaning that I will be studying this topic under the direction of a professor, in this case being Dr. Denhere.
+I have provided a list of topics that I wish to explore in a syllabus
+Dr. Denhere likes to approach independent studies from a theoretical and applied sense. Meaning, I will learn the theory of the different algorithms, and then figure out a way to apply them onto a dataset.
+There is no definitive textbook for this course. Instead I and Dr. Denhere search for materials that we think best demonstrates the topic at hand.
+I have created a Reading Page to keep track of the different reading materials.
+I like to type of the content I learn from different sources. A notes page is created to keep track of the content discussed each meeting.
+Reinforcement learning evaluates the actions taken rather than accepting $instructions$ of the correct actions. This creates the need for active exploration.
+This chapter of the book goes over a simplified version of the reinforcement learning problem, that does not involve learning to act in more than one situation. This is called a nonassociative setting.
+In summation, the type of problem we are about to study is a nonassociative, evaluative feedback problem that is a simplified version of the $k$-armed bandit problem.
+Consider the following learning problem. You are faced repeatedly with a choice among $k$ different options/actions. After each choice you receive a numerical reward chosen from a stationary probability distribution that depends on the action you selected.
+Your objective (if you choose to accept it) is to maximize the expected total reward over some time period. Let's say $1000$ time steps.
+This is called the $k$-armed bandit problem because it's an analogy of a slot machine. Slot machines are nick-named the "one-armed bandit", and the goal here is to play the slot machine that has the greatest value return.
+We want to figure out which slot machine produces the greatest value. Therefore, we want to be able to estimate the value of a slot machine as close to the actual value as possible.
+If we maintain estimates of the action values, then at any time step there is at least one action whose estimated value is the greatest. We call these greedy actions. When you select one of these actions we say that you are exploiting your current knowledge of the values of the actions.
+If you instead select a non-greedy action, then you are exploring, because this enables you to better improve your estimate of the non-greedy action's value.
+Uncertainty is such that at least one of the other actions probably better than the greedy action, you just don't know which one yet.
+In this section, we will look at simple balancing methods in how to gain the greatest reward through exploration and exploitation.
+We begin by looking more closely at some simple methods for estimating the values of actions and for using said estimates to make action selection decisions.
+One natural way to estimate this is by averaging the rewards actually received +$$ +Qt(a) = \frac{\sum{i = 1}^{t - 1}Ri * \mathbb{I}{Ai = 1}}{\sum{i = 1}^{t - 1}\mathbb{I}_{Ai = 1}} +$$ +where $\mathbb{I}{predicate}$ denotes the random variable that is 1 if the predicate is true and 0 if it is not. If the denominator is zero (we have not experienced the reward), then we assume some default value such as zero.
+This is where you choose greedily all the time. +$$ +A_t = argmax_a(Q_t(a)) +$$
+This is where we choose greedily most of the time, except for a small probability $\epsilon$. Where instead of selecting greedily, we select randomly from among all the actions with equal probability.
+The advantage of $\epsilon$-greedy over greedy methods depends on the task. With noisier rewards it takes more exploration to find the optimal action, and $\epsilon$-greedy methods should fare better relative to the greedy method. However, if the reward variances were zero, then the greedy method would know the true value of each action after trying it once.
+Suppose the bandit task were non-stationary, that is, the true values of actions changed over time. In this case exploration is needed to make sure one of the non-greedy actions has not changed to become better than the greedy one.
+There is a way to update averages using small constant computations rather than storing the the numerators and denominators separate.
+Note the derivation for the update formula +$$ +\begin{align} +Q{n + 1} &= \frac{1}{n}\sum{i = 1}^n{R_i} \ +&= \frac{1}{n}(Rn + \sum{i = 1}^{n - 1}{R_i}) \ +&= \frac{1}{n}(Rn + (n - 1)\frac{1}{n-1}\sum{i = 1}^{n - 1}{R_i}) \ +&= \frac{1}{n}{R_n + (n - 1)Q_n} \ +&= \frac{1}{n}(R_n + nQ_n - Q_n) \ +&= Q_n + \frac{1}{n}(R_n - Q_n) \tag{2.3} +\end{align} +$$ +With formula 2.3, the implementation requires memory of only $Q_n$ and $n$.
+This update rule is a form that occurs frequently throughout the book. The general form is +$$ +NewEstimate = OldEstimate + StepSize(Target - OldEstimate) +$$
+As noted earlier, we often encounter problems that are nonstationary, in such cases it makes sense to give more weight to recent rewards than to long-past rewards. One of the most popular ways to do this is to use a constant value for the $StepSize$ parameter. We modify formula 2.3 to be +$$ +\begin{align} +Q_{n + 1} &= Q_n + \alpha(R_n - Q_n) \ +&= \alpha R_n + Q_n - \alpha Q_n \ +&= \alpha R_n + (1 - \alpha)Q_n \ +&= \alpha Rn + (1 - \alpha)(\alpha R{n - 1} + (1-\alpha)Q_{n - 1}) \ +&= \alpha Rn + (1 - \alpha)(\alpha R{n - 1} + (1-\alpha)(\alpha R{n - 2} + (1 - a)Q{n - 2})) \ +&= \alpha Rn + (1-\alpha)\alpha R{n - 1} + (1-\alpha)^2\alpha R_{n - 2} + \dots + (1-\alpha)^nQ_1 \ +&= (1-\alpha)^nQ1 + \sum{i = 1}^n{\alpha(1-\alpha)^{n - i}R_i} +\end{align} +$$ +This is a weighted average since the summation of all the weights equal one. Note here that the farther away a value is from the current time, the more times $(1-\alpha)$ gets multiplied by itself. Hence making it less influential. This is sometimes called an exponential recency-weighted average.
+Sometimes it is convenient to vary the step-size parameter from step to step. We can denote $\alpha_n(a)$ to be a function that determines the step-size parameter after the $n$th selection of action $a$. As noted before $\alpha_n(a) = \frac{1}{n}$ results in the sample average method which is guaranteed to converge to the truth action values assuming a large number of trials.
+A well known result in stochastic approximation theory gives us the following conditions to assure convergence with probability 1: +$$ +\sum_{n = 1}^\infty{\alphan(a) = \infty} \and \sum{n = 1}^{\infty}{\alpha_n^2(a) \lt \infty} +$$ +The first condition is required to guarantee that the steps are large enough to overcome any initial conditions or random fluctuations. The second condition guarantees that eventually the steps become small enough to assure convergence.
+Note: Both convergence conditions are met for the sample-average case but not for the constant step-size parameter. The latter condition is violated in the constant parameter case. This is desirable since if the rewards are changing then we don't want it to converge to any one parameter.
+The methods discussed so far are biased by their initial estimates. Another downside is that these values are another set of parameters that must be chosen by the user. Though these initial values can be used as a simple way to encourage exploration.
+Let's say you set an initial estimate that is wildly optimistic. Whichever actions are initially selected, the reward is less than the starting estimates. Therefore, the learner switches to other actions, being disappointed with the rewards it was receiving.
+The result of this is that all actions are tried several times before their values converge. It even does this if the algorithm is set to choose greedily most of the time!
+
This simple trick is quite effective for stationary problems. Not so much for nonstationary problems since the drive for exploration only happens at the beginning. If the task changes, creating a renewed need for exploration, this method would not catch it.
+Exploration is needed because there is always uncertainty about the accuracy of the action-value estimates. The greedy actions are those that look best at the present but some other options may actually be better. Let's choose options that have potential for being optimal, taking into account how close their estimates are to being maximal and the uncertainties in those estimates. +$$ +A_t = argmax_a{(Q_t(a) + c\sqrt{\frac{ln(t)}{N_t(a)}})} +$$ +where $N_t(a)$ denotes the number of times that $a$ has been selected prior to time $t$ and $c > 0$ controls the degree of exploration.
+So far, we've only considered nonassociative tasks, where there is no need to associate different actions with different situations. However, in a general reinforcement learning task there is more than one situation and the goal is to learn a policy: a mapping from situations to the actions that are best in those situations.
+For sake of continuing our example, let us suppose that there are several different $k$-armed bandit tasks, and that on each step you confront one of these chosen at random. To you, this would appear as a single, nonstationry $k$-armed bandit task whose true action values change randomly from step to step. You could try using one of the previous methods, but unless the true action values change slowly, these methods will not work very well.
+Now suppose, that when a bandit task is selected for you, you are given some clue about its identity. Now you can learn a policy association each task, singled by the clue, with the best action to take when facing that task.
+This is an example of an associative search task, so called because it involves both trial-and-error learning to search for the best actions, and association of these actions with situations in which they are best. Nowadays they're called contextual bandits in literature.
+If actions are allowed to affect the next situation as well as the reward, then we have the full reinforcement learning problem. This will be presented in the next chapter of the book with its ramifications appearing throughout the rest of the book.
+
Dynamic programming refers to a collection of algorithms that can be used to compute optimal policies given a perfect model of the environment as a Markov decision process (MDP).
+Classic DP algorithms are of limited utility due to their assumption of a perfect model and their great computational expense.
+Let's assume that the environment is a finite MDP. We assume its state, action, and reward sets, $\mathcal{S}, \mathcal{A}, \mathcal{R}$ are finite, and that its dynamics are given by a set of probabilities $p(s^\prime, r | s , a)$.
+The key idea of dynamic programming, and of reinforcement learning is the use of value functions to organize and structure the search for good policies. In this chapter, we show how dynamic programming can be used to compute the value functions defined in Chapter 3. We can easily obtain optimal policies once we have found the optimal value functions which satisfy the Bellman optimality equations.
+First we consider how to compute the state-value function for an arbitrary policy. The existence and uniqueness of a state-value function for an arbitrary policy are guaranteed as long as either the discount rate is less than one or eventual termination is guaranteed from all states under the given policy.
+Consider a sequence of approximate value functions. The initial approximation, $v0$, is chosen arbitrarily (except that the terminal state must be given a value of zero) and each successive approximation is obtained by using the Bellman equation for $v\pi$ as an update rule: +$$ +v{k + 1} = \sum{a}{\pi(a |s)\sum_{s^\prime, r}{p(s^\prime,r|s,a)[r + \gamma v_k(s^\prime)]}} +$$ +This algorithm is called iterative policy evaluation.
+To produce each successive approximation, $v_{k + 1}$ from $v_k$, iterative policy evaluation applies the same operation to each state $s$: it replaces the old value of $s$ with a new value obtained from the old values of the successor states of $s$, and the expected immediate rewards, along all the one-step transitions possible under the policy being evaluated.
+Iterative Policy Evaluation
+
Input π, the policy to be evaluated
+Initialize an array V(s) = 0, for all s ∈ S+
+Repeat
+ ∆ ← 0
+ For each s ∈ S:
+ v ← V(s)
+ V(s) ← ∑_a π(a|s) ∑_s′,r p(s′,r|s,a)[r+γV(s′)]
+ ∆ ← max(∆,|v−V(s)|)
+until ∆ < θ (a small positive number)
+Output V ≈ v_πConsider a grid world where the top left and bottom right squares are the terminal state. Now consider that every other square, produces a reward of -1, and the available actions on each state is {up, down, left, right} as long as that action keeps the agent on the grid. Suppose our agent follows the equiprobable random policy.
+
+One reason for computing the value function for a policy is to help find better policies. Suppose we have determined the value function $v\pi$ for an arbitrary deterministic policy $\pi$. For some state $s$ we would like to know whether or not we should change the policy to determinatically chose another action. We know how good it is to follow the current policy from $s$, that is $v\pi(s)$, but would it be better or worse to change to the new policy?
+One way to answer this question is to consider selecting $a$ in $s$ and thereafter follow the existing policy, $\pi$. The key criterion is whether the this produces a value greater than or less than $v_\pi(s)$. If it is greater, then one would expect it to be better still to select $a$ every time $s$ is encountered, and that the new policy would in fact be a better one overall.
+That this is true is a special case of a general result called the policy improvement theorem. Let $\pi$ and $\pi^\prime$ be any pair of deterministic policies such that for all $s \in \mathcal{S}$. +$$ +q\pi(s, \pi^\prime(s)) \ge v\pi(s) +$$ +So far we have seen how, given a policy and its value function, we can easily evaluate a change in the policy at a single state to a particular action. It is a natural extension to consider changes at all states and to all possible actions, selecting at each state the action that appears best according to $q\pi(s, a)$. In other words, to consider the new greedy policy, $\pi^\prime$, given by: +$$ +\pi^\prime = argmax (q\pi(s, a)) +$$ +So far in this section we have considered the case of deterministic policies. We will not go through the details, but in fact all the ideas of this section extend easily to stochastic policies.
+Once a policy, $\pi$, has been improved using $v\pi$ to yield a better policy, $\pi^\prime$, we can then compute $v{\pi^\prime}$ and improve it again to yield an even better $\pi^{\prime\prime}$. We can thus obtain a sequence of monotonically improving policies and value functions.
+Each policy is guaranteed to be a strict improvement over the previous one (unless its already optimal). Since a finite MDP has only a finite number of policies, this process must converge to an optimal policy and optimal value function in a finite number of iterations.
+This way of finding an optimal policy is called policy iteration.
+Algorithm
+1. Initialization
+ V(s) ∈ R and π(s) ∈ A(s) arbitrarily for all s ∈ S
+2. Policy Evaluation
+ Repeat
+ ∆ ← 0
+ For each s∈S:
+ v ← V(s)
+ V(s) ← ∑_{s′,r} p(s′,r|s,π(s))[r + γV(s′)]
+ ∆ ← max(∆,|v − V(s)|)
+ until ∆ < θ (a small positive number)
+3. Policy Improvement
+ policy-stable ← true
+ For each s ∈ S:
+ old-action ← π(s)
+ π(s) ← arg max_a ∑_{s′,r} p(s′,r|s,a)[r + γV(s′)]
+ If old-action != π(s), then policy-stable ← false
+ If policy-stable, then stop and return V ≈ v_∗,
+ and π ≈ π_∗; else go to 2One drawback to policy iteration is that each of its iterations involve policy evaluation, which may itself be a protracted iterative computation requiring multiple sweeps through the state set. If policy evaluation is done iteratively, then convergence exactly to $v_\pi$ occurs only in the limit. Must we wait for exact convergence, or can we stop short of that?
+In fact, the policy evaluation step of policy iteration can be truncated in several ways without losing the convergence guarantee of policy iteration. One important special case is when policy evaluation is stopped after one sweep. This algorithm is called value iteration.
+Another way of understanding value iteration is by reference to the Bellman optimality equation. Note that value iteration is obtained simply by turning the Bellman optimality equation into an update rule. Also note how value iteration is identical to the policy evaluation update except that it requires the maximum to be taken over all actions.
+Finally, let us consider how value iteration terminates. Like policy evaluation, value iteration formally requires an infinite number of iterations to converge exactly. In practice, we stop once the value function changes by only a small amount.
+Initialize array V arbitrarily (e.g., V(s) = 0 for all
+s ∈ S+)
+
+Repeat
+ ∆ ← 0
+ For each s ∈ S:
+ v ← V(s)
+ V(s) ← max_a∑_{s′,r} p(s′,r|s,a)[r + γV(s′)]
+ ∆ ← max(∆, |v − V(s)|)
+until ∆ < θ (a small positive number)
+
+Output a deterministic policy, π ≈ π_∗, such that
+ π(s) = arg max_a ∑_{s′,r} p(s′,r|s,a)[r + γV(s′)]Value iteration effectively combines, in each of its sweeps, one sweep of policy evaluation and one sweep of policy improvement. Faster convergence is often achieved by interposing multiple policy evaluation sweeps between each policy improvement sweep.
+Asynchronous DP algorithms are in-place DP algorithms that are not organized in terms of systematic sweeps of the state set. These algorithms update the values of states in any order whatsoever, using whatever value of other states happen to be available.
+To converge correctly, however, an asynchronous algorithm must continue to update the value of all the states: it can't ignore any state after some point in the computation.
+Policy iteration consists of two simultaneous, iterating processes, one making the value function consistent with the current policy (policy evaluation) and the other making the policy greedy with respect to the current value function (policy improvement).
+We use the term generalized policy iteration (GPI) to competing and cooperating. They compete in the sense that they pull in opposing directions. Making the policy greedy with respect to the value function typically makes the value function incorrect for the changed policy. Making the value function consistent with the policy typically causes the policy to be greedy. In the long run, however, the two processes interact to find a single joint solution.
+Recall that this course is based on the book --
+Reinforcement Learning: An Introduction by Richard S. Sutton and Andrew G. Barto
+These notes really serve as talking points for the overall concepts described in the chapter and are not meant to stand for themselves. Check out the book for more complete thoughts :)
+Reinforcement Learning is learning what to do -- how to map situations to actions -- so as to maximize a numerical reward signal. There are two characteristics, trial-and-error search and delayed reward, that are the two most important distinguishing features of reinforcement learning.
+Markov decision processes are intended to include just these three aspects: sensation, action, and goal(s).
+Reinforcement learning is different than the following categories
+One of the challenges that arise in reinforcement learning is the trade-off between exploration and exploitation. The dilemma is that neither exploration nor exploitation can be pursued exclusively without failing at the task.
+Another key feature of reinforcement learning is that it explicitly considers the whole problem of a goal-directed agent interacting with an uncertain environment. This is different than supervised learning since they're concerned with finding the best classifier/regression without explicitly specifying how such an ability would finally be useful.
+A complete, interactive, goal-seeking agent can also be a component of a larger behaving system. A simple example is an agent that monitors the charge level of a robot's battery and sends commands to the robot's control architecture. This agent's environment is the rest of the robot together with the robot's environment.
+A policy defines the learning agent's way of behaving at a given time
+A reward signal defines the goal in a reinforcement learning problem. The agent's sole objective is to maximize the total reward it receives over the long run
+A value function specifies what is good in the long run. Roughly speaking, the value of a state is the total amount of reward an agent can expect to accumulate over the future, starting from that state. Without rewards there could be no value,s and the only purpose of estimating values is to achieve more reward. We seek actions that bring about states of highest value.
+Unfortunately, it is much harder to determine values than it is to determine rewards. The most important component of almost all reinforcement learning algorithms we consider is a method for efficiently estimating values.
+Look at Tic-Tac-Toe example
+Most of the time in a reinforcement learning algorithm, we move greedily, selecting the move that leads to the state with greatest value. Occasionally, however, we select randomly from amoung the other moves instead. These are called exploratory moves because they cause us to experience states that we might otherwise never see.
+Summary: Reinforcement learning is learning by an agent from direct interaction wit its environment, without relying on exemplary supervision or complete models of the environment.
+Monte Carlo methods do not assume complete knowledge of the environment. They require only experience which is a sample sequence of states, actions, and rewards from actual or simulated interaction with an environment.
+Monte Carlo methods are ways of solving the reinforcement learning problem based on averaging sample returns. To ensure that well-defined returns are available, we define Monte Carlo methods only for episodic tasks. Only on the completion of an episode are value estimates and policies changed.
+Monte Carlo methods sample and average returns for each state-action pair is much like the bandit methods explored earlier. The main difference is that there are now multiple states, each acting like a different bandit problems and the problems are interrelated. Due to all the action selections undergoing learning, the problem becomes nonstationary from the point of view of the earlier state.
+Recall that the value of a state is the expected return -- expected cumulative future discounted reward - starting from that state. One way to do it is to estimate it from experience by averaging the returns observed after visits to that state.
+Each occurrence of state $s$ in an episode is called a visit to $s$. The first-visit MC method estimates $v_\pi(s)$ as the average of the returns following first visits to $s$, whereas the every-visit MC method averages the returns following all visits to $s$. These two Monte Carlo methods are very similar but have slightly different theoretical properties.
+First-visit MC prediction
+Initialize:
+ π ← policy to be evaluated
+ V ← an arbitrary state-value function
+ Returns(s) ← an empty list, for all s ∈ S
+
+Repeat forever:
+ Generate an episode using π
+ For each state s appearing in the episode:
+ G ← the return that follows the first occurrence of
+s
+ Append G to Returns(s)
+ V(s) ← average(Returns(s))If a model is not available then it is particularly useful to estimate action values rather than state values. With a model, state values alone are sufficient to define a policy. Without a model, however, state values alone are not sufficient. One must explicitly estimate the value of each action in order for the values to be useful in suggesting a policy.
+The only complication is that many state-action pairs may never be visited. If $\pi$ is a deterministic policy, then in following $\pi$ one will observe returns only for one of the actions from each state. With no returns to average, the Monte Carlo estimates of the other actions will not improve with experience. This is a serious problem because the purpose of learning action values is to help in choosing among the actions available in each state.
+This is the general problem of maintaining exploration. For policy evaluation to work for action values, we must assure continual exploration. One way to do this is by specifying that the episodes start in a state-action pair, and that each pair has a nonzero probability of being selected as the start. We call this the assumption of exploring starts.
+We made two unlikely assumptions above in order to easily obtain this guarantee of convergence for the Monte Carlo method. One was that the episodes have exploring starts, and the other was that policy evaluation could be done with an infinite number of episodes.
+Monte Carlo Exploring Starts
+Initialize, for all s ∈ S, a ∈ A(s):
+ Q(s,a) ← arbitrary
+ π(s) ← arbitrary
+ Returns(s,a) ← empty list
+
+Repeat forever:
+ Choose S_0 ∈ S and A_0 ∈ A(S_0) s.t. all pairs have probability > 0
+ Generate an episode starting from S_0,A_0, following
+π
+ For each pair s,a appearing in the episode:
+ G ← the return that follows the first occurrence of
+s,a
+ Append G to Returns(s,a)
+ Q(s,a) ← average(Returns(s,a))
+ For each s in the episode:
+ π(s) ← arg max_a Q(s,a)The only general way to ensure that actions are selected infinitely often is for the agent to continue to select them. There are two approaches ensuring this, resulting in what we call on-policy methods and off-policy methods.
+On-policy methods attempt to evaluate or improve the policy that is used to make decisions, whereas off-policy methods evaluate or improve a policy different from that used to generate the data.
+In on-policy control methods the policy is generally soft, meaning that $\pi(a|s)$ for all $a \in \mathcal{A}(s)$. The on-policy methods in this section uses $\epsilon$-greedy policies, meaning that most of the time they choose an action that has maximal estimated action value, but with probability $\epsilon$ they instead select an action at random.
+On-policy first-visit MC control (for $\epsilon$-soft policies)
+Initialize, for all $s \in \mathcal{S}$, $a \in \mathcal{A}(s)$:
+$Q(s, a)$ ← arbitrary
+$Returns(s,a)$ ← empty list
+$\pi(a|s)$ ← an arbitrary $\epsilon$-soft policy
+Repeat forever:
+(a) Generate an episode using $\pi$
+(b) For each pair $s,a$ appearing in the episoe
+ $G$ ← the return that follows the first occurance of s, a
+ Append $G$ to $Returns(s,a)$
+ $Q(s, a)$ ← average($Returns(s,a)$)
+(c) For each $s$ in the episode:
+ $A^*$ ← argmax$_a Q(s,a)$ (with ties broken arbitrarily)
+ For all $a \in \mathcal{A}(s)$:
+ $\pi(a|s)$ ← $\begin{cases} 1 - \epsilon + \epsilon / |\mathcal{A}(s)| & a = A^ \ \epsilon / | \mathcal{A}(s)| & a \neq A^ \end{cases}$
+Markov Decision processes are a classical formalization of sequential decision making, where actions influence not just immediate rewards, but also subsequent situations, or states, and through those future rewards. Thus MDPs involve delayed reward and the need to trade-off immediate and delayed reward. Whereas in bandit problems we estimated the value of $q*(a)$ of each action $a$, in MDPs we estimate the value of $q*(s, a)$ of each action $a$ in state $s$.
+MDPs are a mathematically idealized form of the reinforcement learning problem. As we will see in artificial intelligence, there is often a tension between breadth of applicability and mathematical tractability. This chapter will introduce this tension and discuss some of the trade-offs and challenges that it implies.
+The learner and decision maker is called the agent. The thing it interacts with is called the environment. These interact continually, the agent selecting actions and the environment responding to these actions and presenting new situations to the agent.
+The environment also give rise to rewards, special numerical values that the agent seeks to maximize over time through its choice of actions.
+
+To make the future paragraphs clearer, a Markov Decision Process is a discrete time stochastic control process. It provides a mathematical framework for modeling decision making in situations where outcomes are partly random and partly under the control of the decision maker.
+In a finite MDP, the set of states, actions, and rewards all have a finite number of elements. In this case, the random variables R_t and St have a well defined discrete probability distribution dependent only on the preceding state and action. +$$ +p(s^\prime | s,a) = \sum{r \in \mathcal{R}}{p(s^\prime, r|s, a)} +$$ +Breaking down the above formula, it's just an instantiation of the law of total probability. If you partition the probabilistic space by the reward, if you sum up each partition you should get the overall space. This formula has a special name: the state-transition probability.
+From this we can compute the expected rewards for each state-action pair by multiplying each reward with the probability of getting said reward and summing it all up. +$$ +r(s, a) = \sum{r \in \mathcal{R}}{r}\sum{s^\prime \in \mathcal{S}}{p(s^\prime, r|s,a)} +$$ +The expected reward for a state-action-next-state triple is +$$ +r(s, a, s^\prime) = \sum_{r \in \mathcal{R}}{r\frac{p(s^\prime, r|s,a)}{p(s^\prime|s,a)}} +$$ +I wasn't able to piece together this function in my head like the others. Currently I imagine it as if we took the formula before and turned the universe of discourse from the universal set to the state where $s^\prime$ occurred.
+The MDP framework is abstract and flexible and can be applied to many different problems in many different ways. For example, the time steps need not refer to fixed intervals of real time; they can refer to arbitrary successive states of decision making and acting.
+In particular, the boundary between agent and environment is typically not the same as the physical boundary of a robot's or animals' body. Usually, the boundary is drawn closer to the agent than that. For example, the motors and mechanical linkages of a robot and its sensing hardware should usually be considered parts of the environment rather than parts of the agent.
+The general rule we follow is that anything that cannot be changed arbitrarily by the agent is considered to be outside of it and thus part of its environment. We do not assume that everything in the environment is unknown to the agent. For example, the agent often knows quite a bit about how its rewards are computed as a function of its actions and the states in which they are taken. But we always consider the reward computation to be external to the agent because it defines the task facing the agent and thus must be beyond its ability to change arbitrarily. The agent-environment boundary represents the limit of the agent's absolute control, not of its knowledge.
+This framework breaks down whatever the agent is trying to achieve to three signals passing back and forth between an agent and its environment: one signal to represent the choices made by the agent, one signal to represent the basis on which the choices are made (the states), and one signal to define the agent's goal (the rewards).
+Recall that the agent makes a decision at times determined by external events. At each such time the robot decides whether it should
+(1) Actively search for a can
+(2) Remain stationary and wait for someone to bring it a can
+(3) Go back to home base to recharge its battery
+Suppose the environment works as follows: the best way to find cans is to actively search for them, but this runs down the robot's battery, whereas waiting does not. Whenever the robot is searching the possibility exists that its battery will become depleted. In this case, the robot must shut down and wait to be rescued (producing a low reward.)
+The agent makes its decisions solely as a function of the energy level of the battery, It can distinguish two levels, high and low, so that the state set is $\mathcal{S} = {high, low}$. Let us call the possible decisions -- the agent's actions -- wait, search, and recharge. When the energy level is high, recharging would always be foolish, so we do not include it in the action set for this state. The agent's action sets are +$$ +\begin{align} +\mathcal{A}(high) &= {search, wait} \ +\mathcal{A}(low) &= {search, wait, recharge} +\end{align} +$$ +If the energy level is high, then a period of active search can always be completed without a risk of depleting the battery. A period of searching that begins with a high energy level leaves the energy level high with a probability of $\alpha$ and reduces it to low with a probability of $(1 - \alpha)$. On the other hand, a period of searching undertaken when the energy level is low leaves it low with a probability of $\beta$ and depletes the battery with a probability of $(1 - \beta)$. In the latter case, the robot must be rescued, and the batter is then recharged back to high.
+Each can collected by the robot counts as a unit reward, whereas a reward of $-3$ occurs whenever the robot has to be rescued. Let $r{search}$ and $r{wait}$ with $r{search } > r{wait}$, respectively denote the expected number of cans the robot will collect while searching and waiting. Finally, to keep things simple, suppose that no cans can be collected ruing a run home for recharging and that no cans can be collected on a step in which the battery is depleted.
+| $s$ | +$a$ | +$s^\prime$ | +$p(s^\prime | +s, a)$ | +$r(s | +
|---|---|---|---|---|---|
| high | +search | +high | +$\alpha$ | +$r_{search}$ | +|
| high | +search | +low | +$(1-\alpha)$ | +$r_{search}$ | +|
| low | +search | +high | +$(1 - \beta)$ | +-3 | +|
| low | +search | +low | +$\beta$ | +$r_{search}$ | +|
| high | +wait | +high | +1 | +$r_{wait}$ | +|
| high | +wait | +low | +0 | +$r_{wait}$ | +|
| low | +wait | +high | +0 | +$r_{wait}$ | +|
| low | +wait | +low | +1 | +$r_{wait}$ | +|
| low | +recharge | +high | +1 | +0 | +|
| low | +recharge | +low | +0 | +0 | +
A transition graph is a useful way to summarize the dynamics of a finite MDP. There are two kinds of nodes: state nodes and action nodes. There is a state node for each possible state and an action node for reach state-action pair. Starting in state $s$ and taking action $a$ moves you along the line from state node $s$ to action node $(s, a)$. The the environment responds with a transition ot the next states node via one of the arrows leaving action node $(s, a)$.
+
+The reward hypothesis is that all of what we mean by goals and purposes can be well thought of as the maximization of the expected value of the cumulative sum of a received scalar signal called the reward.
+Although formulating goals in terms of reward signals might at first appear limiting, in practice it has proved to be flexible and widely applicable. The best way to see this is to consider the examples of how it has been, or could be used. For example:
+It is critical that the rewards we set up truly indicate what we want accomplished. In particular, the reward signal is not the place to impart to the agent prior knowledge about how to achieve what we want it to do. For example, a chess playing agent should only be rewarded for actually winning, not for achieving subgoals such as taking its opponent's pieces. If achieving these sort of subgoals were rewarded, then the agent might find a way to achieve them without achieving the real goal. The reward signal is your way of communicating what you want it to achieve, not how you want it achieved.
+In general, we seek to maximize the expected return, where the return is defined as some specific function of the reward sequence. In the simplest case, the return is the sum of the rewards: +$$ +Gt = R{t + 1} + R{t + 2} + R{t + 3} + \dots + R_{T} +$$ +where $T$ is the final time step. This approach makes sense in applications in which there is a natural notion of a final time step. That is when the agent-environment interaction breaks naturally into subsequences or episodes, such as plays of a game, trips through a maze, or any sort of repeated interaction.
+Each episode ends in a special state called the terminal state, followed by a reset to the standard starting state or to a sample from a standard distribution of starting states. Even if you think of episodes as ending in different ways the next episode begins independently of how the previous one ended. Therefore, the episodes can all be considered to ending the same terminal states with different rewards for different outcomes.
+Tasks with episodes of this kind are called episodic tasks. In episodic tasks we sometimes need to distinguish the set of all nonterminal states, denoted $\mathcal{S}$, from the set of all states plus the terminal state, denoted $\mathcal{S}^+$. The time of termination, $T$, is a random variable that normally varies from episode to episode.
+On the other hand, in many cases the agent-environment interaction goes on continually without limit. We call these continuing tasks. The return formulation above is problematic for continuing tasks because the final time step would be $T = \infty$, and the return which we are trying to maximize, could itself easily be infinite. The additional concept that we need is that of discounting. According to this approach, the agent tries to select actions so that the sum of the discounted rewards it receives over the future is maximized. In particular, it chooses $A_t$ to maximize the expected discounted return. +$$ +Gt= \sum{k = 0}^\infty{\gamma^kR_{t+k+1}} +$$ +where $\gamma$ that is a parameter between $0$ and $1$ is called the discount rate.
+The discount rate determines the present value of future rewards: a reward received $k$ time steps in the future is worth only $\gamma^{k - 1}$ time what it would be worth if it were received immediately. If $\gamma < 1$, the infinite sum has a finite value as long as the reward sequence is bounded.
+If $\gamma = 0$, the agent is "myopic" in being concerned only with maximizing immediate rewards. But in general, acting to maximize immediate reward can reduce access to future rewards so that the return is reduced.
+As $\gamma$ approaches 1, the return objective takes future rewards into account more strongly; the agent becomes more farsighted.
+The objective in this task is to apply forces to a car moving along a track so as to keep a pole hinged to the cart from falling over.
+A failure is said to occur if the pole falls past a given angle from the vertical or if the cart runs off the track.
+
+The reward can be a $+1$ for every time step on which failure did not occur. In this case, successful balancing would mean a return of infinity.
+The reward can be $-1$ on each failure and zero all other times. The return at each time would then be related to $-\gamma^k$ where $k$ is the number of steps before failure.
+Either case the return is maximized by keeping the pole balanced for as long as possible.
+Almost all reinforcement learning algorithms involve estimating value functions which estimate what future rewards can be expected. Of course the rewards that the agent can expect to receive is dependent on the actions it will take. Accordingly, value functions are defined with respect to particular ways of acting, called policies.
+Formally, a policy is a mapping from states to probabilities of selecting each possible action. The value of a state s under a policy \pi, denoted $v{\pi}(s)$ is the expected return when starting in $s$ and following $\pi$ thereafter. For MDPs we can define $v{\pi}$ as
+$$ +v{\pi}(s) = \mathbb{E}{\pi}[G_t | St = s] = \mathbb{E}{\pi}[\sum{k = 0}^\infty{\gamma^kR{t+k+1} | S_t = s}] +$$
+We call this function the state-value function for policy $\pi$. Similarly, we define the value of taking action $a$ in state $s$ under a policy $\pi$, denoted as $q_\pi(s,a)$ as the expected return starting from $s$, taking the action $a$, and thereafter following the policy $\pi$. Succinctly, this is called the action-value function for policy $\pi$.
+For some kinds of tasks we are interested in,optimal policies can be generated only with extreme computational cost. A critical aspect of the problemfacing the agent is always teh computational power available to it, in particular, the amount of computation it can perform in a single time step.
+The memory available is also an important constraint. A large amount of memory is often required to build up approximations of value functions, policies, and models. In the case of large state sets, functions must be approximated using some sort of more compact parameterized function representation.
+This presents us with unique oppurtunities for achieving useful approximations. For example, in approximating optimal behavior, there may be many states that the agent faces with such a low probability that selecting suboptimal actions for them has little impact on the amount of reward the agent receives.
+The online nature of reinforcement learning which makes it possible to approximate optimal policies in ways that put more effort into learning to make good decisions for frequently encountered states at the expense of infrequent ones is the key property that distinguishes reinforcement learning from other approaches to approximately solve MDPs.
+Let us summarize the elements of the reinforcement learning problem.
+Reinforcement learning is about learning from interaction how to behave in order to achieve a goal. The reinforcement learning agent and its environment interact over a sequence of discrete time steps.
+The actions are the choices made by the agent; the states are the basis for making the choice; and the rewards are the basis for evaluating the choices.
+Everything inside the agent is completely known and controllable by the agent; everything outside is incompletely controllable but may or may not be completely known.
+A policy is a stochastic rule by which the agent selects actions as a function of states.
+When the reinforcement learning setup described above is formulated with well defined transition probabilities it constitutes a Markov Decision Process (MDP)
+The return is the function of future rewards that the agent seeks to maximize. It has several different definitions depending on the nature of the task and whether one wishes to discount delayed reward.
+A policy's value functions assign to each state, or state-action pair, the expected return from that state, or state-action pair, given that the agent uses the policy. The optimal value functions assign to each state, or state-action pair, the largest expected return achievable by any policy. A policy whose value unctions are optimal is an optimal policy.
+Even if the agent has a complete and accurate environment model, the agent is typically unable to perform enough computation per time step to fully use it. The memory available is also an important constraint. Memory may be required to build up accurate approximations of value functions, policies, and models. In most cases of practical interest there are far more states that could possibly be entries in a table, and approximations must be made.
+Chapter 2: Multi-armed Bandits
+Chapter 3: Markov Decision Processes
+ + +Chapter 1: What is Reinforcement Learning?
+Chapter 2: Multi-armed Bandits
+Chpater 3: Finite Markov Decision Processes Part 1
+Chapter 3: Finite Markov Decision Processes Part 2
+[No Readings] Playing around with Multi-armed Bandits Code
+Lost track of readings around this time period :(
+The goal of this independent study is to gain an introduction to the topic of Reinforcement Learning.
+As such the majority of the semester will be following the textbook to gain an introduction to the topic, and the last part applying it to some problems.
+The majority of the content of this independent study will come from the textbook. This is meant to lessen the burden on the both us of as I already experimented with curating my own content.
+The textbook also includes examples throughout the text to immediately apply what's learned.
+Richard S. Sutton and Andrew G. Barto, "Reinforcement Learning: An Introduction" http://incompleteideas.net/book/bookdraft2017nov5.pdf
+Discussions and notes will be kept track of and published on my tilda space as time and energy permits. This is for easy reference and since it's nice to write down what you learn.
+In this section we will get ourselves familiar with the topics that are commonly discussed in Reinforcement learning problems.
+In this section we will learn the different vocab terms such as:
+This is a type of reinforcement learning problem that is commonly studied and well documented. This helps form an environment for which the agent can operate within. Possible subtopics include:
+Dynamic Programming refers to a collection of algorithms that can be used to compute optimal policies given an environment. Subtopics that we are going over is:
+Now we move onto not having complete knowledge of the environment. This will go into estimating value functions and discovering optimal policies. Possible subtopics include:
+Temporal-Difference learning is a combination of Monte Carlo ideas and Dynamic Programming. This can lead to methods learning directly from raw experience without knowledge of an environment. Subtopics will include:
+Reinforcement learning is the art of analyzing situations and mapping them to actions in order to maximize a numerical reward signal.
+In this independent study, I as well as Dr. Stephen Davies, will explore the Reinforcement Learning problem and its subproblems. We will go over the bandit problem, markov decision processes, and discover how best to translate a problem in order to make decisions.
+I have provided a list of topics that I wish to explore in a syllabus
+In order to spend more time learning, I decided to follow a textbook this time.
+Reinforcement Learning: An Introduction
+By Richard S. Sutton and Andrew G. Barto
+ +The notes for this course, is going to be an extreemly summarized version of the textbook. There will also be notes on whatever side tangents Dr. Davies and I explore.
+ +I wrote a small little quirky/funny report describing the bandit problem. Great for learning about the common considerations for Reinforcement Learning problems.
+ +Code will occasionally be written to solidify the learning material and to act as aids for more exploration.
+ +Specifically, if you want to see agents I've created to solve some OpenAI environments, take a look at this specific folder in the Github Repository
+ +Deep Reinforcement Learning: With Dr. Zacharski I focused more on a particular instance of Reinforcement Learning where deep neural networks are used. During this time, I built out a Reinforcement Learning library written in PyTorch. This library helps me have a test bed for trying out different algorithms and attempts to create my own.
+ + +Reinforcement Learning: Currently studying the fundamentals of reinforcement learning with Dr. Davies. We went over the fundamentals such as value functions, policy functions, how we can describe our environment as a markov decision processes, etc.
+ + +Programming Languages: Studying the design of programming languages. So far I made an implementation of the SLOTH programming language, experimenting with what I want my own programming language to be syntatically and paradigm wise.
+ +Before this study, I worked though a book called "Build your own Lisp" and my implementation of a lisp like language is below.
+ +Competitive Programming: Studying algorithms and data structures necessary for competitive programming. Attending ACM ICPC in November with a team of two other students.
+Worked on an independent study on the topic of Cluster Analysis. This is where you try to group similar observations without knowing what the labels are. +I am studying under the guidance of Dr. Melody Denhere, the link below gives you more of a description of the project along with my course notes.
+ +For the two projects below, I worked on Quantum Research in a physics lab with a fellow student Hannah Killian and an advisor Dr. Hai Nguyen. I mostly assisted with the software support for the project and assisted in the mathematics in whatever way I can.
+Modeling Population Dynamics of Incoherent and Coherent Excitation
+Coherent Control of Atomic Population Using the Genetic Algorithm
+In order to circumvent the frustrations I had with simulation code taking a while, I applied and received funding to build out a Beowulf cluster for the Physics department. Dr. Maia Magrakvilidze was the advisor for this project.
+High Performance Cluster for Research and Education Report (nicknamed LUNA-C)
+ +Below is a list of courses I've taken in the past for credit.
+Deep Reinforcement Learning (CPSC 491)
+Real Analysis I (MATH 471)
+Abstract Algebra I (MATH 431)
+Operating Systems (CPSC 405)
+Digital Storytelling (CPSC 106)
+Reinforcement Learning (CPSC 491)
+Programming Languages (CPSC 491)
+Competitive Programming (CPSC 491)
+Multivariate Statistics (STAT 461)
+Foundations of Advance Mathematics (MATH 330)
+Applications of Databases (CPSC 350)
+Foundations for Data Science (Data 219)
+Theory of Computation (CPSC 326)
+Data Structures and Algorithms (CPSC 340)
+Differential Equations (MATH 312)
+Cluster Analysis (STAT 491)
+Computer Systems and Architecture (CPSC 305)
+Introduction to Data Science (DATA 101)
+Linear Algebra (MATH 300)
+Methods in Mathematical Physics (PHYS 317)
+Probability & Statistical Inference I (STAT 381)
+Undergraduate Research in Physics (URES 197J)
+Computer Systems and Architecture (CPSC 305)
+Introduction to Data Science (DATA 101)
+Linear Algebra (MATH 300)
+Methods in Mathematical Physics (PHYS 317)
+Probability & Statistical Inference I (STAT 381)
+Undergraduate Research in Physics (URES 197J)
+Computer Systems and Architecture (CPSC 305)
+Introduction to Data Science (DATA 101)
+Linear Algebra (MATH 300)
+Methods in Mathematical Physics (PHYS 317)
+Probability & Statistical Inference I (STAT 381)
+Undergraduate Research in Physics (URES 197J)
+Intro to Discrete Mathematics (CPSC 125A)
+Object-Oriented Analysis & Design (CPSC 240)
+Statistical Methods (MATH 280)
+University Physics II, with Lab (PHYS 106)
+Undergraduate Research in Physics (URES 197J)
+Software Development Tools (CPSC 225)
+Numbers Rule Your World (FSEM 100M1)
+Intro to Statistics (MATH 200)
+American Music (MUHL 156)
+University Physics I, with Lab (PHYS 105)
+Computer Programming and Problem Solving (CPSC 220)
+Intro to Human Geography (GEOG 102)
+American History to 1865 (HIST 131)
+American History since 1865 (HIST 132)
+Calculus I (MATH 121)
+Calculus II (MATH 122)
+Calculus III (MATH 223)
+General Psychology (PSYCH 100)
+