st: Dirichlet preference probability

[Marie-Eve St-Onge <marie-eve.st-onge@outlook.com>]

Dear all,
I'm interested in estimating conditional probability parameters from a Dirichlet distribution of 4 brands. Additionally, my program will simulate 10,000 values via MCMC, given prior values provided.

Precisely, I want to simulate the probability that a company identified by column alpha1 beats alpha2 & gets> 30% of the preference. Additionally, I will calibrate the simulation with the size of N for each survey given by the "size" column. Also, the rows in the "prior" matrix is sorted by the order the preference was surveyed, so the values displayed in the first row reflect the more recent considerations than the last one.

I'm frustrated because I though would obtain this by using syntax like this: scalar `beta1' = `size'*invgammap(`alpha1'`alpha2'`alpha3'`alpha4', 1 -`alpha1' 1 -`alpha2' 1 -`alpha3' 1 -`alpha4',uniform())/100+1
Obviously, it didn't work. Any thoughts.

** this is where I'm stuck
capture program drop dirichlet
program define dirichlet, rclass
args size alpha1 alpha2 alpha3 alpha4
tempname beta1 beta2 beta3 beta4
scalar `beta1' = `size'*invgammap(`alpha1',uniform())
scalar `beta2' = `size'*invgammap(`alpha2',uniform())
scalar `beta3' = `size'*invgammap(`alpha3',uniform())
scalar `beta4' = `size'*invgammap(`alpha4',uniform())
return scalar beta1 = `beta1'
return scalar beta2 = `beta2'
return scalar beta3 = `beta3'
return scalar beta4 = `beta4'
end

** Example of prior data
clear
set obs 5
gen alpha1 = .29+(.34-.29)*runiform()
gen alpha2 = .28+(.33-.28)*runiform()
gen alpha3 = .16+(.20-.16)*runiform()
gen alpha4 = .10+(.13-.10)*runiform()
gen err = 1-(alpha1 + alpha2 + alpha3 + alpha4)
gen size = 1000+int((2500-1000)*runiform())

simulate beta1=r(beta1) beta2=r(beta2) beta3=r(beta3) beta4=r(beta4),reps(10000): dirichlet size alpha1 alpha2 alpha3 alpha4 		 	   		  
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/