Как элегантно построить моделирование методом Монте-Карло для уровня достоверности распределения Бернулли

#python #pandas #scipy #statistics #montecarlo

Вопрос:

Задача домашнего задания состоит в том, чтобы, используя список n, p для каждой комбинации n и p, использовать моделирование методом Монте-Карло для оценки фактического уровня достоверности, которого достигает каждый метод генерации интервалов, и проверить, достигается ли номинальный уровень достоверности каждым методом по сетке значений для n и p

Поэтому для каждой пары n, p я должен сгенерировать N копий, чтобы выполнить так называемое моделирование по методу Монте-Карло.

Ниже приведена моя попытка,

 import numpy as np import pandas as pd import matplotlib.pyplot as plt import matplotlib as mpl from scipy import stats as st def ci_prop(  x,  level=0.95,  method="Normal" ):  """  Input:  x : A 1-dimensional NumPy array or compatible sequence type (list, tuple).   confidence level : float, optional.   method: str, optional   Output  -------  lwr  upr  width   """  # check input type  try:  x = np.asarray(x) # or np.array() as instructed.  except TypeError:  print("Could not convert x to type ndarray.")   # check that x is bool or 0/1  if x.dtype is np.dtype('bool'):  pass  elif not np.logical_or(x == 0, x == 1).all():  raise TypeError("x should be dtype('bool') or all 0's and 1's.")   # check method  assert method in ["Normal", "CP", "Jeffrey", "AC"]   # determine the length  n = x.size   # compute estimate  if method == 'AC':  z = st.norm.ppf(1 - (1 - level) / 2)  n = (n   z ** 2)  est = (np.sum(x)   z ** 2 / 2) / n  else:  est = np.mean(x)   # warn for small sample size with "Normal" method  if method == 'Normal' and (n * min(est, 1 - est)) lt; 12:  warn(Warning(  "Normal approximation may be incorrect for n * min(p, 1-p) lt; 12."  ))   # compute bounds for Normal and AC methods  if method in ['Normal', 'AC']:  se = np.sqrt(est * (1 - est) / n)  z = st.norm.ppf(1 - (1 - level) / 2)  lwr, upr = est - z * se, est   z * se  width = upr - lwr   # compute bounds for CP method  if method == 'CP':  alpha = 1 - level  s = np.sum(x)  lwr = beta.ppf(alpha / 2, s, n - s   1)  upr = beta.ppf(1 - alpha / 2, s   1, n - s)  width = upr - lwr    # compute bounds for Jeffrey method  if method == 'Jeffrey':  alpha = 1 - level  s = np.sum(x)  lwr = beta.ppf(alpha / 2, s   0.5, n - s   0.5)  upr = beta.ppf(1 - alpha / 2, s   0.5, n - s   0.5)  width = upr - lwr    # prepare return values    return lwr, upr, width  n = np.linspace(40, 60, 6).astype("int") p = np.linspace(0.7, 0.9, 6) x = len(n) y = len(p) N_rep = 200 # for example, I set N = 200  Data_Set = [] for i in range(x):  for j in range(y):  A = np.random.binomial(1, p[j], n[i])  res = ci_prop(A, 0.8, method = "AC")  Data_Set.append(res)  Data_Set  

Я хочу получить результат, для каждого метода («AC», «Нормальный» и т. Д.) Существует 1-d массив, Такой, что каждый элемент этого массива соответствует пропорции, в которой 200(например, пусть N = 200) реплицировали доверительный интервал, включая истинное среднее значение (p) для соответствующей пары (n, p)

Я попытался использовать циклы, но потерпел неудачу, так как не знаю, как позволить ему генерировать 200 раз. Кроме того, выполнение вложенных циклов for кажется слишком неудобным и неэффективным. Мне интересно, есть ли более эффективный и элегантный способ добиться этого? Любая помощь или подсказка приветствуются.

Комментарии:

1. Как правило, легче получить помощь, если вы зададите один конкретный вопрос с минимальным фрагментом кода, сфокусированным на том, с чем вам нужна помощь, а не на всем задании. В любом случае, для циклов часто помогает упаковать цикл внутри функции, если вы понимаете, что я имею в виду. Или понимание списка, например ci_prop(sample(i,j), 0.8, method="AC") for i in range(X) for j in range(y)] .