rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Greg Francis
# PSY 626
# 17 November 2025

verbose=TRUE
DesiredPower=0.9

# Prior of standardized effect size

dPrior = rcauchy(10000, location=0, scale=sqrt(2)/2)
hist(dPrior[abs(dPrior)<20], breaks=20)
				
# Bayes Factor 

library(BayesFactor)

EvidenceThreshold=3

for(n in seq(8,1000)){
	CountSufficientEvidence=0
	numReps=1000
	# Generate simulated data sets
	for(rep in seq(1, numReps)){
		# Pick a random standardized effect size from prior
		d= sample(dPrior, 1)
		Group1<- rnorm(n, mean=0, sd=1)
		Group2<- rnorm(n, mean=d, sd=1)
		
		test = ttestBF(x = Group1, y = Group2)
		
		bf = extractBF(test, logbf=FALSE, onlybf=TRUE)
		
		if(bf>EvidenceThreshold){
			CountSufficientEvidence = CountSufficientEvidence +1
		}
	
	}
	powerEstimate = CountSufficientEvidence/numReps
	if(verbose){
		cat('With n=', n, ' Bayesian power for threshold based on prior', EvidenceThreshold,' is ', powerEstimate, '\n')
	}
	
	if(powerEstimate>= DesiredPower){
		cat('Bayes Factor, smallest sample size with desired power based on prior is ', n, ' with power of ', powerEstimate, '\n')
		break
	}
}