-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathheli_report.Rmd
227 lines (163 loc) · 14.1 KB
/
heli_report.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
```{r echo=FALSE}
# This is the configuration section. There are two bits of configuration you'll
# Paste the publishing URL from your copy of the data spreadsheet below.
# PUT THE URL BETWEEN THE DOUBLE QUOTES ON LINE 6
csv_url<-"https://docs.google.com/spreadsheet/pub?key=0Avq8qkMnj5AOdGlTLWFVbVZWamx3VWxPZk1pdm9oVFE&single=true&gid=0&output=csv"
# Once you're collected some data, change these to the heli ids you want to compare. The ids are
# from the first column of your data spreadsheet.
heli_id_one<-"6x9.5"
heli_id_two<-"5x8"
##################################################################################
# You should not have to modify anything below this point.
##################################################################################
```
Paper Helicopter Optimization Project
========================================================
This is an R Markdown document. Markdown is a simple formatting syntax for authoring web pages (click the **MD** toolbar button for help on Markdown). When you click the **Knit HTML** button a web page will be generated that includes both content as well as the output of any embedded R code chunks within the document.
Looking at your data
--------------------
```{r include=FALSE}
# This chunk loads libraries that are needed below. include=FALSE means that neither the code nor the output of the code will be included in the output document.
library(RCurl)
library(ggplot2)
library(plyr)
library(pwr)
```
```{r echo=FALSE}
helidata<-read.csv(textConnection(getURL(csv_url)))
helidata<-transform(helidata,heli_id=factor(heli_id))
heli_means<-ddply(helidata,.(heli_id,rotor_width,rotor_length),summarize,avg_flight = mean(flight_time))
```
This plot is a histogram of flight times for each helicopter. The "binwidth", or width of each bar is 0.1 seconds.
```{r fig.width=7, fig.height=4}
ggplot(helidata,aes(flight_time, group=heli_id, fill=heli_id))+
geom_histogram(binwidth=0.1, position="dodge")
```
This plot shows a smoothed estimate of the distribution of flight time for each helicopter.
```{r fig.width=7, fig.height=4}
ggplot(helidata,aes(flight_time, group=heli_id, fill=heli_id))+
geom_density(alpha=1/2)
```
Phase 1: Is one really better?
-------------------------------
Once you've tested a few designs, you'll almost certainly be in the situation where one of the designs seems better in that it has a longer average flight duration, but because you have taken so few flights, you can't be sure that it really is better (maybe it's worse, but random variation has worked against, making it look better when more data would demonstrate that it really is worse).
```{r include=FALSE}
heli_one_data<-with(helidata, flight_time[heli_id==heli_id_one])
heli_two_data<-with(helidata, flight_time[heli_id==heli_id_two])
if (mean(heli_one_data) > mean(heli_two_data)) {
heli_id_long <- heli_id_one
heli_id_short <- heli_id_two
heli_long_data <- heli_one_data
heli_short_data <- heli_two_data
} else {
heli_short_data <- heli_one_data
heli_long_data <- heli_two_data
heli_id_short <- heli_id_one
heli_id_long <- heli_id_two
}
```
```{r echo=FALSE}
heli_long_n<-length(heli_long_data)
heli_short_n<-length(heli_short_data)
heli_long_mean<-ifelse(heli_long_n>0,mean(heli_long_data),NA)
heli_short_mean<-ifelse(heli_short_n>0,mean(heli_short_data),NA)
```
`r if (heli_long_n == 0) { "**WARNING: The id you entered for heli_id_long on line 11 does not exist in your data. The results below will probably be nonsensical and full of errors!** " }`
`r if (heli_short_n == 0) { print("**WARNING: The id you entered for heli_id_short on line 12 does not exist in your data. The results below will probably be nonsensical and full of errors!** ") }`
**Helicopter "`r heli_id_long`" has an average flight time of `r heli_long_mean` based on `r heli_long_n` flights. Helicopter "`r heli_id_short`" has an average flight time of `r heli_short_mean` based on `r heli_short_n` flights.** You think that heli "`r heli_id_long`" is better, but is the result statistically significant, or could it have happened just due to random chance?
Well, if there is no difference in the two helicopters, then the labels "`r heli_id_long`" and "`r heli_id_short`" are meaningless, and we could randomly shuffle the labels around and get results similar to our actual observations. If we randomly shuffle the labels around many times and consistently get results different from our observations, then we know that the labels *do* have meaning, and heli "`r heli_id_long`" really does have longer flights on average.
```{r}
# Extract a subset of the original data containing just the results from the two helis of interest
testdata<-subset(helidata, heli_id %in% c(heli_id_long, heli_id_short))
# Record which rows correspond to which helis
a_rows<-which(testdata$heli_id==heli_id_long)
b_rows<-which(testdata$heli_id==heli_id_short)
# sim_one() randomly shuffles the data and computes the difference in the average flight times
sim_one<- function() {
shuffled<-sample(testdata$flight_time)
return (mean(shuffled[a_rows])-mean(shuffled[b_rows]))
}
N <- 1000
# This line runs sim_one() N times and gathers the results of all the trials
simulated_differences<-replicate(N,sim_one())
```
```{r echo=FALSE}
if (any(is.nan(simulated_differences))) {
simulated_differences<-NA
}
```
Now that we have a `r N` simulated reshufflings, we can plot a histogram of the simulated results. The histogram shows gives you a sense for what you might have observed if there was no difference between your two helis -- that is, if the labels in the first column of your spreadsheet were random. Your actual observed difference between the two helicopters is marked with a red vertical line. If the red line is in the middle of the distribution, then your real data isn't very different from what could have occurred just by chance. If the red line is at the extreme edge of the distribution, then it's very unlikely that what you observed could have happened by chance, and you can be confident than heli "`r heli_id_long`" really has longer flights on average than heli "`r heli_id_short`".
```{r echo=FALSE, fig.height=4}
if (!any(is.na(simulated_differences))) {
qplot(simulated_differences,geom="histogram",binwidth=0.02)+
geom_vline(x=heli_long_mean-heli_short_mean,color="red")+
xlab("Difference in average flight time")
} else {
print("No data to plot!")
}
```
`r ct<-sum(simulated_differences>=(heli_long_mean-heli_short_mean)); ct` out of `r N` random shufflings resulted in a difference as large or larger than your actual observations. Traditionally, a statistician would say a result is statistically significant if it could occur at random less than 5%, or `r round(N*0.05)` times out of `r N`, so your result `r if (ct < round(N*0.05)) { "is" } else { "is not"}` statistically significant.
### An analytic test
If you're interested in another statistical test known as Student's t-test you can read this part. Otherwise you can safely skip it at no loss.
The test we conducted above is known as a permutation test, because we took many random permutations of the data. There is another test we could have run to decide whether the difference in average flights times was statistically significant, Student's t-test. If you took statistics class in school, you probably learned the t-test. I strongly prefer the permutation test. One reason is that it's more flexible (it works with any summary statistic, not just the mean), but the main reason to prefer the permutation test is that it's so beautifully intuitive.
That said, it's good to know Student's t-test because it is so widely known and used, and also because it is computationally more efficient (and just as powerful if your data adhere to its assumptions). I cannot possibly explain how the t-test works here, full development requires a semester long class with lots of calculus. I will just demonstrate running the test in R.
The code below does a Student's t-test to determine whether helicopter "`r heli_id_long`" has a longer average flight time and helicopter "`r heli_id_short`". The output might look all scary, but the only thing you really care about is the p-value. For this test, the p-value is the probability that the average flight time of helicopter "`r heli_id_long`" is longer than helicopter "`r heli_id_long`" just due to random variation. If the p-value is very small (traditionally, less than 0.05), then the there is a very low probability that the result is just due to chance, and you can be confident that helicopter "`r heli_id_long`" is really better than helicopter "`r heli_id_short`". In other words, you can say that your results are statistically significant.
```{r}
test_result<-t.test(heli_long_data, heli_short_data, alternative="greater")
test_result
```
The p-value from your data is `r test_result$p.value`, which `r if (test_result$p.value < 0.05) { "is" } else { "is not" }` less than 0.05, so you `r if (test_result$p.value < 0.05) { "can" } else { "can not" }` be confident that helicopter "`r heli_id_long`" is better than "`r heli_id_short`". If you're not confident, more data might help you decide.
### How much more data will I need?
If the results of your tests above are not significant, you might wonder how much more data you'll need to achieve a statistically significant result. This is a area of statistics known as power analysis.
The larger the difference in the average flight times, the fewer test flights you need to be confident in your results. The code below does a power analysis based on your data so far and predicts the number of flights you'll need to get a statistically significant result. If your average flight times are very close, you might need a **lot** of flights to be sure that one is better than the other.
```{r}
d<-abs(heli_long_mean-heli_short_mean)/sd(c(heli_long_data,heli_short_data))
pwr<-pwr.t.test(d=d, sig.level=0.05, power=0.95, alternative="greater")
pwr
```
If we assume that your current estimates of the average flights times are the *true* average flight times, it looks like you'll need to fly each helicopter `r ceiling(pwr$n)` times to be 95% sure you'll get a significant result.
Phase 2: Optimization: Modeling the response surface
----------------------------------------------------
If we call the flight time $F$, and the rotor length and width $l$, and $w$ respectively, then the simplest possible model for predicting flight times is:
$F(l,w) \sim \beta_1 l + \beta_2 w$
Where $\beta_1$, and $\beta_2$ are unknown coefficients to be fit from your observed data. The call to *lm* below fits this model from your data.
```{r}
model<-lm(flight_time ~ rotor_length + rotor_width, data=helidata)
summary(model)
```
The coefficients of this simple model suggest that adding every 1/2 centimeter added to the length of your rotor changes the average flight time by `r model$coefficients[2]` seconds, and each 1/2 centimeter of width changes the average flight time by `r model$coefficients[3]` seconds.
We can ask the model for predictions of every combination of rotor length and width in 0.5 cm increments.
```{r}
to_predict<-expand.grid(rotor_width=seq(0.5,6,0.5),
rotor_length=seq(0.5,9.5,0.5))
to_predict$predicted_flight_time<-predict(model, to_predict)
max_row<-to_predict[which.max(to_predict$predicted_flight_time),]
```
Below is a heat map showing predicted flight time for each combination of width and length. In optimization, this is known as the response surface. Brighter colors indicate longer predicted flight times. The average flight times for the combinations you've tested are shows as orange circles. The size of the circle shows the average of your observed flights for that combination.
```{r fig.height=5}
ggplot(to_predict,aes(rotor_length,rotor_width))+
geom_tile(aes(fill=predicted_flight_time))+
geom_point(aes(size=avg_flight), data=heli_means, color="orange")+
geom_text(aes(label=heli_id), data=heli_means,hjust=0,vjust=1)
```
Based on this model and the data you've gathered so far, the optimal helicopter has a rotor width of `r max_row$rotor_width`, a rotor length of `r max_row$rotor_length`, and will fly for `r max_row$predicted_flight_time` seconds.
Under this simple model, the predicted response surface will always be a plane, and the model will always predict that the optimal helicopter is in one corner of the space. This can't be right, and is an artifact of the too-simple mathematical form of the model. Still, this model does have some value. In particular, suggests that `r if (model$coefficients[2] > model$coefficients[3]) { "length" } else { "width" }` is the more important variable since it has a larger coefficient.
A better model for this problem is:
$F(l,w) \sim \beta_1 l + \beta_2 w + \beta_3 l^2 + \beta_4 w^2$
```{r}
model2<-lm(flight_time ~ rotor_length + rotor_width + I(rotor_length^2) + I(rotor_width^2), data=helidata)
summary(model2)
to_predict2<-expand.grid(rotor_width=seq(0.5,6,0.5),
rotor_length=seq(0.5,9.5,0.5))
to_predict2$predicted_flight_time<-predict(model2, to_predict2)
max_row2<-to_predict2[which.max(to_predict2$predicted_flight_time),]
```
Fitting this model requires more data, sampled from more combinations of length and width (which is why I started with the simpler model). If you haven't collected much data yet, you may get NA coefficients from the fit, you may get a warning about a "rank-deficient fit", and the model predictions might be non-sensical (like the optimal helicopter has the skinniest, stubbiest possible wings). If this is happening, just get more data and things will probably start to make sense.
This model allows for some curvature of the response surface. The best-fit response surface for your current data is plotted below.
```{r fig.height=5}
ggplot(to_predict2,aes(rotor_length,rotor_width))+
geom_tile(aes(fill=predicted_flight_time))+
geom_point(aes(size=avg_flight), data=heli_means, color="orange")+
geom_text(aes(label=heli_id), data=heli_means,hjust=0,vjust=1)
```
Based on this model and the data you've gathered so far, the optimal helicopter has a rotor width of `r max_row2$rotor_width`, a rotor length of `r max_row2$rotor_length`, and will fly for `r max_row2$predicted_flight_time` seconds.