using Turing, StatsFuns
using Distributions, StatsPlots

using FillArrays

begin #creating synthetic data
	βs = [1 2 3]'
	x = Matrix([1:10 1:2:20 1:3:30])
	X = (x .- mean(x,dims=1)) ./ std(x,dims=1)
	t = [x for x=1:0.5:5.5]
	s = rand([1,2],10)
	σ= [2 3]
	rand_draw1 = randn(10)

	y = x*βs .+ σ[s]
end

@model function JointDurationStateModel(
	DeviationFromExpectedDuration,
	ConclusionStatus,
	SnapshotState,
	CurrentDuration,
)
	# get dimensions
	n,k = size(SnapshotState)
	
	#hyperpriors priors

	#Heirarchal parameters
	#β ~ MvNormal(Fill(0,k),2)
	η ~ MvNormal(Fill(0,k),2)

	#Direct Priors
	#σ_DFED ~ filldist(Exponential(1),2) #TODO: check implication of this form
	
	#model
	#μ = SnapshotState * β
	p = StatsFuns.logistic.(SnapshotState * η)
	
	#estimate ConclusionStatus model
	ConclusionStatus .~ Bernoulli(p)
	#Estimate DFED model
	#=
	for i in eachindex(ConclusionStatus)
		DeviationFromExpectedDuration ~ Normal(
			μ[ConclusionStatus[i]],
			σ_DFED[ConclusionStatus[i]]
		)
	end
	=#
end

model = JointDurationStateModel(y,s,X,t)
prior = JointDurationStateModel(fill(missing,size(y)),fill(missing,size(s)),X,t)

chain = sample(model,NUTS(0.85),2000)



plot(chain)