#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)]
.