Assignment 6

Question 1

lpYear <- function(year){
  if(year %% 4 == 0 && year %% 100 !=0 || year %% 400 == 0){
    print("This is  a leap year!")}
  else 
  {print("This is NOT a leap year!")}
}

I was born in 1996. Let’s check if that year was a leap year!

lpYear(1996)

## [1] "This is  a leap year!"

lifespan<- function(year,yob){
  count <-0
  while(year > yob){
    year <- year-1
    if(year %% 4 == 0 && year %% 100 !=0 || year %% 400 == 0){
      count <-count+1
    } else count <- count
  }
  print(paste("Congratulations!You have lived through:",count, "leap years"))
}

I was born in 1996 and the current year is 2016. Let’s check how many leap years I’ve lived through!

lifespan(2016,1996)

## [1] "Congratulations!You have lived through: 5 leap years"

Question 2

Part A

q <- seq(-1, 1, 0.02)
## I need to remove zero from the list, because we don't want to use it! 
q <- q[ q != 0 ]


x <-sample(q,1)
y <-sample(q,1)
count <- 1

## let c is the  area of circle
## the analogy is if x^2 + y^2 <= 1

count <- 0
c <- 0
q <- seq(-1, 1, 0.01)
q <- q[ q != 0 ]

while(count < 100){ #this is to say we want to replicate 10,000 times
  if (sample(q,1)^2 + sample(q,1)^2 <= 1) {
    c <- c +1
    count <- count +1}
    else {count <- count +1
    }
}

## this is a proportion of points that falls within a unit circle:

p <- (c/count)

## this is a plot of points, I will name is gg
## creating a new data frame -- which contains two variables which is a data frame
## for ggplot
gg <- data.frame(sample(q,100, replace = TRUE))
gg$"gg2" <- (sample(q,100, replace = TRUE))
gg$"Area" <- 0
colnames(gg) <- c("X", "Y","Area")

## "a" is for a point that falls within a circle 
## "b" is for a point that does NOT falls within a circle 
oo <- 0

for (i in 1:nrow(gg)) {
  if(gg[i,1]^2 + gg[i,2]^2 <= 1){
    gg[i,3] <- "Inside Circle"
  } else (gg[i,3] <- "Outside Circle")
}

tt <- ggplot(gg, aes( x = X, y = Y))+ geom_point(size = 3,aes(col = Area)) 

tt + theme_minimal()

Part B

appc <-function(n, plotit = TRUE){
  q <- seq(-1, 1, 0.01)
  ## I need to remove zero from the list, because we don't want to use it! 
  q <- q[ q != 0 ]
  
  gg <- data.frame(sample(q,n, replace = TRUE)) #x
  gg$"gg2" <- (sample(q,n, replace = TRUE)) #y
  gg$"Area" <- 0
  colnames(gg) <- c("X", "Y","Area")
  
  ## "a" is for a point that falls within a circle 
  ## "b" is for a point that does NOT falls within a circle 
  oo <- 0
  prop<-0
  for (i in 1:nrow(gg)) {
    if(gg[i,1]^2 + gg[i,2]^2 <= 1){
      gg[i,3] <- "Inside Circle"
      prop <- prop + 1
    } else (gg[i,3] <- "Outside Circle")
  }
  # gg <-data.frame(gg)
  if (plotit == TRUE) {
    tt <- ggplot(gg, aes( x = X, y = Y))+ 
      geom_point(size = 3,aes(col = Area)) 
    tt+ theme_minimal()
    plot(tt)
  }
  area <- (prop/n)
  return(area)
}

Let’s test the function with n = 200 and plotit = TRUE!

appc(200, TRUE)

## [1] 0.755

Part C Make an R expression and evaluate it a specified number of times. You can think of it as a very simplified for loop in which exactly the same code gets evaluated each time through the loop. Use this function to approximate the area of the circle 100 times, first for n=50 and then for n=500 and storing the results in two different vectors.

a50<-replicate(100,appc(50,FALSE))

a500<-replicate( 100, appc(500, FALSE))

Part D Make two histograms showing the results from part (c). Add a vertical line to the plot to indicate the true area. Give your plot appropriate axis names and titles and a nice readable theme.

Sample_Size_50 <-a50 *4
Sample_Size_100<- a500 *4
com <-data.frame(Sample_Size_50,Sample_Size_100)
com2 <- gather(com,key = size, value = area, Sample_Size_50, Sample_Size_100)

d_d <- ggplot(com2, aes(x =area)) + geom_bar(size = 3.5, fill = c("brown"))

d_d+ theme_minimal()+facet_grid(.~size) + ggtitle("Distribution of Area")+
  geom_vline(aes(xintercept = 3.141592))

Question 3

Part A Write a function that implements Newton’s method for this particular choice of function f. It should have arguments for the initial guess and the tolerance for convergence (currently set to 0.00001 in my code). Pay special attention to whether either of these should have default values. The function should return a single number and not plot anything.

## x is a inital guess we needed to make
## Define your fz and fpz
## Such as using this form fz <- function(x) {x^2+20} 
## and fpz <- function(x) {2*x}

new <- function(a, fz, fpz){
  
  while(abs(fz(a)) > 0.00001){
    a <- a - (fz(a)/fpz(a))
  }
  return(a) ## this value a will make your fp(a) = 0
}

Part B

## This is for fx
fz <- function(x){
  x^3 + 2*x^2 - 7
}

## This is for fpx
fpz <- function(x){
  3*(x^2) + 4*x
}

Part C

 #Rewrite your function from (a) so that it also takes arguments for the function
# whose root you want to find and its first derivative. You may find it helpful to
# first use your functions from (b) and make sure the function still works, then
# modify the arguments to the function.

root_for_fz <- new(0.1,fz,fpz)

root_for_fz

## [1] 1.428818

Part D

Choose a different function f and write R functions for it and the corresponding derivative. Make a plot as I did in Newton.R to verify that the algorithm found a root of your function. (Note: There are some cases in which Newton’s method will fail to converge, which you will experience as the while loop never ending. If that happens to you, just choose another initial guess and/or function.)

This is the same function with a different input function

fq <-function(x) {x^3 +7*x^7 -20}
fpq <- function(x){3*x^2 + 49*x^6}

new2 <- function(a,fz, fpz){
  while(abs(fz(a)) > 0.00001){
    a <- a - (fz(a)/fpz(a))
    value <-fz(a)
    df=data.frame(x =a,y =value)
  }
  return(df)
}

new2(20,fq,fpq)

##          x           y
## 1 1.148794 6.63823e-07

## HE is going to be a set of data that I made a guess for fq
## I guess the zero is betwee 0 to 1.5
he <-seq(0, 1.5, 0.01)

## data set I want for fq
he2 <-data.frame(he,fq(he))

gg <- ggplot(he2, aes(x=he, y = fq.he.)) +
  geom_point() +
  ylim(-40,40) +
  theme_minimal()

gg + geom_vline(xintercept = 1.148794) + geom_hline(yintercept = 0)

## Warning: Removed 15 rows containing missing values (geom_point).

Assignment 6

Daniel Alonzo

March 2, 2016

Question 1

Question 2

Question 3