Point estimates¶
We are trying to reveal some parameter of unknown distribution that we sample from. In this problem we are sampling from distribution $U(0,\theta)$. To imagine a story about this: we play a game where the amount of time we get to play is random, generated by a device set up to generate a uniform random variable from interval $[0,\theta]$ for some value of $\theta$. How to decide the value of $\theta$ from several observations?
A discrete version of this is the "German tank problem". https://en.wikipedia.org/wiki/German_tank_problem
from scipy import stats
from math import sqrt
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
The next function generates $n$ numbers from the distribution $U(0,\theta)$. Here we know $\theta$, in practice it would be an unknown parameter!!!
def sampleU(n, theta):
"""
Returns a sample of num independent values from U(0,theta).
"""
return stats.uniform.rvs(loc=0,scale=theta,size=n)
We present few variants of estimate.
The first one is obtained by the moment method, the second by the maximum likelihood.
(Can you find a better one? The best possible solution is not yet discovered.)
def my_estimate(x):
return 2*x.mean()
def my_estimate2(x):
return x.max()
def my_estimate3(x):
n = len(x)
return x.max()*(n+1)/n
Now: how to test which estimate is better?
We compute the bias and the mean squared error.
Observe the results and think about what it means.
Note the difference between two numbers: $n$ is the size of a sample from which the statistian decides about his/her estimate $\hat\theta$. OTOH, num is the number of repetitions we try this thing, to verify how good are the estimates.
That is, we use num to estimate the true expectation, while $n$ what allow us to compute, e.g., the sample mean.
def test_my_method(theta, sample_method, estimate_method, n=5, num=10**4, image=True, eps=0.1):
if theta is None: theta = 130
L = []
D = []
for _ in range(num):
data = sample_method(n,theta) # one experiment -- perform sample_size measurements, predict what it implies
D.append(data)
L.append(estimate_method(data))
estimates = np.array(L) # theta se striskou
MSE = np.mean((estimates-theta)**2)
bias = estimates.mean()-theta
p = np.mean(abs(estimates-theta)>eps)
print(f"Mean: {estimates.mean():.3f}")
print(f"Bias: {bias:.3f}")
print(f"Mean squared error: {MSE:.3f}")
print(f"Var+Bias^2: {estimates.var()+bias**2:.3f}, ")
print(f"P(error > {eps}) = {p}")
if image:
sns.set_style('whitegrid')
sns.histplot(estimates)
Testing one of the methods.
We can see that $MSE = Var+Bias^2$
test_my_method(130, sampleU, my_estimate, n=100, num=10000, eps=5, image=True)
Mean: 129.975 Bias: -0.025 Mean squared error: 54.634 Var+Bias^2: 54.634, P(error > 5) = 0.5014
Here we can see that with increasing $n$ probability of error $>10$ decreases (rather fast). (The method is consistent.)
Also the MSE decreases.
OTOH, bias is small already for $n=5$ -- the method is unbiased.
test_my_method(130, sampleU, my_estimate, n=5, num=10**4, image=False, eps=10); print()
test_my_method(130, sampleU, my_estimate, n=50, num=10**4, image=False, eps=10); print()
test_my_method(130, sampleU, my_estimate, n=500, num=10**4, image=False, eps=10)
Mean: 130.374 Bias: 0.374 Mean squared error: 1096.876 Var+Bias^2: 1096.876, P(error > 10) = 0.7685 Mean: 130.030 Bias: 0.030 Mean squared error: 111.129 Var+Bias^2: 111.129, P(error > 10) = 0.3485 Mean: 130.010 Bias: 0.010 Mean squared error: 10.929 Var+Bias^2: 10.929, P(error > 10) = 0.0013
This method is not unbiased.
In fact, the histogram of $\hat\theta$ shows, that we always undershoot the true value of $\theta$.
test_my_method(130, sampleU, my_estimate2, sample_size=100, num=10000, eps=5)
Mean: 128.721 Bias: -1.279 Mean squared error: 3.296 Var+Bias^2: 3.296, P(error > 5) = 0.0204
Third method. Note that this is apparently unbiased, also we sometimes overshoot, sometimes undershoot.
The MSE is much smaller than in the other methods.
test_my_method(130, sampleU, my_estimate3, sample_size=100, num=10000, eps=5)
Mean: 130.011 Bias: 0.011 Mean squared error: 1.673 Var+Bias^2: 1.673, P(error > 5) = 0.0076
def sampleGeom(n, secretp):
"""
Returns a sample of num independent values from Geom(secretp).
"""
return stats.geom.rvs(secretp,size=n)
def sampleExp(n, secretp):
"""
Returns a sample of num independent values from Exp(secretp).
"""
return stats.expon.rvs(scale=secretp,size=n)
a = sampleExp(10000,3)
a.mean()
2.975340338068114
def my_estimate4(x):
return 1/x.mean()
def my_estimate5(x):
return x.mean()
test_my_method(1/6, sampleGeom, my_estimate4, n=50, num=10000, eps=0.05)
Mean: 0.170 Bias: 0.003 Mean squared error: 0.000 Var+Bias^2: 0.000, P(error > 0.05) = 0.0279
test_my_method(1/6, sampleGeom, my_estimate4, n=50, num=10000, eps=0.05)
Mean: 0.169 Bias: 0.003 Mean squared error: 0.000 Var+Bias^2: 0.000, P(error > 0.05) = 0.0277
test_my_method(3, sampleExp, my_estimate5, n=5, num=10**5, eps=0.5)
Mean: 3.000 Bias: 0.000 Mean squared error: 1.795 Var+Bias^2: 1.795, P(error > 0.5) = 0.71342
