confidence.m

(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(* ::Input::Initialization:: *)
Off[General::compat];
Needs["Combinatorica`"](*For RandomKSubset*)


(* ::Input::Initialization:: *)
ClearAll[generateRandomSubsample];
generateRandomSubsample::usage="generateRandomSubsample[ssSize,groupIDs,dataArray] generates a subsample of a given size from a data array.";
generateRandomSubsample[ssSize_,groupIDs_,dataArray_]:=Module[{totalGroups,groups,qualifiedIndexes},
totalGroups=Union[Flatten[groupIDs,1]];
groups=RandomKSubset[totalGroups,ssSize];

WriteString["stdout"," ",groups];(*TCHRONIS for verbose results*)

(*Get the indexes of the dataArray rows that correspond to selected groups.*)
qualifiedIndexes=Position[groupIDs,_?(Intersection[groups,#]!={}&),{1},Heads->False];
(*We want to extract certain rows,so we need to transpose before and after extracting the rows we want,due to the matrix layout.*)Transpose[Extract[(*Transpose[*)dataArray(*]*),qualifiedIndexes]]
(*TCHRONIS for precomputed*)];


(* ::Input::Initialization:: *)
ClearAll[pointIdentifiedCR];
Options[pointIdentifiedCR]={progressUpdate->0,confidenceLevel->.95,asymptotics->nests,subsampleMonitor->Null,symmetric->False};
pointIdentifiedCR::usage="pointIdentifiedCR[ssSize,numSubsamples,pointEstimate,args,groupIDs,dataArray,method,permuteinvariant,options] generates a confidence region estimate using subsampling.\r
Parameters:
ssSize - The size of each subsample to be estimated.
numSubsamples -The number of subsamples to use in estimating the confidence region.
pointEstimate - The point estimate to build the confidence region around (typically the output of pairwiseMSE).
objFunc - The objective function used in pairwiseMSE.
args - A list of unique symbols used in pairwiseMSE.
groupIDs - A data map used to generate the subsamples.
dataArray - The dataArray parameter used in pairwiseMSE.
options - An optional parameter specifying options. Available options are:
	progressUpdate - How often to print progress (0 to disable).Default=0. 
	confidenceLevel - The confidence level of the region.Default=.95.
	asymptotics - Type of asymptotics to use (nests or coalitions).Default=nests.
	subsampleMonitor - An expression to evaluate for each subsample.Default=Null.
	symmetric - True or False.If True,the confidence region will be symmetric.Default=False.";

pointIdentifiedCR[ssSize_,numSubsamples_,pointEstimate_,args_,groupIDs_,dataArray_,method_,permuteinvariant_,options___?OptionQ]:=Module[{progress,confLevel,asymp,ssDataArray,estimates,estimate,alpha,cr,sym},{progress,confLevel,asymp,sym}={progressUpdate,confidenceLevel,asymptotics,symmetric}/.Flatten[{options,Options[pointIdentifiedCR]}];

(*This block sets variables that are slightly different for each of the two asymptotics.subNormalization is the standardization multiplier for the subsamples,fullNormalization is the multiplier for the construction of the final confidence interval from all of the subsamples.*)

Switch[asymp,
nests,
	subNormalization=(ssSize)^(1/3);
	fullNormalization=(Length[Union[Flatten[groupIDs,1]]])^(1/3);
	nextRandomSubsample:=generateRandomSubsample[ssSize,groupIDs,dataArray];,
coalitions,
	subNormalization=(ssSize)^(1/2);
	fullNormalization=(Length[Union[Flatten[groupIDs,1]]])^(1/2);(*(numCoalitions[matchIdxMtx])^(1/2);*)
	nextRandomSubsample:=generateRandomSubsample[ssSize,groupIDs,dataArray];
];

(*TCHRONIS for precomputed*)

estimates=Table[0,{numSubsamples}];(*List of standardized subsample estimates.*)
ssEstimates=Table[0,{numSubsamples}];(*List of raw subsample estimates.*)

Do[ssDataArray=Transpose@nextRandomSubsample;
ssEstimate=maximize[ssDataArray,Length[ssDataArray[[1]]],method,permuteinvariant,False][[2]];
	ssEstimates[[i]]=ssEstimate;
	estimates[[i]]=subNormalization (ssEstimate-pointEstimate);
If[progress>0&&Mod[i,progress]==0,WriteString["stdout","\n","Iterations completed:",ToString[i],"\n"](*TCHRONIS*)(*Print["Iterations completed: "<>ToString[i]]*)];
Block[{},subsampleMonitor/.Flatten[{options,Options[pointIdentifiedCR]}]];,
{i,1,numSubsamples}];

For[i=1,i<=Length[args],i=i+1,Histogram[estimates[[All,i]]]];
alpha=1-confLevel;

(*TCHRONIS both symmetric and assymetric to appear in the results*)
{
{{"Symmetric case",(*For the symmetric case,we want to add and subtract the 1-alpha'th quantile from the point estimate.We take the Abs here for simplicity:tn*Abs[x-y]\[Equal]Abs[tn*(x-y)]*)
cr=Table[pointEstimate[[i]]-Reverse[{-1,1}*Quantile[Abs[estimates[[All,i]]],{1-alpha,1-alpha}]]/fullNormalization,{i,1,Length[args]}];
Thread[header[[4+1;;-2]]->cr]
},
{"Asymmetric case",
(*For the asymmetric case we separately take the alpha/2 and 1-alpha/2 quantiles and subtract them.Keep in mind that since estimates has its mean subtracted,only in freakishly unlikely cases will these two have the same sign.This is not true for the symmetric case.*)
cr=Table[pointEstimate[[i]]-Reverse[Quantile[estimates[[All,i]],{alpha/2,1-alpha/2}]]/fullNormalization,{i,1,Length[args]}];
Thread[header[[4+1;;-2]]->cr]
}
},estimates}

(*
If[sym\[Equal]True,(*For the symmetric case,we want to add and subtract the 1-alpha'th quantile from the point estimate.We take the Abs here for simplicity:tn*Abs[x-y]\[Equal]Abs[tn*(x-y)]*)cr=Table[pointEstimate[[i]]-Reverse[{-1,1}*Quantile[Abs[estimates[[All,i]]],{1-alpha,1-alpha}]]/fullNormalization,{i,1,Length[args]}],(*Else*)(*For the asymmetric case we separately take the alpha/2 and 1-alpha/2 quantiles and subtract them.Keep in mind that since estimates has its mean subtracted,only in freakishly unlikely cases will these two have the same sign.This is not true for the symmetric case.*)cr=Table[pointEstimate[[i]]-Reverse[Quantile[estimates[[All,i]],{alpha/2,1-alpha/2}]]/fullNormalization,{i,1,Length[args]}]];
{cr,estimates}
*)

]


(* ::Input::Initialization:: *)
ClearAll[setIdentifiedCR];
Options[setIdentifiedCR]={progressUpdate->0,confidenceLevel->.95,asymptotics->nests,subsampleMonitor->Null,symmetric->False};

setIdentifiedCR[ssSize_,numSubsamples_,pointEstimate_,objFunc_,args_,groupIDs_,dataArray_,method_,permuteinvariant_,options___?OptionQ]:=
Module[{progress,confLevel,ineq,asymp,totalGroups,subsamples,estimates,alpha,dhat,sym},{progress,confLevel,asymp,sym}={progressUpdate,confidenceLevel,asymptotics,symmetric}/.Flatten[{options,Options[setIdentifiedCR]}];
totalGroups=Length[Union[Flatten[groupIDs,1]]];

qwp[x__?NumericQ]:=objFunc[dataArray,x];
(*This block sets variables that are slightly different for each of the two asymptotics.subNormalization is the standardization multiplier for the subsamples,fullNormalization is the multiplier for the construction of the final confidence interval from all of the subsamples.*)
Switch[asymp,nests,subNormalization=ssSize^(-1/3);
fullNormalization=totalGroups^(-1/3);
subsamples=Table[generateRandomSubsample[ssSize,groupIDs,dataArray],{numSubsamples}];,coalitions,subNormalization=1;
fullNormalization=1;
subsamples=Table[generateRandomSubsample[ssSize,groupIDs,dataArray],{numSubsamples}];];
(*This generates[numSubsamples] subsample estimates.They will be reused throughout the rest of the procedure.*)estimates=Table[0,{numSubsamples}];
Do[Print[i];
estimates[[i]]=maximize[subsamples[[i]],Length[subsamples[[i]][[1,1]]],False,method,permuteinvariant][[2]];
If[progress>0&&Mod[i,progress]==0,Print["Iterations completed: "<>ToString[i]]];
Print[estimates[[i]]],{i,1,numSubsamples}];
(*This is the "dhat" quantile function from Shaikh's paper.It returns the 1-alpha'th quantile of the differences of the objective function between x and the estimate value,standardized.*)dhat[x__?NumericQ]:=Module[{objDiff,coverageStatistics},objDiff=Table[(-objFunc[subsamples[[i]],x]- -objFunc[subsamples[[i]],Sequence@@(estimates[[i]])]),{i,Length[subsamples]}];
coverageStatistics=subNormalization*objDiff;
Quantile[coverageStatistics,1-alpha]];
alpha=1-confLevel;
(*For each dimension (length[args]),this does a one-dimensional minimization and a maximization from a point known to be in the confidence region,so that we can find the dimension-by-dimension extents of the region.*)initialRegion=MapIndexed[List[#1,pointEstimate[[First[#2]]]-2,pointEstimate[[First[#2]]]+2]&,args];
Print["Initial Region = ",initialRegion];
cr={};
Print["argsno=",Length[args]];

(*
Print[
{args[[1]],fullNormalization*(-objFunc[dataArray,x1,x2,x3,x4,x5]- -qwp[Sequence@@pointEstimate])\[LessEqual]dhat[Sequence@@args]}

];
*)

Do[Print["#",i];

clow=args[[i]]/.NMinimize[{args[[i]],fullNormalization*(-objFunc[dataArray,Sequence@@args]- -qwp[Sequence@@pointEstimate])<=dhat[Sequence@@args]},initialRegion,StepMonitor:>Print["Minimize candidate value (clow): "<>ToString[args[[i]]]],Method->"DifferentialEvolution",AccuracyGoal->1,PrecisionGoal->1][[2]];


chigh=args[[i]]/.NMaximize[{args[[i]],fullNormalization*(-objFunc[dataArray,Sequence@@args]- -qwp[Sequence@@pointEstimate])<=dhat[Sequence@@args]},initialRegion,StepMonitor:>Print["Maximize candidate value (chigh): "<>ToString[args[[i]]]],Method->"DifferentialEvolution",AccuracyGoal->1,PrecisionGoal->1][[2]];
AppendTo[cr,{clow,chigh}];


,{i,1,Length[args]}];


If[sym,(*We take the maximum distance from the point estimate to either end,and then make the CR have that distance from the point estimate.*)dist=Abs[cr-pointEstimate];
deltas=Map[Max,dist];
pointEstimate+Table[{-1,1},{Length[pointEstimate]}]*deltas,(*Else*)cr]]


(* ::Input::Initialization:: *)
ClearAll[sampleSetIdentifiedCR];
Options[sampleSetIdentifiedCR]={progressUpdate->0,confidenceLevel->.95,asymptotics->nests,subsampleMonitor->Null,samplingVariance->20};

sampleSetIdentifiedCR[ssSize_,numSubsamples_,numSamplePoints_,pointEstimate_,objFunc_,args_,groupIDs_,dataArray_,method_,permuteinvariant_,options___?OptionQ]:=

Module[{progress,confLevel,ineq,asymp,totalGroups,subsamples,estimates,alpha,dhat,svar},{progress,confLevel,asymp,svar}={progressUpdate,confidenceLevel,asymptotics,samplingVariance}/.Flatten[{options,Options[sampleSetIdentifiedCR]}];

totalGroups=Length[Union[Flatten[groupIDs,1]]];

qwp[x__?NumericQ]:=objFunc[dataArray,x];

(*This block sets variables that are slightly different for each of the two asymptotics.subNormalization is the standardization multiplier for the subsamples,fullNormalization is the multiplier for the construction of the final confidence interval from all of the subsamples.*)
Switch[asymp,
nests,
	subNormalization=ssSize^(-1/3);
	fullNormalization=totalGroups^(-1/3);
	subsamples=Table[generateRandomSubsample[ssSize,groupIDs,dataArray],{numSubsamples}];,
coalitions,
	subNormalization=1;
	fullNormalization=1;
	subsamples=Table[generateRandomSubsample[ssSize,groupIDs,dataArray],{numSubsamples}];];

(*This generates[numSubsamples] subsample estimates.They will be reused throughout the rest of the procedure.*)
estimates=Table[0,{numSubsamples}];

Do[estimates[[i]]=(*args/.pairwiseMSE[objFunc,subsamples[[i]],args,options][[2]];*)
maximize[subsamples[[i]],Length[subsamples[[i]][[1,1]]],False,method,permuteinvariant][[2]];
If[progress>0&&Mod[i,progress]==0,
	Print["Iterations completed: "<>ToString[i]]];,
{i,1,numSubsamples}];

(*This is the "dhat" quantile function from Shaikh's paper.It returns the 1-alpha'th quantile of the differences of the objective function between x and the estimate value,standardized.*)

dhat[x__?NumericQ]:=
	Module[{objDiff,coverageStatistics},
		objDiff=Table[(-objFunc[subsamples[[i]],x]
				   - -objFunc[subsamples[[i]],Sequence@@(estimates[[i]])]),
					{i,Length[subsamples]}];
		coverageStatistics=subNormalization*objDiff;
		Quantile[coverageStatistics,1-alpha]
	];

(*This alternate version of dhat is worth remembering.dhat[x__?NumericQ]:=Module[{coverageStatistics},coverageStatistics=subNormalization*MapIndexed[(-objFunc[#1,x]--objFunc[#1,Sequence@@(estimates[[First[#2]]])])&,subsamples];
Quantile[coverageStatistics,1-alpha]];*)

isInCR[point_]:=fullNormalization*(-qwp[Sequence@@point]- -qwp[Sequence@@pointEstimate])<=dhat[Sequence@@point];

alpha=1-confLevel;

cr={};

Do[pt=RandomVariate[MultinormalDistribution[pointEstimate,svar*IdentityMatrix[Length[args]]]];
If[isInCR[pt],
AppendTo[cr,pt];
],
{i,1,numSamplePoints}];
cr

]