Untitled
unknown
python
2 months ago
7.7 kB
3
Indexable
import math
from math import comb, gamma
from itertools import product
import numpy as np
from scipy.stats import norm, beta as beta_dist
from scipy.special import beta as BetaFunc # Beta(x,y) = Γ(x)Γ(y)/Γ(x+y)
# --- ثابت‌های مسأله (معادل globalها در کد R) ---
Z0 = 0.25
Z1 = 0.25
Z2 = 0.75
a = 20
b = 4.22
D0 = 50
D1 = 950
Y = 500
W = 1100
delta = 0 # برای pnorm(delta)
k = 3
lambda_ = 0.002 # (lambda را نمی‌توان نام متغیر در پایتون گذاشت)
# تابع A بر اساس کد R:
# ------------------------------------------------
# A(x) در Banerjee سری‌ای را تا زمانی که اختلاف گام‌ها کمتر از 1e-6 شود جمع می‌کند.
def A(xx):
aw = [0.0] # متناظر با aw<-0 در R ولی با لیست برای نگهداری مقادیر تجمعی
vw = 0 # شمارنده نمای
s = 0.001 # مقدار اولیه برای شروع حلقه
# تا زمانی که مقدار گام اخیر از 1e-6 بیشتر باشد ادامه می‌دهیم
while s >= 1e-6:
# ((1+vw)^(1/k)) * (xx^vw)
cw = ((1 + vw)**(1 / k)) * (xx**vw)
new_val = aw[-1] + cw
s = new_val - aw[-1] # اختلاف گام
aw.append(new_val)
vw += 1
return aw[-1]
# تابع اصلی co که معادل co<-function(x){...} در R است
# x ورودی‌ای به شکل [m,q,h,L]
# ------------------------------------------------
def co(x):
# استخراج پارامترها از بردار x
m = int(math.floor(x[0]))
q = int(math.floor(x[1]))
h = x[2]
L = x[3]
n = m * q # کل نمونه در هر سیکل
alpha = 0.0
beta_ = 0.0 # نام آن را beta_ می‌گذاریم تا با تابع beta تداخلی نداشته باشد
F_theta = norm.cdf(delta) # pnorm(delta)
# Fi(m,i) و delta0(m) در کد R:
# -----------------------------------------------------------------
# Fi(m,i) = m * choose(m-1,i-1)* beta(i, m-i+1)* pbeta(0.5, i, m-i+1)
# delta0(m)=1 - (4/m)* sum((Fi(m,1..m)-0.5)^2)
# اینها برای محاسبه‌ی ucl استفاده می‌شوند.
def Fi(m_, i_):
# choose(m-1, i-1)
ch = comb(m_ - 1, i_ - 1)
# beta(i, m-i+1) => تابع BetaFunc(i, m-i+1)
bt = BetaFunc(i_, m_ - i_ + 1)
# pbeta(0.5, i, m-i+1)
cdf_05 = beta_dist.cdf(0.5, i_, m_ - i_ + 1)
return m_ * ch * bt * cdf_05
def delta0(m_):
vals = []
for i_ in range(1, m_ + 1):
vals.append(Fi(m_, i_) - 0.5)
return 1 - (4 / m_) * sum([v**2 for v in vals])
# ساختن تمام ترکیبات R از طریق itertools.product(range(9), repeat=m)
# معادل با iterpc(9,m, replace=TRUE, ordered=TRUE) و سپس getall(I) - 1.
# در R: getall(I) اعداد 1..9 می‌داد که منهای 1 => 0..8
# در اینجا مستقیماً range(9) => 0..8
R_list = list(product(range(9), repeat=m))
# محاسبه ucl
d0 = delta0(m) # همان delta0(m)
ucl = int(math.floor((n / 2 + (L * np.sqrt(n * d0 / 4)) + n) / 2))
# در کد R اگر ucl>n باشد، EA=100000 مقداردهی می‌شود اما به شکل return انجام نمی‌دهد.
# برای سادگی در پایتون مستقیم برمی‌گردیم:
if ucl > n:
return 1e5 # هزینه بالا برای طرح نامعتبر
# محاسبه alpha با حلقه روی c=0..ucl
# "alpha" = 1 - Sum_{sum(r)=c}(...) => در انتها alpha=1-alpha
cdf_half = [beta_dist.cdf(0.5, j, m - j + 1) for j in range(1, m + 1)]
alpha_sum = 0.0
for c_ in range(ucl + 1):
# فیلتر آنهایی که جمع المانهایشان c_ است
subset_c = [r for r in R_list if sum(r) == c_]
for r_ in subset_c:
# محاسبه حاصلضرب بر اساس j از 1..m
# اما در پایتون ایندکس r_ از 0..m-1 است
pr = 1.0
for j in range(m):
# choose(q, r_[j])
ch = comb(q, r_[j])
# pbeta(0.5, j+1, m-j) = cdf_half[j]
p_half = cdf_half[j]
pr *= ch * (p_half**(q - r_[j])) * ((1 - p_half)**(r_[j]))
alpha_sum += pr
alpha = 1.0 - alpha_sum
# محاسبه beta (beta_) با حلقه مشابه، اما با F_theta
cdf_Ftheta = [beta_dist.cdf(F_theta, j, m - j + 1) for j in range(1, m + 1)]
for y_ in range(ucl + 1):
subset_y = [r for r in R_list if sum(r) == y_]
for r_ in subset_y:
pr = 1.0
for j in range(m):
ch = comb(q, r_[j])
p_ftheta = cdf_Ftheta[j]
pr *= ch * (p_ftheta**(q - r_[j])) * ((1 - p_ftheta)**(r_[j]))
beta_ += pr
# حالا p در Banerjee:
# p = 1 - exp(-lambda * h^k)
p = 1.0 - math.exp(-lambda_ * (h**k))
# محاسبه A(1-p) و A(beta_)
Ap = A(1 - p)
Abeta_ = A(beta_)
# محاسبه ET بر اساس فرمول‌های داخل کد
# ------------------------------------------------
# ET = Z1 + Z2
# + alpha * Z0 * ((1-p)/p)
# + (h*p*Ap)
# + [ beta * h * p * ( (p*Ap) - (1-beta)*Abeta ) ] / (1 - p - beta)
try:
ET = (Z1 + Z2
+ alpha * Z0 * ((1 - p) / p)
+ (h * p * Ap)
+ (beta_ * h * p * ((p * Ap) - (1 - beta_) * Abeta_)) / (1 - p - beta_))
except ZeroDivisionError:
# اگر 1 - p - beta_ صفر یا خیلی کوچک باشد
return 1e5
# محاسبه EC
# ------------------------------------------------
# EC = (a + b*n)* [ (beta/(1-beta)) + 1/p ]
# + alpha*Y*((1-p)/p)
# + (D0 - D1)* ( ((1/lambda)^(1/k)) * gamma(1 + 1/k) - h*(1-p)*p*Ap )
# + D0*h*p*(1-p)*Ap
# + [ (beta*h*D1*p*( p*Ap - (1-beta)*Abeta )) / (1-p-beta) ]
# + (D1*h*(p^2)*Ap)
# + W
try:
EC = ((a + b * n) * ((beta_ / (1 - beta_)) + 1.0 / p)
+ alpha * Y * ((1 - p) / p)
+ (D0 - D1) * (((1.0 / lambda_)**(1.0 / k)) * gamma(1.0 + 1.0 / k)
- h * (1 - p) * p * Ap)
+ (D0 * h * p * (1 - p) * Ap)
+ ((beta_ * h * D1 * p * ((p * Ap) - (1 - beta_) * Abeta_))
/ (1 - p - beta_))
+ (D1 * h * (p**2) * Ap)
+ W)
except ZeroDivisionError:
return 1e5
EA = EC / ET
return EA
# --- نمونه فراخوانی مشابه با کد R: ---
if __name__ == "__main__":
# چند نمونه از ورودی‌ها که در انتهای کد R هم آمده بود:
print("co([2,1,4.17,1]) =", co([2,1,4.17,1]))
print("co([3,1,3.25,1]) =", co([3,1,3.25,1]))
print("co([4,1,3.2,1]) =", co([4,1,3.2,1]))
print("co([2,2,3.13,1]) =", co([2,2,3.13,1]))
print("co([2,3,30,1]) =", co([2,3,30,1]))
# بهینه‌سازی با L-BFGS-B در محدوده خاص:
from scipy.optimize import minimize
# حد پایین و بالای هر پارامتر: [m, q, h, L]
bounds = [(2,4), # m بین 2 و 4
(1,6), # q بین 1 و 6
(1,5), # h بین 1 و 5
(1,2)] # L بین 1 و 2
x0 = [2,1,1,1] # شروع
res = minimize(co, x0, method="L-BFGS-B", bounds=bounds)
print("\nResult of L-BFGS-B optimization:")
print("x* =", res.x)
print("objective =", res.fun)
print("status:", res.message)
Editor is loading...
Leave a Comment