######################################################################
##DDH.txt: Save this file as  DDH.txt                                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read DDH.txt<Enter>                                #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: 
 
print(`Created:  Feb. 6, 2016`):
print(` This is DDH.txt `):
print(`to investigate probability models for Detecting Discrimination in Hiring`):
print(`using an urn model`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):


print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):


ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  PowerG1, PowerP1, PowerSD1  `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  PowerGcurve, PowerPcurve, PowerSDcurve, RejRegG, RejRegP, RejRegSD, Story `):
 print(` `):




elif nops([args])=1 and op(1,[args])=PowerGcurve then
print(`PowerGcurve(N,P): The Power Curve for `):
print(`the Government methodoogy`):
print(` inputs `):
print(` (i)N: the number of people hired `):
print( ` (ii)P: the proportion of qualified minority members in the applicant pool `):
print(` Outputs: the power curve indicating the probability of detecting bias if it is actually there with `):
print(` given bias ratio (P-p)/P from 0 to 1. Try: `):
print(` PowerGcurve(25,0.5); `):
print(` PowerGcurve(16,0.43); `):


elif nops([args])=1 and op(1,[args])=PowerPcurve then
print(`PowerPcurve(N,P,alpha): The Power Curve for `):
print(`the Princeton methodoogy`):
print(` inputs `):
print(` (i) N: the number of people hired `):
print( ` (ii) P: the proportion of qualified minority members in the applicant pool `):
print(` (iii) alpha: the confidence level `):
print(` Outputs: the power curve indicating the probability of detecting bias if it is actually there with `):
print(` given bias ratio (P-p)/P from 0 to 1. Try: `):
print(` PowerPcurve(25,0.5, 1/40); `):
print(` PowerPcurve(16,0.43, 1/40); `):


elif nops([args])=1 and op(1,[args])=PowerSDcurve then
print(`PowerSDcurve(N,P,a): The Power Curve for `):
print(`the Government methodoogy`):
print(` inputs `):
print(` (i) N: the number of people hired `):
print( ` (ii) P: the proportion of qualified minority members in the applicant pool `):
print( ` (iii) a: a real number indicating how many standard deviations `):
print(` Outputs: the power curve indicating the probability of detecting bias if it is actually there, using the constant standard deviation method `):
print(` given bias ratio (P-p)/P from 0 to 1. Try: `):
print(` PowerSDcurve(25,0.5,1); `):
print(` PowerSDcurve(16,0.43,1); `):

elif nops([args])=1 and op(1,[args])=PowerP1 then
print(`PowerP1(N,P,alpha,p): The probability of correctly detecting bias according`):
print(` to the Princeton methodoogy `):
print(` inputs `):
print(` (i)N: the number of people hired `):
print(` (ii)P: the proportion of qualified minority members in the applicant pool `):
print(` (iii) alpha: the confidence level `):
print(` (iv): the actual (biased, less than P) probability of hiring a minority `):
print(` Outputs: the probabilty of detecting bias using the Princeton methodlogy. Try: `):
print(` PowerP1(25,0.5, 1/40,.3); `):
print(`PowerP1(16,0.43, 1/40,.3);`):

elif nops([args])=1 and op(1,[args])=PowerG1 then
print(`PowerG1(N,P,p): The probability of correctly detecting bias according`):
print(` to the Princeton methodoogy `):
print(` inputs `):
print(` (i)N: the number of people hired `):
print(` (ii)P: the proportion of qualified minority members in the applicant pool `):
print(` (iii): the actual (biased, less than P) probability of hiring a minority `):
print(` Outputs: the probabilty of detecting bias using the Princeton methodlogy. Try: `):
print(` PowerG1(25,0.5,.3); `):
print(`PowerG1(16,0.43,.3);`):



elif nops([args])=1 and op(1,[args])=PowerSD1 then
print(`PowerSD1(N,P,p,a): The probability of correctly detecting bias accordinge if it exists if the prob. of bias  is p.`):
print(`using the constant standard-deviation methodology`):
print(` (i)N: the number of people hired `):
print(` (ii)P: the proportion of qualified minority members in the applicant pool `):
print(` (iii): the actual (biased, less than P) probability of hiring a minority `):
print(` (iv) a real numer a such that the min is N*P-a*sqrt(P*(1-P)) `):
print(` Outputs: the probabilty of detecting bias using the Contant Standard Deviation methodlogy. Try: `):
print(` PowerSD1(25,0.5,.3,1); `):
print(`PowerSD1(16,0.43,.3,1);`):


elif nops([args])=1 and op(1,[args])=RejRegG then
print(`RejRegG(N,P): Detecting Discrimination according to the Government methodology`):
print(` (i) N: the number of people hired `):
print(` (ii) P: the proportion of qualified minority members in the applicant pool `):
print(` Outputs: the number, let's call it k, such that if the department hires `):
print(` k or less people then  one can deduce that there is discrimination  `):
print(` RejRegG(5,0.1); `):
print(` RejRegG(25,0.5); `):
print(` RejRegG(16,0.43); `):

elif nops([args])=1 and op(1,[args])=RejRegP then
print(`RejRegP(N,P,alpha): inputs`):
print( ` (i) N: the number of people hired `):
print(` (ii) P: the proportion of qualified minority members in the applicant pool `):
print(` (iii) alpha: the confidence level `):
print(` Outputs: the number, let's call it k, such that if the department hires `):
print(` k or less people then  one can deduce that there is discrimination with `):
print(` confidence level alpha, i.e. if the department has no discrimination the `):
print(` probability that hires k or less people is small, namely alpha/2 `):
print(` (it is a one-tailed test). For example, try: `):
print(` RejRegP(5,0.1, 1/40); `):
print(` RejRegP(25,0.5, 1/40); `):
print(` RejRegP(16,0.43, 1/40); `) :


elif nops([args])=1 and op(1,[args])=RejRegSD then
print(`RejRegSD(N,P,a): Detecting Discrimination according to the  Constant Standard-Deviation`):
print(` methodology where one has to be off by at most a standard deviations`):
print(` (i)N: the number of people hired `):
print(` (ii)P: the proportion of qualified minority members in the applicant pool `):
print(` (iii) The number of standard deviations `):
print(` Outputs: the number, let's call it k, such that if the department hires `):
print(` k or less people then  one can deduce that there is discrimination  `):
print(` RejRegSD(25,0.5,1); `):
print(` RejRegSD(16,0.43,1); `):

elif nops([args])=1 and op(1,[args])=Story then
print(` Story(): comparative study with many scenarios where N goes from 5 to 100 incremented by 10 and `):
print(`p goes from 0.1 to 0.5 incremented by 0.1 and for the Princeton method alpha with 1/40 and 1/100`):
print(`and for the Constant Deviation method a=0.3,0.6,0.9`):


else
print(`There is no ezra for`,args):
fi:
 
end:


#RejRegP(N,P,alpha): Detecting Discrimination according to the Princeton methodology
#inputs
#(i)N: the number of people hired
#(ii)P: the proportion of qualified minority members in the applicant pool
#(iii) alpha: the confidence level
#Outputs: the number, let's call it k, such that if the department hires
#k or less people then  one can deduce that there is discrimination with
#confidence level alpha, i.e. if the department has no discrimination the
#probability that hires k or less people is small, namely alpha/2
#(it is a one-tailed test). For example, try:
#RejRegP(5,0.1, 1/40);
#RejRegP(25,0.5, 1/40);
#RejRegP(16,0.43, 1/40);
RejRegP:=proc(N,P,alpha) local k,i:
for k from 0 to N while add(binomial(N,i)*P^i*(1-P)^(N-i),i=0..k)<=alpha do od:
max(k-1,0):
end:

#RejRegG(N,P): Detecting Discrimination according to the Government methodology
#(i)N: the number of people hired
#(ii)P: the proportion of qualified minority members in the applicant pool
#Outputs: the number, let's call it k, such that if the department hires
#k or less people then  one can deduce that there is discrimination 
#RejRegG(5,0.1);
#RejRegG(25,0.5);
#RejRegG(16,0.43);
RejRegG:=proc(N,P):
if type(N*P,integer) then
 RETURN(N*P-1):
else
 trunc(N*P):
fi:
end:


#PowerP1(N,P,alpha,p): The probability of correctly detecting bias according
#to the Princeton methodoogy
#inputs
#(i)N: the number of people hired
#(ii)P: the proportion of qualified minority members in the applicant pool
#(iii) alpha: the confidence level
#(iv): the actual (biased, less than P) probability of hiring a minority
#Outputs: the probabilty of detecting bias using the Princeton methodlogy. Try:
#PowerP1(25,0.5, 1/40,.3);
#PowerP1(16,0.43, 1/40,.3);
PowerP1:=proc(N,P,alpha,p) local k,i:
k:=RejRegP(N,P,alpha):

add(binomial(N,i)*p^i*(1-p)^(N-i), i=0..k):

end:



#PowerG1(N,P,p): The probability of correctly detecting bias according
#to the Government methodoogy
#inputs
#(i)N: the number of people hired
#(ii)P: the proportion of qualified minority members in the applicant pool
#(iii): the actual (biased, less than P) probability of hiring a minority
#Outputs: the probabilty of detecting bias using the Government methodlogy. Try:
#PowerG1(25,0.5, .3);
#PowerG1(16,0.43,.3);
PowerG1:=proc(N,P,p) local k,i:
k:=RejRegG(N,P):

add(binomial(N,i)*p^i*(1-p)^(N-i), i=0..k):

end:





#PowerPcurve(N,P,alpha): The Power Curve for 
#the Princeton methodoogy
#inputs
#(i)N: the number of people hired
#(ii)P: the proportion of qualified minority members in the applicant pool
#(iii) alpha: the confidence level
#Outputs: the power curve indicating the probability of detecting bias if it is actually there with
#given bias ratio (P-p)/P from 0 to 1. Try:
#PowerPcurve(25,0.5, 1/40);
#PowerPcurve(16,0.43, 1/40);
PowerPcurve:=proc(N,P,alpha) local x:

plot([seq([x*0.01,PowerP1(N,P,alpha,P*(1-x*0.01))],x=0..100)]):

end:

#PowerGcurve(N,P): The Power Curve for 
#the Government methodoogy
#inputs
#(i)N: the number of people hired
#(ii)P: the proportion of qualified minority members in the applicant pool
#Outputs: the power curve indicating the probability of detecting bias if it is actually there with
#given bias ratio (P-p)/P from 0 to 1. Try:
#PowerGcurve(25,0.5);
#PowerGcurve(16,0.43);
PowerGcurve:=proc(N,P) local x:

plot([seq([x*0.01,PowerG1(N,P,P*(1-x*0.01))],x=0..100)]):

end:




#RejRegSD(N,P,a): Detecting Discrimination according to the Constant Standard-Deviation methodology where one has to be off by at most a standard deviations
#(i)N: the number of people hired
#(ii)P: the proportion of qualified minority members in the applicant pool
#(iii) The number of standard deviations
#Outputs: the number, let's call it k, such that if the department hires
#k or less people then  one can deduce that there is discrimination 
#RejRegSD(25,0.5,1);
#RejRegSD(16,0.43,1);
RejRegSD:=proc(N,P,a) local sig,k:
sig:=sqrt(N*P*(1-P)):
k:=trunc(N*P-a*sig):
end:

#PowerSD1(N,P,p,a): The probability of correctly detecting bias according
#to the Government methodoogy
#inputs
#(i)N: the number of people hired
#(ii)P: the proportion of qualified minority members in the applicant pool
#(iii): the actual (biased, less than P) probability of hiring a minority
#(iv): a real number a
#Outputs: the probabilty of detecting bias using the Constant Standard-Deviation methodlogy. Try:
#PowerSD1(25,0.5, .3,1);
#PowerSD1(16,0.43,.3,a);
PowerSD1:=proc(N,P,p,a) local k,i:
k:=RejRegSD(N,P,a):

add(binomial(N,i)*p^i*(1-p)^(N-i), i=0..k):

end:

#PowerSDcurve(N,P,a): The Power Curve for 
#the Constant Standard Deviation methodoogy
#inputs
#(i)N: the number of people hired
#(ii)P: the proportion of qualified minority members in the applicant pool
#(iii) a real number a
#Outputs: the power curve indicating the probability of detecting bias if it is actually there with
#given bias ratio (P-p)/P from 0 to 1. Try:
#PowerSDcurve(25,0.5,1);
#PowerSDcurve(16,0.43,1);
PowerSDcurve:=proc(N,P,a) local x:

plot([seq([x*0.01,PowerSD1(N,P,P*(1-x*0.01),a)],x=0..100)]):

end:

#Story(): comparative study with many scenarios where N goes from 5 to 100 incremented by 10 and
#p goes from 0.1 to 0.5 incremented by 0.1 and for the Princeton method alpha with 1/40 and 1/100
#and for the Constant Deviation method a=0.3,0.6,0.9

Story:=proc() local N,i,P,k,i1:

for N from 5 by 10 to 55 do
 for i from 1 to 5 do
   P:=0.1*i:

print(``):
 print(`When the number of positions is`, N):
 print(``):
 print(`and the ratio of minority is`, P):
print(``):

print(`The Princeton method would detect bias, with confidence  level 1/40, for any number of minority hires less-than-or equal to`):
print(``):
k:=RejRegP(N,P,1/40):
print( k):
print(``):

print(`Hence the probability that an innocent, unbiased department, would be labelled as biased is`):
print(``):
print(add(binomial(N,i1)*P^i1*(1-P)^(N-i1),i1=0..k)):
print(``):

print(`The corresponding power curve is`):
print(``):
print(PowerPcurve(N,P,1/40)):
print(``):


k:=RejRegP(N,P,1/100):
print(`and  with confidence  level 1/100, for any number of minority hires less-than-or equal to`):
print(``):
print( k):
print(``):
print(`Hence the probability that an innocent, unbiased department, would be labelled as biased is`):
print(``):
print(add(binomial(N,i1)*P^i1*(1-P)^(N-i1),i1=0..k)):
print(``):

print(`The corresponding power curve is`):
print(``):
print(PowerPcurve(N,P,1/100)):
print(``):

print(``):
print(`The Government method would need at most`):
print(``):
k:=RejRegG(N,P):
print( k):
print(``):
print(`Hence the probability that an innocent, unbiased department, would be labelled as biased is`):
print(``):
print(add(binomial(N,i1)*P^i1*(1-P)^(N-i1),i1=0..k)):
print(``):

print(`The corresponding power curve is`):
print(``):
print(PowerGcurve(N,P)):
print(``):

print(`The Constant Standard Deviation method with a=0.5, would need `):
print(``):
k:=RejRegSD(N,P,0.5):
print( k):

print(`Hence the probability that an innocent, unbiased department, would be labelled as biased is`):
print(``):
print(add(binomial(N,i1)*P^i1*(1-P)^(N-i1),i1=0..k)):
print(``):

print(``):
print(`The corresponding power curve is`):
print(``):
print(PowerSDcurve(N,P,0.3)):
print(``):


print(`The Constant Standard Deviation method with a=1, would need `):
print(``):
k:=RejRegSD(N,P,1):
print( k):

print(`Hence the probability that an innocent, unbiased department, would be labelled as biased is`):
print(``):
print(add(binomial(N,i1)*P^i1*(1-P)^(N-i1),i1=0..k)):
print(``):

print(``):
print(`The corresponding power curve is`):
print(``):
print(PowerSDcurve(N,P,0.6)):
print(``):


print(`The Constant Standard Deviation method with a=1.5, would need `):
print(``):
k:=RejRegSD(N,P,1.5):
print( k):

print(`Hence the probability that an innocent, unbiased department, would be labelled as biased is`):
print(``):
print(add(binomial(N,i1)*P^i1*(1-P)^(N-i1),i1=0..k)):
print(``):

print(``):
print(`The corresponding power curve is`):
print(``):
print(PowerSDcurve(N,P,0.9)):
print(``):







od:
od:

end: