<jupyter_start><jupyter_text>M市場のデータのみ使って同じことができるかを確かめる結論：<jupyter_code>import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt
% matplotlib inline
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from numpy.linalg import inv
from numpy.random import *
data = pd.read_csv("robust_norm_censored_4.csv")
# M市場についてしか観察できていいない時を考える。
# 実現した参入社数はnumからの経験分布からドロー
M = 10000
T = 1000
seed(1234)
m = data.shape[0]
draw = randint(1,T+1,m) 
data["draw"] = draw
data["realized_entry"] = data.apply(lambda x: 0 if(x["draw"] <= x["num0"]) else 1 if(x["draw"] <= x["num1"]) else 2, axis = 1)<jupyter_output><empty_output><jupyter_text>polynomialは入れない方が良い推定値が得られてるので省略data["popdist1"] = data["Pop"] * data["Dist1"]data["popdist2"]= data["Pop"] * data["Dist2"]data["dist1dist2"]= data["Dist2"] * data["Dist1"]data["pop_sq"] = data["Pop"]**2data["dist1_sq"] = data["Dist1"]**2data["dist2_sq"] = data["Dist2"]**2data["pop_tr"] = data["Pop"]**3data["dist1_tr"] = data["Dist1"]**3data["dist2_tr"] = data["Dist2"]**3data["pop_sqdist1"] = data["Pop"]**2 * data["Dist1"]data["pop_sqdist2"] = data["Pop"]**2 * data["Dist2"]data["dist1_sqpop"] = data["Pop"] * data["Dist1"]**2data["dist1_sqdist2"] = data["Dist2"] * data["Dist1"]**2data["dist2_sqpop"] = data["Pop"] * data["Dist2"]**2data["dist2_sqdist1"] = data["Dist1"] * data["Dist2"]**2 realized entryを使って予測する<jupyter_code>from sklearn import linear_model
logreg = linear_model.LogisticRegression(C=1e5,
                         penalty='l1')
data.head()
Y = data.realized_entry
X = data.drop(["Unnamed: 0", "diff", "num0", "num1", "num2", "draw", "realized_entry"],axis = 1)
logreg.fit(X, Y)
data["diff_est"] = logreg.predict_proba(X)[:, 0] - logreg.predict_proba(X)[:, 2]
d = data["diff"].values
d_est = data["diff_est"].values
plt.scatter(d, d_est)<jupyter_output><empty_output><jupyter_text>d_estを使ってrobustを出してみる<jupyter_code>import numdifftools as nd
def robust(x):
    xb1 = x[0] * pop + x[1] * dist1
    xb2 = x[0] * pop + x[1] * dist2
    r = (diff + s.cdf(xb1 + x[2])*s.cdf(xb2 + x[2]) - s.cdf(-xb1)*s.cdf(-xb2))**2
    return -sum(r)

def NR(initial, fun, tol, maxit):
    theta = initial
    for i in range(0,maxit):
        G = nd.Gradient(fun)(theta)
        H = nd.Hessian(fun)(theta)
        del_theta = inv(H).dot(G)
        theta = theta - del_theta
        if np.linalg.norm(del_theta) < tol:
            print("NR終わり")
            break
    return theta

s = st.norm()
pop = data[["Pop"]].values[:, 0]
dist1 = data[["Dist1"]].values[:, 0]
dist2 = data[["Dist2"]].values[:, 0]
diff = data[["diff_est"]].values[:, 0]
result_rob = NR([0.8,-0.5,-0.25], robust, 0.0000005, 100000)
# そこそこいい。
# あとはサンプルサイズの問題な気がする。
result_rob<jupyter_output><empty_output><jupyter_text>ということでサンプルサイズを増やしてやってみる。summary_noteでやってみた。しかしこれはバイアスがかかってるな  他の予測は？deepはダメでした。<jupyter_code>from sklearn.neural_network import MLPClassifier
clf = MLPClassifier(solver='lbfgs', alpha=1e-5, random_state=1)
clf.fit(X, Y) 
data["diff_est2"] = clf.predict_proba(X)[:, 0] - clf.predict_proba(X)[:, 2]
d = data["diff"].values
d_est = data["diff_est2"].values
plt.scatter(d, d_est)
diff = data[["diff_est2"]].values[:, 0]
result_rob = NR([0.8,-0.5,-0.25], robust, 0.0000005, 100000)
result_rob<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19530773918342473, 'fraction_numerical': 0.04692260816575259, 'mean_word_length': 5.783057851239669, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '323060', 'license': 'MIT', 'n_tokens_mistral': 1538, 'n_tokens_neox': 1464, 'n_words': 222, 'path': 'Mmarketcase.ipynb', 'repo_name': 'keiikegami/nashtest', 'size': '3282'}