-
Notifications
You must be signed in to change notification settings - Fork 4
/
confidence.m
132 lines (107 loc) · 7.34 KB
/
confidence.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
(* ::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:: *)
Echo["Loaded confidence.m"];
(* ::Input::Initialization:: *)
Off[General::compat];
Needs["Combinatorica`"](*For RandomKSubset*)
(* ::Input::Initialization:: *)
ClearAll[generateRandomGroups];
generateRandomGroups::usage="generateRandomGroups[ssSize,groupIDs] generates groups sample.";
generateRandomGroups[ssSize_,groupIDs_]:=Module[{totalGroups,groups},
totalGroups=Union[Flatten[groupIDs,1]];
groups=RandomKSubset[totalGroups,ssSize];
WriteString["stdout"," ",groups];(*TCHRONIS for verbose results*)
groups
]
(* ::Input::Initialization:: *)
ClearAll[generateRandomSubsample];
generateRandomSubsample::usage="generateRandomSubsample[ssSize,groupIDs,dataArray] generates a subsample of a given size from a data array.";
generateRandomSubsample[groups_,groupIDs_,dataArray_]:=Module[{totalGroups,qualifiedIndexes},
(*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,useSavedGroups->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_:True,options___?OptionQ]:=Module[{progress,confLevel,asymp,ssDataArray,estimates,estimate,alpha,cr,sym,useSavedG,nextRandomGroup,nextRandomSubsample},{progress,confLevel,asymp,sym,useSavedG}={progressUpdate,confidenceLevel,asymptotics,symmetric,useSavedGroups}/.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);
,
coalitions,
subNormalization=(ssSize)^(1/2);
fullNormalization=(Length[Union[Flatten[groupIDs,1]]])^(1/2);(*(numCoalitions[matchIdxMtx])^(1/2);*)
];
(*TCHRONIS for precomputed*)
estimates=Table[0,{numSubsamples}];(*List of standardized subsample estimates.*)
ssEstimates=Table[0,{numSubsamples}];(*List of raw subsample estimates.*)
(*randomGroups global variable*)
If[useSavedG,
randomGroups = savedGroups,
randomGroups =Table[generateRandomGroups[ssSize,groupIDs],{numSubsamples}]
];
Do[
nextRandomGroup=randomGroups[[i]];
nextRandomSubsample =generateRandomSubsample[nextRandomGroup,groupIDs,dataArray];
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}
*)
];
Information[pointIdentifiedCR,LongForm->False]