/*This ado file generates random data according to the specified distribution.

Author --Jacob Orchard

v1.2
Update 7/29/2016*/

program rgb2, eclass
version 13.0
	if replay() {
		display "Replay not implemented"
	}
	else {
	args varname a b p q 
		confirm name `varname'
		confirm number `a' 
		confirm number `b'
		confirm number `p'
		confirm number `q'
		
		if `p' <=0{
			di as error "Parameter p must be positive"
		}
		if `q' <=0{
			di as error "Parameter q must be positive"
		}
		if 	`a' <= 0{
			di as error "Parameter a must be positive"
		}
		if `b' <= 0{
			di as error "Parameter b must be positive"
		}
		
		quietly mata: rgb2("`varname'",`a',`b',`p',`q')
	}

end

//Mata functions

version 13
mata:
	function rgb2(string myvar, scalar a, scalar b, scalar p, scalar q)
	{
	nobs = st_nobs()
	base = runiform(nobs,1)
	newvar = J(nobs,1,0)
	paravec = a,b,p,q
	
	if (nobs <1000){
	
		for (i=1; i <=nobs; i++){
			newvar[i] = estimategb2(paravec,base[i])
			}
		}
	else{
		cdflist = interpolategb2(paravec)
		
		for (i=1; i <=nobs; i++){
			newvar[i] = nearest_gb2(cdflist,base[i])
			}
	
	}
	
	(void) st_addvar("double", myvar)
	st_store(.,myvar,newvar[.,1])
	}
	end
	
mata:
	function estimategb2(paravec,unifval)
	{
	
		a = paravec[1]
		b = paravec[2]
		p     = paravec[3]
		q     = paravec[4]
		
		S = optimize_init()
		optimize_init_which(S,"min")
		optimize_init_technique(S,"dfp nr")
		optimize_init_evaluator(S,&closestgb2())
		
		//Expected Value (starting value)
		start = b*((exp(lngamma(p+(1/a)))*exp(lngamma(q-(1/a))))/( exp(lngamma(p))*exp(lngamma(q))))
		optimize_init_params(S,start)
		optimize_init_argument(S,1,paravec)
		optimize_init_argument(S,2,unifval)
		_optimize(S)
	    ans    = optimize_result_params(S)
		return(ans)
	
	}
	
end

mata:
	function interpolategb2(paravec)
	{
		cdflist = J(1000,1,0)
		
		for (i=1; i <=1000; i++){
		cdflist[i] = estimategb2(paravec,(i/1000))
		}
	return(cdflist)
    }
end

mata:
	function nearest_gb2(cdflist,unifval)
	{
		nx = 1000*unifval
		intx = floor(nx)
		remx = nx - intx
		next = intx+1
		
		if (intx ==0) {
			closegb2 = cdflist[1]*remx
		}
		else{
			closegb2 = cdflist[intx]*(1-remx) + cdflist[next]*remx
		}
		return(closegb2)
	}
end

mata: 
	function closestgb2(real scalar todo, real rowvector p, paravec,unifval, v,g,H)
	{
	
	v = abs(gb2_cdf(p,paravec)-unifval)
	
	}
	end
	
mata:
function gb2_cdf(matrix y, paravec)
	{
	ones = J(1,cols(y),1)
	a = paravec[1]
	b = paravec[2]
	p = paravec[3]
	q = paravec[4]
	zu = ((y:/b):^(a)):/(ones+(y:/b):^(a))
	F =  ibeta(p,q,zu)
	return(F)
	}
end