Untitled

import math
from itertools import product
from collections import defaultdict
from scipy.stats import norm, beta

def powerRSS(m, n, Delta, UCL):
    # محاسبه‌ی (F_theta = pnorm(Delta)) در R
    # در پایتون با norm.cdf از scipy.stats:
    F_theta = norm.cdf(Delta)
    
    # k = n/m در R. در صورتی‌که n بخش‌پذیر بر m باشد، k عدد صحیح خواهد بود
    # برای اطمینان از به‌دست آمدن عدد صحیح (در صورت درست بودن ورودی):
    k = n // m
    
    # معادل iterpc(9, m, replace=TRUE, ordered=TRUE) در R، 
    # تمام ترکیب‌های (در واقع پرموته با تکرار) طول m از اعداد {0,...,8}:
    # در R: R = getall(I) - 1 => یعنی همه‌ی عناصر را یک واحد کم می‌کند.
    # پس در پایتون باید ابتدا همه‌ی ترکیت‌ها را از 0..8 بگیریم و بعد 1 واحد کم کنیم.
    # اما طبق کد اصلی، در نهایت اعداد در R از {0..8} تبدیل می‌شوند به {-1..7}.
    all_sequences = list(product(range(9), repeat=m))
    R = []
    for seq in all_sequences:
        # هر عضو را یک واحد کم می‌کنیم
        R.append([x - 1 for x in seq])  # حالا هر x در بازه -1..7
    
    # بردارهای احتمال بتا برای pbeta(0.5, 1:m, m-(1:m)+1)
    # در R: pbeta(0.5, shape1=j, shape2=m-j+1) برای j از 1 تا m
    p_half = [
        beta.cdf(0.5, j, m - j + 1)
        for j in range(1, m+1)
    ]
    
    # بردارهای احتمال بتا برای pbeta(F_theta, 1:m, m-(1:m)+1)
    p_delta = [
        beta.cdf(F_theta, j, m - j + 1)
        for j in range(1, m+1)
    ]
    
    # تعریف تابع comb0 برای حالتی که ممکن است r < 0 یا r > k باشد:
    # در R، choose(k, x) در این شرایط مقدار 0 برمی‌گرداند.
    def comb0(n, r):
        if r < 0 or r > n:
            return 0
        return math.comb(n, r)
    
    # تابعی برای محاسبه‌ی حاصل‌ضرب 
    #    \prod_j choose(k, r_j) * [p[j]^(k - r_j)] * [(1 - p[j])^(r_j)]
    # که در حلقه‌های مختلف استفاده می‌کنیم.
    def product_term(r_row, p):
        prod_val = 1.0
        for j in range(m):
            x = r_row[j]
            cval = comb0(k, x)
            if cval == 0:
                return 0
            prob = p[j]
            prod_val *= cval * (prob ** (k - x)) * ((1 - prob) ** x)
            # اگر در میانه‌ی راه صفر شد، همان صفر را برمی‌گردانیم
            if prod_val == 0:
                return 0
        return prod_val
    
    # طبق فرمول در کد R:
    # U = (UCL + n) / 2
    # اگر U > n باشد، SUM = 0 برمی‌گردانیم و SUM1 = 0
    SUM = 0
    U = (UCL + n) / 2.0
    if U > n:
        return [SUM, 0]
    
    c_ = math.ceil(U)
    
    # برای تسریع عملیات، ابتدا فهرستی از ردیف‌های R را بر اساس مجموع عناصرشان گروه‌بندی می‌کنیم
    sum_map = defaultdict(list)
    for row in R:
        row_sum = sum(row)
        sum_map[row_sum].append(row)
    
    # حلقه‌ی اول: for(w in c:n) و جمع زدن روی p_half
    for w in range(c_, n+1):
        if w in sum_map:
            for r_row in sum_map[w]:
                SUM += product_term(r_row, p_half)
    
    # حلقه‌ی دوم: for(y in 0:(c-1)) و جمع زدن روی p_delta
    SUM1 = 0
    for y in range(0, c_):
        if y in sum_map:
            for r_row in sum_map[y]:
                SUM1 += product_term(r_row, p_delta)
    
    return [SUM, SUM1]

# مثال تست شبیه به R:
res = powerRSS(m=2, n=8, Delta=0, UCL=6)
print("powerRSS(2,8,0,6) in Python => ", res)