from math import sqrt, log, exp
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

gs = np.arange(0,2,0.01)

# Function u(g)
def u(g):
  return g/(1+g)

# Plot
plt.plot(gs, map(lambda g: u(g), gs))
plt.xlabel('g')
plt.ylabel('u(g)')
plt.title('Gun access and use function u(g)')
def stop(g, s, p):
  return 1-(1-s*u(g))**p

"""
Explanation: A well-regulated militia, being necessary to the security of a free state, the right of the people to keep and bear arms shall not be infringed
Second amendment to the US constitution
Who Has Gun Control Right?
<br/>
<font color='red'>This is work in progress - please disregard this draft and check back later for the final version</font>
<br />
By Gilles Pirio October 29, 2015
Gun control is a bitterly debatted topic in today's America. Both sides come with their own certainties -- often unproven and ideological -- about the effect of gun control on society (or the lack thereof). 
The proponent of gun control claims that more guns always induces more violence and death. The opponents, like the NRA, advocate for more people to carry guns to protect themselves. Shooting(!) matches aside, both sides simply refuse to talk to each other.
Unlawful killings must be prevented. But to do so, an understanding of what is really happening is necessary. Even though it is rather morbid, the following models tries to investigate how gun ownership affects the number of unlawful killings in a given country.
Modeling
Let's try to model the amount of victim of unlawful killings by firearms in a given country (or state) in a given year. Even though they understandably take a lot of space in the mainstream media, unlawful killings by security forces like the police are considered negligeable. 
In any given year, some individuals will infortunately try to kill one (like a husband shooting down his ex-wife) or more than one (like the horrific events that happened in US schools in the past years). We note the amount of people a killer intends to eliminate p. The number of victims for a given p, V(p) can be written as:
$$ V(p) = P \, e(p) \, u(g) \, d(p,g) $$
Where d(p,g) is the expected number of death during a killing event:
$$ d(p,g) = p \, r(p,g) $$
P is the total population. e(p) is the propability for a given person in the general population would take action and commit murder and attempt to kill p persons. g is the gun ownership rate, u(g) is the rate of people who have access to guns (own one, or can easily get one) and know how to use it. Finally, r(p,g) is the fatality rate during a event.
Estimating u(g) is a difficult task. It is pretty clear the function is concave. We will use the simple following formula, that we will also plot:
$$ u(g) =\frac{g}{1+g} $$
End of explanation
"""


# Define r and d functions
def r(p, g, s):
  return (1-s*c*u(g))**p
def d(p, g, s):
  return p*r(p, g, s)

"""
Explanation: This matches the US data since in the US there is about 1 gun for every inhabitant (g = 1.0) while the number of people owning a gun is 37%. Even though our value for u(1.0) is 0.5 so above 0.37, it is consistent since it is possible for somebody not owing a gun to easily find one - and learn how to use one. As expected, for small g, we have:
$$ u(g) \approx g $$
We now need to model r(p,g). Let s the probability for someone with a gun present during an event to stop the killer and prevent any killing. Obviously there is an approximation here as the killer could be partially successful. However taking this in consideration would not change the outcome much (TODO: confirm this). We will also need c, the rate of people who have a gun with them all the time(open carry or conceal carry). We will model a situation where 50% of people who own a gun carry it outside their home and a situation where only 10% do (c= 0.5 and c= 0.1, respectively).
The probability for the killer to be stopped becomes:
$$ 1-(1-s\,c\,u(g))^p $$
And therefore:
$$ r(p,g) = (1-s\,c\,u(g))^p $$
$$ d(p,g) = u(g) \, p \, (1-s\,c\,u(g))^p $$
End of explanation
"""


c, s, p = 0.5, 0.5, 1
plt.plot(gs, map(lambda g: 100*u(g)*d(p, g, s), gs), label= "p=%i, c=%.1f" % (p, c))
p = 10
plt.plot(gs, map(lambda g: 100*u(g)*d(p, g, s), gs), label= "p=%i, c=%.1f" % (p, c))
c = 0.1
plt.plot(gs, map(lambda g: 100*u(g)*d(p, g, s), gs), label= "p=%i, c=%.1f" % (p, c))
plt.xlabel('g')
plt.ylabel('Expected deaths')
plt.title('Expected death for assuming s=0.5')
plt.legend(loc= 'upper left')

"""
Explanation: Let's put it all together and plot the number of expected deaths assuming there are 100 people in the a given country who are willing to take the life of 1 (or 10) people:
End of explanation
"""


s = 0.5
p = 1
c = 0.5
ps = range(1,101)
plt.plot(gs, map(lambda g: sum(map(lambda p: 1000/p*u(g)*d(p, g, s), ps)), gs), label= "p=%i, c=%.1f" % (p, c))
c = 0.1
plt.plot(gs, map(lambda g: sum(map(lambda p: 1000/p*u(g)*d(p, g, s), ps)), gs), label= "p=%i, c=%.1f" % (p, c))
plt.xlabel('g')
plt.ylabel('Deaths (arbitrary scale)')
plt.title('Total number of deaths')
plt.legend(loc='best')

"""
Explanation: The number of expected deaths is an increasing function of g for p = 1 but not for p = 10. This reflects the two competing trends: 
* For small g, the limiting factor is availability of a gun and the number of death will be proportional to p*g. The less gun, the less death.
* For higher g, the main limiting factor will come from the fact that the probability that someone could twart the murders using his/her own gun is also higher. The NRA, therefore has a point -- in certain specific circonstances, more guns may mean less death. 
Now let's look at the total number of unlawful killings V, for all p:
$$ V = \sum_{p=1} V(p) = P \, \sum_{p=1} e(p)\,u(g)\,d(p,g) $$
e(p) is not known but it is very likely that e(p) decreases quickly with p -- most of the murders target only 1 person. If we assume that:
$$ e(p) = \frac{1000}{p}  $$
The total death toll can be shown below:
End of explanation
"""


s = 0.5
p = 1
c = 0.5
ps = range(1,101)
plt.plot(gs, map(lambda g: sum(map(lambda p: 1000/p**2*u(g)*d(p, g, s), ps)), gs))
plt.xlabel('g')
plt.ylabel('Deaths (arbitrary scale)')
plt.title('Total number of deaths')

"""
Explanation: In that case, there is definitely a 2 different domains:
* For small g, reducing the amount of guns in circulation will reduce the total deaths
* For g above 0.15, the reverse is true in certain conditions (if the amount of guns in circulation is high)
If we now assume that:
$$ e(p) = \frac{1000}{p^2}  $$
The total death toll can be shown below:
End of explanation
"""


