#!/usr/bin/env python
# coding: utf-8

# $
# \begin{equation*}
# \begin{aligned}
# & \underset{x}{\text{maximize}}
# & & f(x) \\
# & \text{subject to}
# & & g_i(x) \leq b_i, \; i = 1, \ldots, m.\\
# \end{aligned}
# \end{equation*}
# $

# In[1]:


import scipy.stats


# In[80]:


X = np.linspace(0, 1, 200)
plt.figure(figsize=(6, 3))
plt.plot(X, scipy.stats.beta.pdf(X, 150, 150), label='A: 試行300回')
plt.plot(X, scipy.stats.beta.pdf(X, 6, 4), label='B: 試行10回')
plt.legend(loc='upper left')
plt.xlabel('p')
plt.ylabel('尤度')
#plt.title('パラメータpの尤度\n Beta(成功回数, 失敗回数)')


# In[89]:


X = np.linspace(0, 1, 200)
plt.figure(figsize=(6, 3))
plt.plot(X, scipy.stats.beta.pdf(X, 300, 300), label='A: 試行600回')
plt.plot(X, scipy.stats.beta.pdf(X, 700, 400), label='B: 試行1100回')
plt.legend(loc='upper left')
plt.xlabel('p')
plt.ylabel('尤度')
#plt.title('パラメータpの尤度\n Beta(成功回数, 失敗回数)')


# $
# y = w_0 + w_1x_1 + w_2x_2 \cdots + \epsilon
# $

# In[4]:


from sklearn.datasets import load_boston
boston = load_boston()


# In[5]:


print(boston.DESCR)


# In[6]:


boston.feature_names


# In[14]:


import pandas as pd
X = pd.DataFrame(boston.data, columns=boston.feature_names)
X.head()


# In[15]:


y = pd.Series(boston.target)


# In[70]:


import statsmodels.api as sm
model = sm.OLS(y, X[['CRIM', 'ZN', 'CHAS', 'NOX', 'RM', 'DIS', 'RAD', 'TAX']])
result = model.fit()
print(result.summary2())


# In[23]:


import pulp


# $\displaystyle
# \begin{equation*}
# \begin{aligned}
# & \underset{x,y,z}{\text{maximize}}
# & 1800x + 2500y + 3000z \\
# & \text{subject to}
# & 6x + 10y + 15z \leq 12000 \\
# & & x \leq 1000; y \leq 900;z \leq 500 
# \end{aligned}
# \end{equation*}
# $

# In[24]:


# Variable
products = [pulp.LpVariable(f'{v}', cat=pulp.LpInteger, lowBound=0) for v in ['X', 'Y', 'Z']]
# Target
problem = pulp.LpProblem(sense=pulp.LpMaximize)
problem += pulp.lpDot(products, [1800, 2500, 3000])
# Constraints
problem += pulp.lpDot(products, [6, 10, 15]) <= 200 * 60
problem += products[0] <= 1000
problem += products[1] <= 900
problem += products[2] <= 500

status = problem.solve()
print(pulp.LpStatus[status])
print([f"{product.name}: {product.value()}" for product in products])


# In[25]:


1000 * 1800 + 600 * 2500


# In[ ]: