gin commit from MBP-S.local

New files: 76
Modified files: 1

README.md

@@ -1,3 +1,21 @@
-# predictive_coding_in_asd
+---
+title: "Data and codes for the study of Predictive coding in ASD"
+author: "Z. Shi, L. Theisinger, F. Allenmark, R. Pistorius, H. Müller, C. Falter-Wagner"
+output: html_notebook
+---
+
+# Folder Structure
+
+1. `/experiments`: Experimental codes and instructions
+
+This sub-folder contains Matlab codes and instructions for the duration reproduction task. The sequences of the duration reproductions are stored in the sub-folder `/experiments/seqs`. Those sequences were used for matched participants. 
+
+2. `/data`: raw data files
+
+- `rawdata.csv`: Raw reproduction trials
+- `parinfo.csv`: Participant information, including AQ, EQ, SQ, IQ etc. 
+
+3. `/figures`: store figures used in the paper. 
+
+4. `analysis-notebook.rmd`: the main analysis code and report
 
-Predictive coding in ASD: reduced adaptability of priors to environmental changes

analysis-notebook.Rmd

@@ -0,0 +1,675 @@
+---
+title: "Data Analysis of Predictive coding in ASD"
+author: "Z. Shi, L. Theisinger, F. Allenmark, R. Pistorius, H. Müller, C. Falter-Wagner"
+output: html_notebook
+---
+
+# Data Structure
+
+1. `/experiments`: Experimental codes and instructions
+
+This sub-folder contains Matlab codes and instructions for the duration reproduction task. The sequences of the duration reproductions are stored in the sub-folder `/experiments/seqs`. Those sequences were used for matched participants. 
+
+2. `/data`: raw data files
+
+- `rawdata.csv`: Raw reproduction trials
+- `outliers.csv`: those produced almost flat reproduction
+
+3. `/figures`: store figures used in the paper. 
+
+# Data Analysis
+
+## 1. load raw data
+```{r packages, message=FALSE, warning=FALSE, include=FALSE}
+# load packages
+library(tidyverse)
+library(ez)
+library(cowplot)
+library(BayesFactor) # in case need Bayes factor analysis
+library(bayestestR)
+library(rstatix) # using tidyverse friendly statistics
+library(GGally)
+
+# ---- read data and preparation -----
+rawdata = read_csv('./data/rawdata.csv') 
+rawdata$group = toupper(rawdata$group)
+# flag
+saveFig = FALSE
+```
+
+Show the raw trial structure: 
+```{r raw trial structure}
+glimpse(rawdata)
+
+```
+In the raw trials, there are several important columns which are relevant for  further analyses. 
+
+- `Duration`: The test durations generated by computer. 'pdur' is the actual presented durations by the computer. There were some fluctuations, but within 5 ms (within 1 refresh frame). 
+- `sub`, `group`, `sequence` are anonymized subject number, group and the duration sequence used in the experiment. 
+- `Reproduction`, `rep_err`: the reproduced durations and the reproduction errors compared to the given duration. 
+- `itd`, `preDuration`: the inter-trial difference in  a given trial, and the duration used in the previous trial. These were used in the sequential effect analysis. 
+
+### Participants
+
+```{r participants}
+# Total participants
+print(length (unique(rawdata$sub)))
+# total sequences (pairs)
+print(length(unique(rawdata$sequence)))
+```
+
+## The sampled durations and sequences
+
+Let's first illustrate the sequence and the distribution of the sample durations. 
+
+```{r sequence}
+# ---- illustrate one sequence -----
+fig11 = ggplot(rawdata, aes(Duration)) + geom_histogram(binwidth = 0.1, fill = I("white"), col = I('black')) + 
+  theme_classic() + xlab('Duration (secs)')
+
+# a typical sequence
+sub1 = rawdata %>% filter(sub == 'ara27')
+fig12 = ggplot(sub1, aes(trlNo, Duration, color = Volatility)) + geom_line() + 
+  xlab("Trial Sequence") + theme_classic()+
+  theme(legend.position = 'top')
+
+# histogram of the typical sequence
+
+fig13 = ggplot(sub1, aes(y=Duration, color = Volatility)) + geom_density() + theme_classic() + theme(legend.position = 'top')
+fig13
+
+fig1 = plot_grid(fig12,fig13, rel_widths =c(3,1))
+fig1 
+if (saveFig){
+  ggsave("figures/fig_sequence.png", fig1, width=4, height=3.5)
+  ggsave("figures/fig_sequence.pdf", fig1, width=4, height=3.5)
+}
+```
+
+
+## Explorative data analysis and outlier detection
+
+We first estimate two key signatures - the central tendency index (ci) and the sequential dependence index (si). 
+
+```{r linear_models}
+# using error for regression, -slope is the central tendency index. 
+ct_model <- function(df){
+  lm(rep_err ~ Duration, data = df)
+}
+
+# sequential effect using the duration from the previous trial
+seq_model <- function(df){
+  lm(rep_err ~ preDuration, data = df)
+}
+
+# exclude extreme trials (beyond [1/3D, 3D])
+vdata = rawdata %>% 
+  filter(Reproduction > Duration/3, Reproduction < 3*Duration) 
+
+# calculate slope for the central tendency as well as the sequential dependence. 
+slopes_ct <- vdata %>% group_by(sub, sequence, Volatility, group, Order) %>%
+  nest()  %>%  # nested data
+  mutate(model = map(data, ct_model)) %>%  # linear regression
+  mutate(slope = map(model, broom::tidy)) %>%  # get estimates out
+  unnest(slope, .drop = TRUE) %>% # remove raw data
+  select(-std.error,-statistic, -p.value) %>%  # remove unnecessary columns
+  spread(term, estimate) %>%   # spread estimates
+  rename(intercept_ct = `(Intercept)`, slope = Duration) %>%
+  select(-data, -model)
+
+#sequential effect
+slopes_seq <- vdata %>% filter(!is.na(preDuration)) %>%
+  group_by(sub, sequence, Volatility, group, Order) %>%
+  nest()  %>%  # nested data
+  mutate(model = map(data, seq_model)) %>%  # linear regression
+  mutate(slope = map(model, broom::tidy)) %>%  # get estimates out
+  unnest(slope, .drop = TRUE) %>% # remove raw data
+  select(-std.error,-statistic, -p.value) %>%  # remove unnecessary columns
+  spread(term, estimate) %>%   # spread estimates
+  rename(intercept_seq = `(Intercept)`, seq_slope = preDuration) %>% 
+  select(-data, -model)
+
+# merge two tables together
+slopes_all = left_join(slopes_ct, slopes_seq, 
+                       by = c("sub", "sequence", "Volatility","group","Order")) %>% 
+  mutate(ci = -slope, si = seq_slope)  #central tendency index
+
+# estimated general biases (over-/under-estimates)
+
+# individual mean interval (middle point): most of them were around 1 by design
+mInterval  = vdata %>% group_by(sub) %>%  summarise(mDur = mean(Duration))
+
+# join the mean interval, and estimate the general bias
+slopes = slopes_all %>% left_join(., mInterval, by = c('sub')) %>% 
+  mutate(gBias = (intercept_ct + slope * mDur)*1000)
+
+# change factor order for plotting
+slopes$Volatility = factor(slopes$Volatility)
+slopes$Order = factor(slopes$Order)
+slopes$group = factor(slopes$group)
+
+head(slopes)
+```
+
+Let's plot the histograms of the central tendency index and the sequential dependence index, which show some outliers. We visualize the outliers using 3-sigma rule. The 3-sigma rule only rule out 0.3% of the population if we assume the population is normal. 
+
+```{r}
+ci_3sig = mean(slopes$ci) + c(-1,1)*3*sd(slopes$ci)
+si_3sig = mean(slopes$si) + c(-1,1)*3*sd(slopes$si)
+
+hist1 = ggplot(slopes, aes(x = ci, y = ..density..)) + 
+  geom_histogram(colour = 1, fill = 'white', bins = 20)  + 
+  geom_vline(xintercept = ci_3sig, linetype = 'dashed', color = 'red') +
+  theme_classic() + xlab('Central Tendency Index')
+hist2 = ggplot(slopes, aes(x = si, y = ..density..)) + 
+  geom_histogram(colour = 1, fill = 'white', bins = 20)  + 
+    geom_vline(xintercept = ci_3sig, linetype = 'dashed', color = 'red') +
+  theme_classic() + xlab('Sequential Dependence Index')
+
+hist_fig = plot_grid(hist1,hist2, nrow = 2)
+hist_fig
+if (saveFig){
+  ggsave("figures/hist_fig.png", hist_fig, width=4, height=3.5)
+  ggsave("figures/hist_fig.pdf", fig1, width=4, height=3.5)
+}
+
+```
+
+Note, the above outlier detection method is agnostic to the groups and conditions!
+
+Let's find out those outliers and exclude by sequences for further analyses (given that ASD and TD groups were paired). 
+
+```{r}
+slopes %>% ungroup() %>% filter(ci > ci_3sig[2] | ci < ci_3sig[1] | si > si_3sig[2] | si < si_3sig[1]) %>% select(sequence) -> outlier_seq
+outlier_seq$sequence
+```
+
+Now let's visualize the outliers together with their matched pairs. 
+
+```{r outliers}
+# ---- plot  outliers  ----
+fig_outlier = rawdata %>% ungroup() %>% filter(sequence %in% outlier_seq$sequence) %>% 
+  filter(Reproduction > Duration/3, Reproduction < 3*Duration) %>% 
+  group_by(sub, group,sequence, Volatility, Duration) %>%
+  summarise(mRep = mean(Reproduction), sdr=sd(Reproduction), n = n(), 
+            se = sd(Reproduction)/sqrt(n-1)) %>% filter(n>5) %>% #approx. for linear regression without averaging first
+  ggplot(aes(Duration, mRep, color = group, 
+             group = interaction(group, Volatility), 
+             shape = Volatility, linetype = Volatility)) + 
+  geom_point()  + geom_smooth(method = 'lm', se = FALSE) + 
+  facet_wrap(~sequence, ncol = 4) +
+  theme_classic() + theme(legend.pos = 'bottom') +
+  geom_abline(slope = 1, linetype = 3) + 
+  xlab('Duration (Secs)') + ylab(' Reproduction (Secs)') + 
+  theme(strip.background = element_blank(), strip.text.x = element_blank())
+fig_outlier
+if (saveFig){
+  ggsave("figures/fig_outlier.png", fig_outlier, width=5, height=5)
+  ggsave("figures/fig_outlier.pdf", fig_outlier, width=5, height=5)
+}
+```
+Mean slopes and related statistics:
+
+```{r}
+slopes %>% filter(sequence %in% outlier_seq$sequence) %>% arrange(group, Volatility)
+vslopes = slopes %>% filter(!(sequence %in% outlier_seq$sequence))
+
+# outliers
+oslopes = slopes %>% filter(sequence %in% outlier_seq$sequence)
+
+# ANOVA
+ezANOVA(data = oslopes, dv = ci, 
+                 wid = sub, 
+                 within = Volatility,
+                 between = group)
+
+```
+### Outliers and sequential dependence
+
+Visualize the strong sequential dependence. 
+
+```{r}
+fig_sdep31 = rawdata %>% filter(sequence %in% c(31) ) %>% # the extreme sequential dependence pair
+  mutate_at("preDuration", round, 1) %>%
+  group_by(sequence, group, preDuration) %>% 
+  summarise(rep_err = mean(rep_err)) %>% 
+  ggplot(aes(preDuration, rep_err, color = group) ) + 
+    geom_point() + geom_smooth(method = 'lm') + theme_classic() + 
+    geom_abline(slope = 0, linetype = 2) + 
+#  facet_wrap(~sequence) +
+    xlab('Duration trial n-1') + ylab('Repr. Error (Secs)') + 
+  theme(strip.background = element_blank(), strip.text.x = element_blank(), legend.position = 'bottom') 
+fig_sdep31
+
+fig_o = plot_grid(fig_outlier, fig_sdep31, nrow = 1, labels = c('a','b'), rel_widths = c(2.4,1))
+if (saveFig){
+  ggsave("figures/fig_o.png", fig_o, width=7, height=3)
+  ggsave("figures/fig_o.pdf", fig_o, width=7, height=3)
+  
+}
+
+
+```
+
+## Typical reproduction performance 
+
+And here are the two typical participants produced errors from sequence 12:
+```{r typical participants}
+# ---- reproduction figure - an example -----
+# mean analysis
+mrep = rawdata %>% ungroup() %>% filter(!(sequence %in% outlier_seq$sequence)) %>% 
+  filter(Reproduction > Duration/3, Reproduction < 3*Duration) %>% 
+  group_by(sub, group,Order, sequence, Volatility, Duration) %>%
+  summarise(mRep = mean(Reproduction), sdr=sd(Reproduction), n = n(), 
+            se = sd(Reproduction)/sqrt(n-1)) %>% filter(n>3)
+
+# individual examples
+# asd
+fig_asd = mrep%>% 
+  filter( sequence == '11', group == 'ASD') %>% 
+  ggplot(aes(Duration, mRep, color = Volatility, group = Volatility, shape = Volatility)) + 
+  geom_point(size = 2) + 
+  #geom_line(aes(linetype = Volatility)) + 
+  geom_smooth(method = 'lm', aes(fill = Volatility), se = FALSE) +
+  geom_errorbar(aes(ymin = mRep - se, ymax = mRep +se), width = 0.05) + 
+  theme_classic() + theme(strip.background = element_blank()) +
+  geom_abline(slope = 1, linetype = 2) + 
+  theme(legend.position = 'none', plot.margin = unit(c(0,0,0,0),'cm')) +
+  ylab('') + xlab('')
+
+#td
+fig_td = mrep%>% 
+  filter( sequence == '11', group == 'TD') %>% 
+  ggplot(aes(Duration, mRep, color = Volatility, group = Volatility, shape = Volatility)) + 
+  geom_point(size = 2) + 
+  geom_smooth(method = 'lm', aes(fill = Volatility), se = FALSE) +
+  geom_errorbar(aes(ymin = mRep - se, ymax = mRep +se), width = 0.05) + 
+  theme_classic() + theme(strip.background = element_blank()) +
+  geom_abline(slope = 1, linetype = 2) + 
+  theme(legend.position = 'none', plot.margin = unit(c(0,0,0,0),'cm')) +
+  ylab('') + xlab('')
+
+# group mean reproductions
+
+fig_repd_gr = mrep%>% 
+  ggplot(aes(Duration, mRep, color = Volatility, group = Volatility, shape = Volatility)) + 
+  geom_point(alpha = 0.2) + 
+  geom_smooth(method = 'lm', se = FALSE, aes(fill = Volatility)) +
+  theme_classic() + theme(strip.background = element_blank()) +
+  geom_abline(slope = 1, linetype = 2) + 
+  facet_wrap(~group) + theme(legend.position = c(0.1,0.85)) +
+  xlab('Duration (Secs)') + ylab('Reproduction (Secs)')
+
+# plot individual example as inset. 
+
+fig_mrep = ggdraw() + draw_plot(fig_repd_gr) + 
+  draw_plot(fig_asd, x = 0.3, y = 0.13, width  = .2, height = .3) + 
+  draw_plot(fig_td, x = 0.75, y = 0.13, width  = .2, height = .3) 
+  
+fig_mrep
+if(saveFig){
+  ggsave("figures/fig_reproduction.png", fig_mrep, width=7, height=3.5)
+  ggsave("figures/fig_reproduction.pdf", fig_mrep, width=7, height=3.5)
+}
+```
+
+As we can see from the patterns above. The ASD participant produced relative flat errors, while the TD participant showed a strong central tendency effect (shorts being overestimated and longs being underestimated). 
+
+
+### Visulize CTE and sequential dependence
+
+Let's visualize the biases (central tendency and serial dependence) for two groups (ASD vs. TD)
+
+```{r plot_biases}
+pd = position_dodge(width = 0.05)
+
+# plot CTI and SI together
+fig_biases = vslopes%>% 
+  group_by(group, Volatility, Order) %>% 
+  summarise(msi = mean(si), n = n(), se_si = sd(si)/sqrt(n),
+            mci = mean(ci), se_ci = sd(ci)/sqrt(n)) %>%
+  ggplot(aes(msi, mci, color = group, shape = Volatility, 
+              group=interaction(Order, group))) + 
+  geom_hline(yintercept = 0, color = 'gray', linetype = 'dashed') + 
+  geom_vline(xintercept = 0, color = 'gray', linetype = 'dashed') + 
+  geom_point(size = 2) + 
+  geom_line(aes(linetype = Order)) + 
+  geom_errorbar(aes(ymin = mci - se_ci, ymax = mci + se_ci), width = 0.01 ) + 
+  geom_errorbarh(aes(xmin = msi - se_si, xmax = msi + se_si), height = 0.01) + 
+  xlab('Serial Dependence') + ylab('Central Tendency') + 
+  scale_y_continuous(labels = scales::percent) + 
+  scale_x_continuous(labels = scales::percent) + 
+  guides(color = guide_legend(title = 'Group'), 
+         linetype = guide_legend(title = 'Volatility Order'),
+         ) +
+  theme_classic() +
+  theme(legend.position = 'bottom')
+fig_biases
+if (saveFig){
+  ggsave("figures/fig_biases.png", fig_biases, width=5, height=4)
+  ggsave("figures/fig_biases.pdf", fig_biases, width=5, height=4)
+}
+
+
+```
+
+By visual inspection, individuals with ASD exhibited less central tendency relative to their matched TD controls (the red lines below the cyan lines in the above figure), while the local serial dependence was relatively comparable between two groups except in the low volatility condition when that condition started first. Interestingly, the central tendency and serial dependence were similar for both groups when the high-volatility session started first (the solid lines), while they differed when the low-volatility session started first (the dashed lines). 
+
+Given that the main difference was shown in CTE. We further plot the mean CTEs. 
+
+```{r}
+#pd = position_dodge(width = 0.5)
+# separate plots for appendix
+fig_cti = vslopes %>% 
+  group_by(group, Volatility, Order) %>% 
+  summarise(mci = mean(ci), n = n(), se = sd(ci)/sqrt(n)) %>%
+  ggplot(aes(Volatility, mci, shape = Order, color = Order, group = Order)) + 
+  geom_line() + geom_point(size= 2) + 
+  #geom_bar(stat = 'identity', position = pd, width = 0.5) + 
+  facet_wrap(~group)+
+  geom_errorbar(aes(ymin = mci - se, ymax = mci + se), width = 0.2) + 
+  theme_classic() + theme(legend.position = 'bottom', strip.background = element_blank()) +
+  xlab('Volatility') + ylab('CTI') + 
+  scale_y_continuous(labels = scales::percent) + 
+  guides(color = guide_legend(title = 'Order'), fill = guide_legend(title = 'Order'))
+fig_cti
+
+```
+Let's do the sequential dependence as well. 
+
+```{r}
+pd = position_dodge(width = 0.5)
+# separate plots for appendix
+fig_sdi = vslopes %>% 
+  group_by(group, Volatility, Order) %>% 
+  summarise(msi = mean(si), n = n(), se = sd(si)/sqrt(n)) %>%
+  ggplot(aes(Volatility, msi, shape = Order, color = Order, group = Order)) + 
+  geom_line() + geom_point(size = 2) +
+  #geom_bar(stat = 'identity', position = pd, width = 0.5) + 
+  facet_wrap(~group)+
+  geom_errorbar(aes(ymin = msi - se, ymax = msi + se), width = 0.2) + 
+  theme_classic() + theme(legend.position = 'bottom', strip.background = element_blank()) +
+  xlab('Volatility') + ylab('SDI') + 
+  scale_y_continuous(labels = scales::percent) + 
+  guides(color = guide_legend(title = 'Order'), fill = guide_legend(title = 'Order'))
+fig_sdi
+fig3 = plot_grid(fig_cti, fig_sdi, nrow = 1, labels = c("a","b"))
+fig3
+if (saveFig){
+  ggsave('figures/fig3.pdf',fig3, width = 7, height = 3.5)
+  ggsave('figures/fig3.png',fig3, width = 7, height = 3.5)
+}
+
+```
+
+### Statistics
+
+1. Central tendency effect
+
+Average CTIs
+```{r}
+vslopes %>% group_by(group) %>% summarise(mcti = mean(ci))
+vslopes %>% group_by(Volatility) %>% summarise(mcti = mean(ci))
+vslopes %>% group_by(Order) %>% summarise(mcti = mean(ci))
+# calculate mean elevated CTIs between HV First vs. LV First
+vslopes %>% group_by(group, Volatility, Order) %>% summarise(mcti = mean(ci)) %>%
+  pivot_wider(names_from = Order, values_from = mcti) %>% 
+  mutate(dCTI = `HV First` - `LV First`) %>% 
+  group_by(group) %>% summarise(md = mean(dCTI))
+```
+
+
+```{r ANOVAs_cti}
+vslopes = as.data.frame(vslopes) # required by rstatix
+# ---- central tendency index----
+# repeated measures ANOVA on central tendency index
+anova1 = anova_test(data = vslopes, dv = ci, 
+                 wid = sub, 
+                 within = Volatility,
+                 between = c(group, Order))
+anova1
+
+## separate for Volatility order
+anova1a = anova_test(data = vslopes %>% filter(group == 'ASD'), 
+                  dv = ci, 
+                  wid = sub, 
+                  within = Volatility,
+                  between = Order)
+anova1a
+
+anova1b = anova_test(data = vslopes  %>% filter(group == 'TD'), 
+                  dv = ci, 
+                  wid = sub, 
+                  within = Volatility,
+                  between = Order)
+anova1b
+
+```
+
+
+
+2. Bayes factor analysis for pair-wise comparison
+
+The above analysis showed that the main difference came from the central tendency. Here we further get the central tendency bias and Bayes factor analyses:
+```{r cti_bayes}
+
+# --- Bayes t-tests
+bftest = function(df){
+  df = as.data.frame(df)
+  # get the means
+  rdf = df%>% summarise(mci = mean(ci), 
+                        se_ci = sd(ci)/sqrt(n()),
+                        msi = mean(si), 
+                        se_si = sd(si)/sqrt(n()))
+  rdf$ci_bf = ttestBF(df$ci, mu = 0) %>% extractBF() %>% .$bf
+  rdf$si_bf = ttestBF(df$si, mu = 0) %>% extractBF() %>% .$bf
+  return(rdf)
+}
+
+vslopes %>% group_by(group, Volatility, Order) %>% nest() %>%
+  mutate(bf = map(data, bftest)) %>% unnest(bf, .drop = TRUE)
+
+# group comparison for the low-vol-first, high-vol session
+vslopes %>% filter(Volatility == 'High Vola.', Order == 'LV First') %>% as.data.frame() ->v11
+t_test(data = v11, formula = ci ~ group) 
+ttestBF(data = v11, formula = ci ~ group)
+
+vslopes %>% filter(Volatility == 'Low Vola.', Order == 'LV First') %>% as.data.frame() ->v12
+t_test(data = v12, formula = ci ~ group) 
+ttestBF(data = v12, formula = ci ~ group)
+```
+
+
+3. ANOVA analyses for the serial dependence indices:
+
+```{r ANOVAs_sdi}
+# ---- Serial dependence index----
+anova2 = ezANOVA(data = vslopes, dv = si, 
+                 wid = sub, 
+                 within = .(Volatility),
+                 between = .(group, Order))
+anova2$ANOVA
+
+# Bayes factors
+bf = anovaBF(si ~ Volatility + group + Order, data = vslopes, whichRandom = "sub")
+bayesfactor_inclusion(bf) #bayes inclusion values
+
+## separate for session order
+anova2a = ezANOVA(data = vslopes %>% filter(Order == 'LV First'), 
+                  dv = si, 
+                  wid = sub, 
+                  within = .(Volatility),
+                  between = .(group))
+anova2a$ANOVA
+
+anova2b = ezANOVA(data = vslopes %>% filter(Order == 'HV First'), 
+                  dv = si, 
+                  wid = sub, 
+                  within = .(Volatility),
+                  between = .(group))
+anova2b$ANOVA
+
+```
+
+
+
+4. General biases
+
+Additionally we examined the general biases:
+```{r gBias}
+# ---- descriptive of general bias ----
+vslopes %>% group_by(group) %>% summarise(mg = mean(gBias), se = sd(gBias)/sqrt(n())) 
+
+# ANOVA shows no differences
+anova2 = ezANOVA(data = vslopes, dv = gBias, 
+                 wid = sub, 
+                 within = .(Volatility),
+                 between = .(group, Order))
+anova2$ANOVA
+
+bf = anovaBF(gBias ~  group + Volatility + Order, 
+             data = as.data.frame(vslopes), whichRandom = "sub")
+bayesfactor_inclusion(bf) #bayes inclusion values
+
+
+```
+We then visualize the general biases
+```{r plotgBias}
+# --- general bias
+fig_bias = vslopes %>% group_by(group, Volatility, Order) %>%
+  summarise(mbias = mean(gBias), n = n(), se = sd(gBias)/sqrt(n)) %>%
+  ggplot(aes(Volatility, mbias, color = Order,linetype = group, group = interaction(Order, group), shape = group)) + 
+  geom_point(size = 3, position = pd) + 
+  geom_line(position = pd) +
+  geom_errorbar(aes(ymin = mbias - se, ymax = mbias + se), width = 0.2, position = pd) + 
+  theme_classic()  + theme(legend.position = 'bottom') + xlab('Volatility') + 
+  ylab('Mean overestimation (ms)') 
+fig_bias
+```
+
+There was no difference in two groups in general biases, although both groups were positive overestimated. 
+
+Next we examined the reproduced variability:
+```{r rep_var}
+# ---- Reproduction variability -----
+msds_r <- mrep  %>% 
+  filter(n>5) %>% 
+  group_by(group, Order, Volatility,  sub) %>%
+  summarise(msd = mean(sdr), n=n(), msd_se = sd(sdr)/sqrt(n))
+  
+
+sd_ANOVA <- ezANOVA(msds_r, dv=msd, wid=sub, between=.(group,Order), within=Volatility)
+sd_ANOVA$ANOVA
+
+mmsds_r <- msds_r %>% summarize(mmsd=mean(msd*1000), n = n(), se = sd(msd*1000)/sqrt(n))
+
+pd = position_dodge(width = 0.5)
+fig_sd <- ggplot(mmsds_r, aes(Volatility, mmsd, shape = Order, color = Order, group = Order)) + 
+  geom_line() + geom_point() + 
+  #geom_bar(stat = 'identity', position = pd, width = 0.5) + 
+  geom_errorbar(aes(ymin = mmsd - se, ymax = mmsd + se), width = 0.3) +
+  facet_wrap(~group) +
+  ylab('Mean Standard Deviation (ms)') + xlab('Volatility') + theme_classic()  
+#  coord_cartesian(ylim = c(140, 200)) 
+fig_sd
+
+if(saveFig){
+  ggsave('figures/fig_sd.png', fig_sd, width=5, height=4)
+  ggsave('figures/fig_sd.pdf', fig_sd, width=5, height=4)
+}
+
+
+
+```
+
+## Correlation analysis
+
+Next we import individual participant ratings and do correlation analyses. 
+
+```{r parinfo}
+pinfo = read.csv('data/parinfo.csv')
+
+# join with ctis, sdis
+res = left_join(vslopes, pinfo, by = c('group','sequence'))
+res$group = as.factor(res$group)
+
+# scatter plot to see any relations
+fig_corr = ggduo(res, c('AQ','EQ','BDI'), c('ci','si'), mapping = aes(color = group, shape = Order),
+      types = list(continuous = 'smooth_lm', se = FALSE), 
+      showStrips = FALSE, legend = 3,
+      columnLabelsY = c('CTI','SDI')) + 
+  theme_classic() + theme(legend.position = 'bottom', strip.background = element_blank()) 
+if(saveFig){
+  ggsave('figures/fig_corr.png', fig_corr, width=5, height=4)
+  ggsave('figures/fig_corr.pdf', fig_corr, width=5, height=4)
+}
+
+fig_corr
+```
+ Get linear regression out
+ 
+```{r}
+ci_AQ <- function(df){
+  lm(ci ~ AQ, data = df)
+}
+ci_EQ <- function(df){
+  lm(ci ~ EQ, data = df)
+}
+ci_BDI <- function(df){
+  lm(ci ~ BDI, data = df)
+}
+
+si_AQ <- function(df){
+  lm(si ~ AQ, data = df)
+}
+si_EQ <- function(df){
+  lm(si ~ EQ, data = df)
+}
+si_BDI <- function(df){
+  lm(si ~ BDI, data = df)
+}
+
+# regression
+reg_res <- res %>% group_by(Volatility, group, Order) %>%
+  nest()  %>%  # nested data
+  mutate(ci_AQ = map(data, ci_AQ), ci_EQ = map(data, ci_EQ),
+         ci_BDI = map(data, ci_BDI), 
+         si_AQ = map(data, si_AQ), si_EQ = map(data, si_EQ),
+         si_BDI = map(data, si_BDI)) %>%
+  mutate(mci_AQ = map(ci_AQ, broom::tidy), 
+         mci_EQ = map(ci_EQ, broom::tidy),
+         mci_BDI = map(ci_BDI, broom::tidy),
+         msi_AQ = map(si_AQ, broom::tidy),
+         msi_EQ = map(si_EQ, broom::tidy),
+         msi_BDI = map(si_BDI, broom::tidy)  ) %>%
+  unnest(cols = c(mci_AQ, mci_EQ, mci_BDI, msi_AQ, msi_EQ, msi_BDI), names_repair = 'unique' , .drop = TRUE) %>%
+  select(-ci_AQ, -ci_BDI, -ci_EQ, -si_AQ, -si_EQ, -si_BDI, -data, 
+         -starts_with('st')) %>% 
+  filter(term == 'AQ') # remove intercept
+
+```
+
+Given there is no significant of slopes (after correction), but some difference in intercepts, we do general linear regression, without separating groups. 
+
+```{r}
+reg_res2 <- res %>% ungroup() %>% group_by(Order, Volatility) %>%
+  nest()  %>%  # nested data
+  mutate(ci_AQ = map(data, ci_AQ), ci_EQ = map(data, ci_EQ),
+         ci_BDI = map(data, ci_BDI), 
+         si_AQ = map(data, si_AQ), si_EQ = map(data, si_EQ),
+         si_BDI = map(data, si_BDI)) %>%
+  mutate(mci_AQ = map(ci_AQ, broom::tidy), 
+         mci_EQ = map(ci_EQ, broom::tidy),
+         mci_BDI = map(ci_BDI, broom::tidy),
+         msi_AQ = map(si_AQ, broom::tidy),
+         msi_EQ = map(si_EQ, broom::tidy),
+         msi_BDI = map(si_BDI, broom::tidy)  ) %>%
+  unnest(cols = c(mci_AQ, mci_EQ, mci_BDI, msi_AQ, msi_EQ, msi_BDI), names_repair = 'unique' , .drop = TRUE) %>%
+  select(-ci_AQ, -ci_BDI, -ci_EQ, -si_AQ, -si_EQ, -si_BDI, -data, 
+         -starts_with('st')) %>% 
+  filter(term == 'AQ') # remove intercept
+```
+Again, there is no significant of slopes. 
+
+
+

data/parinfo.csv

@@ -0,0 +1,65 @@
+AQ,EQ,SQ,BDI,WST_IQ,group,ID,sequence
+29,46,40,13,80,ASD,A4CB,1
+30,30,8,15,95,ASD,83A6,2
+31,61,21,3,110,ASD,A71D,4
+35,16,46,9,110,ASD,8AE1,5
+33,21,28,18,114,ASD,79B8,6
+45,17,70,13,104,ASD,5741,7
+31,39,19,2,82,ASD,5664,8
+45,25,67,25,118,ASD,31A3,9
+22,36,32,9,107,ASD,81BD,10
+46,13,44,3,118,ASD,34AB,11
+38,29,29,20,110,ASD,422E,12
+43,25,32,9,101,ASD,EEAB,13
+42,15,24,1,101,ASD,FD8F,14
+36,32,36,10,122,ASD,11D9,15
+44,10,49,6,114,ASD,B77F,16
+27,46,34,5,110,ASD,F7B4,17
+27,44,28,3,114,ASD,DFAC,19
+38,31,9,4,99,ASD,2F2D,20
+47,13,40,8,107,ASD,9E75,21
+41,28,26,6,118,ASD,DA48,22
+37,27,35,11,95,ASD,D4D6,23
+42,20,37,7,104,ASD,E785,25
+24,34,26,7,107,ASD,6138,26
+38,22,22,23,122,ASD,EBC4,27
+45,16,52,8,122,ASD,AA04,28
+44,14,34,0,107,ASD,D680,29
+47,22,67,4,118,ASD,40A9,31
+28,49,44,17,110,ASD,3443,32
+39,25,17,9,110,ASD,E6E4,33
+37,13,38,29,110,ASD,8B5B,34
+37,22,39,35,101,ASD,711F,36
+26,22,35,10,97,ASD,ECBF,37
+13,64,18,0,101,TD,C930,1
+11,55,21,1,99,TD,0DFF,2
+5,75,8,0,107,TD,C892,4
+17,69,28,0,99,TD,5249,5
+12,65,16,4,107,TD,81B7,6
+18,57,38,2,133,TD,F19C,7
+13,55,22,3,110,TD,0CC9,8
+32,11,24,5,104,TD,6823,9
+11,68,41,3,110,TD,F4BD,10
+13,46,29,6,97,TD,0400,11
+23,37,20,3,118,TD,E08D,12
+24,45,51,33,133,TD,8F93,13
+24,41,19,27,84,TD,89EA,14
+12,42,22,4,107,TD,5FDC,15
+17,50,23,8,104,TD,6FFF,16
+9,53,31,1,122,TD,5181,17
+11,65,37,0,104,TD,09A6,18
+14,64,19,3,110,TD,63BF,19
+23,35,35,4,110,TD,5EAB,20
+10,61,24,4,114,TD,9419,21
+9,51,25,6,97,TD,DA2B,22
+20,66,21,2,122,TD,5A9C,23
+12,40,39,18,90,TD,0231,25
+15,46,20,7,99,TD,F630,26
+7,58,28,8,122,TD,5C3B,27
+21,61,18,8,122,TD,3169,28
+21,42,43,1,118,TD,0193,31
+19,24,10,0,77,TD,EAA8,32
+18,54,43,2,133,TD,665E,33
+21,32,18,0,118,TD,0476,34
+19,31,29,2,118,TD,B628,36
+11,49,32,3,101,TD,4CD4,37

experiments/CAudio.m

@@ -0,0 +1,185 @@
+classdef CAudio < handle
+    % a class managing audio interface
+    % last modified by Strongway: genNoise with randn, normalize noise
+    properties
+        freq;
+        channels;
+        pins; 
+        latency; 
+        pHandle;
+        hwName; 
+        hwType;
+    end
+    
+    methods
+        function obj = CAudio(varargin)
+            %constructor: frequency, channels, latency
+            p = inputParser;
+            p.addParamValue('freq',96000,@isnumeric);
+            p.addParamValue('channels',2,@isnumeric);
+            p.addParamValue('latency', 0.001, @(x) x>0);
+            p.addParamValue('pins',[1 2],@isnumeric); %define the pin number of the channels
+            p.addParamValue('hardware','sndCard', @(x) any(strcmpi(x,{'sndCard','Aout','Ain'})));
+            p.parse(varargin{:});
+            
+            obj.freq = p.Results.freq;
+            obj.channels = p.Results.channels;
+            obj.latency = p.Results.latency;
+            obj.hwName = p.Results.hardware; 
+            obj.pins = p.Results.pins;
+            switch obj.hwName
+                case {'sndCard'}
+                    obj.hwType = 1;
+                    try
+                        InitializePsychSound(1);
+                        obj.pHandle = PsychPortAudio('Open', [],[],2,obj.freq,obj.channels);
+                        PsychPortAudio('LatencyBias', obj.pHandle, obj.latency);
+                        obj.latency = PsychPortAudio('LatencyBias', obj.pHandle);
+                    catch ME
+                        disp(ME.message);
+                    end
+                case {'Aout','Ain'}
+                    obj.hwType = 2;
+                    obj.pHandle = analogoutput('nidaq',obj.hwName);
+                    achan = addchannel(obj.pHandle,obj.pins);
+                    obj.channels = size(achan,1);
+                    set(obj.pHandle,'SampleRate',obj.freq);
+                    set(obj.pHandle,'TriggerType','manual');
+                otherwise
+                    obj.hwType = 1;
+            end
+        end
+        
+        function prepare(obj, data)
+            switch obj.hwType 
+                case 1
+                    PsychPortAudio('FillBuffer', obj.pHandle, data);
+                case 2
+                    putdata(obj.pHandle,data');
+                    start(obj.pHandle);
+            end
+            
+        end
+        
+        function onsetTime = present(obj,startTime)
+            if nargin < 2
+                startTime = 0;
+            end
+            switch obj.hwType
+                case 1
+                    PsychPortAudio('Start', obj.pHandle, 1, startTime, 0);
+                    active = 0;
+                    while active == 0
+                        s = PsychPortAudio('GetStatus',obj.pHandle);
+                        active = s.Active;
+                    end
+                    onsetTime = s.StartTime; %get true onset time
+                case 2        
+                    if startTime > 0
+                        waitTime = max(0,startTime - GetSecs);
+                        WaitSecs(waitTime);
+                    end
+                    trigger(obj.pHandle);
+                    onsetTime = getSecs;
+            end
+        end
+        
+        function stop(obj,afterFinish)
+            if nargin < 2
+                afterFinish = 0; %immediately stop
+            end
+            switch obj.hwType
+                case 1
+                    PsychPortAudio('Stop', obj.pHandle, afterFinish);
+                case 2
+                    stop(obj.pHandle);
+            end
+        end
+        function open(obj)
+            obj.pHandle = PsychPortAudio('Open', [],[],2,obj.freq,obj.channels);
+            PsychPortAudio('LatencyBias', obj.pHandle, obj.latency);
+            obj.latency = PsychPortAudio('LatencyBias', obj.pHandle);
+        end
+        function close(obj)
+            switch obj.hwType
+                case 1
+                    PsychPortAudio('Stop', obj.pHandle);
+                    % Close the audio device:
+                    PsychPortAudio('Close', obj.pHandle);
+            end    
+        end
+        
+        function tone = genTone(obj, freq, duration)
+            t = MakeBeep(freq, duration, obj.freq);
+            tone = repmat(t(1:end-1),obj.channels(1),1);
+        end
+        
+        function tone = genNoise(obj,duration)
+            % change to randn (normal distribution) - Oct.14. 2013
+            num_samp = obj.freq*duration;
+            t = randn(1,num_samp);
+            tone = repmat(t,obj.channels(1),1);
+            % rescale 
+            tone = tone./3; % 3 sigma
+        end
+        
+        function [tone, pn] = genPinkNoise(obj, f0, duration, dbc_per_hz, num_taps)
+            % based on code phase_noise by Jeff Schenck 11/21/95
+            % f0 reference frequency (must be in Hz.)
+            % dbc_per_hz power per hertz relative to carrier at ref. freq.
+            % num_taps number of filter taps in AR 1/f filter
+            % (optional; default = 100)
+
+            % pn phase-modulated 1/f process
+            % theta 1/f process (before phase modulation)
+
+            % Check input.
+
+            if dbc_per_hz >= 0
+                error('Power per Hz. must be negative.');
+            elseif f0 <= 0
+                error('Reference frequency must be positive.');
+            end
+
+            if nargin < 5
+                num_taps = 100;
+            end
+
+            num_samp = obj.freq * duration;
+            % Generate white noise. Apply gain for desired dBc/Hz. Warn user
+            % if gain is too large (gain thresholds have been chosen somewhat
+            % arbitrarily -- needs work).
+
+            gain = sqrt(2*pi * f0 * 10^(dbc_per_hz/10));
+            wn = gain * randn(1,num_samp);
+
+            fprintf('Gain applied to white noise = %f.\n', gain);
+            if gain >= 1
+                fprintf('WARNING: Narrowband approximation no longer valid.\n');
+            elseif gain >= .5
+                fprintf('WARNING: Narrowband approximation on the verge of collapse.\n');
+            end
+
+
+            % Generate 1/f AR filter and apply to white noise to produce 1/f
+            % noise.
+
+            a = zeros(1,num_taps);
+            a(1) = 1;
+            for ii = 2:num_taps
+                a(ii) = (ii - 2.5) * a(ii-1) / (ii-1);
+            end
+            theta = filter(1,a,wn);
+
+            tone = repmat(theta,obj.channels(1),1);
+            % Phase modulate.
+
+            pn = exp(i*theta);
+
+            % normalize
+            sd3 = std(tone(1,:))*3; % 3 sigma
+            tone = tone./sd3; 
+        end
+        
+    end
+end

experiments/CDisplay.m

@@ -0,0 +1,307 @@
+classdef CDisplay < handle
+    % a display class, manage the most common display functions, such as
+    % display text, fixation, flip, degree of visual angle. It also stores
+    % display informations, such as resolution, center x,y, ifi ...
+    % methods:
+    %   1. constructor
+    %       obj=CDisplay('monitorSize',size,'viewDistance',viewDistance,'bgColor',bgColor,'Color',color,'fontSize',fz,'skipSync',bSync);
+    %   2. display text
+    %       dispText(txt[,flip,clearBackground]);
+    %   3. close screen
+    %       close();
+    %   4. calculate visual angle (deg) to pixels
+    %       deg2pix(degree)
+    %   5. display fixation
+    %       dispFixation(size [,type, flip, clearBackground]);
+    %   6. flip: bring back buffer to front (display)
+    %       flip(obj [, clearBackground])    
+    % Screen class of PTB
+    % Last modify: 7.7.2011
+    % Created by: Z. Shi, shi@lmu.de
+    
+    % 14.07.2011 add createShape function to create simple shapes
+    % 07.07.2014 add debug option - fullWindow, 
+
+    properties
+        wnd = -1;   %window handle
+        bSkipSync = 0;
+        fullWindow = 1; % full screen
+        ifi = -1;   %inter-frame interval (refresh rate)
+        cx = -1;    % center x
+        cy = -1;    % center y
+        pdeg;       % pixels per degree
+        bgColor = 0;     % background color
+        color = 255;     % front color
+        fontSize = 14;
+        lineWidth = 60;
+        lineSpace = 1.5;
+        inch = 20;
+        viewDistance = 57;
+        nItem = 0;
+        items; 
+    end
+    
+    methods
+        function obj = CDisplay(varargin) 
+        	% constructure
+            % inch,viewDistance, bgColor, color, fontsize
+            p = inputParser;
+            p.addParamValue('monitorSize',20,@isnumeric);
+            p.addParamValue('viewDistance',57,@isnumeric);
+            p.addParamValue('bgColor',[0 0 0],@isnumeric);
+            p.addParamValue('Color',[255 255 255],@isnumeric);
+            p.addParamValue('fontSize',14,@isnumeric);
+            p.addParamValue('lineWidth',60,@isnumeric);
+            p.addParamValue('lineSpace',1.5,@isnumeric);
+            p.addParamValue('skipSync',0,@isnumeric);
+            p.addParamValue('fullWindow',1,@isnumeric);
+            p.parse(varargin{:});
+            
+            %init screens
+            obj.inch = p.Results.monitorSize;
+            obj.viewDistance = p.Results.viewDistance;
+            obj.bgColor = p.Results.bgColor;
+            obj.color = p.Results.Color;
+            obj.fontSize = p.Results.fontSize;
+            obj.lineWidth = p.Results.lineWidth;
+            obj.lineSpace = p.Results.lineSpace;
+            obj.bSkipSync = p.Results.skipSync;
+            obj.fullWindow = p.Results.fullWindow;
+           try
+                InitializeMatlabOpenGL;
+                AssertOpenGL;
+                priority = MaxPriority('KbCheck');
+                Priority(priority);
+                HideCursor; 
+                if obj.bSkipSync
+                    Screen('Preference','SkipSyncTests',1);
+                else
+                    Screen('Preference','SkipSyncTests',0);
+                end
+
+                screens=Screen('Screens');
+                screenNumber=max(screens);
+                if obj.fullWindow
+                 [obj.wnd wsize] = Screen('OpenWindow',screenNumber); % Open On Screen window, mainWnd
+                else
+                    [obj.wnd wsize] = Screen('OpenWindow',screenNumber,[],[0 0 600 400]); 
+                end
+                obj.ifi=Screen('GetFlipInterval', obj.wnd); 
+                
+                obj.cx = wsize(3)/2; %center x
+                obj.cy = wsize(4)/2; %center y
+                pix = obj.inch*2.54/sqrt(1+9/16)/wsize(3);  % calculate one pixel in cm
+                obj.pdeg = round(2*tan((1/2)*pi/180) * obj.viewDistance / pix); 
+            catch ME
+                Screen('CloseAll');
+                Priority(0);
+                disp(ME.message);
+                disp('error in initial display');
+            end
+        end % end of constructor
+
+        function nframes = sec2frames(obj,secs)
+            %convert seconds to number of frames
+            nframes = round(secs/obj.ifi);
+        end
+        
+        function dispText(obj,txt, flip,clearBackground)
+        % dispText
+            try
+                if nargin < 3  
+                    flip = 1;
+                    clearBackground = 1;
+                end
+                if nargin == 3
+                    clearBackground = 1;
+                end
+                if clearBackground 
+                	Screen('FillRect',obj.wnd,obj.bgColor);
+                end
+                Screen('TextSize',obj.wnd,obj.fontSize);
+                DrawFormattedText(obj.wnd,txt,'center','center',obj.color,obj.lineWidth,[],[],obj.lineSpace);
+                if flip
+                    Screen('Flip', obj.wnd);
+                end  
+            catch ME
+                Screen('CloseAll');
+                Priority(0);
+                disp(ME.message);
+            end
+        end % end of dispText method
+        
+        function itemIndex = createItem(obj,itemData)
+            %create items (texture, image, etc.)
+             itemIndex= Screen('MakeTexture', obj.wnd, itemData);
+             obj.nItem = obj.nItem + 1;
+             obj.items(obj.nItem) = itemIndex;
+        end
+        
+        function itemIndex = createShape(obj,name,x,y,varargin)
+            %create simple shape 
+            %be sure to do this before the trial starts
+            p = inputParser;
+            p.addRequired('name', @(x) any(strcmpi(x,{'rectangle','circle'})));
+            p.addRequired('x',@(x) x>0);
+            p.addRequired('y',@(x) x>0);
+            p.addParamValue('fill',1,@isnumeric);
+            p.addParamValue('border',0.1,@(x) x>0);
+            p.addParamValue('bgColor',obj.bgColor,@isnumeric);
+            p.addParamValue('color',obj.color,@isnumeric);
+            p.parse(name,x,y,varargin{:});
+            
+            xp = round(p.Results.x * obj.pdeg)/2;
+            yp = round(p.Results.y * obj.pdeg)/2; %convert to pixels
+            bp = round(p.Results.border * obj.pdeg);
+            bc = p.Results.bgColor;
+            fc = p.Results.color;
+            if length(bc) == 1
+                bc = repmat(bc,3,1);
+            end
+            if length(fc) == 1
+                fc = repmat(fc,3,1);
+            end
+            data = zeros(xp*2,yp*2,3); %store in rgb format
+            switch p.Results.name
+                case {'circle'}
+                    if p.Results.fill == 1 %fill
+                        for ix = 1:xp*2
+                            for iy = 1:yp*2
+                                if (ix-xp)*(ix-xp)/xp/xp + (iy-yp)*(iy-yp)/yp/yp < 1 
+                                    data(ix,iy,:) = fc;
+                                else
+                                    data(ix,iy,:) = bc;
+                                end
+                            end
+                        end
+                    else % frame
+                        for ix = 1:xp*2
+                            for iy = 1:yp*2
+                                if (ix-xp)*(ix-xp)/xp/xp + (iy-yp)*(iy-yp)/yp/yp < 1 && ...
+                                        (ix-xp)*(ix-xp)/(xp-bp)/(xp-bp) + (iy-yp)*(iy-yp)/(yp-bp)/(yp-bp) >= 1 
+                                    data(ix,iy,:) = fc;
+                                else
+                                    data(ix,iy,:) = bc;
+                                end
+                            end
+                        end
+                    end
+                case {'rectangle'}
+                    if p.Results.fill == 1 %fill
+                        for ix = 1:xp*2
+                            for iy = 1:yp*2
+                                    data(ix,iy,:) = fc;
+                            end
+                        end
+                    else %frame
+                        for ix = 1:xp*2
+                            for iy = 1:yp*2
+                                if abs(ix-xp)>xp-bp || abs(iy-yp)>yp-bp
+                                    data(ix,iy,:) = fc;
+                                else
+                                    data(ix,iy,:) = bc;
+                                end
+                            end
+                        end
+                    end
+            end %end of switch
+            %create texture
+            itemIndex= Screen('MakeTexture', obj.wnd, data);
+            obj.nItem = obj.nItem + 1;
+            obj.items(obj.nItem) = itemIndex;
+            
+        end
+        function dispItems(obj, xys, itemIndex, itemSizes,rotations, flip)
+            %disp items at xys (in visual angle, center 0,0)
+            if nargin < 6
+                flip = 1;
+            end
+            if nargin < 5
+                rotations = [];
+            end
+            if nargin < 4
+                itemSizes = [obj.cx / obj.pdeg, obj.cy/obj.pdeg]*2;
+            end
+            destRects = zeros(4, length(itemIndex));
+            for iObj = 1: length(itemIndex)
+                if size(itemSizes,1) == 1
+                    itemRect = [0 0 itemSizes*obj.pdeg];
+                else
+                    itemRect = [0 0 itemSizes(iObj,:)*obj.pdeg];
+                end
+                rect = CenterRectOnPoint(itemRect,obj.pdeg*xys(iObj,1)+obj.cx, obj.pdeg*xys(iObj,2)+obj.cy);
+                destRects(:,iObj) = rect';
+            end
+            Screen('DrawTextures',obj.wnd,itemIndex,[],destRects,rotations);
+            if flip
+                Screen('Flip', obj.wnd);
+            end
+        end
+        
+        function close(obj)
+        %% dispClose: close display
+            %close display and delete Screen
+            Screen('CloseAll');
+            ShowCursor;
+            Priority(0);
+            obj.wnd = -1;
+        end % end of closeDisp
+        
+        function pixs=deg2pix(obj, degree) 
+        %% deg2pix: calculate degree to pixels
+            screenWidth = obj.inch*2.54/sqrt(1+9/16);  % calculate screen width in cm
+            pix=screenWidth/obj.cx/2;  %calculates the size of a pixel in cm 
+            pixs = round(2*tan((degree/2)*pi/180) * obj.viewDistance / pix); 
+        end %end of deg2pix
+        
+        function deg = pix2deg(obj, pixels)
+            deg = 1/obj.pdeg*pixels;
+        end
+        
+        function dispFixation(obj,sz,type,flip,clearBackground)
+        %% dispCrossFix: display cross fixation
+            if nargin < 3
+                type = 1;
+                flip = 1;
+                clearBackground = 1;
+            end
+            if nargin == 3
+                flip = 1;
+                clearBackground = 1;
+            end
+            if nargin == 4
+                clearBackground = 1;
+            end
+            if clearBackground
+                Screen('FillRect',obj.wnd, obj.bgColor);
+            end
+            if type == 1
+                rect = CenterRectOnPoint([0 0 2 sz],obj.cx,obj.cy);
+                Screen('FillRect',obj.wnd,obj.color,rect);
+                rect = CenterRectOnPoint([0 0 sz 2],obj.cx,obj.cy);
+                Screen('FillRect',obj.wnd,obj.color, rect);
+            else
+                rect = CenterRectOnPoint([0 0 sz sz],obj.cx,obj.cy);
+                Screen('FillOval',obj.wnd,obj.color,rect);
+            end
+            if flip
+                Screen('Flip',obj.wnd);
+            end
+        end %end of dispCrossFix
+                        
+        function [vbl visonset] = flip(obj,clearBackground, when)
+        %% flip: put back buffer to screen
+            if nargin < 3
+                when = 0;
+            end
+            if nargin < 2
+                clearBackground = 0;
+            end
+            if clearBackground
+                Screen('FillRect',obj.wnd,obj.bgColor);
+            end
+            [vbl visonset] = Screen('Flip',obj.wnd, when);
+        end
+        
+    end
+end

experiments/CExp.m

@@ -0,0 +1,301 @@
+classdef CExp < handle
+% class of experiment
+% store experiment data, structure, procedure etc
+% 10 July, 2011
+% Created by: Z. Shi, shi@lmu.de
+% ver. 0.2
+% 21 June, 2012
+%   add option: Continue. If will continue run experiment with given subject name
+%   add additional Data variable
+%
+%constructor: two ways to initialize: 
+    % constant method: 
+    %   p1, repetition, 
+    %   p2, factors
+    % adapative method:
+    %   p1, alpha 
+    %   p2, beta
+    %optional parameters: blockReptition, blockFactors;  
+    % adaptive method:maxTrials,maxRespRev, gamma, lambda
+
+    properties
+        sName = 'test'; %subject information
+        sAge = 20;
+        sGender = 'f';
+        sPara; %additional parameter
+        eType = 1; % 1 - constant stimuli; 2- bayesian adaptive method        
+        % type 1
+        seq;    % sequences of the trials
+        resp;   % corresponding responses
+        aData;  % store additional data, can be anything...
+        % only for type 2
+        beta;
+        alpha;
+        gamma;
+        lambda;
+        maxTrls; %maximum trial number
+        maxRev; % maximum reversal number
+        % index
+        curTrl; %current trial index
+        curIntensity;
+        practice;
+        % other information
+        startTime;
+        endTime;
+     end
+    
+    methods
+        function obj = CExp(p1,p2,varargin)
+            % constant method: 
+            %   p1, repetition, 
+            %   p2, factors
+            % adapative method:
+            %   p1, alpha 
+            %   p2, beta
+            %optional parameters: blockRepetition, blockFactors;  
+            % adaptive method:maxTrials,maxRespRev, gamma, lambda
+            p = inputParser;
+            p.addOptional('eType','c', @(x) any(strcmpi(x,{'c','a'})));
+            p.addParamValue('blockRepetition',1,@(x) x>0);
+            p.addParamValue('blockFactors', [], @isnumeric);
+            p.addParamValue('maxTrials', 40, @isnumeric);
+            p.addParamValue('maxRespRev',8, @isnumeric);
+            p.addParamValue('gamma',0.025,@(x) min(x)>=0);
+            p.addParamValue('lambda',0.025, @(x) x>=0);
+            p.addParamValue('practice',5, @(x) x>=0);
+            p.parse(varargin{:});
+            
+            
+            if p.Results.eType == 'c' %constant method
+                obj.eType = 1;
+                obj.curTrl = 1;
+                obj.seq = obj.genTrials(p1,p2,p.Results.blockRepetition,p.Results.blockFactors);
+                obj.resp = [];
+                obj.maxTrls = length(obj.seq);
+                
+            else % adaptive method
+                obj.alpha = p1;
+                obj.beta = p2;
+                obj.gamma = p.Results.gamma;
+                obj.lambda = p.Results.lambda;
+                obj.maxTrls = p.Results.maxTrials;
+                obj.maxRev = p.Results.maxRespRev;
+                obj.practice = p.Results.practice;
+                obj.eType = 2;
+                obj.curTrl = 1;
+                obj.seq = zeros(obj.maxTrls,2);
+                obj.resp = [];
+            end
+        end
+        
+        function obj = guessThreshold(obj,x)
+            obj.seq(1,1:2) = [x 0];
+            obj.curIntensity = x;
+        end
+        
+        function obj = setResp(obj,resp)
+            obj.resp(obj.curTrl,1:size(resp,2)) = resp;
+            obj.curTrl = obj.curTrl+1;
+        end
+        
+        function curSeq = getCondition(obj)
+            curSeq = obj.seq(obj.curTrl,:);
+        end
+        
+        function obj = updateThreshold(obj, p_target)
+            % logistic function
+            % p = Logistic(x, alpha, beta, gamma, lambda)
+            %       alpha - threshold
+            %       beta - slope
+            %       gamma - chance performance / fa
+            %       lambda - lapsing rate
+            % Based on MLP toolbox
+            % updated intensity and fa are stored in obj.seq
+            if obj.curTrl < obj.practice
+                % practice
+                ra = randperm(length(obj.alpha));
+                obj.seq(obj.curTrl,1)= obj.alpha(ra(1));
+                obj.seq(obj.curTrl,2) = obj.gamma(1);
+            else
+                ll=zeros(length(obj.alpha), 1);
+                x = obj.seq(obj.practice:max(obj.curTrl-1,1),1);
+                responses = obj.resp(obj.practice:max(obj.curTrl-1,1));
+
+                % calculate the likelihood of each psychometric function
+                for i=1:length(obj.alpha)
+                    for j=1:length(obj.gamma)
+                        ll(i, j)=CalculateLikelihood(obj,x, responses, obj.alpha(i), obj.gamma(j));
+                    end
+                end
+
+                % find the most likely psychometric function
+                [i, j]=find(ll==max(max(ll)));
+                if length(i)+length(j) > 2
+                    i = i(1);
+                    j = j(1);
+                end;
+                % calculate the level of the stimulus at p_target performance
+                % using inverse logistic function
+                obj.curIntensity = obj.alpha(i)-(1/obj.beta)*log(((1-obj.lambda-obj.gamma(j))./(p_target-obj.gamma(j)))-1);
+                obj.seq(obj.curTrl,1)= obj.curIntensity;
+                obj.seq(obj.curTrl,2) = obj.gamma(j);
+            end
+        end
+
+        function bFinish = canStop(obj)
+            if obj.eType == 2
+                bFinish = 0;
+                revnum = sum(abs(diff(obj.resp(obj.practice:obj.curTrl-1,1))));
+                if revnum > obj.maxRev
+                    bFinish = 1;
+                end
+                if obj.curTrl > obj.maxTrls
+                    bFinish = 1;
+                end
+            else
+                disp('This is only available for adaptive method');
+            end
+        end
+        
+        function ll=CalculateLikelihood(obj,x, responses, alpha, gamma)
+            if obj.eType == 2
+                warning off
+                ll = 0;
+                % calculate logistic probablity
+                p=gamma+((1-obj.lambda-gamma).*(1./(1+exp(obj.beta.*(alpha-x)))));
+
+                ll = sum(log(p(responses ==1)))+ sum(log(1-p(responses == 0)));
+            else
+                disp('This is only available for adaptive method');
+            end
+        end
+        
+        function h = plotLogistic(obj, a, g)
+            %plot a logistic function for specify parameter, only available
+            %for adaptive method
+            if obj.eType == 2
+                x = obj.alpha;
+                if nargin < 2
+                    a = median(obj.alpha);
+                    g = median(obj.gamma);
+                end
+                y = g+((1-obj.lambda-g).*(1./(1+exp(obj.beta.*(a-x)))));
+                h = figure; 
+                plot(x,y);
+            else
+                disp('This is only available for adaptive method');
+                h = -1;
+            end
+        end
+        
+        function trials=genTrials(obj,withinblock_repetitions, withinblock_factors, betweenblock_repetitions, betweenblock_factors)
+            % syntax: genTrials(withinblock_repetition, betweenblock_repetition, withinblock_factors, [ betweenblock_factors])
+            % eg: genTrials(2, [2 3],10); generate two factors (2 levels and 3 levels resp. ) with within block repetition 2 and between block repetition 10
+            % coded by: strongway
+            rand('seed',sum(clock*100));
+            if  nargin < 3
+                    error('Incorrect number of arguments for genTrials function');
+            elseif nargin == 3
+                    betweenblock_repetitions = 1;
+                    betweenblock_factors  = [];
+            elseif nargin == 4
+                % when there's only betweenblock_repetitions, which is equal to swap
+                % between block repetitions and factors -strongway 14. Dec. 2006
+                    betweenblock_factors = betweenblock_repetitions;
+                    betweenblock_repetitions  = 1;
+            end
+
+            trials = [];
+            block_design = [];
+            numblock = betweenblock_repetitions;
+            if ~isempty(betweenblock_factors)
+                block_design = fullfact(betweenblock_factors);
+                block_design = repmat(block_design, betweenblock_repetitions,1);
+                idxb = randperm(length(block_design));
+                block_design = block_design(idxb,:);
+                numblock = length(block_design);
+            end
+
+            for iblock = 1:numblock
+                %generate within block trials
+                inblock_trials = fullfact(withinblock_factors);
+                inblock_trials = repmat(inblock_trials,withinblock_repetitions,1);
+                idx=randperm(size(inblock_trials,1));
+                inblock_trials = inblock_trials(idx,:);
+                if ~isempty(block_design)
+                    %add between factors
+                    blockwise_factors = repmat(block_design(iblock,:),length(inblock_trials),1);
+                    inblock_trials = [inblock_trials blockwise_factors];
+                end
+                trials = [trials; inblock_trials];
+            end
+        end
+
+        function obj = subInfo(obj,sp1,sp2)
+            % Get Subject information, sp1 and sp2 for additional parameters caption and parameters.
+            promptParameters = {'Subject Name', 'Age', 'Gender (F or M?)'};
+            defaultParameters = {'test', '20','F'};
+            if nargin == 2
+                promptParameters = [promptParameters, sp1];
+                defaultParameters = {'test', '20','F','[]'};
+            end
+            if nargin == 3
+                promptParameters = [promptParameters, sp1];
+                defaultParameters = {'test', '20','F',sp2};
+            end
+            sub = inputdlg(promptParameters, 'Subject Info  ', 1, defaultParameters); 
+            obj.sName = sub{1};
+            obj.sAge = eval(sub{2});
+            obj.sGender = sub{3};
+            if nargin >= 2
+                obj.sPara = eval(sub{4});
+            end
+            
+            obj.startTime = now;
+            % if it is continue, read the data
+            if nargin == 3 && strcmp(sp1,'continue') 
+                if obj.sPara == 1
+                    subfile = ['data' filesep obj.sName '.mat'];
+                    if ~exist(subfile)
+                        error('No such data file');
+                    else
+                        load(subfile); % trials, expInfo, seq, resp
+                        %restore the data
+                        obj.seq = seq;
+                        obj.resp = resp;
+                        obj.aData = aData;
+                        obj.curTrl = length(obj.resp)+1;
+                        clear trials expInfo aData seq resp; 
+                    end
+                else % new subject, check if there is already a file
+                    subfile = ['data' filesep obj.sName '.mat'];
+                    if exist(subfile)
+                        error('Data file already exist');
+                    end
+                end
+            end
+        end
+        
+        function saveData(obj,filename)
+            if nargin<2
+                filename = obj.sName;
+            end
+            obj.endTime = now;
+            if ~isdir('data')
+                mkdir('data'); %create a directory for storing files
+            end
+            seq = obj.seq;
+            resp = obj.resp;
+            trials = [seq(1:length(obj.resp),:), resp];
+            expInfo.name = obj.sName;
+            expInfo.age = obj.sAge;
+            expInfo.sex = obj.sGender;
+            expInfo.para = obj.sPara;
+            expInfo.time = [obj.startTime, obj.endTime];
+            aData = obj.aData;
+            save(['data' filesep filename],'trials','expInfo','aData','seq','resp');
+            
+        end
+        
+    end
+end

experiments/CInput.m

@@ -0,0 +1,270 @@
+classdef CInput < handle
+    % A class of peripheral input device
+    % Input device:
+    % 1. keyboard: CInput('k',response values, {KeyName List}). Keys are
+    % mapped into response values. For example, 
+    %   kb = CInput('k',[1 2],{'LeftArrow','RightArrow'); % Left arrow 1,   right arrow 2
+    %   by default: kb=CInput('k'); % Left and Right Arrow
+    % 2. Mouse: CInput('m', [1 2]);% left/right mouse button to 1 ,2
+    % 3. Parallel port key: CInput('p', [1 2]);
+    %       note: parallel keypad encoding depends on pins, so be checked
+    %       before use this method
+    % 4. NI input: CInput('n',defaultValues, pins, portname);
+    %       e.g. nk = CInput('n', [0 0], [8 9],'AOut');
+    %           Card 'AOut', pin 8 and 9, default values [0 0]
+    % Methods:
+    % 1. wait(); %wait for input
+    % 2. [kinput ktime x y] = response(bWait); %by default, wait for response. 
+    %   when bWait == 0, it samples reponse values only. 
+    % 3. releaseTime = keyRelease();
+    %      return releaseTime
+    % last modify: 15. March, 2012
+    % First created: 2nd April, 2001
+    % created by: Z. Shi, shi@lmu.de
+    properties
+        esckey = 41;
+        key = [80 79];
+        kval = [1 2];
+        ktype = 1;
+        wantStop = 0;
+        dio ;
+    end
+    
+    methods 
+        function obj = CInput(keytype,keyvalues,keys,portname)
+            try
+                KbName('UnifyKeyNames');
+                obj.esckey = KbName('ESCAPE');
+                switch lower(keytype)
+                    case {'keyboard','k'}
+                        obj.ktype = 1;
+                        if nargin < 3
+                            obj.key = [KbName('LeftArrow') KbName('RightArrow')];
+                        end
+                        if nargin < 2
+                            keyvalues = [1 2];
+                        end
+                        if nargin == 3
+                            for k = 1:length(keys)
+                                obj.key(k) = KbName(keys{k});
+                            end
+                        end
+                            obj.kval= keyvalues;
+                    case {'mouse','m'}
+                        obj.ktype = 2;
+                        if nargin < 2
+                            keyvalues = [1 2];
+                        end
+                        obj.kval = keyvalues;
+
+                    case {'parallel','p'}
+                        obj.ktype = 3;
+                        % todo
+                        obj.key = keys; % actually parallel port pins
+                        obj.kval = keyvalues;
+                        obj.dio = digitalio('parallel', 'LPT1');
+                        addline(obj.dio,obj.key,'in');
+
+                    case {'n','nidaq'}
+                        obj.ktype = 4;
+                        obj.key = keys; % pins
+                        obj.kval = keyvalues; % default key values
+                        obj.dio = digitalio('nidaq',portname);
+                        addline(obj.dio,obj.key,'in');
+                        
+                    otherwise % default as keyboard
+                        obj.ktype = 1;
+                        obj.key(1) = KbName('LeftArrow');
+                        obj.key(2) = KbName('RightArrow');
+                        obj.kval = [1 0];
+                end
+            catch ME
+                disp(ME.message);
+            end
+        end
+        
+        function wait(obj)
+            % halt program for key input 
+            switch obj.ktype
+                case 1 %keyboard
+                    while KbCheck; end %clear buffer
+                    KbWait;
+                case 2 % mouse
+                    GetClicks;
+                case 3 % parallel port keypad
+                    % need to test for different parallel port pad
+                    %here assume: initial status: 1111 ...
+                    k = getvalue(obj.dio);
+                    while sum(k) == length(k)
+                        k = getvalue(obj.dio);
+                    end
+                case 4 %ni daq
+                    k = getvalue(obj.dio);
+                    while sum(xor(obj.kval,k)) == 0 % sum(k) == length(k)
+                        k = getvalue(obj.dio);
+                    end                    
+            end
+        end
+        
+        %
+        function [kinput ktime x y] = response(obj,bWaitPress)
+            % wait for response (by default, it waits for response)
+            % if bWaitPress == 0, it samples response values
+            % return input keys, time, and x, y
+            x = -1;
+            y = -1;
+            if nargin < 2
+                bWaitPress = 1;
+            end
+            validKey = 0;
+            kinput = -1;
+            ktime = GetSecs;
+            switch obj.ktype
+                case 1 %keyboard
+                    while 1
+                        [ keyIsDown, seconds, keyCode ] = KbCheck;
+                        if keyIsDown
+                            if keyCode(obj.esckey)
+                                kinput = -1;
+                                ktime = GetSecs;
+                                obj.wantStop = 1;
+                                break;
+                            end
+                            
+                            for k = 1: length(obj.key)
+                                if keyCode(obj.key(k))
+                                    kinput = obj.kval(k);
+                                    ktime = GetSecs;
+                                    validKey = 1;
+                                end
+                            end
+                            if validKey
+                                break;
+                            end
+                            while KbCheck; end %clear buffer
+                        end
+                        if ~bWaitPress
+                            break;
+                        end
+                    end
+                case 2 % mouse
+                    while 1
+                      [x,y,buttons] = GetMouse;
+                      if buttons(1)
+                          kinput = obj.kval(1);
+                          ktime = GetSecs;
+                          break;
+                      end
+                      if length(buttons)>=2 & buttons(end) 
+                          kinput = obj.kval(2);
+                          ktime = GetSecs;
+                          break;
+                      end
+                      % stop running
+                      [ keyIsDown, seconds, keyCode ] = KbCheck;
+                        if keyIsDown && keyCode(obj.esckey)
+                            kinput = -1;
+                            ktime = GetSecs;
+                            obj.wantStop = 1;
+                            break;
+                        end
+                        if ~bWaitPress
+                            break;
+                        end
+                    end 
+                case 3 % parallel port keypad
+                    k = getvalue(obj.dio);
+%                    while sum(k) == length(k)
+                    while sum(xor(obj.kval,k)) == 0
+                        k = getvalue(obj.dio);
+                        if ~bWaitPress
+                            break;
+                        end
+                      % stop running
+                      [ keyIsDown, seconds, keyCode ] = KbCheck;
+                        if keyIsDown && keyCode(obj.esckey)
+                            kinput = -1;
+                            ktime = GetSecs;
+                            obj.wantStop = 1;
+                            break;
+                        end
+                    end
+%                    kidx = find(k == 0);
+                    ktime = GetSecs;
+%                    kinput = obj.kval(kidx(1));
+                    kinput = sum(find(xor(obj.kval,k)));
+                case 4 %NI
+                    k = getvalue(obj.dio);
+                    while sum(xor(obj.kval,k)) == 0 %sum(k) == length(k)
+                        k = getvalue(obj.dio);
+                        if ~bWaitPress
+                            break;
+                        end
+                    end
+                  % stop running
+                  [ keyIsDown, seconds, keyCode ] = KbCheck;
+                    if keyIsDown && keyCode(obj.esckey)
+                        kinput = -1;
+                        ktime = GetSecs;
+                        obj.wantStop = 1;
+                    end
+                    %kidx = find(k == 0);
+                    ktime = GetSecs;
+                    kinput = sum(find(xor(obj.kval,k)));
+                    %kinput = obj.kval(kidx(1));
+            end % end of switch
+
+            
+        end % end of function response
+ 
+        function releaseTime = keyRelease(obj)
+         %return release time  when key is released,    
+            switch obj.ktype
+                case 1 %keyboard
+                    while 1
+                        [ keyIsDown, seconds, keyCode ] = KbCheck;
+                        if ~keyIsDown
+                                break;
+                        end
+                    end
+                case 2 % mouse
+                    while 1
+                      [x,y,buttons] = GetMouse;
+                      if sum(buttons) == 0
+                          break;
+                      end
+                    end
+                case 3 % parallel port keypad
+                    k = getvalue(obj.dio);
+                    while sum(k) < length(k)
+                        k = getvalue(obj.dio);
+                    end
+                case 4
+                    %keys are too sensitive, we need to apply filter mechanisms
+                    noise_ticks = 0;
+                    while 1
+                        k = getvalue(obj.dio);
+                        if sum(xor(obj.kval,k))==0 % key releases
+                            noise_ticks = noise_ticks +1;
+                        else 
+                            noise_ticks = 0;
+                        end
+                        if noise_ticks >10 % 10 ms key releases
+                            break;
+                        end
+                        WaitSecs(0.001);
+                    end
+            end % end of switch
+            releaseTime = GetSecs-0.01; %- 10 ms
+            
+            % stop running
+            [ keyIsDown, seconds, keyCode ] = KbCheck;
+            if keyIsDown && keyCode(obj.esckey)
+                obj.wantStop = 1;
+            end
+            
+        end %end of function keyResease
+        
+        
+    end
+end

experiments/ReadMe.txt

@@ -0,0 +1,4 @@
+Date: 25.07.2018
+
+This is the experimental folder, which contains the necessary matlab files for running the experiment. 
+The main program called 'main.m'. Run this file will start the experiment. 

experiments/data/data_ana.R

@@ -0,0 +1,22 @@
+library(tidyverse)
+library(R.matlab)
+
+res = readMat('test_RS04.mat')
+
+dat = as.tibble(res$trials)
+names(dat) = c('Duration', 'session','trlNo','blkNo','dur1','pdur','production','vrep','Reproduction')
+
+
+mrep =  dat %>% group_by(session, Duration) %>%
+  summarise(mRep = mean(Reproduction), n = n(), se = sd(Reproduction)/sqrt(n-1)) 
+
+mrep$session = factor(mrep$session, labels = c("RWalk", "Random"))
+
+mrep%>%
+  ggplot(aes(Duration, mRep, color = session, group = session)) + 
+  geom_point() + geom_line() + 
+  geom_errorbar(aes(ymin = mRep - se, ymax = mRep +se), width = 0.05) + 
+  theme_classic() + 
+  geom_abline(slope = 1)
+
+ggplot(dat, aes(trlNo, Duration)) + geom_line() + geom_point()

experiments/main.m

@@ -0,0 +1,197 @@
+ function main
+% This experiment is to examine the difference between ASP and TN groups
+% using reproduction and central tendency effect. 
+% coded by Strongway (shi@lmu.de)
+% date: 25th July., 2018
+
+% experiment related parameters
+para.viewDistance = 57; % viewing distance 57 cm
+para.monitor = 22; % monitor size
+para.fntSize = 24; % font size
+para.bkColor = [128,128,128]; % background color
+para.fColorCircle = [128,128,0]; % foreground color
+para.fColorRectangle = [0,200,200]; % white color
+para.fColorW = [200];
+para.green = [0, 192, 0];
+para.red = [192, 0, 0];
+
+para.xPosition = 0; % 0 degree above
+para.yPosition = 0; % 0 degree above
+para.iti = [0.8, 1]; % range of inter-trial interval
+para.iprp = 0.5; % inter production-reproduction interval
+para.xyFeedbackArray = [-2,0; -1, 0; 0, 0; 1, 0; 2,  0]; % location of feedback array
+para.fbRange = [-100, -0.3; -0.3, -0.1; -0.1, 0.1; 0.1, 0.3; 0.3, 100]; %feedback range with respect to the reproduction error
+para.withFeedback = true;
+para.vSize = 10; % 10 degree
+para.vFeedDivid = 10;
+
+nTrlsBlk = 25;
+numHalfTrial = 250;
+%[para.nDurations(1,:) ,  para.nDurations(2,:)] = genDurSeq(numHalfTrial, 1);
+%[withblkFactors, nDurs] = size(para.nDurations);
+%sessionType = withblkFactors;% 1: short matched to circle and long matched to square
+%inBlkRep = 1; %inblock repetition
+try
+%    kb = CInput('k',[1], {'downArrow'});
+    kb = CInput('m'); % replace with mouse
+    myexp = CExp(1, numHalfTrial, 'blockFactors',  2,...
+        'blockRepetition',1);
+    % acquire subject information
+    myexp.subInfo('Sequence', '1');
+    % load sequence w1 from subfolder 'seqs'
+    load(['seqs', filesep, 'seq', num2str(myexp.sPara)],'w1');
+   % Dur, session, ntrl, nblock
+    %production duration, session (1= predictive sequence, 2 = non predictive sequence) 
+   myexp.seq( myexp.seq(:, 2) == 1, 1) = w1;
+   myexp.seq( myexp.seq(:, 2) == 2, 1) = w1(randperm( length(w1) ) ); % randomized sequence
+    % set second column ad block number, 4th colum as block information
+    myexp.seq(:,3) = 1: myexp.maxTrls; % trial sequence no. 
+    myexp.seq(:,4) = floor((myexp.seq(:,3)-1)/nTrlsBlk) + 1;
+    
+    v = CDisplay('bgColor',para.bkColor,'fontSize',para.fntSize,'monitorSize',para.monitor,...
+        'viewDistance',para.viewDistance,'fullWindow',1, 'skipSync',1);
+    HideCursor;
+    
+    % create stimuli
+    para.vObjCircle = v.createShape('circle',para.vSize,para.vSize,'color',para.fColorCircle);
+    para.vObjRectangle = v.createShape('rectangle',para.vSize,para.vSize,'color',para.fColorRectangle);
+
+    % create feedback stimuli
+    vGreenDisk = v.createShape('circle', para.vSize/para.vFeedDivid, para.vSize/para.vFeedDivid, 'color',para.green);
+    vRedDisk = v.createShape('circle', para.vSize/para.vFeedDivid, para.vSize/para.vFeedDivid, 'color',para.red);
+    vDiskFrame = v.createShape('circle',para.vSize/para.vFeedDivid, para.vSize/para.vFeedDivid, 'color',para.fColorW,'fill',0);
+    para.vFullFrames = [vDiskFrame, vDiskFrame, vDiskFrame, vDiskFrame, vDiskFrame];
+    para.vFullDisks = [vRedDisk, vGreenDisk, vGreenDisk, vGreenDisk, vRedDisk];
+    %initialize text
+    infoText = init_text;
+    
+    % start instruction
+    v.dispText(infoText.instruction);
+    kb.wait;
+    WaitSecs(2);
+    
+    for iTrl = 1:myexp.maxTrls
+        %get current condition
+        cond = myexp.getCondition;  % duration, range
+        curDuration = cond(1); % current standard
+        
+        %start trial presentation
+        results = trialPresentation(v,kb, cond, curDuration, para);
+        %store results
+        myexp.setResp(results);%
+        % debugging
+        if kb.wantStop
+            break;
+        end 
+        % ITI 0.8-1 seconds
+        WaitSecs(para.iti(1) + (para.iti(2)-para.iti(1))*rand);
+        if  mod(iTrl,nTrlsBlk) == 0 
+             v.dispText(infoText.startBlock);
+             kb.wait;
+        end
+        
+        if  iTrl == numHalfTrial
+            myexp.saveData;   %save data
+            v.dispText(infoText.sessionBreak);
+            kb.wait;
+            WaitSecs(2);
+            kb.wait;
+        end
+        
+    end
+    myexp.saveData;   %save data
+    v.dispText(infoText.thankyou);
+    kb.wait;
+    v.close;
+    ShowCursor; 
+  catch ME
+    % debugging
+    v.close;
+  
+    disp(ME.message);
+    disp(ME.stack);
+    for i=1:length(ME.stack)
+        disp(ME.stack(i).name);
+        disp(ME.stack(i).line);
+    end
+    v.close;
+    ShowCursor;
+end
+end
+
+
+function results = trialPresentation(v, kb, cond, curDuration, para)
+
+        v.dispFixation(20);
+        WaitSecs(0.500); %at least 500 ms
+
+        % initiated by key pressing, and measured by key release
+        v.dispFixation(5,2); % change fixation from cross to a dot 
+        [key, keyInitTime] = kb.response;
+        dispObj = para.vObjCircle;
+        v.dispItems([para.xPosition, para.yPosition],  dispObj,[para.vSize para.vSize],0,0); % draw texture    
+        [vbl, vInitTime] = v.flip;
+        WaitSecs(curDuration - 0.005); % 5 ms earlier, so make sure to clear next frame on time.
+        [vbl, vStopTime] = v.flip(1); 
+  
+        keyReleaseTime = kb.keyRelease;
+        v.flip(1); %clear screen
+
+        phyDuration = vStopTime - vInitTime;        %visual duration
+        proDuration = keyReleaseTime - keyInitTime; %key production
+
+        % reproduction
+        WaitSecs(para.iprp); %wait at least 250 ms
+        v.dispFixation(5,2);
+        
+        [key, keyInitTime] = kb.response;
+        v.dispItems([para.xPosition, para.yPosition],  dispObj,[para.vSize para.vSize],0,0);
+        [vbl, vInitTime] = v.flip;
+        keyReleaseTime = kb.keyRelease;
+
+        [vbl, vStopTime]=v.flip(1);
+        repDuration = keyReleaseTime - keyInitTime; % key reproduction
+        repVDuration = vStopTime - vInitTime; % visual reproduction
+       
+        %store results
+        results = [curDuration, phyDuration, proDuration, repVDuration, repDuration];
+   
+        if para.withFeedback
+            % present a feedback display
+            feedbackDisplay = para.vFullFrames;
+            delta = (repDuration - phyDuration)/phyDuration;
+            % find the range of the error
+            cIdx = para.fbRange > delta; % column index of left and right boundary
+            idx = find(xor(cIdx(:,1),cIdx(:,2)));
+            feedbackDisplay(idx(1)) = para.vFullDisks(idx(1));
+
+            WaitSecs(0.25); % wait 250 ms
+            v.dispItems(para.xyFeedbackArray, feedbackDisplay,[para.vSize/para.vFeedDivid para.vSize/para.vFeedDivid]); % draw texture 
+            WaitSecs(0.500); % display the feedback for 500 ms
+            v.flip;  
+        end
+        
+end
+
+function infoText = init_text
+    % specify experimental text
+    infoText.instruction = [ 'Erklärung \n\n',...
+        'Die folgende Sitzung starten Sie durch Drücken der linken Maustaste. ', ...
+        'Bitte fixieren Sie das Kreuz in der Mitte des Bildschirms. ' ...
+        'Nach einer sehr kurzen Zeit wird das Kreuz zu einem Punkt. Nun kann das Experiment starten. ', ...
+        'Halten Sie nun die Maustaste solange gedrückt wie Sie den gelben Kreis sehen. Bitte achten Sie genau auf die Dauer. ', ...
+        'Der Kreis verschwindet automatisch. Lassen Sie die Maustaste wieder los sobald der Kreis verschwindet. ' ...
+        'Nachdem Sie nun wieder den kleinen Punkt sehen sollen Sie die Dauer des gelben Kreises reproduzieren. ' ...
+        'Hierfür halten Sie die Maustaste solange gedrückt wie Sie die Dauer eingeschätzt haben. ' ...
+        'Der gelbe Kreis erscheint dabei als Hilfsmittel. \n'];
+
+    infoText.blockInfo = 'Bitte machen Sie eine Pause. Wenn Sie bereit sind drücken Sie die Maustaste. \n ';
+    infoText.endTrial = 'Bitte lassen Sie die Maustaste los';
+    infoText.production = '+';
+    infoText.reproduction = '+';
+    infoText.sessionBreak = ' Die erste Sitzung ist beendet! Bitte öffnen Sie die Tür und machen eine kurze Pause! \n';
+    infoText.startBlock = ' Bitte machen Sie eine Pause! \n\n Bitte drücken Sie die Maustaste um den Block zu starten.\n';
+    infoText.goingon = ' Bitte machen Sie eine Pause! \n\n Bitte drücken Sie die Maustaste um den Block fortzufahren.';
+    infoText.thankyou = 'Das Experiment ist beendet! \nVielen herzlichen Dank!';
+
+end

experiments/seqs/genDurSeq.m

@@ -0,0 +1,23 @@
+function [w1, w2] = genDurSeq(nTrl, plotSign)
+
+% first generate the random walk sequence
+w1 = cumsum(randn(1, nTrl));
+w1 = (w1 - mean(w1))/std(w1); % normalize
+% appr. scale to 0.5 to 2
+w1 = w1/(max(w1)-min(w1))*1.6 + 1;
+
+% discretize to 100 ms
+w1 = round(w1*10)/10;
+
+% another sequence with randomization
+w2 = w1(randperm( length(w1) ) );
+
+if plotSign == 1
+    figure(); hold on;  subplot(2,2,1);plot(w1)
+    subplot(2,2,2);hist(w1);
+    subplot(2,2,3);plot(w2)
+    subplot(2,2,4);hist(w2);
+end
+% 
+%estHighHurst = estimate_hurst_exponent(w1)
+%estLowHurst = estimate_hurst_exponent(w2)

experiments/seqs/saveSeq.m

@@ -0,0 +1,18 @@
+function saveSeq()
+% to keey the ASD and control group equal, we save the randomwalk sequence
+% first
+for i= 1:60
+    w1 = genDurSeq(250,1);
+    save(['seq', num2str(i)], 'w1');
+end
+
+%%
+w = [];
+for i=1:27
+   load(['seq' num2str(i)]);
+   figure;
+   plot(w1);
+   w = [w, w1];
+end
+figure;
+hist(w);

volatility-reproduction.Rproj

@@ -0,0 +1,13 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: XeLaTeX