# 3.5 编程练习-最大似然估计
### 最大似然估计
这一节我们演示如何利用最大似然估计(maximum likelihood estimation)的方法来估计一组数据背后的高斯分布
首先我们生成数据：
{% tabs %}
{% tab title="Python" %}
```python
from scipy.stats import norm
import matplotlib.pyplot as plt
sigma = 3 # 高斯分布的均值
u = 5 # 均值
normobj = norm(loc=u, scale=sigma) # 定义一个高斯分布
data = normobj.rvs(size=1000) # 从高斯分布中采样1000个数据
```
{% endtab %}
{% endtabs %}
随后对数据数据进行简单可视化：
{% tabs %}
{% tab title="Python" %}
```python
plt.hist(data)
plt.ylabel('count')
plt.xlabel('data value')
```
{% endtab %}
{% endtabs %}
> 
现在我们假定分布的均值未知，利用数据来估计均值的值。
首先写出负对数似然函数(negative log-likelihood function)：
{% tabs %}
{% tab title="Python" %}
```python
def negloglikeli(params):
# data 和sigma是已知的，params(即均值)是需要传入的
likeli = norm.pdf(data, loc=params, scale=sigma)
loglikeli = np.log(likeli).sum()
# 对总的log likeli值取负，最大化loglikeli就等于最小化负loglikeli
nll = -loglikeli
return nll
```
{% endtab %}
{% endtabs %}
我们将不同的均值传入定义的函数，看其得出的负对数似然值，可以看到当mu=5时，负对数似然函数取到最小值。
{% tabs %}
{% tab title="Python" %}
# 输入不同的mu值
mu = np.arange(0,10)
nll = []
# 储存不同的mu值对应的负对数似然函数值
for i in mu:
nll.append(negloglikeli(i))
#可视化
plt.plot(nll)
plt.ylabel('negative log-likeli')
plt.xlabel('mu')
{% endtab %}
{% endtabs %}
> 
或者我们也可以直接使用scipy 中的minimize函数，来直接找到使得负对数似然函数取到最小的mu值。
{% tabs %}
{% tab title="Python" %}
```python
from scipy.optimize import minimize
res = minimize(fun=negloglikeli, x0=(3,))
print('The estimated mean of the distribution is', res.x[0])
```
{% endtab %}
{% endtabs %}
> The estimated mean of the distribution is 4.9821100894980646
## 用最大似然估计解决线性回归问题
假定$$y =0.5x+3$$，我们生成$$x$$与对应的$$y$$并进行可视化。
{% tabs %}
{% tab title="Python" %}
```python
import numpy as np
a = 0.5 # 系数
b = 3 # 截距
noiseVar = 1 # noise的标准差
x = np.arange(1, 100, 0.1) #生成x值
y = a * x + b + noiseVar* np.random.randn(x.size) #代入公式生成对应的y
#可视化
plt.plot(x, y, 'o', color='r')
```
{% endtab %}
{% endtabs %}
> 
我们先复习一下上一章讲的最小二乘法求解线性回归问题：
{% tabs %}
{% tab title="Python" %}
```python
# 构建loss函数
def loss(params):
# params是个(2,)的数组, 第一个是斜率，第二个是截距
a = params[0]
b = params[1]
y_pred = a * x + b
return ((y-y_pred)**2).sum() # 计算squared error作为loss函数
from scipy.optimize import minimize
# minimize会自动寻找让loss取得最小值的两个参数
res = minimize(fun=loss, x0=(1, 1))
print('Linear coefficient is', res.x[0])
print('Intercept is', res.x[1])
```
{% endtab %}
{% endtabs %}
> Linear coefficient is 0.5008892586858943
>
> Intercept is 2.943209971344448
下面用最大似然估计来求解线性回归问题：
{% tabs %}
{% tab title="Python" %}
```python
from scipy.optimize import minimize
from scipy.stats import norm
eps = np.finfo(float).eps
rv = norm() # 创建一个正态分布的对象
def lossfun2(params):
# params是个(2,)的数组, 第一个是斜率，第二个是截距
a = params[0]
b = params[1]
y_pred = a * x + b
pp = norm.pdf(y, loc=y_pred, scale=noiseVar)
pp = 0.999*pp + eps # 防止概率过小结果变成inf
negloglikeli = -np.log(pp).sum()
return negloglikeli # 计算negative log likelihood作为loss函数
res = minimize(fun=lossfun2, x0=(1, 1))
print('Linear coefficient is', res.x[0])
print('Intercept is', res.x[1])
```
{% endtab %}
{% endtabs %}
> Linear coefficient is 0.5008892577373644
>
> Intercept is 2.943210034713618
可以看到，两者的结果非常接近，实际上两者在数学上是完全等价的。
---
# Agent Instructions: Querying This Documentation
If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.
Perform an HTTP GET request on the current page URL with the `ask` query parameter:
```
GET https://ruyuanzhang.gitbook.io/compmodcogpsy/di-san-zhang-gailtui-duan-he-bei-ye-si-li-lun/3.5-bian-cheng-lian-xi-zui-da-si-ran-gu-ji.md?ask=
```
The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.
Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.