# Scott Merrill
# Sept 19th Assignment 3 example

# I will discuss this remove list (memory) statement 
#	during R code discussion
# 	In brief - include as part of your template! 
rm(list = ls(all=TRUE))


# What is my population
prob.space = matrix(ncol = 2, nrow = 6, data = 0)
possibilities = c("Red1(Red2)","Red2(Red1)","Black1(Black2)",
			"Black2(Black1)","Red3(Black3)","Black3(Red3)")


population = c(TRUE,TRUE,TRUE,TRUE,FALSE,FALSE)

# making a visual table

headings = c("possible.draws","Other.side.same")
prob.space[,1] = possibilities
prob.space[,2] = population
prob.space2 = list(location = headings,composition = prob.space)
prob.space2


### 10 draws (repeated 100 times)
# Create a vector of to iterate over
n = numeric(100)

# Create a blank vector of to be filled by means
m <- numeric(length(n))

# reminder: length() returns the length of a vector


# The For Loop

# i is an index to iterate over
# the values of i are from the vector 1:length(n)

# sampling from the three card deck 10 times
for (i in 1:length(n)) {
	m[i] <- sum(sample(population, size = 10,replace=TRUE))
}


histogram(m,main = "number of times that the opposite side was the same color with 10 draws",
xlab = c("Probability Space: 
	Red1 with Red2 on flip side, 
	Red2 with Red1 on flip side, 
	etc.:



",possibilities))

###### 1000 draws


# Create a vector of sample sizes you want to iterate over
n = numeric(1000)

# Create a blank vector of means
m <- rep(NA, times = length(n))

# length() returns the length of a vector

# What is my population

possibilities = c("R1(R2)","R2(R1)","B1(B2)","B2(B1)","R3(B3)","B3(R3)")
population = c(TRUE,TRUE,TRUE,TRUE,FALSE,FALSE)

# The For Loop

# i is an index to iterate over
# the values of i are from the vector 1:length(n)

for (i in 1:length(n)) {
	m[i] <- sum(sample(population, size = 1000,replace=TRUE))
}


histogram(m,nint=40,main = "number of times that the opposite side was the same color with 1000 draws",
xlab = c("Probability Space: 
	Red1 with Red2 on flip side, 
	Red2 with Red1 on flip side, 
	etc.:



",possibilities))



### 100000 draws (repeated 1000 times)
# Create a vector of to iterate over
n = numeric(1000)

# Create a blank vector of to be filled by means
m <- numeric(length(n))

# reminder: length() returns the length of a vector


# The For Loop

# i is an index to iterate over
# the values of i are from the vector 1:length(n)

# sampling from the three card deck 10 times
for (i in 1:length(n)) {
	m[i] <- sum(sample(population, size = 100000,replace=TRUE))
}


histogram(m/1000,nint=50,main = "Percentage of time that the opposite side was the same color with 100000 draws",
xlab = c("Probability Space: 
	Red1 with Red2 on flip side, 
	Red2 with Red1 on flip side, 
	etc.:



",possibilities))

##########################################################
# Scott Merrill
# variation and distributions revisited
# For September 19th class

rm(list =ls(all= TRUE))


require(lattice)

# sampling 50 times from a normal distribution with a mean of 50 and a sd of 10
histogram(rnorm(100000,mean = 50, sd = 10),nint = 200,main = "100000 samples from a normal distribution 
with a mean of 50 and standard deviation of 10: 
The assumed distribution of blueberry fruit in bushels per acre")

# realistic number of samples = 50 yield samples / year
yield = rnorm(50, mean = 50, sd = 10)
histogram(yield)


yield2 = rnorm(50, mean = 50, sd = 10)

######## looking at inoculation effect 

histogram(rnorm(100000,mean = 2, sd = 0.7),nint = 200,main = "100000 samples from a normal distribution 
with a mean of 2 and standard deviation of 0.7: 
The assumed distribution of increased yield in plants with 100% rootmass inoculation in bushels per acre")
# sampling 50 times from a normal distribution with a mean of 2 and a sd of .7
inoc = rnorm(50, mean = 2, sd = .7)


#### setting 1 row and 1 column of plots per graphic device display
par(mfrow = c(1,1)) 

#### Plotting simulated data
plot (yield,main = "Fruit Yield Variation Only",xlab = "sample number", col = 2,cex = .5, lwd = 5)
plot (inoc,main = "Inoculation Variation Only", xlab = "sample number", col = 3,cex = .5, lwd = 5)
y.and.inoc = yield2 + inoc

#### Plot of combined data
plot (y.and.inoc,main = "Treatment Plots in Blue, 
Control plots in red", xlab = "sample number", ylab = "yield" ,col = 4,cex = .5, lwd = 5)
points(yield, col = 2,cex = .5, lwd = 5)

# t-test for comparing two populations
t.test(x=yield,y=y.and.inoc)


#### looping the simulation 9 times using a sample number of 100
# setting graphic display to be 3 rows * 3 columns
par(mfrow = c(3,3)) 
	
for (plot.count in 1:9) { 
	
	# sampling 100 times from a normal distribution with a mean of 50 and a sd of 10
	yield = rnorm(100, mean = 50, sd = 10)
	yield2 = rnorm(100, mean = 50, sd = 10)

	# sampling 100 times from a normal distribution with a mean of 2 and a sd of .7
	inoc = rnorm(100, mean = 2, sd = .7)

	y.and.inoc = yield2 + inoc

	plot (y.and.inoc,main = "Treatment plots in Blue, 
	Control plots in red", xlab = "sample number", ylab = "yield" ,col = 4,cex = .5, lwd = 5)
	points(yield, col = 2,cex = .5, lwd = 5)
	
	print(t.test(x=yield,y=y.and.inoc))
}


#### Increasing the number of samples
### Same test but sampling 200 times instead of 100

# sampling 200 times from a normal distribution with a mean of 50 and a sd of 10
yield = rnorm(200, mean = 50, sd = 10)
yield2 = rnorm(200, mean = 50, sd = 10)

# sampling 200 times from a normal distribution with a mean of 2 and a sd of .7
inoc = rnorm(200, mean = 2, sd = .7)

# setting 1 row and 1 column of plots per graphic device display
par(mfrow = c(1,1)) 

plot (yield,main = "Fruit Yield Variation Only",xlab = "sample number", col = 2,cex = .5, lwd = 5)
plot (inoc,main = "Inoculation Variation Only", xlab = "sample number", col = 3,cex = .5, lwd = 5)
y.and.inoc = yield2 + inoc

plot (y.and.inoc,main = "Treatment Plots in Blue, 
Control plots in red", xlab = "sample number", ylab = "yield" ,col = 4,cex = .5, lwd = 5)
points(yield, col = 2,cex = .5, lwd = 5)

# setting 2 row and 1 column of plots per graphic device display
#par(mfrow = c(2,1)) 

#histogram(yield,nint=10)
#histogram(y.and.inoc,nint=10, col = 2)
#mean(yield)
#mean(y.and.inoc)
t.test(x=yield,y=y.and.inoc)


###### Decreasing the variability in the system
# Can you reduce (via management etc.) the variability in your system
# Changed Standard deviation to 5 bushels per acre
# setting graphic display to be 3*3
par(mfrow = c(3,3)) 
	
for (plot.count in 1:9) { 
	
	# sampling 100 times from a normal distribution with a mean of 50 and a sd of 5
	yield = rnorm(100, mean = 50, sd = 5)
	yield2 = rnorm(100, mean = 50, sd = 5)

	# sampling 100 times from a normal distribution with a mean of 2 and a sd of .7
	inoc = rnorm(100, mean = 2, sd = .7)

	y.and.inoc = yield2 + inoc

	plot (y.and.inoc,main = "Treatment plots in Blue, 
	Control plots in red", xlab = "sample number", ylab = "yield" ,col = 4,cex = .5, lwd = 5)
	points(yield, col = 2,cex = .5, lwd = 5)
	
	print(t.test(x=yield,y=y.and.inoc))
}

# histograms of values with means in title
histogram(yield,nint=10,main = c("yield with 0% inoc mean = ",mean(yield)))
histogram(y.and.inoc,nint=10, col = 2, main = c("yield with 100% inoc mean = ",mean(y.and.inoc)))
histogram(inoc,nint=10, col = 3,main = c("yield added with 100% inoc mean = ",mean(inoc)))




##############################################################################################
# Drunken walk
##############
# Original air date December 17, 1989
# 8312 nights
chair = 0
drunk = matrix(nrow = 8312,ncol = 1,data = 0)
for (nights in 1:8312) {
	for (n in 1:100) {
		stumble = sample(c(-1,1),1)
	 	chair = chair + stumble
	}
	drunk[nights,1] = chair
	chair = 0
}

histogram(drunk[,1],nint = 40)

########################################################################################

#sampling from other probability distribution functions

# Poisson
histogram(rpois(100000, .5),nint=30)
histogram(rpois(100000, 1),nint=30)
histogram(rpois(100000, 2),nint=30)
histogram(rpois(100000, 4),nint=30)
histogram(rpois(100000, 10),nint=30)

# Binomial
histogram(rbinom(100000,10, .5),nint=10)
histogram(rbinom(100000,10,.7),nint=10)
histogram(rbinom(100000,10,.9),nint=7)

# gamma
# rgamma
histogram(rgamma(100000,10, .5),nint=10)
histogram(rbinom(100000,10,.7),nint=10)
histogram(rbinom(100000,10,.9),nint=7)


# Talk about the list function
ls()


# give an example of the sample function
sample(1:6,10,replace = TRUE)