Estimating class birth year

Functions

elementcount[list_] :=
Module[{elements,i},
elements = Union[list];
Table[
{elements[[i]],Count[list,elements[[i]]]},
{i,1,Length[elements]}
]
];

Import SSA CA names

names =
Import[
"/Users/dmeyer/Teaching/111fall11/names/CA.txt",
"CSV"
];

Length[names]
names[[1]]
Length[names[[1]]]

Format SSA CA names

Do[
totalpop[year,"F"] = 0;
totalpop[year,"M"] = 0,
{year,1910,2010}
];
namecount[year_,sex_,name_]=0;
counts[year_,sex_]={};
Fnames1990 = {};
Mnames1990 = {};
Do[
{namestring,count} = names[[i]];
namelength = StringLength[namestring]-8;
year = ToExpression[StringTake[namestring,{5,8}]];
sex =StringTake[namestring,{4}];
AppendTo[counts[year,sex],count];
totalpop[year,sex] = totalpop[year,sex] + count;
name =  StringTake[namestring,-namelength];
namecount[year,sex,name]=count;
If[year == 1990,
If[sex == "F",
AppendTo[Fnames1990,{totalpop[year,sex],name}],
AppendTo[Mnames1990,{totalpop[year,sex],name}]
]
],
{i,1,Length[names]}
];
(*
Fnames1990
Mnames1990
*)

Select random names for 1990

Fsimnames = {};
Do[
f = RandomInteger[{1,totalpop[1990,"F"]}];
n = 1;
While[f>Fnames1990[[n,1]],n = n + 1];
AppendTo[Fsimnames,Fnames1990[[n,2]]],
{Fcount,1,16}
];
Fsimnames = Sort[Fsimnames];
Fsimnamescounts = elementcount[Fsimnames] 

Msimnames = {};
Do[
m = RandomInteger[{1,totalpop[1990,"M"]}];
n = 1;
While[m>Mnames1990[[n,1]],n = n + 1];
AppendTo[Msimnames,Mnames1990[[n,2]]],
{Mcount,1,36}
];
Msimnames = Sort[Msimnames];
Msimnamescounts = elementcount[Msimnames]

Look up simulated name counts in SSA CA names for year range

y1 = 1988;y2 = 1996;
y1 = 1981;y2 = 2001;
y1 = 1970;y2 = 2010;

Frarenamepositions = {};
Mrarenamepositions = {};
Do[
Fsimcounts = 
Table[
namecount[year,"F",Fsimnamescounts[[i,1]]],
{i,1,Length[Fsimnamescounts]}
];
Fzeropositions = Position[Fsimcounts,0];
Print[year," ",Fzeropositions];
AppendTo[Frarenamepositions,Fzeropositions];
Msimcounts = 
Table[
namecount[year,"M",Msimnamescounts[[i,1]]],
{i,1,Length[Msimnamescounts]}
];
Mzeropositions = Position[Msimcounts,0];
Print[year," ",Mzeropositions];
AppendTo[Mrarenamepositions,Mzeropositions],
{year,y1,y2}
];

Frarenamepositions = Union[Flatten[Frarenamepositions]]
Mrarenamepositions = Union[Flatten[Mrarenamepositions]]

Remove names with 0 counts in SSA CA names in year range

Fcommonsimnamescounts = {};
Do[
If[!MemberQ[Frarenamepositions,i],
AppendTo[Fcommonsimnamescounts,Fsimnamescounts[[i]]]
],
{i,1,Length[Fsimnamescounts]}
];
Fcommonsimnamescounts

Mcommonsimnamescounts = {};
Do[
If[!MemberQ[Mrarenamepositions,i],
AppendTo[Mcommonsimnamescounts,Msimnamescounts[[i]]]
],
{i,1,Length[Msimnamescounts]}
];
Mcommonsimnamescounts

Clear[Fcounts];Clear[Mcounts];
Do[
Fcounts[year] = 
Table[
namecount[year,"F",Fcommonsimnamescounts[[i,1]]],
{i,1,Length[Fcommonsimnamescounts]}
];
Print[year," ",Position[Fcounts[year],0]];
Mcounts[year] = 
Table[
namecount[year,"M",Mcommonsimnamescounts[[i,1]]],
{i,1,Length[Mcommonsimnamescounts]}
];
Print[year," ",Position[Mcounts[year],0]],
{year,1910,2010}
]

Compute Pr(year|simulated names) and plot

prob1[year_] := 
Which[
y1≤year≤y2,1/(y2-y1+1),
True,0
];
prob2[year_] :=
Module[{norm,y},
norm = Sum[E^(-(y-1990)^2/(2*4)),{y,y1,y2}];
Which[
y1≤year≤y2,E^(-(year-1990)^2/(2*4))/norm,
True,0
]
];

probsimNgY[year_] := 
(Product[
Binomial[Fcounts[year][[i]],Fcommonsimnamescounts[[i,2]]],
{i,1,Length[Fcommonsimnamescounts]}
]/
Binomial[
totalpop[year,"F"],
Sum[Fcommonsimnamescounts[[i,2]],{i,1,Length[Fcommonsimnamescounts]}]
]
)*
(Product[
Binomial[Mcounts[year][[i]],Mcommonsimnamescounts[[i,2]]],
{i,1,Length[Mcommonsimnamescounts]}
]/
Binomial[
totalpop[year,"M"],
Sum[Mcommonsimnamescounts[[i,2]],{i,1,Length[Mcommonsimnamescounts]}]
]
);

probYgsimN1[year_] :=
probsimNgY[year]*prob1[year]/
Sum[probsimNgY[y]*prob1[y],{y,y1,y2}];

plotprior1 =
Table[{{y},prob1[y]},{y,y1,y2}];
plotdata1 =
Table[{{y},N[probYgsimN1[y]]},{y,y1,y2}]
maxy = y1; maxp = plotdata1[[1,2]];
Do[
If[plotdata1[[i,2]]> maxp,
{{maxy},maxp} = plotdata1[[i]]
],
{i,2,Length[plotdata1]}
];
Print[maxy," ",maxp];

DateListPlot[
{plotprior1,plotdata1},
PlotStyle->{Directive[Blue,Thick],Directive[Red,Thick]},
FrameTicksStyle->Directive[18],
Joined->True,
PlotRange->All
]

probYgsimN2[year_] :=
probsimNgY[year]*prob2[year]/
Sum[probsimNgY[y]*prob2[y],{y,y1,y2}];

plotprior2 =
Table[{{y},prob2[y]},{y,y1,y2}];
plotdata2 =
Table[{{y},N[probYgsimN2[y]]},{y,y1,y2}]

DateListPlot[
{plotprior2,plotdata2},
PlotStyle->{Directive[Blue,Thick],Directive[Red,Thick]},
FrameTicksStyle->Directive[18],
Joined->True,
PlotRange->All
]

Run multiple simulations

Nsims = 10000;
yearcount[year_]=0;
Do[
Fsimnames = {};
Do[
f = RandomInteger[{1,totalpop[1990,"F"]}];
n = 1;
While[f>Fnames1990[[n,1]],n = n + 1];
AppendTo[Fsimnames,Fnames1990[[n,2]]],
{Fcount,1,16}
];
Fsimnames = Sort[Fsimnames];
Fsimnamescounts = elementcount[Fsimnames]; 
Msimnames = {};
Do[
m = RandomInteger[{1,totalpop[1990,"M"]}];
n = 1;
While[m>Mnames1990[[n,1]],n = n + 1];
AppendTo[Msimnames,Mnames1990[[n,2]]],
{Mcount,1,36}
];
Msimnames = Sort[Msimnames];
Msimnamescounts = elementcount[Msimnames];
Frarenamepositions = {};
Mrarenamepositions = {};
Do[
Fsimcounts = 
Table[
namecount[year,"F",Fsimnamescounts[[i,1]]],
{i,1,Length[Fsimnamescounts]}
];
Fzeropositions = Position[Fsimcounts,0];
AppendTo[Frarenamepositions,Fzeropositions];
Msimcounts = 
Table[
namecount[year,"M",Msimnamescounts[[i,1]]],
{i,1,Length[Msimnamescounts]}
];
Mzeropositions = Position[Msimcounts,0];
AppendTo[Mrarenamepositions,Mzeropositions],
{year,y1,y2}
];
Frarenamepositions = Union[Flatten[Frarenamepositions]];
Mrarenamepositions = Union[Flatten[Mrarenamepositions]];
Fcommonsimnamescounts = {};
Do[
If[!MemberQ[Frarenamepositions,i],
AppendTo[Fcommonsimnamescounts,Fsimnamescounts[[i]]]
],
{i,1,Length[Fsimnamescounts]}
];
Mcommonsimnamescounts = {};
Do[
If[!MemberQ[Mrarenamepositions,i],
AppendTo[Mcommonsimnamescounts,Msimnamescounts[[i]]]
],
{i,1,Length[Msimnamescounts]}
];
Clear[Fcounts];Clear[Mcounts];
Do[
Fcounts[year] = 
Table[
namecount[year,"F",Fcommonsimnamescounts[[i,1]]],
{i,1,Length[Fcommonsimnamescounts]}
];
Mcounts[year] = 
Table[
namecount[year,"M",Mcommonsimnamescounts[[i,1]]],
{i,1,Length[Mcommonsimnamescounts]}
],
{year,y1,y2}
]; 
probsimNgY[year_] := 
(Product[
Binomial[Fcounts[year][[i]],Fcommonsimnamescounts[[i,2]]],
{i,1,Length[Fcommonsimnamescounts]}
]/
Binomial[
totalpop[year,"F"],
Sum[Fcommonsimnamescounts[[i,2]],{i,1,Length[Fcommonsimnamescounts]}]
]
)*
(Product[
Binomial[Mcounts[year][[i]],Mcommonsimnamescounts[[i,2]]],
{i,1,Length[Mcommonsimnamescounts]}
]/
Binomial[
totalpop[year,"M"],
Sum[Mcommonsimnamescounts[[i,2]],{i,1,Length[Mcommonsimnamescounts]}]
]
);
probYgsimN1[year_] :=
probsimNgY[year]*prob1[year]/
Sum[probsimNgY[y]*prob1[y],{y,y1,y2}];
plotdata1 =
Table[{{y},N[probYgsimN1[y]]},{y,y1,y2}];
maxy = y1; maxp = plotdata1[[1,2]];
Do[
If[plotdata1[[i,2]]> maxp,
{{maxy},maxp} = plotdata1[[i]]
],
{i,2,Length[plotdata1]}
];
yearcount[maxy] = yearcount[maxy]+1,
{count,1,Nsims}
];

plotdata =
Table[{{y},yearcount[y]/Nsims},{y,y1,y2}];

DateListPlot[
plotdata,
PlotStyle->{Red,Thick},
FrameTicksStyle->Directive[18],
Joined->True,
PlotRange->All
]