From 41188ac2bab6ae872dc93f7c81448008d83a73d3 Mon Sep 17 00:00:00 2001 From: Luis Baroja Date: Mon, 27 Mar 2023 17:41:19 -0700 Subject: [PATCH] Prior over sd instead tau --- GenerativeMatching/generative_matching.bug | 17 ++++---- GenerativeMatching/plot_simulations.R | 6 +++ GenerativeMatching/posterior_quantiles.pdf | Bin 0 -> 5862 bytes GenerativeMatching/test_matching.R | 46 +++++++++++++-------- 4 files changed, 43 insertions(+), 26 deletions(-) create mode 100644 GenerativeMatching/plot_simulations.R create mode 100644 GenerativeMatching/posterior_quantiles.pdf diff --git a/GenerativeMatching/generative_matching.bug b/GenerativeMatching/generative_matching.bug index 4b3d103..85f0e8f 100644 --- a/GenerativeMatching/generative_matching.bug +++ b/GenerativeMatching/generative_matching.bug @@ -1,12 +1,13 @@ model{ alpha ~ dnorm(mean_alpha_prior,1/(sd_alpha_prior)^2) beta ~ dnorm(mean_beta_prior,1/(sd_beta_prior)^2) - tau ~ dgamma(shape_tau_prior,rate_tau_prior)T(0.01,) - #tau ~ dlnorm(shape_tau_prior,rate_tau_prior) - for(i in 1:n_obs){ - lambda_Br[i] ~ dlnorm(alpha/2+beta*log(Wr[i])/2,tau) - lambda_Bl[i] ~ dlnorm(-alpha/2-beta*log(Wl[i])/2,tau) - Br[i] ~ dpois(lambda_Br[i]) - Bl[i] ~ dpois(lambda_Bl[i]) - } + tau ~ dgamma(shape_tau_prior,rate_tau_prior) + sd <- 1/pow(tau,1/2) + #tau ~ dlnorm(shape_tau_prior,rate_tau_prior) + for(i in 1:n_obs){ + lambda_Br[i] ~ dlnorm(alpha/2+beta*log(Wr[i])/2,tau) + lambda_Bl[i] ~ dlnorm(-alpha/2-beta*log(Wl[i])/2,tau) + Br[i] ~ dpois(lambda_Br[i]) + Bl[i] ~ dpois(lambda_Bl[i]) + } } diff --git a/GenerativeMatching/plot_simulations.R b/GenerativeMatching/plot_simulations.R new file mode 100644 index 0000000..e1a5d40 --- /dev/null +++ b/GenerativeMatching/plot_simulations.R @@ -0,0 +1,6 @@ +pdf('posterior_quantiles.pdf',width=4,height=12) +layout(1:3) +hist(quantiles_alpha,breaks=60) +hist(quantiles_beta,breaks=60) +hist(quantiles_tau,breaks=60) +dev.off() diff --git a/GenerativeMatching/posterior_quantiles.pdf b/GenerativeMatching/posterior_quantiles.pdf new file mode 100644 index 0000000000000000000000000000000000000000..c4da6895f1e73226190683b819573d789357032b GIT binary patch literal 5862 zcmb7Ic|4T++ZK_%Y)O%atc5XS-^y+f*>^IH!7ygVEcUVQBxEgXgh(o6%T^>=vx^qX zLDsTlYx&J|PUoE8Iq&;^KJWZ7pLwqDy?md0{|%4G75*yhRuZygbh33 z7{P!vfVt+usH(~c(!wEN1SAHn4I>}`F>PgtG(<*PMj8T^lL5;qh%%?9SO$+SBL~Gkw8u}!C_p8a0HHWDk}w%f{?MFy)YC=oblpCMg-c0atp|x zw*f#F{#XP6(t;6SD2y8bGJ&}v@PG_Czzp~+LE&EslsHN>4uK~85f0MApb6whJfQe9 z9b|^UV~997xt{Dl>QMrIjpH%`?M841zzRyz07w^!B9Pwzq>CaCM;ighxF9I^#uLah z^J4VNsl8>0;Zm87mwKTV@{RA;8d_%`6ZtY@5Tk*_+Q54c(q6wP%Q8QlF@Rlq02--U#_evScIOu|`SZ#Z z6q;O7oIjrW;~N4ur=N$qVbUGx>@S%}ImMJ4NH+!$znGG?=9+F#T?ui_53$I9S6(r& zm|3@R(P8>)@b{Vg)`3<_w+-L-`I1*i>(~j%jmFrqPI-^Wdj~a%0d49PhrUrP+(oZM zmc!>%TU+W6g1q+^S}#&-#|eH~9px7Mv^+WRX`FYowQ0Onyy*M1?GE+#=Jzw3A#bAY zYgu&D@;yE%bGWZuA{An?JUSMxZBO)ZHz}9$mjhCGY&L2tDs(a%^jVedT{R&x-twMj zCM@(9M21(#N2PCSK^BrKS**CT8oEht8bo)yw1Vbb(%TvFbi$WvyyAqCrPthQ4lB`T zY}v+j5%&m}k!S(`V6fR8xi8n1X#7AL27n?@*+u(l+SZ zqn1vI?sq{Eu+Dg&Y~!+sf~-3^;Jz3SWtF^6MCYZ+NMpS!-oUHX%n!e1m zYrHfdX?1wsv&eWWwE)U9Rcq6|A71o?)Gt>97wIdP^3CUz@XR3}nhm{UAr1-0!tKTM zU6Y7*!E4JC(W%WdW-qTn?{77}T=BT=F^}D?Fp4eP(eDCcXb+Md=G(PKUc{ZbrZsBZ zFL&<|^{uqH7kl|U z6n8Z80;0@Yl^R|ws}9{;&pmw7bFZt`5~NZk*;abSLvbH>QK{a*j1xa;TPMr=EMxa9 z6soZJVVm3gH`arZ`T6v^O^4J~*380mj#nD@c1`mFuwDrZY4l#Vo6dTKX6nV*dCkwx z&pNR@G+D*n=ndfVS4@DL8-40j?E*|GkY4cmISy#3 zY@iEuq1P!jmF`6KXdUjgX!U3f?zIf{=+0#Fx7AIlDao+==4+4NvcLT9vHF}nHpXz( zf;v0shH!HCz72JD&L2Nu{eeekBNmUcpHB;oPO1k)`8LB9ednW)>kXq74BW-_?2@rcTJ?Z&o}Wg%%IEOa|*dh&1j<=9?-G+O858*R5$DT*Z`_v?MWK{SbYO zCi=~*_=uiav~#`Mk03u^>#3A*i^=D-FwX3g60Ys=A!>A+dzGon{lOifeHL+ooPKZ* zgX4Rh*1M4~^*iq&&$58HU_<&TuDzIrh`fj`s1f#T?PXhxR@PQ9*Rnz zI=G9}i))pGR<7U3q)X@+sq{%7&DBOz>8Fk#fzH!4Ji&NIABpd6Yk#D_LNA1)aTu0% zHdK8W7qc>_y&6g|v}v79$s0THG0izDA+Me>XU`J2jmM+ZhV>uL6<8AYbUXMC9mWAe zy1Ob5inv>!HIPn=XZ#5Ho=uleh4VRs8y=Q3* zNySPWy!wv>21?&?}K`oE)iRZj!m+az{`=FD|Ym>%hm+=(3Ph(>kB; zgyxi}@#T&YtaZujjB{-ITOt5^moAW+8g?LKdR`yNE_jU zgd-?LLB?by0{TNwK;)F<|3^JRxznF|A~(a=!hG^%%s{h3vq|cS__RS!cQa}zOEWY@ z>O7d?gh5*}m7u(~1ynyow3&w+?J<8E$b1@;fm4_``br}~kf)?h0LrZUM&T-@B7bl~CcR#daj5w5$^ zo%KZ4s9)NkU(hxjVs#V0DPbS?RULZz`hA`a%!*Kcwp%nS?|EhcpKO0z0lKQ{}W$mX%-oxxGL4t%xc^!MU*s6Pli)u$U2WU|K9U2*1_2DyF9$43Fe{ii4Im4)R$;w5|pXn6F+K)VK%Xm`BqpT+JoFXYGU@gMt6JbRMut58cUC{0JZr{pqlIiuG8B}Xi6wq~^Y zUjD7UuE#k*-iUbZ+*3+mIbJE)M=L0&n(e!v>koY##T>t`~wZOhFI)T zZ)fm_aL(i-+5k6aXl=o%wD3x2(K0$>0h?@uJoXeJ4E{k@@t9OtkMq&#BWBX9d7tQ@)cR&35vcu4am#5CfY)MKXIaHB$Xe^4VhkhInX7=}UV~IvpLb z$>lVrSNvR2vv)EItqSFf%1w`%SJWtu*o!Ce`-nsRu=QCnY zH|TgU{iURQ0Wpq5h6j;Q&GRtH!qh=sMZl<-O#wo4pY0k7DWV$kI_Rl5X(tjB?~1j|7u*q zChc!}oD!X}vsyEz@nVN!ZennY$nq-)^M(=Pi0-oOvLcI?ObP9v8S-uo#yz4&oGU&) z*PzBA&*1tiZ!U4Z!ZWXFmOgGn4qykw0~Je2%PW>G%Oa!lqZi7nMi@DQ30oGSVxk)>ZzfP{8Kchyb7eTVZr_^-9L zw~T*}h$xD9S=a{~dff(o{b=knUBU}SxR-0LucnH>Ql*db!uHr?yX|Vw#6CB7iQz57 zgq|Y*7h8KAM!uW z9I-nx9l8*DPra)xJ&dGLTtKTas&PELFI=ZxO#lmZe)%Q2X^Ul%Y0)<-d@yuyt#J2Q z;`59s5`PlIai%JUkW(-Y(-W^-LmQpD!!u~B1e*nW1UIy^we{}=6mb=C+__zJwac$7 ztqTgTyk`BvDSkMXql4H5pDUWXJzL6sATiy5Rc*gms5o;VIH%*EASE$olhd-d$KENN zBdo{%Iu0Bs8MmQSm_ygUJrsJk*08i5wJx>Jk#`Jp43oO?3bc4ds6ixMRL#K0=v|&_ z_H3S6%BDrzNeiaD(u2~myDi2{MpOA;#MMNa)o3zPvb-&zY$XhbkXo57NKAEab$oSk zPfq5VgS;>!yBtDmGU0Jkr|jgwNN9)lV<*V1xeW80TOta|yU;vbo;1d}VRF-`M|6|4 z)#?;x*Y5b*zQF$B6KBHR#a9QGgXMQC0<2f6r0mpMUSRsdY9niro*`m2QWBt`1}Rs^ zD=Aj5Z{2yc6TD)dz!k64sUTmEDw~RZBQ~K}k6bS)H{v(+6XrIoG`wK&9rTp4$y|ME0<$))OL%Db z*wBv0VbT#{Qd_%BA{y*KuUi`hR=Td0ycuuuZ>n1TJTYulYH4k^i>`9s=nH5)=%*`+ zjwfZ2ln&3z>neN+={sm_&@Z}I^YMA`Oz_9BSQ=7aQeU6sElETrwxXE0>HMs{=e=LZ z@Y@%t3QTET=2_|(p_sEFxA)d|o@>Me-guZ)Lvq`5^g32sRLrsL@J?woeMopMY31z2 z(|$E|U3W|80(&-g`G?i4@D9KEO)@--yTg%hypW%BF`c;l@#Fk3yxA~dBPUJ`UkQ4O z`Ggvr+*%A8$o;sxdddIY&eumw{oF}TE{>~>^Voa%w%~%fV~708bbNFwF}yL2AQi`T zo;-!1+4#fsVe96qyrz@bwV{_^Q)X67$1(4q^#aY^P5Vn$2bGI;c=W)uEuNcHH=j4u z9DjIg_e`+tw?~f4fZD}nHRVsv?GY;Q=|`fM5I?qqpDqdP!#X4TwYSbq>2-;aL_wFF z_WZu;kKN2J?Ywm1hr?m-{G@NI-&cq^HjMy zWj-|x!GDlk@%qt}<;x+{u<>Z)!lkU>rTtUepVvzhujB{Q>_%=K{!4W${V78k|1i&V z&}5$qi6*;e6xsUMlQj~pfyN_$0gNDRq^m1}?3ht3F?&E(4&;SI6Y+qOJjk1fAs|o) zR{|qgRu<%fA)9=7JQ4u8;b1;wjSGVlDN+}X#KDPPt|){b03sk!E@a3H2FGF00O%SH zK?x*h!r*YS3C9R>L6TeHk$3aum8r#-Yt1^`UeJ)AunM6z+~jfJJUZgrGM2U zo55s{_P^t#Axi(kV{|cS}!x4T?OO9NJQCL{#iZ0`S02SF6aR2}S literal 0 HcmV?d00001 diff --git a/GenerativeMatching/test_matching.R b/GenerativeMatching/test_matching.R index a5987b9..509deb7 100644 --- a/GenerativeMatching/test_matching.R +++ b/GenerativeMatching/test_matching.R @@ -9,6 +9,7 @@ posterior_quantile <- function(mean_alpha_prior,sd_alpha_prior, alpha0 <- rnorm(n=1,mean=mean_alpha_prior,sd=sd_alpha_prior) beta0 <- rnorm(n=1,mean=mean_beta_prior,sd=sd_beta_prior) tau0 <- rgamma(n=1,shape=shape_tau_prior,rate=rate_tau_prior) + sd0 <- 1/sqrt(tau0) #tau0 <- rlnorm(n=1,meanlog=shape_tau_prior,sdlog=1/sqrt(rate_tau_prior)) # Sample data 'Br, Bl' from the likelihood given prior sample 'th0' Br <- NA @@ -20,8 +21,8 @@ posterior_quantile <- function(mean_alpha_prior,sd_alpha_prior, for(i in 1:n_sessions){ lambda_Br0[i] <<- rlnorm(n=1,alpha0/2+beta0*log(Wr[i])/2,1/sqrt(tau0)) lambda_Bl0[i] <<- rlnorm(n=1,-alpha0/2-beta0*log(Wl[i])/2,1/sqrt(tau0)) - Br[i] <- rpois(n=1,lambda=lambda_Br0[i]) - Bl[i] <- rpois(n=1,lambda=lambda_Bl0[i]) + Br[i] <- rpois(n=1,lambda=lambda_Br0[i])+1 + Bl[i] <- rpois(n=1,lambda=lambda_Bl0[i])+1 } # Compute a (to-be-tested) posterior sample [th1,th2,...,thL] # given the sample from the likelihood @@ -32,46 +33,55 @@ posterior_quantile <- function(mean_alpha_prior,sd_alpha_prior, 'mean_beta_prior','sd_beta_prior', 'shape_tau_prior','rate_tau_prior', 'n_obs') - unobserved <- c('alpha','beta','tau','lambda_Br','lambda_Bl') + unobserved <- c('alpha','beta','tau','sd','lambda_Br','lambda_Bl') write('model{ alpha ~ dnorm(mean_alpha_prior,1/(sd_alpha_prior)^2) beta ~ dnorm(mean_beta_prior,1/(sd_beta_prior)^2) - tau ~ dgamma(shape_tau_prior,rate_tau_prior)T(0.01,) - #tau ~ dlnorm(shape_tau_prior,rate_tau_prior) - for(i in 1:n_obs){ - lambda_Br[i] ~ dlnorm(alpha/2+beta*log(Wr[i])/2,tau) - lambda_Bl[i] ~ dlnorm(-alpha/2-beta*log(Wl[i])/2,tau) - Br[i] ~ dpois(lambda_Br[i]) - Bl[i] ~ dpois(lambda_Bl[i]) - } + tau ~ dgamma(shape_tau_prior,rate_tau_prior) + sd <- 1/pow(tau,1/2) + #tau ~ dlnorm(shape_tau_prior,rate_tau_prior) + for(i in 1:n_obs){ + lambda_Br[i] ~ dlnorm(alpha/2+beta*log(Wr[i])/2,tau) + lambda_Bl[i] ~ dlnorm(-alpha/2-beta*log(Wl[i])/2,tau) + Br[i] ~ dpois(lambda_Br[i]) + Bl[i] ~ dpois(lambda_Bl[i]) + } }','generative_matching.bug') bayes <- jags(data=observed, parameters.to.save=unobserved, - model.file='generative_matching.bug') + model.file='generative_matching.bug', + inits=list(list(alpha=0,beta=0,tau=1), + list(alpha=0,beta=0,tau=1), + list(alpha=0,beta=0,tau=1)), + n.chains=3) nds<- bayes$BUGSoutput$sims.list alpha <- nds$alpha # Posterior sample beta <- nds$beta # Posterior sample tau <- nds$tau # Posterior sample + sd <- nds$sd # Posterior sample L <- dim(alpha)[1] # Size of the posterior sample # Compute quantile (1/L)*sum(I_{th0>th_m}) q0_alpha <- sum(alpha0>alpha)/L q0_beta <- sum(beta0>beta)/L q0_tau <- sum(tau0>tau)/L - return(list(q0_alpha=q0_alpha,q0_beta=q0_beta,q0_tau=q0_tau)) + q0_sd <- sum(sd0>sd)/L + return(list(q0_alpha=q0_alpha,q0_beta=q0_beta,q0_tau=q0_tau,q0_sd=q0_sd)) } quantiles_alpha <- NA quantiles_beta <- NA quantiles_tau <- NA - for(i in 1:1000){ - iteration <- posterior_quantile(mean_alpha_prior=0,sd_alpha_prior=1/sqrt(0.1), - mean_beta_prior=0,sd_beta_prior=1/sqrt(0.1), - shape_tau_prior=2,rate_tau_prior=0.5, + quantiles_sd <- NA + for(i in 1:100){ + iteration <- posterior_quantile(mean_alpha_prior=0,sd_alpha_prior=1/sqrt(1), + mean_beta_prior=0,sd_beta_prior=1/sqrt(1), + shape_tau_prior=1,rate_tau_prior=1, #shape_tau_prior=5,rate_tau_prior=1, - n_sessions=10) + n_sessions=5) quantiles_alpha[i] <- iteration$q0_alpha quantiles_beta[i] <- iteration$q0_beta quantiles_tau[i] <- iteration$q0_tau + quantiles_sd[i] <- iteration$q0_sd print(i) }