*! 1.7 16 Apr 2017 Austin Nichols add new default for age indicator option
*  1.6 21 Mar 2017 Austin Nichols removed lnmvnormal and replaced with (-.5*(x-m)*invsym(v)*(x-m)')-ln(pi()*det(v))
*  1.5 20 Mar 2017 Austin Nichols removed grouplist and gvar options, added agevar and modevar options
*  1.4 16 Mar 2017 Austin Nichols renumbered groups
*  1.3 13 Mar 2017 Austin Nichols added check for correct average correlations across items
*  1.2 09 Mar 2017 Austin Nichols removed obsolete groups
*  1.2 08 Mar 2017 Austin Nichols added round() and missok options
*  1.1 03 Mar 2017 Austin Nichols changed variable name defaults
*  1.0 17 Oct 2016 Austin Nichols <austinnichols@gmail.com>
prog pfwb, rclass
 version 10.2
 syntax name [if] [in] [, modevar(varname numeric) agevar(varname numeric) gvar(varname numeric) qlist(varlist numeric min=10 max=10) replace INTPoints(int 49) zlimit(real 6) ghq round(real 1) missok skipcorr]
 marksample touse
 if "`replace'"=="" qui g double `namelist'=.
 else {
  cap g double `namelist'=.
  }
 *check rounding value
 if `round'<0 {
  di as err "Rounding value for scores must be positive"
  error 198
  }
 *check mode var
 if "`modevar'"=="" {
  loc modevar "self"
  }
 cap assert inlist(`modevar',0,1) if `touse'
 if _rc {
  di as err "mode variable `modevar' not defined correctly"
  err 459
  }
 *check age var
 if "`agevar'"=="" {
  loc agevar "age18_61"
  }
 cap assert inlist(`agevar',0,1) if `touse'
 if _rc {
  di as err "age variable `agevar' not defined correctly"
  err 459
  }
 *check group var
 tempvar gvar
 qui g byte `gvar'=1+2*`agevar'+`modevar' if `touse'
 cap conf numeric var `gvar'
 if _rc!=0 {
  di as err "Error: group var `gvar' not found"
  error 198
  }
 cap assert inlist(`gvar',1,2,3,4) if `touse'
 if _rc!=0 {
  di as err "Error: group var not in {1,2,3,4}"
  error 459
  }
 *check ghq option
 loc GHQ=cond("`ghq'"!="",1,0)
 if "`ghq'"!="" {
  if "`zlimit'"!="6" {
   di as err "zlimit option is ignored when option ghq is specified"
   }
  }
 *check questions
if "`qlist'"=="" loc qlist "fwb1_exp fwb2_getby fwb3_secure fwb4_concern fwb5_never fwb6_enjoy fwb7_behind fwb8_control fwb9_strain fwb10_left"
if wordcount("`qlist'")!=10 {
 di as err "Must specify 10 variables holding all 10 items in FWB instrument (even if some are everywhere missing)"
 error 198
 }
if "`skipcorr'"=="" {
 tokenize `qlist'
 loc avgcorr=0
 loc ctcorr=0
 loc c 0
 foreach pos in 1 3 6 10 {
  foreach neg in 2 4 5 7 8 9 {
   cap corr ``pos'' ``neg''
   if _rc==0 loc c=r(rho)
   if `c'<. {
    loc avgcorr=`avgcorr'+`c'
    loc ctcorr=`ctcorr'+1
    }
   }
  }
 if `avgcorr'<0 {
  di as err "Reverse-coded items seem to have a negative correlation with other items; check coding of items."
  di as err "If you are sure of your coding of items, use the" as txt " skipcorr " as err "option." _n
  err 459
  }
 }
if "`grouplist'"=="" {
  loc grouplist "1 2 3 4"
  }
tempvar allmiss
g byte `allmiss'=1
 foreach v in `qlist' {
  cap conf numeric var `v'
  if _rc!=0 {
   di as err "Error: question var `v' not found"
   di as err "Every question needs to be present, though all values may be missing."
   error 198
   }
  cap assert inlist(`v',0,1,2,3,4)|mi(`v')
  if _rc!=0 {
   di as err "Error: question var `v' responses not in {0,1,2,3,4}"
   di as err "Every question needs to be coded 0 to 4, or missing."
   error 198
   }
  qui replace `allmiss'=0 if !mi(`v')
  }
di as res "Computing scores... this may take some time."
 * check gsem option
 if "`gsem'"!="" {
 mat g1=(2.1987183,2.5977786,2.0262112,-.2961961,1.6585993,1.9563691,1.9234236,-.9979853,2.5144634,-.5726793,2.1467741,2.3913435,2.1146504,.5220543,1.4970444,.159058,2.8716455,.3471574,1.5260218,1.4633664)
 mat g2=(2.578371,2.5977786,2.0262112,-.2961961,2.2843891,1.9563691,1.9234236,-.9979853,2.5144634,-.5726793,2.1467741,2.3913435,2.1146504,.5220543,1.9950379,.159058,2.8716455,.3471574,1.9723144,1.4633664)
 mat g3=(2.1987183,2.5977786,2.0262112,-.2961961,1.6585993,1.9563691,1.9234236,-.9979853,2.5144634,-.5726793,2.1467741,2.3913435,2.1146504,.5220543,1.5421585,.159058,2.8716455,.3471574,1.9723144,1.4633664)
 mat g4=(2.1987183,2.5977786,2.0262112,-.2961961,1.6585993,1.9563691,1.9234236,-.9979853,2.5144634,-.5726793,2.1467741,2.3913435,2.1146504,.5220543,1.9950379,.159058,2.8716455,.3471574,1.9723144,1.4633664)
 mat h1=(-2.2047837,-1.0882027,.3585419,2.5735247,-1.3208536,-.2740339,.8621188,1.9373029,-2.2144952,-1.059253,.4227569,2.1435156,-2.2209814,-.4245686,1.4446934,3.3603167,-3.0937724,-1.3091939,.4448066,2.4752427,-3.8399637,-2.1489973,-.1686061,1.7348393,-2.0926086,-1.0962606,-.3216849,.3565704,-2.1420109,-1.064027,.3304075,1.6338181,-2.0624651,-1.1108175,.2679782,1.535325,-1.7448388,-.7562154,.3777863,1.4842417,1.0871438,.7192873,3.2854107,0,0,0)
 mat h2=(-2.4179007,-.4122027,2.0406173,4.7737384,-2.2328062,-.5681815,.9966695,2.622802,-2.6647402,-.8492241,1.6940806,4.5478367,-2.2209814,-.4245686,1.4446934,3.3603167,-3.0937724,-1.3091939,.4448066,2.4752427,-3.5076333,-1.2281664,1.5150107,4.3227005,-3.4129491,-1.8072053,-.1601308,1.7784225,-2.8620269,-1.0488892,1.0262574,3.1092984,-3.6906161,-1.5936348,.6377263,3.1938995,-3.2028665,-1.0689037,.9702596,2.9018261,1.3285903,.4652468,.8905083,0,0,0)
 mat h3=(-1.2052772,.2585893,1.9945245,3.7360925,-2.2328062,-.5681815,.9966695,2.622802,-3.1553966,-1.2531268,.9456383,3.2224684,-2.2209814,-.4245686,1.4446934,3.3603167,-3.0937724,-1.3091939,.4448066,2.4752427,-2.9830786,-1.289892,1.2312826,3.4551957,-2.0225894,-1.1311798,.125659,1.4371861,-1.2378236,-.4001816,1.0566471,2.295852,-2.465411,-1.5280219,.0812848,1.8023567,-1.5709306,-.6284243000000001,.8178498,2.2357033,1.3252162,.512746,2.5283897,0,0,0)
 mat h4=(-2.4179007,-.4122027,2.0406173,4.7737384,-2.2328062,-.5681815,.9966695,2.622802,-3.1553966,-1.2531268,.9456383,3.2224684,-2.2209814,-.4245686,1.4446934,3.3603167,-3.0937724,-1.3091939,.4448066,2.4752427,-3.5076333,-1.2281664,1.5150107,4.3227005,-3.4129491,-1.8072053,-.1601308,1.7784225,-2.8620269,-1.0488892,1.0262574,3.1092984,-3.6906161,-1.5936348,.6377263,3.1938995,-3.2028665,-1.0689037,.9702596,2.9018261,1,1,1,0,0,0)
 mat f1=g1,h1
 mat f2=g2,h2
 mat f3=g3,h3
 mat f4=g4,h4
 set more 1
 *loop over groups
 tokenize `qlist'
 foreach i of numlist `grouplist' {
  mat coln f`i'= `"`1':F"' `"`1':P"' `"`2':F"' `"`2':N"' `"`3':F"' `"`3':P"' `"`4':F"' `"`4':N"' `"`5':F"' `"`5':N"' `"`6':F"' `"`6':P"' `"`7':F"' `"`7':N"' `"`8':F"' `"`8':N"' `"`9':F"' `"`9':N"' `"`10':F"' `"`10':P"' `"`1'_cut1:_cons"' `"`1'_cut2:_cons"' `"`1'_cut3:_cons"' `"`1'_cut4:_cons"' `"`2'_cut1:_cons"' `"`2'_cut2:_cons"' `"`2'_cut3:_cons"' `"`2'_cut4:_cons"' `"`3'_cut1:_cons"' `"`3'_cut2:_cons"' `"`3'_cut3:_cons"' `"`3'_cut4:_cons"' `"`4'_cut1:_cons"' `"`4'_cut2:_cons"' `"`4'_cut3:_cons"' `"`4'_cut4:_cons"' `"`5'_cut1:_cons"' `"`5'_cut2:_cons"' `"`5'_cut3:_cons"' `"`5'_cut4:_cons"' `"`6'_cut1:_cons"' `"`6'_cut2:_cons"' `"`6'_cut3:_cons"' `"`6'_cut4:_cons"' `"`7'_cut1:_cons"' `"`7'_cut2:_cons"' `"`7'_cut3:_cons"' `"`7'_cut4:_cons"' `"`8'_cut1:_cons"' `"`8'_cut2:_cons"' `"`8'_cut3:_cons"' `"`8'_cut4:_cons"' `"`9'_cut1:_cons"' `"`9'_cut2:_cons"' `"`9'_cut3:_cons"' `"`9'_cut4:_cons"' `"`10'_cut1:_cons"' `"`10'_cut2:_cons"' `"`10'_cut3:_cons"' `"`10'_cut4:_cons"' `"var(F):_cons"' `"var(P):_cons"' `"var(N):_cons"' `"cov(P,F):_cons"' `"cov(N,F):_cons"' `"cov(N,P):_cons"'
  mat rown f`i'="y1"
  gsem (F -> `1' `2' `3' `4' `5' `6' `7' `8' `9' `10', ologit) (P -> `1' `3' `6' `10' , ologit) (N -> `2' `4' `5' `7' `8' `9', ologit) if grp==`i'*`touse', from(f`i') iter(0) intp(`intpoints')
  tempvar f`i'
  qui predict double `f`i'' if e(sample), latent(F)
  qui replace `namelist'=(`f`i''*15+50) if `f`i''<.
  }
 }
 else {
 * define parameter vectors
 mat g1m=(.6291655,-.0685166,.442397)
 mat g1v=(1.0871438,.7192873,3.2854107)
 mat g1q1a=(2.1987183,2.5977786,0)
 mat g1q1c=(-2.2047837,-1.0882027,.3585419,2.5735247)
 mat g1q2a=(2.0262112,0,-.2961961)
 mat g1q2c=(-1.3208536,-.2740339,.8621188,1.9373029)
 mat g1q3a=(1.6585993,1.9563691,0)
 mat g1q3c=(-2.2144952,-1.059253,.4227569,2.1435156)
 mat g1q4a=(1.9234236,0,-.9979853)
 mat g1q4c=(-2.2209814,-.4245686,1.4446934,3.3603167)
 mat g1q5a=(2.5144634,0,-.5726793)
 mat g1q5c=(-3.0937724,-1.3091939,.4448066,2.4752427)
 mat g1q6a=(2.1467741,2.3913435,0)
 mat g1q6c=(-3.8399637,-2.1489973,-.1686061,1.7348393)
 mat g1q7a=(2.1146504,0,.5220543)
 mat g1q7c=(-2.0926086,-1.0962606,-.3216849,.3565704)
 mat g1q8a=(1.4970444,0,.159058)
 mat g1q8c=(-2.1420109,-1.064027,.3304075,1.6338181)
 mat g1q9a=(2.8716455,0,.3471574)
 mat g1q9c=(-2.0624651,-1.1108175,.2679782,1.535325)
 mat g1q10a=(1.5260218,1.4633664,0)
 mat g1q10c=(-1.7448388,-.7562154,.3777863,1.4842417)
 mat g2m=(.7185538,-.3572805,.2846914)
 mat g2v=(1.3285903,.4652468,.8905083)
 mat g2q1a=(2.578371,2.5977786,0)
 mat g2q1c=(-2.4179007,-.4122027,2.0406173,4.7737384)
 mat g2q2a=(2.0262112,0,-.2961961)
 mat g2q2c=(-2.2328062,-.5681815,.9966695,2.622802)
 mat g2q3a=(2.2843891,1.9563691,0)
 mat g2q3c=(-2.6647402,-.8492241,1.6940806,4.5478367)
 mat g2q4a=(1.9234236,0,-.9979853)
 mat g2q4c=(-2.2209814,-.4245686,1.4446934,3.3603167)
 mat g2q5a=(2.5144634,0,-.5726793)
 mat g2q5c=(-3.0937724,-1.3091939,.4448066,2.4752427)
 mat g2q6a=(2.1467741,2.3913435,0)
 mat g2q6c=(-3.5076333,-1.2281664,1.5150107,4.3227005)
 mat g2q7a=(2.1146504,0,.5220543)
 mat g2q7c=(-3.4129491,-1.8072053,-.1601308,1.7784225)
 mat g2q8a=(1.9950379,0,.159058)
 mat g2q8c=(-2.8620269,-1.0488892,1.0262574,3.1092984)
 mat g2q9a=(2.8716455,0,.3471574)
 mat g2q9c=(-3.6906161,-1.5936348,.6377263,3.1938995)
 mat g2q10a=(1.9723144,1.4633664,0)
 mat g2q10c=(-3.2028665,-1.0689037,.9702596,2.9018261)
 mat g3m=(.5131654,.0516339,.5999753)
 mat g3v=(1.3252162,.512746,2.5283897)
 mat g3q1a=(2.1987183,2.5977786,0)
 mat g3q1c=(-1.2052772,.2585893,1.9945245,3.7360925)
 mat g3q2a=(2.0262112,0,-.2961961)
 mat g3q2c=(-2.2328062,-.5681815,.9966695,2.622802)
 mat g3q3a=(1.6585993,1.9563691,0)
 mat g3q3c=(-3.1553966,-1.2531268,.9456383,3.2224684)
 mat g3q4a=(1.9234236,0,-.9979853)
 mat g3q4c=(-2.2209814,-.4245686,1.4446934,3.3603167)
 mat g3q5a=(2.5144634,0,-.5726793)
 mat g3q5c=(-3.0937724,-1.3091939,.4448066,2.4752427)
 mat g3q6a=(2.1467741,2.3913435,0)
 mat g3q6c=(-2.9830786,-1.289892,1.2312826,3.4551957)
 mat g3q7a=(2.1146504,0,.5220543)
 mat g3q7c=(-2.0225894,-1.1311798,.125659,1.4371861)
 mat g3q8a=(1.5421585,0,.159058)
 mat g3q8c=(-1.2378236,-.4001816,1.0566471,2.295852)
 mat g3q9a=(2.8716455,0,.3471574)
 mat g3q9c=(-2.465411,-1.5280219,.0812848,1.8023567)
 mat g3q10a=(1.9723144,1.4633664,0)
 mat g3q10c=(-1.5709306,-.6284243000000001,.8178498,2.2357033)
 mat g4m=(0,0,0)
 mat g4v=(1,1,1)
 mat g4q1a=(2.1987183,2.5977786,0)
 mat g4q1c=(-2.4179007,-.4122027,2.0406173,4.7737384)
 mat g4q2a=(2.0262112,0,-.2961961)
 mat g4q2c=(-2.2328062,-.5681815,.9966695,2.622802)
 mat g4q3a=(1.6585993,1.9563691,0)
 mat g4q3c=(-3.1553966,-1.2531268,.9456383,3.2224684)
 mat g4q4a=(1.9234236,0,-.9979853)
 mat g4q4c=(-2.2209814,-.4245686,1.4446934,3.3603167)
 mat g4q5a=(2.5144634,0,-.5726793)
 mat g4q5c=(-3.0937724,-1.3091939,.4448066,2.4752427)
 mat g4q6a=(2.1467741,2.3913435,0)
 mat g4q6c=(-3.5076333,-1.2281664,1.5150107,4.3227005)
 mat g4q7a=(2.1146504,0,.5220543)
 mat g4q7c=(-3.4129491,-1.8072053,-.1601308,1.7784225)
 mat g4q8a=(1.9950379,0,.159058)
 mat g4q8c=(-2.8620269,-1.0488892,1.0262574,3.1092984)
 mat g4q9a=(2.8716455,0,.3471574)
 mat g4q9c=(-3.6906161,-1.5936348,.6377263,3.1938995)
 mat g4q10a=(1.9723144,1.4633664,0)
 mat g4q10c=(-3.2028665,-1.0689037,.9702596,2.9018261)
 loc qn 1
 loc qqlist
 foreach q in `qlist' {
  * turn each response into five columns in {0,1} to multiply probability
  tempvar q`qn'c1
  tempvar q`qn'c2
  tempvar q`qn'c3
  tempvar q`qn'c4
  tempvar q`qn'c5
  g byte `q`qn'c1'=(`q'==0)
  g byte `q`qn'c2'=(`q'==1)
  g byte `q`qn'c3'=(`q'==2)
  g byte `q`qn'c4'=(`q'==3)
  g byte `q`qn'c5'=(`q'==4)
  loc qqlist `qqlist' `q`qn'c1' `q`qn'c2' `q`qn'c3' `q`qn'c4' `q`qn'c5'
  loc qn=`qn'+1
  }
 *loop over groups
 foreach g of numlist `grouplist' {
  qui count if `touse' & `gvar'==`g'
  if r(N)>0 {
   timer clear 1
   timer on 1
   mat FWBmn=g`g'm
   mat FWBva=g`g'v
   forv q=1/10 {
    mat FWBq`q'a=g`g'q`q'a
    mat FWBq`q'c=g`g'q`q'c
    }
   tempname tmp`g'
   g byte `tmp`g''=min(`touse',`gvar'==`g')
   mata:ipr_fwb("`namelist'","`qqlist'","`tmp`g''",`intpoints',`zlimit',`GHQ')
   timer off 1
   qui timer list 1
   di as res "Group `g' complete in " r(t1) " seconds."
   }
  }
 ret scalar ghq=`GHQ'
 ret scalar intpoints=`intpoints'
 if "`ghq'"==""  ret scalar zlimit=`zlimit'
 }
 qui replace `namelist'=round(`namelist',`round')
 if "`missok'"=="" qui replace `namelist'=. if `allmiss'==1
 ret scalar round=`round'
 ret local name "`namelist'"
end
version 10.2
mata:
void ipr_fwb(string scalar f, string scalar x, string scalar tousename, real intpt, real zlim, real ghq) {
 st_view(y, ., tokens(f), tousename)
 q = st_data(., tokens(x), tousename)
 mean= st_matrix("FWBmn")
 v = st_matrix("FWBva")
 var= diag(v)
 q1a= st_matrix("FWBq1a")
 q1c= st_matrix("FWBq1c")
 q2a= st_matrix("FWBq2a")
 q2c= st_matrix("FWBq2c")
 q3a= st_matrix("FWBq3a")
 q3c= st_matrix("FWBq3c")
 q4a= st_matrix("FWBq4a")
 q4c= st_matrix("FWBq4c")
 q5a= st_matrix("FWBq5a")
 q5c= st_matrix("FWBq5c")
 q6a= st_matrix("FWBq6a")
 q6c= st_matrix("FWBq6c")
 q7a= st_matrix("FWBq7a")
 q7c= st_matrix("FWBq7c")
 q8a= st_matrix("FWBq8a")
 q8c= st_matrix("FWBq8c")
 q9a= st_matrix("FWBq9a")
 q9c= st_matrix("FWBq9c")
 q10a=st_matrix("FWBq10a")
 q10c=st_matrix("FWBq10c")
 P=J(rows(q),1,0)
 Pf=J(rows(q),1,0)
 Pft=J(rows(q),1,0)
 Z=_gauss_hermite_nodes(intpt)
 for (i=1; i<=intpt; i++) {
  if (ghq==1) t1=Z[1,i]
  else t1=(i-1)*2*zlim/(intpt-1)-(zlim)
  for (j=1; j<=intpt; j++) {
   if (ghq==1) t2=Z[1,j]
   else t2=(j-1)*2*zlim/(intpt-1)-(zlim)
   for (k=1; k<=intpt; k++) {
    if (ghq==1) t3=Z[1,k]
    else t3=(k-1)*2*zlim/(intpt-1)-(zlim)
 theta=(t1,t2,t3)
 a1theta=theta*q1a'
 a2theta=theta*q2a'
 a3theta=theta*q3a'
 a4theta=theta*q4a'
 a5theta=theta*q5a'
 a6theta=theta*q6a'
 a7theta=theta*q7a'
 a8theta=theta*q8a'
 a9theta=theta*q9a'
 a10theta=theta*q10a'
 lp1=(ln(1-invlogit(a1theta-q1c[1]))\ln(invlogit(a1theta-q1c[1])-invlogit(a1theta-q1c[2]))\ln(invlogit(a1theta-q1c[2])-invlogit(a1theta-q1c[3]))\ln(invlogit(a1theta-q1c[3])-invlogit(a1theta-q1c[4]))\ln(invlogit(a1theta-q1c[4])))
 lp2=(ln(1-invlogit(a2theta-q2c[1]))\ln(invlogit(a2theta-q2c[1])-invlogit(a2theta-q2c[2]))\ln(invlogit(a2theta-q2c[2])-invlogit(a2theta-q2c[3]))\ln(invlogit(a2theta-q2c[3])-invlogit(a2theta-q2c[4]))\ln(invlogit(a2theta-q2c[4])))
 lp3=(ln(1-invlogit(a3theta-q3c[1]))\ln(invlogit(a3theta-q3c[1])-invlogit(a3theta-q3c[2]))\ln(invlogit(a3theta-q3c[2])-invlogit(a3theta-q3c[3]))\ln(invlogit(a3theta-q3c[3])-invlogit(a3theta-q3c[4]))\ln(invlogit(a3theta-q3c[4])))
 lp4=(ln(1-invlogit(a4theta-q4c[1]))\ln(invlogit(a4theta-q4c[1])-invlogit(a4theta-q4c[2]))\ln(invlogit(a4theta-q4c[2])-invlogit(a4theta-q4c[3]))\ln(invlogit(a4theta-q4c[3])-invlogit(a4theta-q4c[4]))\ln(invlogit(a4theta-q4c[4])))
 lp5=(ln(1-invlogit(a5theta-q5c[1]))\ln(invlogit(a5theta-q5c[1])-invlogit(a5theta-q5c[2]))\ln(invlogit(a5theta-q5c[2])-invlogit(a5theta-q5c[3]))\ln(invlogit(a5theta-q5c[3])-invlogit(a5theta-q5c[4]))\ln(invlogit(a5theta-q5c[4])))
 lp6=(ln(1-invlogit(a6theta-q6c[1]))\ln(invlogit(a6theta-q6c[1])-invlogit(a6theta-q6c[2]))\ln(invlogit(a6theta-q6c[2])-invlogit(a6theta-q6c[3]))\ln(invlogit(a6theta-q6c[3])-invlogit(a6theta-q6c[4]))\ln(invlogit(a6theta-q6c[4])))
 lp7=(ln(1-invlogit(a7theta-q7c[1]))\ln(invlogit(a7theta-q7c[1])-invlogit(a7theta-q7c[2]))\ln(invlogit(a7theta-q7c[2])-invlogit(a7theta-q7c[3]))\ln(invlogit(a7theta-q7c[3])-invlogit(a7theta-q7c[4]))\ln(invlogit(a7theta-q7c[4])))
 lp8=(ln(1-invlogit(a8theta-q8c[1]))\ln(invlogit(a8theta-q8c[1])-invlogit(a8theta-q8c[2]))\ln(invlogit(a8theta-q8c[2])-invlogit(a8theta-q8c[3]))\ln(invlogit(a8theta-q8c[3])-invlogit(a8theta-q8c[4]))\ln(invlogit(a8theta-q8c[4])))
 lp9=(ln(1-invlogit(a9theta-q9c[1]))\ln(invlogit(a9theta-q9c[1])-invlogit(a9theta-q9c[2]))\ln(invlogit(a9theta-q9c[2])-invlogit(a9theta-q9c[3]))\ln(invlogit(a9theta-q9c[3])-invlogit(a9theta-q9c[4]))\ln(invlogit(a9theta-q9c[4])))
 lp10=(ln(1-invlogit(a10theta-q10c[1]))\ln(invlogit(a10theta-q10c[1])-invlogit(a10theta-q10c[2]))\ln(invlogit(a10theta-q10c[2])-invlogit(a10theta-q10c[3]))\ln(invlogit(a10theta-q10c[3])-invlogit(a10theta-q10c[4]))\ln(invlogit(a10theta-q10c[4])))
 lp=lp1\lp2\lp3\lp4\lp5\lp6\lp7\lp8\lp9\lp10
 P=exp(cross(q',lp):+((-.5*(theta-mean)*invsym(var)*(theta-mean)')-ln(pi()*det(var))))
 Pf=Pf+P
 Pft=Pft+P:*t1
   }
  }
 }
y[.,.]=Pft:/Pf:*15:+50
}
end