fit_uvm_random.R

library(CircStats)

#define a data generator for the uvm model
get_mix_data = function( n , rho , kappa , seed ){
	set.seed(seed)
	to_return = rep(NA,n)
	trial_is_vm = as.logical(rbinom(n,1,rho))
	to_return[trial_is_vm] = rvm(sum(trial_is_vm),pi,kappa)
	to_return[!trial_is_vm] = runif(sum(!trial_is_vm),0,2*pi)
	to_return = to_return-pi #center the data on 0
	return(to_return)
}

#define a modification of the est.kappa function from CircStats where mu is assumed to be zero.
est_kappa = function( x , bias ){
	kappa <- A1inv(mean(cos(x)))
	if (bias == TRUE) {
		kappa.ml <- kappa
		n <- length(x)
		if (kappa.ml < 2) 
			kappa <- max(kappa.ml - 2 * (n * kappa.ml)^-1, 0)
		if (kappa.ml >= 2) 
			kappa <- ((n - 1)^3 * kappa.ml)/(n^3 + n)
	}
	return(kappa)
}

#define a function to fit the uvm model
fit_uvm_random = function( x , iterations , fix_seed=TRUE , show_progress=TRUE , show_time=TRUE ){
	rho_prime = 2:length(x) #rho prime is simply rho*length(x)
	best_seed = NA
	best_rho_prime = NA
	best_kappa = NA
	best_SLL = -Inf
	d_un = dunif(0,-pi,pi)
	start=proc.time()[3]
	if(show_progress){
		pb = txtProgressBar(max=length(rho_prime)-2,style=3)
		progress=1
	}
	for(this_rho_prime in rho_prime){
		n_un = length(x)-this_rho_prime
		kappa_bias_correct = ifelse(this_rho_prime<16,T,F)
		for(j in 1:iterations){
			if(fix_seed){
				set.seed(j)
			}
			vm = sample(
				x = 1:length(x)
				, size = this_rho_prime
			)
			vm_kappa = est_kappa(x[vm],kappa_bias_correct)
			vm_lik = dvm(x[vm],0,vm_kappa)
			un_lik = rep(d_un,n_un)
			this_SLL = sum(log(c(vm_lik,un_lik)))
			if(this_SLL>best_SLL){
				best_seed = j
				best_rho_prime = this_rho_prime
				best_kappa = vm_kappa
				best_SLL = this_SLL
			}
		}
		if(show_progress){
			setTxtProgressBar(pb,progress)
			progress=progress+1
		}
	}
	if(show_progress){
		close(pb)
	}
	if(show_time){
		print(proc.time()[3]-start)
	}
	to_return = list( best_rho=best_rho_prime/length(x) , best_kappa , best_SLL )
	names(to_return) = c( 'best_rho' , 'best_kappa' , 'best_SLL' )
	return(to_return)
}

#generate some data
x = get_mix_data( n=100 , rho=.9 , kappa=10 , seed=1 )

#try to fit the data
this_fit = fit_uvm_random( x=x , iterations=1e3 )

print(this_fit)

#let the seed vary to get a different fit
this_fit = fit_uvm_random( x=x , iterations=1e3 , fix_seed=FALSE )

print(this_fit)