Homework 6
Hannah Shafer
2/23/2022
- Running the sample code
Open libraries
library(ggplot2) # for graphics
library(MASS) # for maximum likelihood estimationRead in data vector and generate fake data
# quick and dirty, a truncated normal distribution to work on the solution set
z <- rnorm(n=3000,mean=0.2)
z <- data.frame(1:3000,z)
names(z) <- list("ID","myVar")
z <- z[z$myVar>0,]
str(z)## 'data.frame': 1737 obs. of 2 variables:
## $ ID : int 1 3 4 5 7 9 12 13 19 20 ...
## $ myVar: num 0.828 0.726 0.65 0.176 0.219 ...# data.frame': 1725 obs. of 2 variables:
# $ ID : int 8 10 11 13 16 17 19 22 24 29 ...
# $ myVar: num 0.276 0.878 1.04 0.431 2.484
summary(z$myVar)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000781 0.337917 0.726945 0.848326 1.215819 3.754458# Min. 1st Qu. Median Mean 3rd Qu.
# 0.000904 0.348635 0.760045 0.877109 1.270632
# Max.
# 4.085390Plot a histogram of the data
p1 <- ggplot(data=z, aes(x=myVar, y=..density..)) +
geom_histogram(color="grey60",fill="cornsilk",size=0.2, bins=200)
print(p1) # prints out histogram
Add an empirical density curve
p1 <- p1 + geom_density(linetype="dotted",size=0.75)
print(p1) # prints histogram from before but with a dotted density curve
Get maximum likelihood parameters for normal
normPars <- fitdistr(z$myVar,"normal") # fit normal distribution to the data
print(normPars)## mean sd
## 0.84832570 0.63270331
## (0.01518099) (0.01073458) # mean sd
# 0.87710945 0.63904730
# (0.01538645) (0.01087986)
str(normPars)## List of 5
## $ estimate: Named num [1:2] 0.848 0.633
## ..- attr(*, "names")= chr [1:2] "mean" "sd"
## $ sd : Named num [1:2] 0.0152 0.0107
## ..- attr(*, "names")= chr [1:2] "mean" "sd"
## $ vcov : num [1:2, 1:2] 0.00023 0 0 0.000115
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "mean" "sd"
## .. ..$ : chr [1:2] "mean" "sd"
## $ n : int 1737
## $ loglik : num -1670
## - attr(*, "class")= chr "fitdistr"# List of 5
# $ estimate: Named num [1:2] 0.877 0.639
# ..- attr(*, "names")= chr [1:2] "mean" "sd"
# $ sd : Named num [1:2] 0.0154 0.0109
# ..- attr(*, "names")= chr [1:2] "mean" "sd"
# $ vcov : num [1:2, 1:2] 0.000237 0 0 0.000118
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : chr [1:2] "mean" "sd"
# .. ..$ : chr [1:2] "mean" "sd"
# $ n : int 1725
# $ loglik : num -1675
# - attr(*, "class")= chr "fitdistr"
normPars$estimate["mean"] # note structure of getting a named attribute## mean
## 0.8483257# mean
# 0.8771095 Plot the normal probability density; involves calling the dnorm function inside ggplot’s stat_function stat_function in help system- “draws a function as a continuous curve” - this would allow us to add a smooth function to any ggplot
meanML <- normPars$estimate["mean"]
sdML <- normPars$estimate["sd"]
xval <- seq(0,max(z$myVar),len=length(z$myVar))
stat <- stat_function(aes(x = xval, y = ..y..), fun = dnorm, colour="red", n = length(z$myVar), args = list(mean = meanML, sd = sdML))
p1 + stat # prints the histogram with empirical density curve from before, but adds the normal density curve (note biased mean)
Plot exponential probability density
expoPars <- fitdistr(z$myVar,"exponential")
rateML <- expoPars$estimate["rate"]
stat2 <- stat_function(aes(x = xval, y = ..y..), fun = dexp, colour="blue", n = length(z$myVar), args = list(rate=rateML))
p1 + stat + stat2 # exponential curve added to plot
Plot uniform probability density
stat3 <- stat_function(aes(x = xval, y = ..y..), fun = dunif, colour="darkgreen", n = length(z$myVar), args = list(min=min(z$myVar), max=max(z$myVar)))
p1 + stat + stat2 + stat3
Plot gamma probability density
gammaPars <- fitdistr(z$myVar,"gamma")## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs producedshapeML <- gammaPars$estimate["shape"]
rateML <- gammaPars$estimate["rate"]
stat4 <- stat_function(aes(x = xval, y = ..y..), fun = dgamma, colour="brown", n = length(z$myVar), args = list(shape=shapeML, rate=rateML))
p1 + stat + stat2 + stat3 + stat4
Plot beta probability density
pSpecial <- ggplot(data=z, aes(x=myVar/(max(myVar + 0.1)), y=..density..)) +
geom_histogram(color="grey60",fill="cornsilk",size=0.2) +
xlim(c(0,1)) +
geom_density(size=0.75,linetype="dotted")
betaPars <- fitdistr(x=z$myVar/max(z$myVar + 0.1),start=list(shape1=1,shape2=2),"beta")## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs producedshape1ML <- betaPars$estimate["shape1"]
shape2ML <- betaPars$estimate["shape2"]
statSpecial <- stat_function(aes(x = xval, y = ..y..), fun = dbeta, colour="orchid", n = length(z$myVar), args = list(shape1=shape1ML,shape2=shape2ML))
pSpecial + statSpecial## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 2 rows containing missing values (geom_bar).
- Import data and run the analysis above on it; data imported and truncated from: https://datadryad.org/stash/dataset/doi:10.5061%2Fdryad.vdncjsxrx
z <- read.table("../Bio381/MyDataFile.csv",header=TRUE,sep=",")
names(z) <- list("ID","pop", "myVar")
z <- z[z$myVar>0,]
str(z)## 'data.frame': 1550 obs. of 3 variables:
## $ ID : chr "ALB-M011A01" "ALB-M011A02" "ALB-M011A03" "ALB-M011A04" ...
## $ pop : chr "ALB" "ALB" "ALB" "ALB" ...
## $ myVar: int 101 102 108 113 100 100 102 102 114 100 ...summary(z$myVar)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 93 100 102 106 114 127str(z)## 'data.frame': 1550 obs. of 3 variables:
## $ ID : chr "ALB-M011A01" "ALB-M011A02" "ALB-M011A03" "ALB-M011A04" ...
## $ pop : chr "ALB" "ALB" "ALB" "ALB" ...
## $ myVar: int 101 102 108 113 100 100 102 102 114 100 ...summary(z)## ID pop myVar
## Length:1550 Length:1550 Min. : 93
## Class :character Class :character 1st Qu.:100
## Mode :character Mode :character Median :102
## Mean :106
## 3rd Qu.:114
## Max. :127- Run the above code on the imported data rather than the “fake data”
p1 <- ggplot(data=z, aes(x=myVar, y=..density..)) +
geom_histogram(color="grey60",fill="cornsilk",size=0.2)
print(p1) # plot a histogram of the data## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
p1 <- p1 + geom_density(linetype="dotted",size=0.75)
print(p1) # add empirical density curve; black dotted line## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
normPars <- fitdistr(z$myVar,"normal")
print(normPars)## mean sd
## 105.9516129 6.8019623
## ( 0.1727700) ( 0.1221668)str(normPars)## List of 5
## $ estimate: Named num [1:2] 106 6.8
## ..- attr(*, "names")= chr [1:2] "mean" "sd"
## $ sd : Named num [1:2] 0.173 0.122
## ..- attr(*, "names")= chr [1:2] "mean" "sd"
## $ vcov : num [1:2, 1:2] 0.0298 0 0 0.0149
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "mean" "sd"
## .. ..$ : chr [1:2] "mean" "sd"
## $ n : int 1550
## $ loglik : num -5171
## - attr(*, "class")= chr "fitdistr"normPars$estimate["mean"] # get maximum likelihood parameters for normal density## mean
## 105.9516meanML <- normPars$estimate["mean"]
sdML <- normPars$estimate["sd"]
xval <- seq(0,max(z$myVar),len=length(z$myVar))
stat <- stat_function(aes(x = xval, y = ..y..), fun = dnorm, colour="red", n = length(z$myVar), args = list(mean = meanML, sd = sdML))
p1 + stat # plot normal probability density; red line## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
expoPars <- fitdistr(z$myVar,"exponential")
rateML <- expoPars$estimate["rate"]
stat2 <- stat_function(aes(x = xval, y = ..y..), fun = dexp, colour="blue", n = length(z$myVar), args = list(rate=rateML))
p1 + stat + stat2 # plot exponential probability density; blue line## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
stat3 <- stat_function(aes(x = xval, y = ..y..), fun = dunif, colour="darkgreen", n = length(z$myVar), args = list(min=min(z$myVar), max=max(z$myVar)))
p1 + stat + stat2 + stat3 # plot uniform probability density; green line## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
gammaPars <- fitdistr(z$myVar,"gamma")
shapeML <- gammaPars$estimate["shape"]
rateML <- gammaPars$estimate["rate"]
stat4 <- stat_function(aes(x = xval, y = ..y..), fun = dgamma, colour="brown", n = length(z$myVar), args = list(shape=shapeML, rate=rateML))
p1 + stat + stat2 + stat3 + stat4 # plot gamma probability density; brown line## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
pSpecial <- ggplot(data=z, aes(x=myVar/(max(myVar + 0.1)), y=..density..)) +
geom_histogram(color="grey60",fill="cornsilk",size=0.2) +
xlim(c(0,1)) +
geom_density(size=0.75,linetype="dotted")
betaPars <- fitdistr(x=z$myVar/max(z$myVar + 0.1),start=list(shape1=1,shape2=2),"beta")## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs producedshape1ML <- betaPars$estimate["shape1"]
shape2ML <- betaPars$estimate["shape2"]
statSpecial <- stat_function(aes(x = xval, y = ..y..), fun = dbeta, colour="orchid", n = length(z$myVar), args = list(shape1=shape1ML,shape2=shape2ML))
pSpecial + statSpecial # plot beta probability density## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 2 rows containing missing values (geom_bar).
The best fitting probability density curve for this set of data is either the normal or gamma - the lines are nearly overlapping, so I think they fit equally well. These two curves matched the high peak of the data. However, neither line is a true match because the data appears to be bimodal, and neither the normal nor gamma curves have more than one peak. This dataset had multiple columns of data, so the above analysis was performed on just the first column. Below I tried another column to see if it fit a normal, uniform, exponential, or gamma distribution rather than bimodal.
z <- read.table("../Bio381/MyDataFile2.csv",header=TRUE,sep=",")
names(z) <- list("ID","pop", "myVar")
z <- z[z$myVar>0,]
str(z)## 'data.frame': 1550 obs. of 3 variables:
## $ ID : chr "ALB-M011A01" "ALB-M011A02" "ALB-M011A03" "ALB-M011A04" ...
## $ pop : chr "ALB" "ALB" "ALB" "ALB" ...
## $ myVar: int 113 114 116 117 100 127 115 115 115 100 ...summary(z$myVar)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 93.0 114.0 116.0 114.2 117.0 135.0str(z)## 'data.frame': 1550 obs. of 3 variables:
## $ ID : chr "ALB-M011A01" "ALB-M011A02" "ALB-M011A03" "ALB-M011A04" ...
## $ pop : chr "ALB" "ALB" "ALB" "ALB" ...
## $ myVar: int 113 114 116 117 100 127 115 115 115 100 ...summary(z)## ID pop myVar
## Length:1550 Length:1550 Min. : 93.0
## Class :character Class :character 1st Qu.:114.0
## Mode :character Mode :character Median :116.0
## Mean :114.2
## 3rd Qu.:117.0
## Max. :135.0p1 <- ggplot(data=z, aes(x=myVar, y=..density..)) +
geom_histogram(color="grey60",fill="cornsilk",size=0.2)
print(p1) # plot a histogram of the data## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
p1 <- p1 + geom_density(linetype="dotted",size=0.75)
print(p1) # add empirical density curve; black dotted line## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
normPars <- fitdistr(z$myVar,"normal")
print(normPars)## mean sd
## 114.1658065 6.3604822
## ( 0.1615564) ( 0.1142376)str(normPars)## List of 5
## $ estimate: Named num [1:2] 114.17 6.36
## ..- attr(*, "names")= chr [1:2] "mean" "sd"
## $ sd : Named num [1:2] 0.162 0.114
## ..- attr(*, "names")= chr [1:2] "mean" "sd"
## $ vcov : num [1:2, 1:2] 0.0261 0 0 0.0131
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "mean" "sd"
## .. ..$ : chr [1:2] "mean" "sd"
## $ n : int 1550
## $ loglik : num -5067
## - attr(*, "class")= chr "fitdistr"normPars$estimate["mean"] # get maximum likelihood parameters for normal density## mean
## 114.1658meanML <- normPars$estimate["mean"]
sdML <- normPars$estimate["sd"]
xval <- seq(0,max(z$myVar),len=length(z$myVar))
stat <- stat_function(aes(x = xval, y = ..y..), fun = dnorm, colour="red", n = length(z$myVar), args = list(mean = meanML, sd = sdML))
p1 + stat # plot normal probability density; red line## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
expoPars <- fitdistr(z$myVar,"exponential")
rateML <- expoPars$estimate["rate"]
stat2 <- stat_function(aes(x = xval, y = ..y..), fun = dexp, colour="blue", n = length(z$myVar), args = list(rate=rateML))
p1 + stat + stat2 # plot exponential probability density; blue line## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
stat3 <- stat_function(aes(x = xval, y = ..y..), fun = dunif, colour="darkgreen", n = length(z$myVar), args = list(min=min(z$myVar), max=max(z$myVar)))
p1 + stat + stat2 + stat3 # plot uniform probability density; green line## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
gammaPars <- fitdistr(z$myVar,"gamma")
shapeML <- gammaPars$estimate["shape"]
rateML <- gammaPars$estimate["rate"]
stat4 <- stat_function(aes(x = xval, y = ..y..), fun = dgamma, colour="brown", n = length(z$myVar), args = list(shape=shapeML, rate=rateML))
p1 + stat + stat2 + stat3 + stat4 # plot gamma probability density; brown line## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
The second column of data was less bimodal, and fit either a normal or gamma curve - the two curves were again nearly identical.
- Simulating a new data set using the maximum likelihood parameters from the above gamma distribution, and comparing histograms
Maximum likelihood estimators of the parameters:
data <- read.table("../Bio381/MyDataFile2.csv",header=TRUE,sep=",")
names(data) <- list("ID","pop", "myVar")
data <- z[z$myVar>0,]
print(data$myVar)## [1] 113 114 116 117 100 127 115 115 115 100 118 116 117 114 117 117 116 108
## [19] 116 120 118 116 115 118 115 115 115 118 116 116 119 115 115 102 116 114
## [37] 115 116 118 115 118 118 117 114 116 115 116 117 115 116 115 134 119 113
## [55] 118 118 118 120 117 117 119 100 115 116 116 117 118 117 119 115 119 117
## [73] 117 113 116 116 117 100 119 119 114 119 127 118 117 102 118 118 119 119
## [91] 100 115 114 118 117 115 113 117 111 115 116 116 102 115 100 115 116 116
## [109] 117 117 127 114 102 104 103 115 100 113 118 118 108 111 100 118 117 119
## [127] 135 116 114 111 117 118 119 114 115 117 111 114 117 102 103 103 103 100
## [145] 116 116 113 118 100 118 127 116 117 118 116 117 117 116 102 102 118 100
## [163] 127 116 118 115 127 117 127 116 114 118 112 105 115 117 115 117 113 119
## [181] 115 116 109 103 118 118 100 101 103 118 116 101 115 117 103 114 118 118
## [199] 115 113 127 120 118 113 116 100 117 102 117 102 116 116 100 100 103 117
## [217] 115 100 118 127 100 118 116 117 109 113 100 117 117 116 117 117 104 115
## [235] 115 102 116 115 113 119 119 119 100 100 104 101 102 116 109 126 102 115
## [253] 115 114 104 114 119 118 117 115 102 117 115 114 117 116 114 116 115 101
## [271] 118 116 118 102 115 127 110 117 118 108 116 116 114 118 116 117 111 116
## [289] 105 118 119 116 116 116 114 113 118 112 115 105 116 119 108 114 115 100
## [307] 114 103 119 102 115 100 120 117 113 118 117 116 116 115 115 117 109 119
## [325] 116 100 117 127 116 117 115 110 115 111 114 115 117 116 117 115 117 117
## [343] 111 117 111 117 119 126 116 102 117 117 115 117 116 115 116 115 116 115
## [361] 117 119 116 118 118 109 116 119 109 100 115 116 115 111 116 115 116 118
## [379] 117 116 116 116 115 114 117 115 117 114 114 115 113 118 119 115 115 115
## [397] 116 116 115 127 129 115 117 100 120 120 116 102 115 102 115 107 117 118
## [415] 118 102 117 117 116 100 119 102 116 117 100 114 100 115 117 102 102 113
## [433] 117 116 100 102 102 118 115 118 115 114 112 115 115 111 118 115 116 118
## [451] 116 114 117 116 113 114 100 113 116 117 114 114 119 117 102 115 115 119
## [469] 118 117 113 118 116 116 117 117 117 111 115 102 113 111 135 115 102 116
## [487] 117 119 117 115 100 117 111 100 115 114 113 109 117 103 113 116 113 102
## [505] 116 117 117 118 118 115 100 116 102 116 115 118 115 118 117 115 115 135
## [523] 117 116 118 116 118 116 116 116 135 101 115 115 118 101 115 118 115 115
## [541] 116 115 114 115 116 115 109 116 116 115 117 135 115 115 116 115 127 104
## [559] 100 118 117 116 116 115 118 116 120 100 109 114 115 102 115 116 135 116
## [577] 135 114 115 120 117 115 115 116 117 115 114 108 117 116 100 118 103 119
## [595] 117 127 117 116 114 118 116 117 120 114 101 100 108 114 114 113 117 118
## [613] 115 101 113 116 116 119 115 117 113 112 113 115 119 117 103 103 117 117
## [631] 116 109 118 116 115 116 118 126 115 127 118 127 109 115 117 114 117 116
## [649] 116 114 115 118 127 114 116 116 118 116 100 118 113 111 118 111 118 102
## [667] 118 115 117 127 116 115 119 116 117 102 118 113 113 102 118 102 117 118
## [685] 117 115 117 102 114 114 111 114 127 118 119 116 110 114 118 118 118 116
## [703] 119 100 114 119 114 116 118 117 115 119 118 116 113 118 116 117 116 113
## [721] 117 117 108 108 113 115 117 118 127 116 117 100 115 117 102 127 115 117
## [739] 115 113 116 116 115 116 117 115 115 117 115 102 116 103 100 100 117 109
## [757] 103 102 118 116 127 127 116 115 123 113 118 100 101 114 118 117 118 117
## [775] 118 117 116 116 114 117 117 118 117 109 102 102 117 116 101 103 115 116
## [793] 115 116 117 115 127 115 100 116 116 118 117 119 103 117 114 115 115 109
## [811] 117 126 114 115 118 115 118 116 118 118 115 117 111 115 117 117 102 114
## [829] 120 100 116 100 114 114 118 118 102 117 108 110 103 114 117 117 111 114
## [847] 100 102 116 118 117 100 117 102 100 116 127 115 117 117 114 117 118 100
## [865] 118 117 102 116 113 100 100 117 116 114 111 113 116 117 115 119 102 114
## [883] 117 101 116 118 118 117 118 127 118 117 113 101 115 127 116 118 116 119
## [901] 117 100 119 117 116 116 116 116 102 115 115 114 117 114 117 100 117 114
## [919] 102 100 115 115 119 117 118 116 115 115 108 135 114 119 101 101 116 117
## [937] 108 127 117 111 115 115 117 109 110 118 115 110 103 116 115 114 116 116
## [955] 119 113 117 119 114 115 116 116 115 114 100 119 114 114 115 119 103 100
## [973] 100 116 115 118 116 116 110 100 117 118 114 103 115 116 115 114 103 110
## [991] 118 113 103 118 118 115 117 115 118 128 115 116 116 100 100 118 116 118
## [1009] 114 114 117 118 112 116 118 108 117 101 103 126 116 115 118 116 117 102
## [1027] 113 112 102 117 117 116 118 100 118 117 108 116 101 118 114 117 118 116
## [1045] 112 117 117 118 118 118 115 117 114 117 111 114 114 117 117 116 117 119
## [1063] 118 102 102 118 113 127 117 118 115 118 117 117 102 111 116 118 118 114
## [1081] 117 115 117 117 119 127 117 117 118 102 113 118 116 127 118 116 126 100
## [1099] 100 127 118 117 115 115 116 114 119 111 118 118 126 113 117 119 116 116
## [1117] 115 114 102 115 117 115 119 117 116 102 115 117 116 118 100 117 115 118
## [1135] 116 117 117 116 116 118 115 100 117 116 100 116 102 116 115 119 118 127
## [1153] 102 102 115 114 117 126 115 100 115 115 118 119 117 102 118 126 115 114
## [1171] 101 116 116 116 111 113 118 116 116 116 102 117 114 115 117 113 116 100
## [1189] 116 117 117 116 113 114 100 115 116 115 115 115 114 115 115 115 114 100
## [1207] 116 116 115 117 108 113 117 114 115 116 113 118 116 116 115 115 100 115
## [1225] 108 114 117 117 114 111 117 117 117 100 117 116 114 117 115 119 109 104
## [1243] 109 109 117 119 115 111 109 111 119 100 115 114 114 118 115 109 119 115
## [1261] 116 102 102 116 115 116 127 117 117 117 117 117 117 116 100 117 117 117
## [1279] 102 115 117 118 115 100 118 117 114 117 116 118 118 118 117 126 116 117
## [1297] 117 117 118 114 116 116 103 116 117 119 116 115 117 116 102 100 117 116
## [1315] 104 115 118 116 116 102 119 114 118 115 115 118 135 114 119 117 117 113
## [1333] 100 113 104 102 117 100 111 117 114 117 100 114 100 100 116 119 118 116
## [1351] 114 119 115 118 115 111 93 115 116 116 120 110 118 118 114 116 119 117
## [1369] 119 115 119 119 118 115 117 115 116 116 116 117 116 113 127 100 113 115
## [1387] 102 118 103 117 118 115 114 117 116 117 110 114 117 103 119 115 117 115
## [1405] 114 102 117 115 116 116 118 117 127 118 117 118 114 100 115 110 115 135
## [1423] 115 135 117 117 115 115 111 116 116 115 118 115 118 117 135 118 113 115
## [1441] 117 117 127 115 116 100 135 116 115 111 117 118 102 116 109 117 117 113
## [1459] 120 115 100 127 118 127 117 115 113 115 117 118 117 115 102 113 126 117
## [1477] 117 119 103 115 116 100 115 116 115 117 115 117 113 100 101 102 114 115
## [1495] 102 101 116 115 100 116 102 117 117 115 116 116 115 117 115 115 104 115
## [1513] 117 117 117 117 115 103 117 114 105 114 115 115 103 116 118 103 115 100
## [1531] 115 116 117 118 100 118 117 100 120 117 101 120 117 101 115 117 116 115
## [1549] 109 116data <- fitdistr(data$myVar, "gamma")
print(data)## shape rate
## 321.967275 2.820173
## ( 11.556479) ( 0.101304) # shape rate
# 321.967275 2.820173
# ( 11.556479) ( 0.101304)Next, use those ML parameters and length of my imported dataset to simulate new data:
sim_data <- rgamma(n=1550, shape=data$estimate["shape"], rate=data$estimate["rate"])
summary(sim_data)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 93.94 109.59 113.96 114.15 118.29 134.62Finally, compare the resulting histograms & probability curves for the respective data sets.
#Simulated data:
data <- read.table("../Bio381/MyDataFile2.csv",header=TRUE,sep=",")
names(data) <- list("ID","pop", "myVar")
data <- z[z$myVar>0,]
print(data$myVar)## [1] 113 114 116 117 100 127 115 115 115 100 118 116 117 114 117 117 116 108
## [19] 116 120 118 116 115 118 115 115 115 118 116 116 119 115 115 102 116 114
## [37] 115 116 118 115 118 118 117 114 116 115 116 117 115 116 115 134 119 113
## [55] 118 118 118 120 117 117 119 100 115 116 116 117 118 117 119 115 119 117
## [73] 117 113 116 116 117 100 119 119 114 119 127 118 117 102 118 118 119 119
## [91] 100 115 114 118 117 115 113 117 111 115 116 116 102 115 100 115 116 116
## [109] 117 117 127 114 102 104 103 115 100 113 118 118 108 111 100 118 117 119
## [127] 135 116 114 111 117 118 119 114 115 117 111 114 117 102 103 103 103 100
## [145] 116 116 113 118 100 118 127 116 117 118 116 117 117 116 102 102 118 100
## [163] 127 116 118 115 127 117 127 116 114 118 112 105 115 117 115 117 113 119
## [181] 115 116 109 103 118 118 100 101 103 118 116 101 115 117 103 114 118 118
## [199] 115 113 127 120 118 113 116 100 117 102 117 102 116 116 100 100 103 117
## [217] 115 100 118 127 100 118 116 117 109 113 100 117 117 116 117 117 104 115
## [235] 115 102 116 115 113 119 119 119 100 100 104 101 102 116 109 126 102 115
## [253] 115 114 104 114 119 118 117 115 102 117 115 114 117 116 114 116 115 101
## [271] 118 116 118 102 115 127 110 117 118 108 116 116 114 118 116 117 111 116
## [289] 105 118 119 116 116 116 114 113 118 112 115 105 116 119 108 114 115 100
## [307] 114 103 119 102 115 100 120 117 113 118 117 116 116 115 115 117 109 119
## [325] 116 100 117 127 116 117 115 110 115 111 114 115 117 116 117 115 117 117
## [343] 111 117 111 117 119 126 116 102 117 117 115 117 116 115 116 115 116 115
## [361] 117 119 116 118 118 109 116 119 109 100 115 116 115 111 116 115 116 118
## [379] 117 116 116 116 115 114 117 115 117 114 114 115 113 118 119 115 115 115
## [397] 116 116 115 127 129 115 117 100 120 120 116 102 115 102 115 107 117 118
## [415] 118 102 117 117 116 100 119 102 116 117 100 114 100 115 117 102 102 113
## [433] 117 116 100 102 102 118 115 118 115 114 112 115 115 111 118 115 116 118
## [451] 116 114 117 116 113 114 100 113 116 117 114 114 119 117 102 115 115 119
## [469] 118 117 113 118 116 116 117 117 117 111 115 102 113 111 135 115 102 116
## [487] 117 119 117 115 100 117 111 100 115 114 113 109 117 103 113 116 113 102
## [505] 116 117 117 118 118 115 100 116 102 116 115 118 115 118 117 115 115 135
## [523] 117 116 118 116 118 116 116 116 135 101 115 115 118 101 115 118 115 115
## [541] 116 115 114 115 116 115 109 116 116 115 117 135 115 115 116 115 127 104
## [559] 100 118 117 116 116 115 118 116 120 100 109 114 115 102 115 116 135 116
## [577] 135 114 115 120 117 115 115 116 117 115 114 108 117 116 100 118 103 119
## [595] 117 127 117 116 114 118 116 117 120 114 101 100 108 114 114 113 117 118
## [613] 115 101 113 116 116 119 115 117 113 112 113 115 119 117 103 103 117 117
## [631] 116 109 118 116 115 116 118 126 115 127 118 127 109 115 117 114 117 116
## [649] 116 114 115 118 127 114 116 116 118 116 100 118 113 111 118 111 118 102
## [667] 118 115 117 127 116 115 119 116 117 102 118 113 113 102 118 102 117 118
## [685] 117 115 117 102 114 114 111 114 127 118 119 116 110 114 118 118 118 116
## [703] 119 100 114 119 114 116 118 117 115 119 118 116 113 118 116 117 116 113
## [721] 117 117 108 108 113 115 117 118 127 116 117 100 115 117 102 127 115 117
## [739] 115 113 116 116 115 116 117 115 115 117 115 102 116 103 100 100 117 109
## [757] 103 102 118 116 127 127 116 115 123 113 118 100 101 114 118 117 118 117
## [775] 118 117 116 116 114 117 117 118 117 109 102 102 117 116 101 103 115 116
## [793] 115 116 117 115 127 115 100 116 116 118 117 119 103 117 114 115 115 109
## [811] 117 126 114 115 118 115 118 116 118 118 115 117 111 115 117 117 102 114
## [829] 120 100 116 100 114 114 118 118 102 117 108 110 103 114 117 117 111 114
## [847] 100 102 116 118 117 100 117 102 100 116 127 115 117 117 114 117 118 100
## [865] 118 117 102 116 113 100 100 117 116 114 111 113 116 117 115 119 102 114
## [883] 117 101 116 118 118 117 118 127 118 117 113 101 115 127 116 118 116 119
## [901] 117 100 119 117 116 116 116 116 102 115 115 114 117 114 117 100 117 114
## [919] 102 100 115 115 119 117 118 116 115 115 108 135 114 119 101 101 116 117
## [937] 108 127 117 111 115 115 117 109 110 118 115 110 103 116 115 114 116 116
## [955] 119 113 117 119 114 115 116 116 115 114 100 119 114 114 115 119 103 100
## [973] 100 116 115 118 116 116 110 100 117 118 114 103 115 116 115 114 103 110
## [991] 118 113 103 118 118 115 117 115 118 128 115 116 116 100 100 118 116 118
## [1009] 114 114 117 118 112 116 118 108 117 101 103 126 116 115 118 116 117 102
## [1027] 113 112 102 117 117 116 118 100 118 117 108 116 101 118 114 117 118 116
## [1045] 112 117 117 118 118 118 115 117 114 117 111 114 114 117 117 116 117 119
## [1063] 118 102 102 118 113 127 117 118 115 118 117 117 102 111 116 118 118 114
## [1081] 117 115 117 117 119 127 117 117 118 102 113 118 116 127 118 116 126 100
## [1099] 100 127 118 117 115 115 116 114 119 111 118 118 126 113 117 119 116 116
## [1117] 115 114 102 115 117 115 119 117 116 102 115 117 116 118 100 117 115 118
## [1135] 116 117 117 116 116 118 115 100 117 116 100 116 102 116 115 119 118 127
## [1153] 102 102 115 114 117 126 115 100 115 115 118 119 117 102 118 126 115 114
## [1171] 101 116 116 116 111 113 118 116 116 116 102 117 114 115 117 113 116 100
## [1189] 116 117 117 116 113 114 100 115 116 115 115 115 114 115 115 115 114 100
## [1207] 116 116 115 117 108 113 117 114 115 116 113 118 116 116 115 115 100 115
## [1225] 108 114 117 117 114 111 117 117 117 100 117 116 114 117 115 119 109 104
## [1243] 109 109 117 119 115 111 109 111 119 100 115 114 114 118 115 109 119 115
## [1261] 116 102 102 116 115 116 127 117 117 117 117 117 117 116 100 117 117 117
## [1279] 102 115 117 118 115 100 118 117 114 117 116 118 118 118 117 126 116 117
## [1297] 117 117 118 114 116 116 103 116 117 119 116 115 117 116 102 100 117 116
## [1315] 104 115 118 116 116 102 119 114 118 115 115 118 135 114 119 117 117 113
## [1333] 100 113 104 102 117 100 111 117 114 117 100 114 100 100 116 119 118 116
## [1351] 114 119 115 118 115 111 93 115 116 116 120 110 118 118 114 116 119 117
## [1369] 119 115 119 119 118 115 117 115 116 116 116 117 116 113 127 100 113 115
## [1387] 102 118 103 117 118 115 114 117 116 117 110 114 117 103 119 115 117 115
## [1405] 114 102 117 115 116 116 118 117 127 118 117 118 114 100 115 110 115 135
## [1423] 115 135 117 117 115 115 111 116 116 115 118 115 118 117 135 118 113 115
## [1441] 117 117 127 115 116 100 135 116 115 111 117 118 102 116 109 117 117 113
## [1459] 120 115 100 127 118 127 117 115 113 115 117 118 117 115 102 113 126 117
## [1477] 117 119 103 115 116 100 115 116 115 117 115 117 113 100 101 102 114 115
## [1495] 102 101 116 115 100 116 102 117 117 115 116 116 115 117 115 115 104 115
## [1513] 117 117 117 117 115 103 117 114 105 114 115 115 103 116 118 103 115 100
## [1531] 115 116 117 118 100 118 117 100 120 117 101 120 117 101 115 117 116 115
## [1549] 109 116gammaPars <- fitdistr(data$myVar, "gamma")
# shape rate
# 321.967275 2.820173
# ( 11.556479) ( 0.101304)
shapeML <- gammaPars$estimate["shape"]
rateML <- gammaPars$estimate["rate"]
sim_data <- data.frame(myVar=rgamma(n=1550, shape=shapeML, rate=rateML))
print(sim_data)## myVar
## 1 120.65311
## 2 111.62617
## 3 104.97734
## 4 103.76994
## 5 109.30203
## 6 100.04502
## 7 121.59968
## 8 111.83596
## 9 121.81386
## 10 118.85732
## 11 108.65200
## 12 118.37263
## 13 115.45388
## 14 118.18821
## 15 118.34363
## 16 120.92640
## 17 119.81436
## 18 110.29291
## 19 111.78549
## 20 127.65958
## 21 113.33367
## 22 124.17942
## 23 106.83113
## 24 116.33113
## 25 113.64823
## 26 116.88586
## 27 113.27281
## 28 115.75474
## 29 111.98699
## 30 106.73495
## 31 120.82234
## 32 109.98257
## 33 117.82288
## 34 117.71537
## 35 106.66268
## 36 119.64596
## 37 116.77091
## 38 120.35263
## 39 116.73426
## 40 101.32288
## 41 119.12131
## 42 120.06824
## 43 101.81483
## 44 110.74363
## 45 111.04718
## 46 126.32579
## 47 109.47957
## 48 109.74457
## 49 115.38255
## 50 112.71352
## 51 103.94115
## 52 120.94450
## 53 108.30152
## 54 114.17327
## 55 110.34601
## 56 106.59145
## 57 111.96088
## 58 100.87631
## 59 107.90881
## 60 109.46415
## 61 124.24800
## 62 111.93420
## 63 117.79107
## 64 113.23104
## 65 116.28443
## 66 117.50757
## 67 106.54467
## 68 113.14593
## 69 108.08240
## 70 108.80087
## 71 110.04003
## 72 114.75545
## 73 123.07056
## 74 118.61783
## 75 121.06454
## 76 107.17208
## 77 115.72460
## 78 111.79111
## 79 104.34428
## 80 113.56725
## 81 119.15067
## 82 110.48234
## 83 101.89913
## 84 106.95034
## 85 117.84524
## 86 113.27174
## 87 118.41225
## 88 116.56984
## 89 114.21607
## 90 112.36474
## 91 111.54822
## 92 116.48679
## 93 117.18350
## 94 129.82237
## 95 109.13018
## 96 118.79766
## 97 118.48480
## 98 111.54777
## 99 125.62830
## 100 105.00543
## 101 111.61627
## 102 111.65559
## 103 116.84433
## 104 114.24914
## 105 118.62933
## 106 120.93481
## 107 110.97384
## 108 125.97400
## 109 108.91304
## 110 119.44895
## 111 119.97242
## 112 110.89698
## 113 116.58769
## 114 117.27074
## 115 116.71989
## 116 119.93343
## 117 105.19690
## 118 117.15749
## 119 116.32219
## 120 113.60255
## 121 114.46559
## 122 118.60884
## 123 121.11671
## 124 111.85970
## 125 127.14122
## 126 119.23651
## 127 120.60902
## 128 115.06196
## 129 112.17094
## 130 122.53056
## 131 110.04095
## 132 117.72942
## 133 114.20366
## 134 121.31731
## 135 121.45561
## 136 107.73881
## 137 122.37439
## 138 110.21868
## 139 113.89878
## 140 107.30209
## 141 121.71812
## 142 116.32972
## 143 114.62295
## 144 109.51927
## 145 128.37981
## 146 114.15315
## 147 117.51836
## 148 106.89053
## 149 120.45028
## 150 115.74632
## 151 120.52408
## 152 115.81107
## 153 111.87747
## 154 120.50771
## 155 124.13922
## 156 115.51181
## 157 108.34690
## 158 118.96801
## 159 118.13350
## 160 120.94622
## 161 116.05594
## 162 121.83363
## 163 123.35581
## 164 107.63958
## 165 118.59533
## 166 110.80161
## 167 119.47677
## 168 105.82644
## 169 115.40075
## 170 113.65631
## 171 109.89732
## 172 108.53349
## 173 123.68949
## 174 121.76326
## 175 117.09581
## 176 114.65580
## 177 119.82267
## 178 121.80989
## 179 108.77819
## 180 113.63252
## 181 113.74964
## 182 113.13861
## 183 111.28331
## 184 125.94300
## 185 111.15053
## 186 114.37976
## 187 103.56412
## 188 103.09202
## 189 116.70945
## 190 112.30553
## 191 108.00108
## 192 109.08703
## 193 117.88741
## 194 116.54690
## 195 115.07312
## 196 108.59938
## 197 112.27111
## 198 115.63991
## 199 120.70417
## 200 114.42087
## 201 113.35342
## 202 124.94708
## 203 113.58178
## 204 102.50988
## 205 122.54437
## 206 123.03685
## 207 115.92287
## 208 111.31890
## 209 116.93211
## 210 118.88105
## 211 121.45268
## 212 105.09374
## 213 123.99219
## 214 109.78509
## 215 109.57345
## 216 112.63878
## 217 125.97081
## 218 118.11949
## 219 120.57243
## 220 119.21199
## 221 108.05821
## 222 106.40818
## 223 110.20128
## 224 118.08316
## 225 103.33396
## 226 110.46569
## 227 114.13566
## 228 111.13020
## 229 114.78352
## 230 99.95078
## 231 111.43867
## 232 118.12713
## 233 112.36849
## 234 117.60498
## 235 106.71608
## 236 120.66478
## 237 100.86013
## 238 120.31520
## 239 111.22587
## 240 112.53928
## 241 116.79850
## 242 119.62728
## 243 126.34236
## 244 110.40765
## 245 100.70277
## 246 110.39759
## 247 116.58342
## 248 123.03244
## 249 113.62736
## 250 98.27796
## 251 111.83802
## 252 103.22871
## 253 118.66463
## 254 124.46535
## 255 113.60729
## 256 118.48873
## 257 126.52366
## 258 119.14163
## 259 112.33390
## 260 108.25885
## 261 122.61528
## 262 114.42948
## 263 112.46351
## 264 107.36232
## 265 114.80760
## 266 109.56375
## 267 113.54977
## 268 109.68553
## 269 117.31499
## 270 113.87877
## 271 110.02325
## 272 112.47105
## 273 109.39365
## 274 118.49179
## 275 113.29025
## 276 116.38304
## 277 118.77190
## 278 105.89436
## 279 99.80010
## 280 114.65587
## 281 107.16027
## 282 110.80971
## 283 121.23848
## 284 126.95782
## 285 111.63897
## 286 102.34778
## 287 111.86574
## 288 115.12989
## 289 121.66300
## 290 124.00992
## 291 121.68802
## 292 111.15915
## 293 110.89140
## 294 108.90808
## 295 129.79024
## 296 107.19791
## 297 103.95938
## 298 127.43810
## 299 120.75281
## 300 118.28250
## 301 100.81299
## 302 123.72091
## 303 107.55625
## 304 107.40690
## 305 112.16039
## 306 115.81825
## 307 114.13387
## 308 110.30513
## 309 113.74964
## 310 111.24565
## 311 106.46262
## 312 120.16015
## 313 118.66097
## 314 110.98721
## 315 127.07480
## 316 119.65602
## 317 115.01675
## 318 111.89265
## 319 123.74399
## 320 118.35742
## 321 114.40319
## 322 113.45524
## 323 106.38889
## 324 124.10506
## 325 118.73410
## 326 104.65092
## 327 125.83740
## 328 110.04914
## 329 116.70236
## 330 120.04451
## 331 115.80084
## 332 116.37686
## 333 113.14594
## 334 117.48841
## 335 108.94647
## 336 109.17817
## 337 107.83450
## 338 116.77634
## 339 108.36722
## 340 120.17410
## 341 119.45182
## 342 121.11013
## 343 110.48138
## 344 125.13119
## 345 123.43294
## 346 111.85352
## 347 125.92957
## 348 114.02165
## 349 114.86099
## 350 108.40400
## 351 113.50834
## 352 117.31683
## 353 118.92797
## 354 108.56518
## 355 118.38775
## 356 108.61596
## 357 107.89851
## 358 117.13815
## 359 112.58546
## 360 116.22959
## 361 117.91129
## 362 129.83575
## 363 121.57186
## 364 114.06026
## 365 109.36222
## 366 104.43932
## 367 122.74783
## 368 120.99498
## 369 102.41259
## 370 122.57199
## 371 112.51407
## 372 106.64353
## 373 117.31797
## 374 112.43846
## 375 105.98414
## 376 119.74025
## 377 116.24902
## 378 120.93633
## 379 119.30561
## 380 96.22491
## 381 109.46877
## 382 118.67939
## 383 110.93793
## 384 115.20400
## 385 102.33390
## 386 114.49552
## 387 116.55941
## 388 114.56397
## 389 110.91732
## 390 133.78924
## 391 116.82412
## 392 123.14798
## 393 114.61838
## 394 119.27425
## 395 119.23208
## 396 121.32355
## 397 119.38821
## 398 112.52145
## 399 110.15539
## 400 116.96115
## 401 107.94363
## 402 116.04190
## 403 113.89395
## 404 112.33049
## 405 121.06223
## 406 119.12156
## 407 115.88052
## 408 121.94541
## 409 112.18052
## 410 113.09908
## 411 122.86066
## 412 112.38817
## 413 123.60145
## 414 118.89651
## 415 113.78932
## 416 112.77905
## 417 114.86420
## 418 114.86502
## 419 116.46226
## 420 106.97340
## 421 118.07347
## 422 111.05874
## 423 118.87943
## 424 119.88668
## 425 103.21430
## 426 110.87580
## 427 110.47326
## 428 108.92896
## 429 130.06997
## 430 114.34459
## 431 120.43701
## 432 114.62637
## 433 126.81655
## 434 118.29842
## 435 111.51286
## 436 110.03707
## 437 110.72717
## 438 117.68008
## 439 110.76394
## 440 103.98462
## 441 100.68700
## 442 121.71624
## 443 120.74513
## 444 123.20601
## 445 118.52454
## 446 104.12707
## 447 102.51193
## 448 117.24321
## 449 105.03880
## 450 119.02457
## 451 114.99143
## 452 104.59774
## 453 114.81590
## 454 115.38200
## 455 113.01282
## 456 111.11562
## 457 113.27675
## 458 102.62509
## 459 109.96509
## 460 114.68693
## 461 116.11532
## 462 101.40577
## 463 106.07118
## 464 116.70342
## 465 117.05983
## 466 97.99271
## 467 109.89243
## 468 109.10361
## 469 118.40155
## 470 126.31873
## 471 121.73565
## 472 115.94661
## 473 117.33583
## 474 112.10482
## 475 127.99756
## 476 117.58258
## 477 110.41411
## 478 120.52988
## 479 121.63799
## 480 108.65226
## 481 116.86707
## 482 111.47709
## 483 110.18063
## 484 109.21178
## 485 112.57435
## 486 112.31997
## 487 117.81780
## 488 102.14385
## 489 122.56182
## 490 115.48460
## 491 103.68983
## 492 110.98642
## 493 112.10363
## 494 113.08594
## 495 103.37426
## 496 119.49163
## 497 111.97680
## 498 112.66620
## 499 122.12623
## 500 110.44235
## 501 120.96261
## 502 114.25422
## 503 111.05697
## 504 127.24278
## 505 113.13365
## 506 110.51703
## 507 119.49684
## 508 111.49310
## 509 104.69958
## 510 124.84247
## 511 122.49593
## 512 113.72686
## 513 117.56671
## 514 110.52743
## 515 107.94695
## 516 101.90257
## 517 129.30958
## 518 108.63523
## 519 114.46252
## 520 116.52849
## 521 122.97830
## 522 126.43438
## 523 113.96500
## 524 114.36272
## 525 118.48576
## 526 120.85132
## 527 127.69343
## 528 114.50711
## 529 117.57472
## 530 109.09581
## 531 109.99790
## 532 121.02436
## 533 120.29288
## 534 111.25902
## 535 104.36435
## 536 113.87356
## 537 112.59878
## 538 119.87261
## 539 102.75154
## 540 118.46097
## 541 105.83186
## 542 121.72101
## 543 119.27546
## 544 114.42219
## 545 122.37844
## 546 106.54789
## 547 111.24286
## 548 117.87161
## 549 116.05600
## 550 122.23772
## 551 111.07465
## 552 114.03112
## 553 117.87395
## 554 122.56973
## 555 110.20420
## 556 99.13556
## 557 101.56923
## 558 113.82342
## 559 116.28310
## 560 112.82828
## 561 106.79919
## 562 111.71176
## 563 115.57318
## 564 111.88126
## 565 117.85068
## 566 117.70882
## 567 107.45835
## 568 105.90429
## 569 119.29875
## 570 108.94405
## 571 113.05517
## 572 110.86839
## 573 109.56455
## 574 107.25021
## 575 102.09514
## 576 125.52907
## 577 119.44660
## 578 114.49830
## 579 115.28589
## 580 122.86242
## 581 118.74323
## 582 117.05170
## 583 119.51635
## 584 117.03127
## 585 107.99071
## 586 111.22500
## 587 124.89363
## 588 123.78959
## 589 119.30360
## 590 117.33408
## 591 106.21866
## 592 102.85284
## 593 116.13957
## 594 109.95220
## 595 118.33955
## 596 115.99519
## 597 115.55646
## 598 116.86267
## 599 113.07138
## 600 117.46301
## 601 109.03311
## 602 107.75486
## 603 103.37785
## 604 116.75335
## 605 118.39727
## 606 118.65455
## 607 118.94770
## 608 115.94005
## 609 112.42419
## 610 114.62161
## 611 119.86581
## 612 115.10730
## 613 123.72748
## 614 109.20105
## 615 115.98336
## 616 107.73026
## 617 111.52061
## 618 108.55478
## 619 108.25242
## 620 114.53423
## 621 111.68888
## 622 110.87523
## 623 113.96194
## 624 120.26529
## 625 106.27456
## 626 112.21822
## 627 123.91783
## 628 105.97087
## 629 116.74484
## 630 115.59109
## 631 120.48515
## 632 123.21707
## 633 120.96013
## 634 101.22473
## 635 126.87051
## 636 107.32773
## 637 111.45471
## 638 117.94704
## 639 112.67173
## 640 127.72043
## 641 114.05801
## 642 122.36572
## 643 120.43763
## 644 125.76562
## 645 116.55418
## 646 118.50038
## 647 109.14212
## 648 121.17105
## 649 120.20786
## 650 115.49805
## 651 115.37650
## 652 122.90115
## 653 116.83800
## 654 126.88627
## 655 118.63060
## 656 112.09024
## 657 120.84014
## 658 112.28792
## 659 114.62161
## 660 114.79907
## 661 125.86647
## 662 127.55258
## 663 116.38845
## 664 123.51763
## 665 101.46287
## 666 115.39749
## 667 108.46403
## 668 107.45025
## 669 119.75600
## 670 112.33650
## 671 107.47334
## 672 117.51164
## 673 124.16357
## 674 116.77540
## 675 110.48364
## 676 117.23981
## 677 125.59026
## 678 106.49418
## 679 113.03801
## 680 115.84969
## 681 109.48361
## 682 102.36455
## 683 113.59280
## 684 96.12437
## 685 128.42602
## 686 125.84034
## 687 109.96822
## 688 114.45556
## 689 123.22247
## 690 116.53660
## 691 119.62153
## 692 117.29462
## 693 112.55062
## 694 101.20515
## 695 108.82934
## 696 105.44592
## 697 109.26730
## 698 114.67069
## 699 124.20232
## 700 123.66079
## 701 121.14226
## 702 105.78272
## 703 121.47771
## 704 120.40494
## 705 122.00588
## 706 100.18015
## 707 114.16659
## 708 110.83178
## 709 120.11424
## 710 108.30964
## 711 113.05645
## 712 116.14343
## 713 119.63702
## 714 117.72433
## 715 116.34668
## 716 117.89278
## 717 114.02154
## 718 115.28813
## 719 117.04725
## 720 122.81961
## 721 118.13652
## 722 118.50667
## 723 120.99700
## 724 112.38653
## 725 110.55860
## 726 114.00673
## 727 109.27968
## 728 111.53129
## 729 122.11942
## 730 120.62727
## 731 119.66793
## 732 108.99965
## 733 113.34031
## 734 113.04793
## 735 119.42822
## 736 114.49316
## 737 108.38074
## 738 107.62191
## 739 116.69430
## 740 105.45842
## 741 105.27183
## 742 127.43247
## 743 120.78481
## 744 102.06368
## 745 115.87950
## 746 102.81359
## 747 121.16002
## 748 111.26094
## 749 124.55771
## 750 111.20620
## 751 113.40028
## 752 112.05318
## 753 116.39517
## 754 115.30661
## 755 117.85140
## 756 104.81887
## 757 118.58848
## 758 115.85942
## 759 112.99237
## 760 114.40448
## 761 118.86732
## 762 118.29130
## 763 114.54700
## 764 106.90991
## 765 108.75181
## 766 123.58844
## 767 114.14979
## 768 114.53923
## 769 114.92256
## 770 118.67924
## 771 126.42803
## 772 125.64062
## 773 119.28471
## 774 116.30126
## 775 116.36028
## 776 113.24238
## 777 113.83599
## 778 109.46854
## 779 103.26540
## 780 113.95285
## 781 114.39977
## 782 113.19189
## 783 106.67449
## 784 114.72514
## 785 96.00080
## 786 116.94007
## 787 126.37578
## 788 119.20380
## 789 102.14877
## 790 111.23354
## 791 124.63526
## 792 108.12637
## 793 118.45390
## 794 126.32164
## 795 118.54211
## 796 108.21735
## 797 132.61531
## 798 115.67399
## 799 100.65525
## 800 113.16189
## 801 116.98172
## 802 103.53070
## 803 115.11568
## 804 120.11335
## 805 109.10751
## 806 110.35873
## 807 106.56219
## 808 118.07198
## 809 112.30941
## 810 116.89600
## 811 117.67879
## 812 117.54305
## 813 126.36835
## 814 118.88078
## 815 115.49444
## 816 106.21276
## 817 122.00681
## 818 115.24530
## 819 115.35131
## 820 123.05692
## 821 113.60563
## 822 108.02822
## 823 110.65783
## 824 96.21999
## 825 107.25754
## 826 119.44237
## 827 108.56032
## 828 121.10737
## 829 109.69577
## 830 127.53977
## 831 107.51343
## 832 107.80397
## 833 110.78190
## 834 119.53646
## 835 101.95686
## 836 108.99195
## 837 111.51785
## 838 127.20894
## 839 112.66662
## 840 108.70958
## 841 115.17368
## 842 122.73318
## 843 102.88378
## 844 108.70029
## 845 121.84143
## 846 126.74446
## 847 106.34174
## 848 107.65631
## 849 107.31127
## 850 119.97418
## 851 109.49515
## 852 98.82029
## 853 99.45553
## 854 114.15930
## 855 106.96730
## 856 123.36217
## 857 104.81154
## 858 113.37558
## 859 109.38071
## 860 112.02864
## 861 114.13581
## 862 108.44751
## 863 116.89962
## 864 118.90887
## 865 112.58988
## 866 104.41930
## 867 131.34067
## 868 117.90940
## 869 120.99204
## 870 114.28088
## 871 119.99232
## 872 110.82095
## 873 122.62258
## 874 123.99419
## 875 113.02862
## 876 108.76185
## 877 125.14874
## 878 113.68704
## 879 115.87625
## 880 100.45000
## 881 102.38865
## 882 109.88964
## 883 119.86717
## 884 108.89240
## 885 114.12890
## 886 115.75607
## 887 104.23633
## 888 102.85351
## 889 111.26697
## 890 116.32628
## 891 108.22443
## 892 113.07349
## 893 110.05897
## 894 111.46333
## 895 116.70315
## 896 107.86795
## 897 114.29878
## 898 118.64883
## 899 104.62646
## 900 115.41094
## 901 122.39055
## 902 123.71998
## 903 131.42209
## 904 115.77547
## 905 130.92730
## 906 105.37761
## 907 111.71724
## 908 112.19791
## 909 118.00278
## 910 110.90503
## 911 114.41610
## 912 115.00116
## 913 121.85785
## 914 121.72160
## 915 121.09707
## 916 122.43731
## 917 111.12217
## 918 109.46068
## 919 105.69710
## 920 110.50737
## 921 114.20986
## 922 120.75889
## 923 108.43723
## 924 113.47047
## 925 120.20248
## 926 118.82153
## 927 108.46038
## 928 107.72397
## 929 115.45475
## 930 114.66816
## 931 104.04150
## 932 109.73763
## 933 113.76633
## 934 111.73118
## 935 115.90323
## 936 117.84858
## 937 113.62479
## 938 116.60700
## 939 114.50965
## 940 112.24964
## 941 115.14963
## 942 120.99246
## 943 118.86237
## 944 103.71462
## 945 132.84244
## 946 132.03663
## 947 101.57013
## 948 112.41812
## 949 112.36712
## 950 104.07612
## 951 111.93690
## 952 123.08568
## 953 115.39565
## 954 125.93261
## 955 108.35385
## 956 119.70240
## 957 109.77706
## 958 116.74795
## 959 101.76995
## 960 121.80677
## 961 119.10560
## 962 116.15633
## 963 111.66990
## 964 111.34686
## 965 110.89061
## 966 110.44174
## 967 113.75092
## 968 109.79657
## 969 111.34030
## 970 116.88790
## 971 120.18138
## 972 124.03506
## 973 124.56035
## 974 120.53974
## 975 116.45955
## 976 107.89795
## 977 111.57146
## 978 110.78269
## 979 122.98981
## 980 111.46111
## 981 118.55060
## 982 132.74137
## 983 113.35316
## 984 113.60161
## 985 113.12467
## 986 117.01820
## 987 108.65456
## 988 115.46747
## 989 116.94805
## 990 103.01360
## 991 108.58562
## 992 104.64691
## 993 116.20362
## 994 108.81216
## 995 120.93100
## 996 118.90987
## 997 118.99991
## 998 103.04464
## 999 120.67592
## 1000 105.98332
## 1001 107.39676
## 1002 99.42732
## 1003 110.02874
## 1004 107.83698
## 1005 116.80460
## 1006 112.53923
## 1007 112.77830
## 1008 117.39734
## 1009 124.18814
## 1010 115.34649
## 1011 108.42574
## 1012 112.53779
## 1013 117.16994
## 1014 121.93180
## 1015 110.35529
## 1016 118.53323
## 1017 126.68058
## 1018 110.82617
## 1019 112.89810
## 1020 121.64189
## 1021 114.81144
## 1022 109.72482
## 1023 121.73654
## 1024 120.93752
## 1025 117.13175
## 1026 113.28760
## 1027 108.72849
## 1028 120.87175
## 1029 107.90905
## 1030 114.42077
## 1031 107.20326
## 1032 123.95712
## 1033 113.83297
## 1034 113.42256
## 1035 103.34093
## 1036 117.31018
## 1037 114.22317
## 1038 114.26984
## 1039 114.53527
## 1040 122.20039
## 1041 119.60713
## 1042 120.65091
## 1043 112.55784
## 1044 107.75633
## 1045 104.59076
## 1046 119.68565
## 1047 115.04017
## 1048 113.43647
## 1049 115.04858
## 1050 111.68772
## 1051 117.34634
## 1052 102.86703
## 1053 106.40782
## 1054 120.21352
## 1055 118.85486
## 1056 119.60239
## 1057 122.42786
## 1058 125.39507
## 1059 109.73039
## 1060 129.24019
## 1061 104.35494
## 1062 121.32563
## 1063 120.51355
## 1064 110.68463
## 1065 109.70743
## 1066 114.31979
## 1067 122.47630
## 1068 109.33660
## 1069 116.63026
## 1070 104.13307
## 1071 114.02190
## 1072 112.96301
## 1073 108.63326
## 1074 115.52846
## 1075 118.12090
## 1076 102.98863
## 1077 109.99656
## 1078 106.20855
## 1079 108.29314
## 1080 101.12784
## 1081 111.50772
## 1082 115.66501
## 1083 115.11895
## 1084 118.02281
## 1085 120.37978
## 1086 122.94257
## 1087 120.37196
## 1088 107.66282
## 1089 122.27264
## 1090 121.80008
## 1091 114.20818
## 1092 115.60338
## 1093 120.07308
## 1094 104.56740
## 1095 106.18902
## 1096 113.08772
## 1097 123.65409
## 1098 111.34184
## 1099 117.71781
## 1100 109.27216
## 1101 116.28959
## 1102 109.83281
## 1103 113.97744
## 1104 123.93122
## 1105 115.91349
## 1106 118.62359
## 1107 116.88809
## 1108 113.31973
## 1109 115.56216
## 1110 123.70220
## 1111 117.49252
## 1112 111.87122
## 1113 122.81465
## 1114 108.13230
## 1115 108.43394
## 1116 110.97207
## 1117 113.56092
## 1118 122.32425
## 1119 115.35797
## 1120 109.94631
## 1121 117.28133
## 1122 114.91506
## 1123 105.48877
## 1124 112.30777
## 1125 108.97363
## 1126 112.51868
## 1127 115.51567
## 1128 118.79914
## 1129 108.32687
## 1130 120.83004
## 1131 118.52322
## 1132 116.91906
## 1133 127.01978
## 1134 107.44997
## 1135 126.22792
## 1136 112.45110
## 1137 123.05301
## 1138 112.51350
## 1139 116.41741
## 1140 111.88860
## 1141 109.63609
## 1142 116.89757
## 1143 122.61168
## 1144 112.74659
## 1145 116.41676
## 1146 106.53811
## 1147 122.07810
## 1148 137.02480
## 1149 117.74585
## 1150 118.39682
## 1151 106.99130
## 1152 116.17625
## 1153 108.09086
## 1154 116.35265
## 1155 120.39358
## 1156 109.74894
## 1157 116.03596
## 1158 99.01670
## 1159 121.79667
## 1160 119.98497
## 1161 104.96032
## 1162 121.01812
## 1163 116.55244
## 1164 107.58994
## 1165 110.01751
## 1166 114.22114
## 1167 112.42751
## 1168 120.35666
## 1169 113.72492
## 1170 113.59464
## 1171 111.29093
## 1172 112.34751
## 1173 113.52420
## 1174 117.20471
## 1175 119.05001
## 1176 113.15087
## 1177 124.54592
## 1178 113.52543
## 1179 106.98187
## 1180 116.51891
## 1181 105.61484
## 1182 111.49334
## 1183 100.43088
## 1184 115.90246
## 1185 116.37439
## 1186 111.01936
## 1187 116.14034
## 1188 104.09496
## 1189 105.97810
## 1190 124.30307
## 1191 108.60817
## 1192 117.62985
## 1193 118.30709
## 1194 122.41795
## 1195 125.77975
## 1196 109.64150
## 1197 125.85815
## 1198 110.92823
## 1199 122.84808
## 1200 104.82501
## 1201 133.62583
## 1202 118.28384
## 1203 116.74694
## 1204 101.81029
## 1205 113.24943
## 1206 105.48394
## 1207 111.30572
## 1208 112.61971
## 1209 109.74968
## 1210 109.79495
## 1211 109.32231
## 1212 109.29901
## 1213 111.38644
## 1214 123.56745
## 1215 117.98234
## 1216 118.43184
## 1217 105.84193
## 1218 109.96097
## 1219 104.56164
## 1220 104.62480
## 1221 109.61287
## 1222 122.65344
## 1223 104.44920
## 1224 108.88404
## 1225 121.61847
## 1226 112.71167
## 1227 126.72238
## 1228 107.09429
## 1229 109.68591
## 1230 101.97364
## 1231 99.25005
## 1232 107.46512
## 1233 125.30189
## 1234 115.40665
## 1235 110.39091
## 1236 113.76895
## 1237 115.29665
## 1238 109.43016
## 1239 110.46057
## 1240 100.85908
## 1241 113.92393
## 1242 117.92551
## 1243 117.91616
## 1244 121.87788
## 1245 114.96081
## 1246 116.07141
## 1247 98.08407
## 1248 112.60613
## 1249 111.33056
## 1250 112.81289
## 1251 128.46114
## 1252 112.55117
## 1253 121.78702
## 1254 117.98642
## 1255 112.05869
## 1256 112.17444
## 1257 125.21526
## 1258 115.07260
## 1259 113.33808
## 1260 113.53907
## 1261 111.75330
## 1262 105.01567
## 1263 113.32399
## 1264 116.14636
## 1265 114.58627
## 1266 121.30034
## 1267 122.03011
## 1268 115.16407
## 1269 115.38944
## 1270 116.37500
## 1271 114.36458
## 1272 111.07091
## 1273 124.99094
## 1274 121.65442
## 1275 110.99704
## 1276 117.03834
## 1277 106.57677
## 1278 117.59494
## 1279 106.84696
## 1280 113.03181
## 1281 117.97754
## 1282 108.61366
## 1283 106.38725
## 1284 107.22261
## 1285 105.93956
## 1286 117.25610
## 1287 123.35486
## 1288 121.84363
## 1289 116.27857
## 1290 110.70243
## 1291 115.68645
## 1292 113.42808
## 1293 115.96382
## 1294 115.46266
## 1295 101.53327
## 1296 101.97036
## 1297 116.91082
## 1298 112.64344
## 1299 115.45398
## 1300 114.49709
## 1301 123.84321
## 1302 122.26496
## 1303 112.36225
## 1304 115.70328
## 1305 108.58859
## 1306 107.23245
## 1307 112.97223
## 1308 124.69150
## 1309 125.16647
## 1310 120.92660
## 1311 108.62908
## 1312 113.62180
## 1313 122.97960
## 1314 114.47642
## 1315 109.19868
## 1316 109.79452
## 1317 109.90556
## 1318 118.44801
## 1319 110.21642
## 1320 105.06789
## 1321 113.39167
## 1322 116.59223
## 1323 116.54769
## 1324 106.38990
## 1325 109.03869
## 1326 110.61229
## 1327 115.79207
## 1328 118.93880
## 1329 114.30905
## 1330 118.11469
## 1331 115.71046
## 1332 112.57247
## 1333 108.07279
## 1334 131.76577
## 1335 119.86693
## 1336 125.41810
## 1337 112.00091
## 1338 107.07893
## 1339 117.41528
## 1340 118.51856
## 1341 115.96032
## 1342 114.88809
## 1343 113.99888
## 1344 115.86342
## 1345 109.30283
## 1346 119.20762
## 1347 108.80144
## 1348 109.77109
## 1349 117.41817
## 1350 122.09066
## 1351 125.86479
## 1352 118.57707
## 1353 119.57007
## 1354 113.08948
## 1355 116.53206
## 1356 121.12793
## 1357 122.09735
## 1358 117.37581
## 1359 113.03700
## 1360 123.31345
## 1361 105.28295
## 1362 108.65872
## 1363 120.38605
## 1364 116.86841
## 1365 118.66964
## 1366 130.47491
## 1367 110.39091
## 1368 111.67616
## 1369 108.17790
## 1370 115.18128
## 1371 102.07697
## 1372 108.73700
## 1373 113.55302
## 1374 115.97214
## 1375 116.80995
## 1376 118.28639
## 1377 118.62764
## 1378 109.52818
## 1379 111.18679
## 1380 115.81685
## 1381 115.66672
## 1382 111.19914
## 1383 127.67876
## 1384 106.32345
## 1385 110.72494
## 1386 110.75158
## 1387 107.12961
## 1388 110.80983
## 1389 119.63880
## 1390 112.09906
## 1391 124.30318
## 1392 122.01000
## 1393 117.17351
## 1394 102.38084
## 1395 122.20356
## 1396 123.68668
## 1397 110.61227
## 1398 118.42915
## 1399 114.05236
## 1400 116.66375
## 1401 106.68414
## 1402 115.03455
## 1403 116.79569
## 1404 107.52210
## 1405 117.11242
## 1406 108.25120
## 1407 114.47757
## 1408 112.65300
## 1409 109.07116
## 1410 106.99030
## 1411 111.18925
## 1412 118.03157
## 1413 114.26045
## 1414 127.29533
## 1415 103.23449
## 1416 112.98530
## 1417 113.27475
## 1418 126.40289
## 1419 118.64690
## 1420 114.72266
## 1421 118.53123
## 1422 120.72132
## 1423 116.24403
## 1424 111.06906
## 1425 122.80490
## 1426 115.83929
## 1427 117.84073
## 1428 105.91110
## 1429 115.59821
## 1430 123.53564
## 1431 109.02291
## 1432 100.67064
## 1433 114.61600
## 1434 113.12137
## 1435 122.47315
## 1436 112.34474
## 1437 104.45614
## 1438 116.56086
## 1439 109.37724
## 1440 106.37446
## 1441 113.58377
## 1442 106.63404
## 1443 116.74376
## 1444 124.48867
## 1445 115.22871
## 1446 118.66825
## 1447 110.41522
## 1448 113.85012
## 1449 114.21785
## 1450 115.70182
## 1451 118.31209
## 1452 118.12072
## 1453 111.08750
## 1454 109.10689
## 1455 119.52703
## 1456 99.04859
## 1457 126.64506
## 1458 110.98295
## 1459 128.02487
## 1460 108.00945
## 1461 122.05885
## 1462 118.17640
## 1463 114.93315
## 1464 113.33343
## 1465 127.82365
## 1466 125.59832
## 1467 109.39209
## 1468 108.25975
## 1469 116.71491
## 1470 116.51490
## 1471 117.11950
## 1472 107.46223
## 1473 99.14449
## 1474 111.78634
## 1475 117.43609
## 1476 121.01825
## 1477 130.33615
## 1478 114.14506
## 1479 113.79779
## 1480 115.01271
## 1481 115.40507
## 1482 116.27281
## 1483 106.99765
## 1484 111.15927
## 1485 114.19564
## 1486 118.03660
## 1487 117.38057
## 1488 107.06130
## 1489 108.62669
## 1490 104.00006
## 1491 107.72538
## 1492 121.02457
## 1493 114.62937
## 1494 119.03128
## 1495 102.64641
## 1496 124.33428
## 1497 118.74612
## 1498 111.53948
## 1499 131.08193
## 1500 103.42940
## 1501 124.79388
## 1502 117.86687
## 1503 118.26596
## 1504 115.51566
## 1505 120.18002
## 1506 110.88156
## 1507 115.71475
## 1508 111.19558
## 1509 114.01921
## 1510 114.86283
## 1511 104.43688
## 1512 123.04082
## 1513 123.24534
## 1514 105.74978
## 1515 108.70097
## 1516 123.27749
## 1517 96.75194
## 1518 119.01621
## 1519 111.99872
## 1520 113.26512
## 1521 110.48632
## 1522 116.77566
## 1523 111.43287
## 1524 115.98805
## 1525 108.75393
## 1526 121.04364
## 1527 125.04929
## 1528 105.39378
## 1529 109.31036
## 1530 107.30189
## 1531 112.97711
## 1532 112.68369
## 1533 120.22364
## 1534 114.74826
## 1535 107.05027
## 1536 115.84235
## 1537 106.04065
## 1538 109.33356
## 1539 127.55640
## 1540 107.11656
## 1541 113.65302
## 1542 116.34802
## 1543 115.14639
## 1544 109.70165
## 1545 112.25580
## 1546 109.05799
## 1547 114.10387
## 1548 113.60222
## 1549 116.69471
## 1550 114.06747sim_plot <- ggplot(data=sim_data, aes(x=myVar, y=..density..)) +
geom_histogram(color="grey60",fill="cornsilk",size=0.2)
print(sim_plot) # plot a histogram of the data## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
sim_curve <- stat_function(aes(x = xval, y = ..y..), fun = dgamma, colour="brown", n = length(sim_data$myVar), args = list(shape=shapeML, rate=rateML))
sim_plot + stat4 # plot gamma probability density; brown line## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#Imported data stuff:
z <- read.table("../Bio381/MyDataFile2.csv",header=TRUE,sep=",")
names(z) <- list("ID","pop", "myVar")
z <- z[z$myVar>0,]
str(z)## 'data.frame': 1550 obs. of 3 variables:
## $ ID : chr "ALB-M011A01" "ALB-M011A02" "ALB-M011A03" "ALB-M011A04" ...
## $ pop : chr "ALB" "ALB" "ALB" "ALB" ...
## $ myVar: int 113 114 116 117 100 127 115 115 115 100 ...summary(z$myVar)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 93.0 114.0 116.0 114.2 117.0 135.0str(z)## 'data.frame': 1550 obs. of 3 variables:
## $ ID : chr "ALB-M011A01" "ALB-M011A02" "ALB-M011A03" "ALB-M011A04" ...
## $ pop : chr "ALB" "ALB" "ALB" "ALB" ...
## $ myVar: int 113 114 116 117 100 127 115 115 115 100 ...summary(z)## ID pop myVar
## Length:1550 Length:1550 Min. : 93.0
## Class :character Class :character 1st Qu.:114.0
## Mode :character Mode :character Median :116.0
## Mean :114.2
## 3rd Qu.:117.0
## Max. :135.0p1 <- ggplot(data=z, aes(x=myVar, y=..density..)) +
geom_histogram(color="grey60",fill="cornsilk",size=0.2)
print(p1) # plot a histogram of the data## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
gammaPars <- fitdistr(z$myVar,"gamma")
shapeML <- gammaPars$estimate["shape"]
rateML <- gammaPars$estimate["rate"]
stat4 <- stat_function(aes(x = xval, y = ..y..), fun = dgamma, colour="brown", n = length(z$myVar), args = list(shape=shapeML, rate=rateML))
p1 + stat4 # plot gamma probability density; brown line## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
The curve seems to match the simulated data better than my imported data. This means that the gamma distribution may not be the best probability distribution to describe the imported data. It’s not as obviously bimodal as the first data set I went through in question 3, but there is a slight peak before the large peak which skews the curve a bit. Therefore, the model does not do a great job of simulating realistic data, but the curve patterns are at least a little bit lined up with what we may expect to see.