rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Weber's Law
# Check predictions of analysis if there is 10% effect for each illusion
# Greg Francis
# PSY 646
# 12 December 2016


# Generate data frame with the various conditions
Subject<-c()
TargetLength <-c()
WingType <-c()
MatchOrientation <-c()
MatchLength<-c()
SDPhysical<-c()
SDPerceptual<-c()

for(i in c(1:50)){
	for(TL in c(100, 150)){
	for(Wings in c("Out", "In", "None")){
	for(MO in c("Vertical", "Horizontal")){	
		Subject<-c(Subject, i)
		TargetLength<-c(TargetLength, TL)
		WingType<-c(WingType, Wings)
		MatchOrientation<-c(MatchOrientation, MO)
		
		# Compute predicted match length for this condition
		predictedMatch = TL
		# Vertical-horizontal illusion (vertical will have smaller match size because it looks bigger than horizontal target)
		if(MO =="Vertical"){
			predictedMatch = predictedMatch * 0.9
		}
		# Outward wings (will have smaller match size because it looks bigger than horizontal target)
		if(Wings =="Out"){
			predictedMatch = predictedMatch * 0.9
		}
		# Inward wings (will have bigger match size because it looks smaller than horizontal target)
		if(Wings == "In"){
			predictedMatch = predictedMatch * 1.1
		}	
		MatchLength <-c(MatchLength, predictedMatch)
		
		# Compute SD using Weber's law for physical 
		SDPhysical <- c(SDPhysical, rnorm(1, 6+0.6725*0.06*predictedMatch, 4.2) )  # Numbers estimated from pilot study
		# Compute SD using Weber's law for perceptual 
		SDPerceptual <-c(SDPerceptual, rnorm(1, 6+0.6725*0.06*TL, 4.2) )  # Numbers estimated from pilot study	
	}
	}
	}
}
		
WLdata =data.frame(Participant=Subject, Wings = WingType, MatchOrientation=MatchOrientation, TargetLength=TargetLength, MeanMatch=MatchLength, SDPhysical=SDPhysical, SDPerceptual=SDPerceptual)		

																			
# load the rethinking library
library(rethinking)
	
# Physical length data generation
PhysicalData1 = data.frame(Participant = WLdata$Participant, SD=WLdata$SDPhysical,  MeanMatch=WLdata$MeanMatch)	
										
PhysicalLengthLinearModel <- ulam(
			alist( SD ~ dnorm(mu, sigma), 
			mu <- a[Participant] + b[Participant]*MeanMatch,
			a[Participant] ~ dnorm(grand_mua, grand_sa),
			grand_mua ~ dnorm(5, 5),
			grand_sa ~ dcauchy(0, 5),
			b[Participant] ~ dnorm(grand_mub, grand_s),
			grand_mub ~ dnorm(5, 5),
			grand_s ~ dcauchy(0, 5),
			sigma ~ dcauchy(0, 5)
			), data= PhysicalData1, iter=10000, warmup=4000, log_lik=TRUE, chains=4, cores=4 ) 

print("Finished physical model on physical data\n")					

PhysicalData2 = data.frame(Participant = WLdata$Participant, SD=WLdata$SDPhysical, TargetLength = WLdata$TargetLength)	
													
PerceptualMatchLinearModel <- ulam(
			alist( SD ~ dnorm(mu, sigma), 
			mu <- a[Participant] + b[Participant]* TargetLength,
			a[Participant] ~ dnorm(grand_mua, grand_sa),
			grand_mua ~ dnorm(5, 5),
			grand_sa ~ dcauchy(0, 5),
			b[Participant] ~ dnorm(grand_mub, grand_s),
			grand_mub ~ dnorm(5, 5),
			grand_s ~ dcauchy(0, 5),
			sigma ~ dcauchy(0, 5)
			), data= PhysicalData2, iter=10000, warmup=4000, log_lik=TRUE, chains=4, cores=4 ) 
			
print("Finished perceptual model on physical data\n")					

print("Comparing physical and perceptual models for data generated with a physical model.\n")
print(compare(PhysicalLengthLinearModel, PerceptualMatchLinearModel))
	
# Perceptual length data generation
PerceptualData1 = data.frame(Participant = WLdata$Participant, SD=WLdata$SDPerceptual, MeanMatch=WLdata$MeanMatch)	

PhysicalLengthLinearModel2 <- ulam(
			alist( SD ~ dnorm(mu, sigma), 
			mu <- a[Participant] + b[Participant]*MeanMatch,
			a[Participant] ~ dnorm(grand_mua, grand_sa),
			grand_mua ~ dnorm(5, 5),
			grand_sa ~ dcauchy(0, 5),
			b[Participant] ~ dnorm(grand_mub, grand_s),
			grand_mub ~ dnorm(5, 5),
			grand_s ~ dcauchy(0, 5),
			sigma ~ dcauchy(0, 5)
			), data= PerceptualData1, iter=10000, warmup=4000, log_lik=TRUE, chains=4, cores=4 ) 

print("Finished physical model on perceptual data\n")					

				
PerceptualData2 = data.frame(Participant = WLdata$Participant, SD=WLdata$SDPerceptual, TargetLength = WLdata$TargetLength)	
										
	
PerceptualMatchLinearModel2 <- ulam(
			alist( SD ~ dnorm(mu, sigma), 
			mu <- a[Participant] + b[Participant]* TargetLength,
			a[Participant] ~ dnorm(grand_mua, grand_sa),
			grand_mua ~ dnorm(5, 5),
			grand_sa ~ dcauchy(0, 5),
			b[Participant] ~ dnorm(grand_mub, grand_s),
			grand_mub ~ dnorm(5, 5),
			grand_s ~ dcauchy(0, 5),
			sigma ~ dcauchy(0, 5)
			), data= PerceptualData2, iter=10000, warmup=4000, log_lik=TRUE, chains=4, cores=4 ) 
			
print("Finished perceptual model on perceptual data\n")					

	

print("Comparing physical and perceptual models for data generated with a perceptual model.\n")
print(compare(PhysicalLengthLinearModel2, PerceptualMatchLinearModel2))