daijiajie1@1: cat("\n____________________________\n")
daijiajie1@1: cat("\n____ MAIN PAPER RESULTS ____\n")
daijiajie1@1: cat("\n____________________________\n\n")
daijiajie1@1: cat("This script assumes that the script data_import.R has been run.\n")
daijiajie1@1: cat("This script generates statistics and figures used in the JASA submission.\n\n")
daijiajie1@1: 
daijiajie1@1: ##############################################################################
daijiajie1@1: # METHOD SECTION OUTPUT
daijiajie1@1: ##############################################################################
daijiajie1@1: 
daijiajie1@1: cat("\n____ METHOD SECTION ____\n\n")
daijiajie1@1: cat("- data for self-assessment table:\n")
daijiajie1@1: print(with(meta[meta$instudy,], summary(data.frame(age, musical.background, sing.ability, sing.experience, choir.experience))))
daijiajie1@1: 
daijiajie1@1: cat("\n____ ACCURACY SECTION ____\n\n")
daijiajie1@1: 
daijiajie1@1: cat("- Interval error significance tests:\n\n")
daijiajie1@1: significance.level <- 0.0001
daijiajie1@1: res <- with(data[data$ispoor == F, ],
daijiajie1@1:     t.test(intervalerror[ expectedinterval == 8 ]))
daijiajie1@1: cat(        "  - Upward minor sixth is flat.\n")
daijiajie1@1: cat(sprintf("    Is the difference significant at significance level %0.4f? %s (p = %g, df=%g, t = %g)\n\n", significance.level, res$p.value < significance.level, res$p.value, res$parameter, res$statistic))
daijiajie1@1: 
daijiajie1@1: res <- with(data[data$ispoor == F, ],
daijiajie1@1:     t.test(intervalerror[ expectedinterval == 12 ]))
daijiajie1@1: cat(        "  - Upward octave is flat.\n")
daijiajie1@1: cat(sprintf("    Is the difference significant at significance level %0.4f? %s (p = %g, df=%g, t = %g)\n\n", significance.level, res$p.value < significance.level, res$p.value, res$parameter, res$statistic))
daijiajie1@1: 
daijiajie1@1: res <- with(data[data$ispoor == F, ],
daijiajie1@1:     t.test(intervalerror[ expectedinterval == 0 ]))
daijiajie1@1: cat(sprintf("  - Prime (0 interval) is sharp by %0.2f.\n", res$estimate))
daijiajie1@1: cat(sprintf("    Is the difference significant at significance level %0.4f? %s (p = %g, df=%g, t = %g)\n\n", significance.level, res$p.value < significance.level, res$p.value, res$parameter, res$statistic))
daijiajie1@1: 
daijiajie1@1: cat("- explaining pitch error by following interval or beginning of phrase\n\n")
daijiajie1@1: is.beginning.of.phrase <- I(data$noteid[!data$ispoor] %in% c(1,7,13,20))
daijiajie1@1: res <- with(data[!data$ispoor, ],
daijiajie1@1:     summary(lm(pitcherror ~ is.beginning.of.phrase + expectedfutureinterval)))
daijiajie1@1: cat(sprintf("  - p: %0.3g, %0.3g\n",res$coefficients[3,4],res$coefficients[2,4]))
daijiajie1@1: cat(sprintf("  - F = %g\n",res$fstatistic[1]))
daijiajie1@1: cat(sprintf("  - DF: %0.0f\n",res$fstatistic[3]))
daijiajie1@1: cat(sprintf("  - Variance explained: %0.1f%s\n",100*res$r.squared,"%"))
daijiajie1@1: cat(sprintf("  - Sharpening in cents per note in interval: %0.1f\n\n",100*res$coefficients[3,1]))
daijiajie1@1: 
daijiajie1@1: res <- with(data[!data$ispoor, ],
daijiajie1@1:     summary(lm(pitcherror ~ is.beginning.of.phrase)))
daijiajie1@1: cat(sprintf("  - Variance explained only by intercept and following interval: %0.1f%s\n",100*res$r.squared,"%"))
daijiajie1@1: cat(sprintf("  - p: %0.3g\n",res$coefficients[2,4]))
daijiajie1@1: cat(sprintf("  - F = %g\n",res$fstatistic[1]))
daijiajie1@1: cat(sprintf("  - DF: %0.0f\n\n",res$fstatistic[3]))
daijiajie1@1: # cat(sprintf("  - Variance explained: %0.1f%s\n\n",100*res$r.squared,"%"))
daijiajie1@1: 
daijiajie1@1: res <- with(data[!data$ispoor, ],
daijiajie1@1:     summary(lm(pitcherror ~ expectedfutureinterval)))
daijiajie1@1: cat(sprintf("  - Variance explained only by intercept and beginning of phrase: %0.1f%s\n",100*res$r.squared,"%"))
daijiajie1@1: cat(sprintf("  - p: %0.3g\n",res$coefficients[2,4]))
daijiajie1@1: cat(sprintf("  - F = %g\n",res$fstatistic[1]))
daijiajie1@1: cat(sprintf("  - DF: %0.0f\n\n",res$fstatistic[3]))
daijiajie1@1: # cat(sprintf("  - Variance explained: %0.1f%s\n\n",100*res$r.squared,"%"))
daijiajie1@1: 
daijiajie1@1: et.ji.test <- with(data[!data$ispoor,], t.test(abs(pitcherror), abs(jipitcherror), paired=T))
daijiajie1@1: etMAPE <- with(data[!data$ispoor,], mean(abs(pitcherror), na.rm=T))
daijiajie1@1: jiMAPE <- with(data[!data$ispoor,], mean(abs(jipitcherror), na.rm=T))
daijiajie1@1: cat("- Error difference between ET and JI error:\n")
daijiajie1@1: cat("  Is it significant? Paired t-test says:")
daijiajie1@1: cat(sprintf("%s (p = %g, df: %0.0f, t = %g)\n\n", et.ji.test$p.value < 0.01, et.ji.test$p.value, et.ji.test$parameter, et.ji.test$statistic), "\n")
daijiajie1@1: cat(sprintf("ET (MAPE: %g) fits better than JI (MAPE: %g), by %g\n\n", etMAPE, jiMAPE, et.ji.test$estimate))
daijiajie1@1: 
daijiajie1@1: ##############################################################################
daijiajie1@1: # RESULTS SECTION OUTPUT
daijiajie1@1: ##############################################################################
daijiajie1@1: cat("\n____ RESULTS SECTION ____\n\n")
daijiajie1@1: cat("Basic metrics statistics.\n")
daijiajie1@1: cat(sprintf("  - MAPE. Mean: %0.3f;  median: %0.3f;  sd: %0.3f\n", mean(agg.mape$mape), median(agg.mape$mape), sd(agg.mape$mape)))
daijiajie1@1: cat(sprintf("  - MAIE. Mean: %0.3f;  median: %0.3f;  sd: %0.3f\n", mean(agg.maie$maie), median(agg.maie$maie), sd(agg.maie$maie)))
daijiajie1@1: # cat(sprintf("  - IB. Mean: %0.3f;  median: %0.3f;  sd: %0.3f\n", mean(agg.ib$ib), median(agg.ib$ib), sd(agg.ib$ib)))
daijiajie1@1: cat(sprintf("  - D. Mean: %0.3f;  median: %0.3f;  sd: %0.3f\n", mean(agg.drifts$diffs), median(agg.drifts$diffs), sd(agg.drifts$diffs)))
daijiajie1@1: cat(sprintf("  - |D|. Mean: %0.3f;  median: %0.3f;  sd: %0.3f\n\n", mean(abs(agg.drifts$diffs)), median(abs(agg.drifts$diffs)), sd(abs(agg.drifts$diffs))))
daijiajie1@1: cat(sprintf("  - lin D. in cents: Mean: %0.3f;  median: %0.3f;  sd: %0.3f\n", 100*mean(linear.drift$lineardrift), 100*median(linear.drift$lineardrift), 100*sd(linear.drift$lineardrift)))
daijiajie1@1: 
daijiajie1@1: cat(sprintf("- Significance of drift:"))
daijiajie1@1: cat(sprintf("  - significant drifts: %i of %i (%0.0f%s)\n", sum(agg.drifts.p$p<0.01), length(agg.drifts.p$p), 100*mean(agg.drifts.p$p<0.01), "%"))
daijiajie1@1: with(agg.drifts.p[agg.drifts$diffs>0,], cat(sprintf("  - number of upward drifts: %i, of which %i significantly so.\n", length(p), sum(p<0.01))))
daijiajie1@1: with(agg.drifts.p[agg.drifts$diffs<0,], cat(sprintf("  - number of downward drifts: %i, of which %i significantly so.\n\n", length(p), sum(p<0.01))))
daijiajie1@1: 
daijiajie1@1: cat(sprintf("- Significance of linear drift"))
daijiajie1@1: cat(sprintf("  - significant drifts: %i of %i (%0.0f%%) at conf. level 0.01, %i of %i (%0.0f%%) at conf. level 0.05 \n", sum(linear.drift$p.value<0.01), length(linear.drift$p.value), 100*mean(linear.drift$p.value<0.01), sum(linear.drift$p.value<0.05), length(linear.drift$p.value), 100*mean(linear.drift$p.value<0.05)))
daijiajie1@1: with(linear.drift[linear.drift$lineardrift>0,], cat(sprintf("  - number of upward drifts: %i, of which %i significantly so at 0.01.\n", length(p.value), sum(p.value<0.01))))
daijiajie1@1: with(linear.drift[linear.drift$lineardrift<0,], cat(sprintf("  - number of downward drifts: %i, of which %i significantly so at 0.01.\n\n", length(p.value), sum(p.value<0.01))))
daijiajie1@1: with(linear.drift[linear.drift$lineardrift>0,], cat(sprintf("  - number of upward drifts: %i, of which %i significantly so at 0.05.\n", length(p.value), sum(p.value<0.05))))
daijiajie1@1: with(linear.drift[linear.drift$lineardrift<0,], cat(sprintf("  - number of downward drifts: %i, of which %i significantly so at 0.05.\n\n", length(p.value), sum(p.value<0.05))))
daijiajie1@1: 
daijiajie1@1: 
daijiajie1@1: cat("- Correlation table (LaTeX syntax)\n")
daijiajie1@1: # cor.df <- data.frame(sg.abl=as.numeric(merged.agg$sing.ability), sg.exp=as.numeric(merged.agg$sing.experience), mus.bg=as.numeric(merged.agg$musical.background), ch.exp=as.numeric(merged.agg$choir.experience), MAIE=merged.agg$maie, MAPE=merged.agg$mape, IB=merged.agg$ib, D=merged.agg$diffs, abs.D=abs(merged.agg$diffs))
daijiajie1@1: cor.df <- data.frame(sg.abl=as.numeric(merged.agg$sing.ability), sg.exp=as.numeric(merged.agg$sing.experience), mus.bg=as.numeric(merged.agg$musical.background), ch.exp=as.numeric(merged.agg$choir.experience), MAIE=merged.agg$maie, MAPE=merged.agg$mape, lin.D=merged.agg$lineardrift, abs.D=abs(merged.agg$diffs), D=merged.agg$diffs)
daijiajie1@1: # cortable(data.frame(sg.abl=as.numeric(merged.agg$sing.ability), mus.bg=as.numeric(merged.agg$musical.background), ch.exp=as.numeric(merged.agg$choir.experience), MAIE=merged.agg$maie, ib=merged.agg$ib, MAPE=merged.agg$mape, D=merged.agg$diffs, abs.D=abs(merged.agg$diffs)))
daijiajie1@1: cortable(cor.df)
daijiajie1@1: cat("\n\n")
daijiajie1@1: 
daijiajie1@1: # testing singer predicting MAPE, MAIE, and drift
daijiajie1@1: cat("- Variance explained by singer for different measures\n\n")
daijiajie1@1: res <- summary(lm(mape~singer, data=agg.mape))
daijiajie1@1: pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
daijiajie1@1: df <- res$fstatistic[2]
daijiajie1@1: cat(sprintf("  - MAPE: %0.2f%%, p-value: %0.5g, df = %g, F = %g\n", res$r.squared*100, pval, df, res$fstatistic[1]))
daijiajie1@1: res <- summary(lm(maie~singer, data=agg.maie))
daijiajie1@1: pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
daijiajie1@1: df <- res$fstatistic[2]
daijiajie1@1: cat(sprintf("  - MAIE: %0.2f%%, p-value: %0.5g, df = %g, F = %g\n", res$r.squared * 100, pval, df, res$fstatistic[1]))
daijiajie1@1: # res <- summary(lm(ib~singer, data=agg.ib))
daijiajie1@1: # pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
daijiajie1@1: # cat(sprintf("  - IB: %0.2f%%, p-value: %0.5g\n", res$r.squared * 100, pval))
daijiajie1@1: res <- summary(lm(lineardrift~singer, data=linear.drift))
daijiajie1@1: pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
daijiajie1@1: df <- res$fstatistic[2]
daijiajie1@1: cat(sprintf("  - lin.D: %0.2f%%, p-value: %0.5g, df = %g, F = %g\n", res$r.squared * 100, pval, df, res$fstatistic[1]))
daijiajie1@1: res <- summary(lm(diffs~singer, data=agg.drifts))
daijiajie1@1: pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
daijiajie1@1: df <- res$fstatistic[2]
daijiajie1@1: cat(sprintf("  - D: %0.2f%%, p-value: %0.5g, df = %g, F = %g\n", res$r.squared * 100, pval, df, res$fstatistic[1]))
daijiajie1@1: res <- summary(lm(abs(diffs)~singer, data=agg.drifts))
daijiajie1@1: pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
daijiajie1@1: df <- res$fstatistic[2]
daijiajie1@1: cat(sprintf("  - |D|: %0.2f%%, p-value: %0.5g, df = %g, F = %g\n\n", res$r.squared * 100, pval, df, res$fstatistic[1]))
daijiajie1@1: 
daijiajie1@1: 
daijiajie1@1: # testing influence of condition
daijiajie1@1: cat("- Testing influence of singing condition (Normal/Masked/Imagined) via Kruskal-Wallis tests (sign. level 0.01):\n\n")
daijiajie1@1: res <- kruskal.test(mape ~ mode, data = agg.mape)
daijiajie1@1: cat(sprintf("  - MAPE significant difference due to condition? %s (p-value: %0.5g, df = %g, chisquared = %g)\n", res$p.value<0.01, res$p.value, res$parameter, res$statistic))
daijiajie1@1: res <- kruskal.test(maie ~ mode, data = agg.maie)
daijiajie1@1: cat(sprintf("  - MAIE significant difference due to condition? %s (p-value: %0.5g, df = %g, chisquared = %g)\n", res$p.value<0.01, res$p.value, res$parameter, res$statistic))
daijiajie1@1: # res <- kruskal.test(ib ~ mode, data = agg.ib)
daijiajie1@1: # cat(sprintf("  - IB significant difference due to condition? %s (p-value: %0.5g)\n", res$p.value<0.01, res$p.value))
daijiajie1@1: res <- kruskal.test(lineardrift ~ mode, data = linear.drift)
daijiajie1@1: cat(sprintf("  - lin. D significant difference due to condition? %s (p-value: %0.5g, df = %g, chisquared = %g)\n", res$p.value<0.01, res$p.value, res$parameter, res$statistic))
daijiajie1@1: res <- kruskal.test(diffs ~ mode, data = agg.drifts)
daijiajie1@1: cat(sprintf("  - D significant difference due to condition? %s (p-value: %0.5g, df = %g, chisquared = %g)\n", res$p.value<0.01, res$p.value, res$parameter, res$statistic))
daijiajie1@1: res <- kruskal.test(abs(diffs) ~ mode, data = agg.drifts)
daijiajie1@1: cat(sprintf("  - |D| significant difference due to condition? %s (p-value: %0.5g, df = %g, chisquared = %g)\n", res$p.value<0.01, res$p.value, res$parameter, res$statistic))
daijiajie1@1: 
daijiajie1@1: # cat("- \"Flat\" singer: \n\n")
daijiajie1@1: # "hl"
daijiajie1@1: # with(agg.drifts.all, cat(sprintf("  - drift in normal condition: %0.2f semitones\n", diffs[singer=="hl" & mode=="normal"])))
daijiajie1@1: # with(linear.drift.all, cat(sprintf("  - linear drift in normal condition: %0.2f cents per note\n", 100 * lineardrift[singer=="hl" & mode=="normal"])))
daijiajie1@1: # with(agg.drifts.all, cat(sprintf("  - drift in masked condition: %0.2f semitones\n", diffs[singer=="hl" & mode=="masked"])))
daijiajie1@1: # with(linear.drift.all, cat(sprintf("  - linear drift in masked condition: %0.2f cents per note\n", 100 * lineardrift[singer=="hl" & mode=="masked"])))
daijiajie1@1: # with(agg.drifts.all, cat(sprintf("  - drift in imagined condition: %0.2f semitones\n", diffs[singer=="hl" & mode=="imagined"])))
daijiajie1@1: # with(linear.drift.all, cat(sprintf("  - linear drift in imagined condition: %0.2f cents per note\n\n", 100 * lineardrift[singer=="hl" & mode=="imagined"])))
daijiajie1@1: 
daijiajie1@1: ##############################################################################
daijiajie1@1: # MODEL SECTION OUTPUT
daijiajie1@1: ##############################################################################
daijiajie1@1: cat("\n____ MODEL SECTION ____\n\n")
daijiajie1@1: 
daijiajie1@1: # notevar <- with(mdr[!mdr$ispoor,], aggregate(I(p-s-tstart), list(notenumber=notenumber), function(x)var(x,na.rm=T))[,2])
daijiajie1@1: # notevar2 <- with(data[!data$ispoor,], aggregate(pitcherror, list(notenumber=notenumber), function(x)var(x,na.rm=T))[,2])
daijiajie1@1: # driftvar <- notevar-notevar2
daijiajie1@1: # driftvar.fit <- lm(driftvar[7:75]~I(7:75))
daijiajie1@1: # sqrt(driftvar.fit$coefficients[2]*50)
daijiajie1@1: 
daijiajie1@1: cat("- existence of 'spread'\n")
daijiajie1@1: var.by.notenumber <- aggregate(data$relativetonic[!data$ispoor], list(data$notenumber[!data$ispoor]), function(x) {var(x, na.rm=T)})
daijiajie1@1: var2.by.notenumber <- aggregate(data$pitcherror[!data$ispoor], list(data$notenumber[!data$ispoor]), function(x) {var(x, na.rm=T)})
daijiajie1@1: counts.by.notenumber <- aggregate(data$relativetonic[!data$ispoor], list(data$notenumber[!data$ispoor]), function(x) {sum(!is.na(x))})
daijiajie1@1: dist.by.notenumber <- aggregate(abs(data$relativetonic[!data$ispoor]), list(data$notenumber[!data$ispoor]), function(x) {mean(x, na.rm=T)})
daijiajie1@1: fit1 <- lm(var.by.notenumber[,2]~var2.by.notenumber[,2], weights=counts.by.notenumber[,2])
daijiajie1@1: fit2 <- lm(fit1$residuals~var.by.notenumber[(burnin+1):75,1], weights=counts.by.notenumber[(burnin+1):75,2])
daijiajie1@1: summary(fit2)
daijiajie1@1: res <- summary(fit2)
daijiajie1@1: pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
daijiajie1@1: cat(sprintf("  - notenumber: %0.2f%%, p-value: %0.5g, df = %g, F = %g\n\n", res$r.squared * 100, pval, res$fstatistic[1], res$fstatistic[3]))
daijiajie1@1: plot(names(fit1$residuals), fit1$residuals)
daijiajie1@1: abline(fit2)
daijiajie1@1: 
daijiajie1@1: intvar <- var(data$intervalerror[!data$ispoor], na.rm=T)
daijiajie1@1: expected.D.var <- 50 * intvar
daijiajie1@1: expected.D.sd <- expected.D.var^0.5
daijiajie1@1: pitch.err.sd <- sd(data$pitcherror[!data$ispoor], na.rm=T)
daijiajie1@1: cat(sprintf("- inadequacy of mu = 0 model:\n\n"))
daijiajie1@1: cat(sprintf("  - interval (error) variance: %0.3f\n", intvar))
daijiajie1@1: cat(sprintf("  - expected variance of D (interval variance summed over 50 notes): %0.3f\n", expected.D.var))
daijiajie1@1: cat(sprintf("  - ... and the corresponding standard deviation: %0.3f\n\n", expected.D.sd))
daijiajie1@1: cat(sprintf("  - In comparison the standard deviation of pitch erros is: %0.3f\n\n", pitch.err.sd))
daijiajie1@1: 
daijiajie1@1: model.mape <- apply(md[!md$ispoor,(1:nmu)+8], 2, function(x) {mean(abs(x),na.rm=T)})
daijiajie1@1: model.mape.singerwise <- aggregate(md[!md$ispoor,(1:nmu)+8], list(md$singer[!md$ispoor]), function(x) {mean(abs(x),na.rm=T)})
daijiajie1@1: model.mape.recordingwise <- aggregate(md[!md$ispoor,(1:nmu)+8], list(md$singer[!md$ispoor],md$mode[!md$ispoor]), function(x) {mean(abs(x),na.rm=T)})
daijiajie1@1: 
daijiajie1@1: recording.mus <- data.frame(singer = model.mape.recordingwise[,1], mu = mus[apply(model.mape.recordingwise[,-(1:2)], 1, which.min)])
daijiajie1@1: # singer.mus <- data.frame(singer = model.mape.singerwise[,1], mu = mus[apply(model.mape.singerwise[,-1], 1, which.min)])
daijiajie1@1: singer.mus <- aggregate(recording.mus[,2], list(recording.mus[,1]),mean)
daijiajie1@1: 
daijiajie1@1: cat(sprintf("- Better parameter for mu, firstly via 'model MAPE':\n\n"))
daijiajie1@1: cat(sprintf("  - the smallest error of %0.2f is achieved with model '%s'\n", min(model.mape), names(which.min(model.mape))))
daijiajie1@1: cat(sprintf("  - error at mu = 0: %0.2f\n", model.mape[1]))
daijiajie1@1: cat(sprintf("  - error at mu = 1: %0.2f\n\n", model.mape[length(model.mape)]))
daijiajie1@1: 
daijiajie1@1: #-------------------------------
daijiajie1@1: mdr <- merge(md, mr[,c(1,2,5)], by=c("singer", "mode"))
daijiajie1@1: # plot(I(p-s-tstart)~factor(notenumber), data=mdr)
daijiajie1@1: 
daijiajie1@1: sig2 <- var(data$pitcherror[!data$ispoor], na.rm=T)
daijiajie1@1: sig2int <- var(data$intervalerror[!data$ispoor], na.rm=T)
daijiajie1@1: n <- 50
daijiajie1@1: # vd <- var(drifts[!drifts$singer %in% poor.singers,4], na.rm=T)
daijiajie1@1: # vd <- var(drifts$diffs[!drifts$singer %in% poor.singers], na.rm=T)
daijiajie1@1: vd <- var(agg.drifts$diffs)
daijiajie1@1: est.mu <- 1-sqrt(vd/(n*sig2))
daijiajie1@1:   
daijiajie1@1: cat("- Estimate according to drift variance formula:\n\n")
daijiajie1@1: cat(sprintf("  - overall pitch error variance: %0.3f\n", sig2))
daijiajie1@1: # cat(sprintf("  - overall interval error variance: %0.3f", sig2int))
daijiajie1@1: cat(sprintf("  - observed variance over 50-note gap: %0.3f\n", vd))
daijiajie1@1: cat(sprintf("  - estimated mu based on these: %0.3f\n\n", est.mu))
daijiajie1@1: 
daijiajie1@1: cat("- Stats for mu by singer:\n")
daijiajie1@1: cat(sprintf("  - range: %0.3f to %0.3f\n", range(singer.mus[,2])[1], range(singer.mus[,2])[2]))
daijiajie1@1: cat(sprintf("  - mean: %0.3f, median: %0.3f, sd: %0.3f\n\n", mean(singer.mus[,2]), median(singer.mus[,2]), sd(singer.mus[,2])))
daijiajie1@1: 
daijiajie1@1: 
daijiajie1@1: # n <- 50
daijiajie1@1: # mu <- 0.916
daijiajie1@1: # sig2 <- var(data$pitcherror[ !data$ispoor ], na.rm=T)
daijiajie1@1: # vv <- var(drifts$diffs[!drifts$singer %in% poor.singers], na.rm=T)
daijiajie1@1: # vv <- var(drifts$diffs[!drifts$singer %in% poor.singers & drifts$noteid ==iNote], na.rm=T)
daijiajie1@1: 
daijiajie1@1: # n <- 50
daijiajie1@1: # sig2 <- var(data$pitcherror[ !data$ispoor ], na.rm=T)
daijiajie1@1: # vv <- var(agg.drifts$diffs)
daijiajie1@1: # est.mu <- 1-sqrt(vv/(n*sig2))
daijiajie1@1: # 
daijiajie1@1: # vv <- var(agg.drifts$diffs)
daijiajie1@1: # for (mu in c(0,0.8,0.9,est.mu,  0.95, 1))
daijiajie1@1: # {
daijiajie1@1: #     dvar <- (1-mu)^2 * n * sig2
daijiajie1@1: #     cat(sprintf("Example: mu = %g, sigma^2 = %g, n = %g\n", mu, sig2, n))
daijiajie1@1: #     cat(sprintf(" ... var = %g, sd = %g\n", dvar, sqrt(dvar)))
daijiajie1@1: # }
daijiajie1@1: # CHECK!!!!!
daijiajie1@1: 
daijiajie1@1: # singermodemu <- c()
daijiajie1@1: # for (singer in singers[!singers %in% poor.singers])
daijiajie1@1: # {
daijiajie1@1: #     for (mode in c("normal","imagined","masked"))
daijiajie1@1: #     {
daijiajie1@1: #         sig2 <- var(data$pitcherror[drifts$singer==singer & drifts$mode==mode], na.rm=T)
daijiajie1@1: #         # sig2int <- var(data$intervalerror[data$singermode==singermode], na.rm=T)
daijiajie1@1: #         n <- 50
daijiajie1@1: #         vd <- var(drifts$diffs[drifts$singer==singer & drifts$mode==mode], na.rm=T)
daijiajie1@1: #         est.mu <- 1-sqrt(vd/(n*sig2))
daijiajie1@1: #         # cat(sprintf("%s: %0.3f\n", singer, est.mu))
daijiajie1@1: #         singermode <- sprintf("%s%s", singer, mode)
daijiajie1@1: #         singermodemu <- rbind(singermodemu, data.frame(singermode=singermode, singer=singer, mode=mode, mu=est.mu))
daijiajie1@1: #     }
daijiajie1@1: # }
daijiajie1@1: # hist(singermodemu$mu, col="gray")
daijiajie1@1: 
daijiajie1@1: # make drift figures
daijiajie1@1: for (singermode in sort(unique(data$singermode)))
daijiajie1@1: {
daijiajie1@1: 	pdf(sprintf("../figures/drift_figure_%s.pdf", singermode), width = 3.5, height = 3.5)
daijiajie1@1: 	makeDriftGraph(data, singermode)
daijiajie1@1: 	dev.off()
daijiajie1@1: }
daijiajie1@1: 
daijiajie1@1: for (singermode in sort(unique(data$singermode)))
daijiajie1@1: {
daijiajie1@1: 	pdf(sprintf("../figures/drift_figure_model0850_%s.pdf", singermode), width = 5.5, height = 4.5)
daijiajie1@1: 	with(data[data$singermode==singermode,], plot(notenumber, newmedian-semitone, pch="+", ylab="pitch", xlab="note"))
daijiajie1@1: 	with(md[md$singermode==singermode,], points(notenumber, p-s-err0850, pch=20, type='o'))
daijiajie1@1: 	dev.off()
daijiajie1@1: }
daijiajie1@1: 
daijiajie1@1: countsByInterval <- with(data[!data$ispoor,], aggregate(data.frame(count=data$intervalerror), list(interval=data$expectedinterval), function(x){sum(!is.na(x))})) 
daijiajie1@1: 
daijiajie1@1: pdf("../figures/interval_error_by_interval_boxplot.pdf", width = 5.5, height = 4)
daijiajie1@1: par(mar=c(2,4,1,.5))
daijiajie1@1: plot(c(1,18), c(-1.6,1.6), type='n', xaxt='n', xlab = '', ylab='', yaxt='n')
daijiajie1@1: abline(h=0)
daijiajie1@1: abline(h=-1, lty=2)
daijiajie1@1: abline(h=1, lty=2)
daijiajie1@1: plot(intervalerror ~ factor(expectedinterval, levels=-5:12), data=data, subset = !ispoor,
daijiajie1@1: 	pch=20, col=rgb(0.9,0.9,0.9), ylim = c(-2,2), outline=F, xaxt='n', yaxt='n', add=T,
daijiajie1@1: 	ylab="interval error (with respect to ET)", xlab="interval in semitones",
daijiajie1@1: 	cex.axis=0.8, lty=1)
daijiajie1@1: posi <- c(1:6, 8, 11, 13, 14, 18)
daijiajie1@1: axis(2, at=c(-1,0,1), cex.axis=0.8)
daijiajie1@1: axis(1, at=posi, labels=posi-6, cex.axis=0.8)
daijiajie1@1: # text(counts$interval + 6, rep(2,11), counts$count, cex=0.7)
daijiajie1@1: mtext(c("n=", countsByInterval$count), at=c(-.1,countsByInterval$interval + 6), cex=0.7)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: countsByNoteID <- with(data[!data$ispoor,], aggregate(data.frame(count=data$intervalerror), list(noteid=data$noteid), function(x){sum(!is.na(x))}))
daijiajie1@1: 
daijiajie1@1: pdf("../figures/interval_error_by_notenumber_boxplot.pdf", width = 7, height = 4)
daijiajie1@1: par(mar=c(2,.5,1,1))	
daijiajie1@1: plot(c(1,25), c(-1.6,1.6), type='n', xaxt='n', xlab = '', ylab='', yaxt='n')
daijiajie1@1: abline(h=0)
daijiajie1@1: abline(h=-1, lty=2)
daijiajie1@1: abline(h=1, lty=2)
daijiajie1@1: plot(intervalerror ~ factor(noteid, levels=1:75), data=data, subset = !ispoor,
daijiajie1@1: 	pch=20, col=rgb(0.9,0.9,0.9), ylim = c(-2,2), outline=F, add=T,
daijiajie1@1: 	ylab="interval error (with respect to ET)", xlab="interval in semitones",
daijiajie1@1: 	cex.axis=0.8, lty=1, yaxt='n', boxwex=0.6)
daijiajie1@1: # posi <- c(1:6, 8, 11, 13, 14, 18)
daijiajie1@1: # axis(1, at=posi, labels=posi-6, cex.axis=0.8)
daijiajie1@1: # text(countsByNoteID$noteid, rep(2,11), countsByNoteID$count, cex=0.7)
daijiajie1@1: mtext(c(countsByNoteID$count), at=c(countsByNoteID$noteid), cex=0.7)
daijiajie1@1: text(1:25, rep(-1.5,11), c("(-5)", gtdata$interval_ST[-25]), cex=0.7)
daijiajie1@1: text(2.7, -1.35 , "intervals in semitones", cex=0.7)
daijiajie1@1: text(1-.5, 1.3, "phrase 1", pos=4)
daijiajie1@1: abline(v=6.5, lty=1, col='gray')
daijiajie1@1: text(7-.5, 1.3, "phrase 2", pos=4)
daijiajie1@1: abline(v=12.5, lty=1, col='gray')
daijiajie1@1: text(13-.5, 1.3, "phrase 3", pos=4)
daijiajie1@1: abline(v=19.5, lty=1, col='gray')
daijiajie1@1: text(20-.5, 1.3, "phrase 4", pos=4)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: pdf("../figures/pitch_error_by_interval_boxplot.pdf", width = 5.5, height = 4)
daijiajie1@1: par(mar=c(2,4,1,.5))
daijiajie1@1: plot(c(1,18), c(-1.6,1.6), type='n', xaxt='n', xlab = '', ylab='', yaxt='n')
daijiajie1@1: abline(h=0)
daijiajie1@1: abline(h=-1, lty=2)
daijiajie1@1: abline(h=1, lty=2)
daijiajie1@1: plot(pitcherror ~ factor(expectedinterval, levels=-5:12), data=data, subset = !ispoor,
daijiajie1@1: 	pch=20, col=rgb(0.9,0.9,0.9), ylim = c(-2,2), outline=F, xaxt='n', yaxt='n', add=T,
daijiajie1@1: 	ylab="pitch error (with respect to ET)", xlab="interval in semitones",
daijiajie1@1: 	cex.axis=0.8, lty=1)
daijiajie1@1: posi <- c(1:6, 8, 11, 13, 14, 18)
daijiajie1@1: axis(2, at=c(-1,0,1), cex.axis=0.8)
daijiajie1@1: axis(1, at=posi, labels=posi-6, cex.axis=0.8)
daijiajie1@1: # text(counts$interval + 6, rep(2,11), counts$count, cex=0.7)
daijiajie1@1: # mtext(c("n=", countsByInterval$count), at=c(-.1,countsByInterval$interval + 6), cex=0.7)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: countsBySemitone <- with(data[!data$ispoor,], aggregate(data.frame(count=data$pitcherror), list(interval=data$semitone), function(x){sum(!is.na(x))})) 
daijiajie1@1: pdf("../figures/pitch_error_by_semitone_boxplot.pdf", width = 4.5, height = 4)
daijiajie1@1: par(mar=c(2,4,1,.5))
daijiajie1@1: plot(c(.5,13.5), c(-1.6,1.6), type='n', xaxt='n', xlab = '', ylab='', yaxt='n')
daijiajie1@1: abline(h=0)
daijiajie1@1: abline(h=-1, lty=2)
daijiajie1@1: abline(h=1, lty=2)
daijiajie1@1: plot(pitcherror ~ factor(semitone, levels=-5:7), data=data, subset = !ispoor,
daijiajie1@1: 	pch=20, col=rgb(0.9,0.9,0.9), ylim = c(-2,2), outline=F, xaxt='n', yaxt='n', add=T,
daijiajie1@1: 	ylab="pitch error (with respect to ET)", xlab="interval in semitones",
daijiajie1@1: 	cex.axis=0.8, lty=1)
daijiajie1@1: posi <- c(1, 3, 5, 6, 8, 10, 11, 13)
daijiajie1@1: axis(2, at=c(-1,0,1), cex.axis=0.8)
daijiajie1@1: axis(1, at=posi, labels=posi-6, cex.axis=0.8)
daijiajie1@1: # text(counts$interval + 6, rep(2,11), counts$count, cex=0.7)
daijiajie1@1: # mtext(c("n=", countsByInterval$count), at=c(-.1,countsByInterval$interval + 6), cex=0.7)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: pdf("../figures/pitch_error_by_notenumber_boxplot.pdf", width = 7, height = 4)
daijiajie1@1: par(mar=c(2,.5,1,1))	
daijiajie1@1: plot(c(1,25), c(-1.6,1.6), type='n', xaxt='n', xlab = '', ylab='', yaxt='n')
daijiajie1@1: abline(h=0)
daijiajie1@1: abline(h=-1, lty=2)
daijiajie1@1: abline(h=1, lty=2)
daijiajie1@1: plot(pitcherror ~ factor(noteid, levels=1:75), data=data, subset = !ispoor,
daijiajie1@1: 	pch=20, col=rgb(0.9,0.9,0.9), ylim = c(-2,2), outline=F, add=T,
daijiajie1@1: 	ylab="pitch error (with respect to ET)", xlab="interval in semitones",
daijiajie1@1: 	cex.axis=0.8, lty=1, yaxt='n', boxwex=0.6)
daijiajie1@1: # posi <- c(1:6, 8, 11, 13, 14, 18)
daijiajie1@1: # axis(1, at=posi, labels=posi-6, cex.axis=0.8)
daijiajie1@1: # text(countsByNoteID$noteid, rep(2,11), countsByNoteID$count, cex=0.7)
daijiajie1@1: # mtext(c(countsByNoteID$count), at=c(countsByNoteID$noteid), cex=0.7)
daijiajie1@1: text(1:25, rep(-1.5,11), gtdata$pitch_ion_cents_et/100, cex=0.7)
daijiajie1@1: text(2.7, -1.35 , "pitch relative to tonic", cex=0.7)
daijiajie1@1: text(1-.5, 1.3, "phrase 1", pos=4)
daijiajie1@1: abline(v=6.5, lty=1, col='gray')
daijiajie1@1: text(7-.5, 1.3, "phrase 2", pos=4)
daijiajie1@1: abline(v=12.5, lty=1, col='gray')
daijiajie1@1: text(13-.5, 1.3, "phrase 3", pos=4)
daijiajie1@1: abline(v=19.5, lty=1, col='gray')
daijiajie1@1: text(20-.5, 1.3, "phrase 4", pos=4)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: # linear fit example plot
daijiajie1@1: pdf("../figures/pitch_error_example.pdf", width = 6, height = 5)
daijiajie1@1: par(mar=c(4,4,.1,.1))
daijiajie1@1: ind <- which(data$singer == "jp" & data$mode == "imagined" & data$run == 1)
daijiajie1@1: plot(data$newmedian[ind], pch = 5,
daijiajie1@1: 	xlab="note number", ylab="pitch",
daijiajie1@1: 	xlim = c(0,25), las = 1)
daijiajie1@1: points((data$newmedian-data$semitone)[ind], pch=18)
daijiajie1@1: fit <- lm(I(newmedian-semitone)~notenumber, data = data[ind,])
daijiajie1@1: for (ii in 1:25)
daijiajie1@1: {
daijiajie1@1: 	lines(c(ii,ii), c(fit$fitted.values[ii], (data$newmedian-data$semitone)[ind[ii]]))
daijiajie1@1: }
daijiajie1@1: abline(a=fit$coefficients[1], b=fit$coefficients[2])
daijiajie1@1: text(0, fit$coefficients[1]+0.3, 0, cex = 0.7)
daijiajie1@1: for (offs in 1:8)
daijiajie1@1: {
daijiajie1@1: 	abline(a=fit$coefficients[1]+offs, b=fit$coefficients[2], lty = 2)
daijiajie1@1: 	abline(a=fit$coefficients[1]-offs, b=fit$coefficients[2], lty = 2)
daijiajie1@1: 	text(0, fit$coefficients[1]+offs+0.3, +offs, cex = 0.7)
daijiajie1@1: 	text(0, fit$coefficients[1]-offs+0.3, -offs, cex = 0.7)
daijiajie1@1: }
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: # drift histogram
daijiajie1@1: pdf("../figures/drift_histogram.pdf", width = 6, height = 4)
daijiajie1@1: par(mar=c(4,4,.3,.5))
daijiajie1@1: hist(agg.drifts$diffs, seq(-1.5,1.5,by=0.2), col="dark gray", border="black", 
daijiajie1@1: 	xlab = "drift", ylab = "number of recordings", main = "")
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: # mean absolute pitch error histogram
daijiajie1@1: pdf("../figures/mape_histogram.pdf", width = 3.3, height = 4)
daijiajie1@1: par(mar=c(4,4,.3,.5))
daijiajie1@1: hist(agg.mape$mape, seq(0,0.4,by=0.05), col="dark gray", border="black", 
daijiajie1@1: 	xlab = "semitones", ylab = "number of recordings", main = "",
daijiajie1@1: 	ylim=c(0,24))
daijiajie1@1: abline(h=5, col = "gray")
daijiajie1@1: abline(h=10, col = "gray")
daijiajie1@1: abline(h=15, col = "gray")
daijiajie1@1: abline(h=20, col = "gray")
daijiajie1@1: hist(agg.mape$mape, seq(0,0.4,by=0.05), col="dark gray", border="black", 
daijiajie1@1: 	xlab = "semitones", ylab = "number of recordings", main = "",
daijiajie1@1: 	ylim=c(0,24), add = T)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: # mean absolute interval error histogram
daijiajie1@1: pdf("../figures/maie_histogram.pdf", width = 3.3, height = 4)
daijiajie1@1: par(mar=c(4,1,.3,.5))
daijiajie1@1: hist(agg.maie$maie, seq(0,.5,by=0.05), col="dark gray", border="black", 
daijiajie1@1: 	xlab = "semitones", ylab = "number of recordings", main = "",
daijiajie1@1: 	ylim=c(0,24), yaxt="n")
daijiajie1@1: abline(h=5, col = "gray")
daijiajie1@1: abline(h=10, col = "gray")
daijiajie1@1: abline(h=15, col = "gray")
daijiajie1@1: abline(h=20, col = "gray")
daijiajie1@1: hist(agg.maie$maie, seq(0,.5,by=0.05), col="dark gray", border="black", 
daijiajie1@1: 	xlab = "semitones", ylab = "number of recordings", main = "",
daijiajie1@1: 	ylim=c(0,24), yaxt="n", add = T)
daijiajie1@1: dev.off()
daijiajie1@1: mean(agg.maie$maie)
daijiajie1@1: median(agg.maie$maie)
daijiajie1@1: sd(agg.maie$maie)
daijiajie1@1: 
daijiajie1@1: # MAR
daijiajie1@1: pdf("../figures/mar_histogram.pdf", width = 3.3, height = 4)
daijiajie1@1: par(mar=c(4,1,.3,.5))
daijiajie1@1: hist(agg.rise$rise, seq(0,.5,by=0.05), col="dark gray", border="black", 
daijiajie1@1: 	xlab = "semitones", ylab = "number of recordings", main = "",
daijiajie1@1: 	ylim=c(0,24), yaxt="n")
daijiajie1@1: abline(h=5, col = "gray")
daijiajie1@1: abline(h=10, col = "gray")
daijiajie1@1: abline(h=15, col = "gray")
daijiajie1@1: abline(h=20, col = "gray")
daijiajie1@1: hist(agg.rise$rise, seq(0,.5,by=0.05), col="dark gray", border="black", 
daijiajie1@1: 	xlab = "semitones", ylab = "number of recordings", main = "",
daijiajie1@1: 	ylim=c(0,24), yaxt="n", add = T)
daijiajie1@1: dev.off()
daijiajie1@1: mean(agg.rise$rise)
daijiajie1@1: median(agg.rise$rise)
daijiajie1@1: sd(agg.rise$rise)
daijiajie1@1: 
daijiajie1@1: # drift
daijiajie1@1: pdf("../figures/d_histogram.pdf", width = 3.3, height = 4)
daijiajie1@1: par(mar=c(4,1,.3,.5))
daijiajie1@1: hist(agg.drifts$diffs, seq(-.6,.6,by=0.1), col="dark gray", border="black", 
daijiajie1@1: 	xlab = "semitones", ylab = "number of recordings", main = "",
daijiajie1@1: 	ylim=c(0,24), yaxt="n")
daijiajie1@1: abline(h=5, col = "gray")
daijiajie1@1: abline(h=10, col = "gray")
daijiajie1@1: abline(h=15, col = "gray")
daijiajie1@1: abline(h=20, col = "gray")
daijiajie1@1: hist(agg.drifts$diffs, seq(-.6,.6,by=0.1), col="dark gray", border="black", 
daijiajie1@1: 	xlab = "semitones", ylab = "number of recordings", main = "",
daijiajie1@1: 	ylim=c(0,24), yaxt="n", add = T)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: # absolute drift
daijiajie1@1: pdf("../figures/abs_d_histogram.pdf", width = 3.3, height = 4)
daijiajie1@1: par(mar=c(4,1,.3,.5))
daijiajie1@1: hist(abs(agg.drifts$diffs), seq(0,.6,by=0.05), col="dark gray", border="black", 
daijiajie1@1: 	xlab = "semitones", ylab = "number of recordings", main = "",
daijiajie1@1: 	ylim=c(0,24), yaxt="n")
daijiajie1@1: abline(h=5, col = "gray")
daijiajie1@1: abline(h=10, col = "gray")
daijiajie1@1: abline(h=15, col = "gray")
daijiajie1@1: abline(h=20, col = "gray")
daijiajie1@1: hist(abs(agg.drifts$diffs), seq(0,.6,by=0.05), col="dark gray", border="black", 
daijiajie1@1: 	xlab = "semitones", ylab = "number of recordings", main = "",
daijiajie1@1: 	ylim=c(0,24), yaxt="n", add = T)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: # linear drift
daijiajie1@1: pdf("../figures/lineard_histogram.pdf", width = 3.3, height = 4)
daijiajie1@1: par(mar=c(4,1,.3,.5))
daijiajie1@1: hist(100*linear.drift$lineardrift, seq(-1.2,1.2,by=0.2), col="dark gray", border="black", 
daijiajie1@1: 	xlab = "cents", ylab = "number of recordings", main = "",
daijiajie1@1: 	ylim=c(0,24), yaxt="n")
daijiajie1@1: abline(h=5, col = "gray")
daijiajie1@1: abline(h=10, col = "gray")
daijiajie1@1: abline(h=15, col = "gray")
daijiajie1@1: abline(h=20, col = "gray")
daijiajie1@1: hist(100*linear.drift$lineardrift, seq(-1.2,1.2,by=0.2), col="dark gray", border="black", 
daijiajie1@1: 	xlab = "cents", ylab = "number of recordings", main = "",
daijiajie1@1: 	ylim=c(0,24), yaxt="n", add = T)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: # mu
daijiajie1@1: pdf("../figures/mu_histogram.pdf", width = 4.3, height = 3)
daijiajie1@1: par(mar=c(4,4,.3,.5))
daijiajie1@1: hist(singer.mus[,2], seq(0.6, 1.0, 0.4/5), col="dark gray", border="black", 
daijiajie1@1: 	xlab = expression(mu), ylab = "number of singers", main = "",
daijiajie1@1: 	add = F)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: 
daijiajie1@1: 
daijiajie1@1: # drift significance plot
daijiajie1@1: pdf("../figures/drift_significance_plot.pdf", width = 6, height = 3.5)
daijiajie1@1: plot(agg.drifts$diffs, agg.drifts.p$p, log="y", col = rgb(0,0,0,0.4), 
daijiajie1@1: 	pch=19, yaxt = "n", ylab = "p-value                  ", 
daijiajie1@1: 	xlab = "drift magnitude in semitones",
daijiajie1@1: 	xlim=c(-1,1))
daijiajie1@1: axis(2, at = 10^(-(1:4)), label = c("0.1","0.01","0.001","0.0001"),las=1)
daijiajie1@1: abline(h=0.01, lty = 2)
daijiajie1@1: abline(v=0.0)
daijiajie1@1: text(0.75,0.0000005,"significant\nupward\ndrift")
daijiajie1@1: text(-0.75,0.0000005,"significant\ndownward\ndrift")
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: # drift significance plot
daijiajie1@1: pdf("../figures/lineardrift_significance_plot.pdf", width = 6, height = 4.8)
daijiajie1@1: plot(linear.drift$lineardrift, linear.drift$p.value, log="y", col = rgb(0,0,0,0.4), 
daijiajie1@1: 	pch=19, yaxt = "n", ylab = "p-value                  ", 
daijiajie1@1: 	xlab = "linear drift per note in semitones")
daijiajie1@1: axis(2, at = 10^(-(1:4)), label = c("0.1","0.01","0.001","0.0001"),las=1)
daijiajie1@1: abline(h=0.01, lty = 2)
daijiajie1@1: abline(v=0.0)
daijiajie1@1: text(0.0035,0.0000005,"significant\nupward\ndrift")
daijiajie1@1: text(-0.0035,0.0000005,"significant\ndownward\ndrift")
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: 
daijiajie1@1: # how many are significantly adrift?
daijiajie1@1: # sum(agg.drifts.p$p<0.01)
daijiajie1@1: # mean(agg.drifts.p$p<0.01)
daijiajie1@1: 
daijiajie1@1: # drift vs mape
daijiajie1@1: pdf("../figures/drift_vs_mape.pdf", width = 5.5, height = 4)
daijiajie1@1: par(mar=c(4,4,.3,.5))
daijiajie1@1: outl <- agg.mape$mape>0.35
daijiajie1@1: fit1 <- lm(agg.mape$mape~agg.drifts$diffs)
daijiajie1@1: fit2 <- lm(agg.mape$mape~agg.drifts$diffs, subset=!outl)
daijiajie1@1: plot(agg.drifts$diffs, agg.mape$mape, type="n", 
daijiajie1@1: 	xlab="drift magnitude in semitones",
daijiajie1@1: 	ylab="mean absolute pitch error")
daijiajie1@1: points(agg.drifts$diffs[outl], agg.mape$mape[outl], pch=19)
daijiajie1@1: points(agg.drifts$diffs[!outl], agg.mape$mape[!outl], pch=21)
daijiajie1@1: abline(fit1, lty=2)
daijiajie1@1: abline(fit2)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: # drift vs maie
daijiajie1@1: pdf("../figures/drift_vs_maie.pdf", width = 5.5, height = 4)
daijiajie1@1: par(mar=c(4,4,.3,.5))
daijiajie1@1: outl <- agg.mape$mape>0.35
daijiajie1@1: fit1 <- lm(agg.maie$maie~agg.drifts$diffs)
daijiajie1@1: fit2 <- lm(agg.maie$maie~agg.drifts$diffs, subset=!outl)
daijiajie1@1: plot(agg.drifts$diffs, agg.maie$maie, type="n", 
daijiajie1@1: 	xlab="drift magnitude in semitones",
daijiajie1@1: 	ylab="mean absolute pitch error")
daijiajie1@1: points(agg.drifts$diffs[outl], agg.maie$maie[outl], pch=19)
daijiajie1@1: points(agg.drifts$diffs[!outl], agg.maie$maie[!outl], pch=21)
daijiajie1@1: abline(fit1, lty=2)
daijiajie1@1: abline(fit2)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: # drift vs IB
daijiajie1@1: pdf("../figures/drift_vs_ib.pdf", width = 5.5, height = 4)
daijiajie1@1: par(mar=c(4,4,.3,.5))
daijiajie1@1: fit1 <- lm(agg.ib$ib~agg.drifts$diffs)
daijiajie1@1: plot(agg.drifts$diffs, agg.ib$ib, type="n", 
daijiajie1@1: 	xlab="drift magnitude in semitones",
daijiajie1@1: 	ylab="mean interval error")
daijiajie1@1: points(agg.drifts$diffs, agg.ib$ib, pch=21)
daijiajie1@1: abline(fit1, lty=2)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: 
daijiajie1@1: # cor(merged.singeragg$maie, merged.singeragg$mape)
daijiajie1@1: 
daijiajie1@1: pdf("../figures/mape_maie_drift.pdf", width = 5, height = 5)
daijiajie1@1: # par(mar=c(1,1,1,1))
daijiajie1@1: plot(data.frame(MAPE=agg.mape$mape, MAIE=agg.maie$maie, Drift=agg.drifts$diffs))
daijiajie1@1: dev.off()
daijiajie1@1: # and their correlations
daijiajie1@1: cor.test(agg.mape$mape, agg.maie$maie)
daijiajie1@1: cor.test(agg.mape$mape, agg.drifts$diffs)
daijiajie1@1: cor.test(agg.maie$maie, agg.drifts$diffs)
daijiajie1@1: 
daijiajie1@1: # corrgram(data.frame(merged.singeragg$maie, merged.singeragg$mape, merged.singeragg$diffs, as.numeric(merged.singeragg$sing.ability), as.numeric(merged.singeragg$musical.background), as.numeric(merged.singeragg$choir.experience)))
daijiajie1@1: # corrgram(data.frame(merged.agg$maie, merged.agg$mape, merged.agg$diffs, as.numeric(merged.agg$sing.ability), as.numeric(merged.agg$musical.background), as.numeric(merged.agg$choir.experience)))
daijiajie1@1: 
daijiajie1@1: # pdf("../figures/corrgram.pdf", width = 5, height = 5)
daijiajie1@1: # corrgram(data.frame(MAIE=merged.agg$maie, ib=merged.agg$ib, MAPE=merged.agg$mape, D=merged.agg$diffs, singing.ability=as.numeric(merged.agg$sing.ability), musical.bg=as.numeric(merged.agg$musical.background), choir.exp=as.numeric(merged.agg$choir.experience)), lower.panel=panel.conf, upper.panel=panel.pie)
daijiajie1@1: # dev.off()
daijiajie1@1: 
daijiajie1@1: # slope figure
daijiajie1@1: pdf("../figures/rise_by_notenumber_boxplot.pdf", width = 5.5, height = 4)
daijiajie1@1: par(mar=c(2,3.8,.1,.1))	
daijiajie1@1: plot(c(1,25), c(-1.6,1.6), type='n', xaxt='n', xlab = '', ylab='', yaxt='n')
daijiajie1@1: abline(h=0)
daijiajie1@1: abline(h=-1, lty=2)
daijiajie1@1: abline(h=1, lty=2)
daijiajie1@1: plot(I(slope * duration * 44100/1024) ~ factor(noteid, levels=1:75), data=data, subset = !ispoor,
daijiajie1@1: 	pch=20, col=rgb(0.9,0.9,0.9), ylim = c(-2,2), outline=F, add=T,
daijiajie1@1: 	ylab="rise (semitones)", xlab="interval in semitones",
daijiajie1@1: 	cex.axis=0.8, lty=1, boxwex=0.6)
daijiajie1@1: # posi <- c(1:6, 8, 11, 13, 14, 18)
daijiajie1@1: # axis(1, at=posi, labels=posi-6, cex.axis=0.8)
daijiajie1@1: # text(countsByNoteID$noteid, rep(2,11), countsByNoteID$count, cex=0.7)
daijiajie1@1: mtext(c(countsByNoteID$count), at=c(countsByNoteID$noteid), cex=0.7)
daijiajie1@1: # text(1:25, rep(-1.5,11), c("(-5)", gtdata$interval_ST[-25]), cex=0.7)
daijiajie1@1: # text(2.7, -1.35 , "intervals in semitones", cex=0.7)
daijiajie1@1: text(1-.5, 1.3, "phrase 1", pos=4)
daijiajie1@1: abline(v=6.5, lty=1, col='gray')
daijiajie1@1: text(7-.5, 1.3, "phrase 2", pos=4)
daijiajie1@1: abline(v=12.5, lty=1, col='gray')
daijiajie1@1: text(13-.5, 1.3, "phrase 3", pos=4)
daijiajie1@1: abline(v=19.5, lty=1, col='gray')
daijiajie1@1: text(20-.5, 1.3, "phrase 4", pos=4)
daijiajie1@1: dev.off()
daijiajie1@1: 
daijiajie1@1: # mu versus model MAPE
daijiajie1@1: pdf("../figures/memory_vs_model_mape.pdf", width = 4.3, height = 3)
daijiajie1@1: par(mar=c(4,4,.3,.5))
daijiajie1@1: plot(mus,
daijiajie1@1:      apply(md[!md$ispoor,(1:nmu)+8], 2, function(x) {mean(abs(x),na.rm=T)}), 
daijiajie1@1:      type='l', 
daijiajie1@1:      pch=20,
daijiajie1@1:      xlab=expression(mu),
daijiajie1@1:      ylab="mean absolute pitch error",
daijiajie1@1:      ylim=c(0.18,0.30))
daijiajie1@1: abline(h=mean(abs(data$pitcherror[!data$ispoor]), na.rm=T), lty=2)
daijiajie1@1: dev.off()