*! 1.1.1 NJC 4 Nov 2006 * 1.1.0 NJC 12 Oct 2006 * 1.0.0 NJC 9 Sept 2006 program nrunsi, sort rclass version 9 syntax anything(id="values in sequence") /// [, nreps(int 10000) saving(str) ] tokenize "`anything'" local nvals : word count `anything' quietly { preserve drop _all set obs `nvals' tempvar values runs freq prob gen `values' = "" forval i = 1/`nvals' { replace `values' = "``i''" in `i' } count if `values' != `values'[_n-1] local obsruns = r(N) gen long `runs' = _n mata: work(`nreps', "`c(seed)'", "`values'", "`freq'") gen double `prob' = `freq' / `nreps' } di _n as txt "observed number of values = " as res `nvals' di as txt "observed number of runs = " as res `obsruns' di _n as txt "simulated number of runs in " as res `nreps' /// as txt " shuffles" _c label var `runs' "# of runs" label var `freq' "frequency" label var `prob' " probability" su `runs' if `freq', meanonly tabdisp `runs' if inrange(`runs', `r(min)', `r(max)'), /// c(`freq' `prob') su `prob' if `runs' < `obsruns', meanonly return scalar pless = r(sum) di _n as txt " P(# of runs < " as res `obsruns' /// as txt ") = " %7.5f as res r(sum) su `prob' if `runs' == `obsruns', meanonly return scalar psame = r(sum) di as txt " P(# of runs = " as res `obsruns' /// as txt ") = " %7.5f as res r(sum) su `prob' if `runs' > `obsruns', meanonly return scalar pmore = r(sum) di as txt " P(# of runs > " as res `obsruns' /// as txt ") = " %7.5f as res r(sum) if `"`saving'"' != "" { quietly { clonevar values = `values' clonevar runs = `runs' clonevar freq = `freq' clonevar prob = `prob' label var prob "probability" sort runs keep values runs freq prob } di save `saving' } end mata: void work(real scalar nreps, string scalar seed, string scalar valuesname, string scalar freqname) { string colvector values real colvector freq real scalar n, runs values = st_sdata(., valuesname) n = rows(values) freq = J(n, 1, 0) uniformseed(seed) for(i = 1; i <= nreps; i++) { _jumble(values) runs = 1 + sum(values[|1\n-1|] :!= values[|2\n|]) freq[runs] = freq[runs] + 1 } st_addvar("long", freqname) st_store(., freqname, freq) } end