ホーム » コードあり (ページ 4)

「コードあり」カテゴリーアーカイブ

2022 Q3(3)

投稿日:2024年1月29日 　最終更新日:2024年8月5日

ポアソン分布のパラメータλがガンマ分布に従うと負の二項分布になる。

コード

数式を使った計算

# 2022 Q3(3)  2024.8.5

import sympy as sp

# 定義
k = sp.symbols('k', integer=True)
alpha, beta = sp.symbols('alpha beta', positive=True)
lambda_var = sp.symbols('lambda', positive=True)

# ポアソン分布の確率質量関数
poisson_pmf = (lambda_var**k * sp.exp(-lambda_var)) / sp.factorial(k)

# ガンマ分布の確率密度関数
gamma_pdf = (beta**alpha / sp.gamma(alpha)) * lambda_var**(alpha - 1) * sp.exp(-beta * lambda_var)

# 周辺分布の計算
marginal_distribution = sp.integrate(poisson_pmf * gamma_pdf, (lambda_var, 0, sp.oo)).simplify()

marginal_distribution

次式と同じなので、負の二項分布となる。

シミュレーションによる計算

# 2022 Q3(3)  2024.8.5

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma, nbinom

# パラメータの設定
alpha = 3.0
beta = 2.0
sample_size = 10000

# ガンマ分布から λ を生成
lambda_samples = np.random.gamma(alpha, 1/beta, sample_size)

# 生成された λ を使ってポアソン分布から X を生成
poisson_samples = [np.random.poisson(lam) for lam in lambda_samples]

# 負の二項分布のパラメータを設定
r = alpha
p = beta / (beta + 1)

# 負の二項分布から理論的な確率質量関数を計算
x = np.arange(0, max(poisson_samples) + 1)
nbinom_pmf = nbinom.pmf(x, r, p)

# ヒストグラムの描画
plt.hist(poisson_samples, bins=np.arange(0, max(poisson_samples) + 1) - 0.5, density=True, alpha=0.75, color='blue', edgecolor='black', rwidth=0.8)

# 理論的な負の二項分布の確率質量関数をプロット
plt.plot(x, nbinom_pmf, 'r', linestyle='-', label='理論的な負の二項分布')

# グラフのタイトルとラベル
plt.title('ポアソン-ガンマ混合分布のシミュレーション結果')
plt.xlabel('値')
plt.ylabel('確率密度')
plt.legend()

# グラフの表示
plt.grid(True)
plt.show()

2022 Q3(2)

投稿日:2024年1月28日 　最終更新日:2024年8月4日

ガンマ分布の期待値と分散を求めました。

コード

数式を使った計算

## 2022 Q3(2)  2024.8.4

import sympy as sp

# 定義
lambda_var = sp.symbols('lambda')
alpha, beta = sp.symbols('alpha beta', positive=True)

# ガンマ分布の確率密度関数
gamma_pdf = (beta**alpha / sp.gamma(alpha)) * lambda_var**(alpha-1) * sp.exp(-beta * lambda_var)

# 期待値 E[Λ] の計算
expected_value = sp.integrate(lambda_var * gamma_pdf, (lambda_var, 0, sp.oo)).simplify()

# E[Λ^2] の計算
expected_value_Lambda2 = sp.integrate(lambda_var**2 * gamma_pdf, (lambda_var, 0, sp.oo)).simplify()

# 分散 V[Λ] の計算
variance = (expected_value_Lambda2 - expected_value**2).simplify()

# 結果を辞書形式で表示
{
    "期待値 E[Λ]": expected_value,
    "分散 V[Λ]": variance
}
{'期待値 E[Λ]': alpha/beta, '分散 V[Λ]': alpha/beta**2}

シミュレーションによる計算

## 2022 Q3(2)  2024.8.4

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma

# パラメータ α と β の設定
alpha = 3.0
beta = 2.0
# サンプルサイズ
sample_size = 10000

# ガンマ分布に従う乱数を生成
samples = np.random.gamma(alpha, 1/beta, sample_size)

# 期待値と分散の計算
expected_value = np.mean(samples)
variance = np.var(samples)

# 結果の表示
print(f"期待値のシミュレーション結果: {expected_value}")
print(f"分散のシミュレーション結果: {variance}")

# ヒストグラムの描画
plt.hist(samples, bins=50, density=True, alpha=0.75, color='blue', edgecolor='black')

# 理論的なガンマ分布の確率密度関数をプロット
x = np.linspace(0, max(samples), 100)
gamma_pdf = gamma.pdf(x, alpha, scale=1/beta)
plt.plot(x, gamma_pdf, 'r', linestyle='-', label='理論的なガンマ分布')

# グラフのタイトルとラベル
plt.title('ガンマ分布のシミュレーション結果')
plt.xlabel('値')
plt.ylabel('確率密度')
plt.legend()

# グラフの表示
plt.grid(True)
plt.show()
期待値のシミュレーション結果: 1.4947527551017439
分散のシミュレーション結果: 0.7452822182213243
# パラメータ α と β の設定
alpha = 3.0
beta = 2.0

# 理論値の計算
theoretical_expected_value = alpha / beta
theoretical_variance = alpha / beta**2

# 結果の表示
print(f"理論値 - 期待値: {theoretical_expected_value}")
print(f"理論値 - 分散: {theoretical_variance}")
理論値 - 期待値: 1.5
理論値 - 分散: 0.75

2022 Q3(1)

投稿日:2024年1月27日 　最終更新日:2024年8月3日

ポアソン分布の期待値と分散の導出。

コード

数式を使った計算

# 2022 Q3(1)  2024.8.3

import sympy as sp

# 定義
k = sp.symbols('k')
lambda_param = sp.symbols('lambda')

# ポアソン分布の確率関数
poisson_pmf = (lambda_param**k * sp.exp(-lambda_param)) / sp.factorial(k)

# 期待値 E[X] の計算
expected_value = sp.summation(k * poisson_pmf, (k, 0, sp.oo)).simplify()

# E[X^2] の計算
expected_value_X2 = sp.summation(k**2 * poisson_pmf, (k, 0, sp.oo)).simplify()

# 分散 V[X] の計算
variance = (expected_value_X2 - expected_value**2).simplify()

# 結果を表示
{
    "期待値 E[X]": expected_value,
    "分散 V[X]": variance
}
{'期待値 E[X]': lambda, '分散 V[X]': lambda}

シミュレーションによる計算

import numpy as np

# パラメータ λ の設定
lambda_param = 5
# シミュレーションの回数
num_simulations = 10000
# 各シミュレーションで生成するサンプルの数
sample_size = 1000

# シミュレーション結果の保存用
expected_values = []
variances = []

for _ in range(num_simulations):
    # ポアソン分布に従う乱数を生成
    samples = np.random.poisson(lambda_param, sample_size)
    # 期待値を計算
    expected_values.append(np.mean(samples))
    # 分散を計算
    variances.append(np.var(samples))

# シミュレーション結果の平均を計算
average_expected_value = np.mean(expected_values)
average_variance = np.mean(variances)

print(f"期待値のシミュレーション結果: {average_expected_value}")
print(f"分散のシミュレーション結果: {average_variance}")

期待値のシミュレーション結果: 4.999602
分散のシミュレーション結果: 4.9992533574

プロット

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import factorial

# パラメータ λ の設定
lambda_param = 5
# サンプルサイズ
sample_size = 10000

# ポアソン分布に従う乱数を生成
samples = np.random.poisson(lambda_param, sample_size)

# ヒストグラムの描画
plt.hist(samples, bins=np.arange(0, max(samples) + 1) - 0.5, density=True, alpha=0.75, color='blue', edgecolor='black')

# 理論的なポアソン分布の確率質量関数をプロット
x = np.arange(0, max(samples) + 1)
poisson_pmf = (lambda_param**x * np.exp(-lambda_param)) / factorial(x)
plt.plot(x, poisson_pmf, 'r', marker='o', linestyle='-', label='理論的なポアソン分布')

# グラフのタイトルとラベル
plt.title('ポアソン分布のシミュレーション結果')
plt.xlabel('値')
plt.ylabel('確率')
plt.legend()

# グラフの表示
plt.grid(True)
plt.show()

2022 Q2(5)

投稿日:2024年1月26日 　最終更新日:2024年8月2日

2変数の条件付き確率を表す図形から相関係数を求める問題をやりました。

コード

二重積分による計算

# 2022 Q2(5)  2024.8.2

import numpy as np
from scipy.integrate import dblquad

# 確率密度関数
def pdf(u, v):
    return 1/3 if (u**2 - 2*u*v + v**2 <= 1) else 0

# 期待値 E[U] を計算
E_U, _ = dblquad(lambda v, u: u * pdf(u, v), -1, 1, lambda u: -1, lambda u: 1)
# 期待値 E[V] を計算
E_V, _ = dblquad(lambda v, u: v * pdf(u, v), -1, 1, lambda u: -1, lambda u: 1)

# 期待値 E[U^2] を計算
E_U2, _ = dblquad(lambda v, u: u**2 * pdf(u, v), -1, 1, lambda u: -1, lambda u: 1)
# 期待値 E[V^2] を計算
E_V2, _ = dblquad(lambda v, u: v**2 * pdf(u, v), -1, 1, lambda u: -1, lambda u: 1)

# 分散 Var(U) を計算
Var_U = E_U2 - E_U**2
# 分散 Var(V) を計算
Var_V = E_V2 - E_V**2

# 共分散 Cov(U, V) を計算
Cov_UV, _ = dblquad(lambda v, u: u*v * pdf(u, v), -1, 1, lambda u: -1, lambda u: 1)

# 相関係数 ρ(U, V) を計算
rho_UV = Cov_UV / np.sqrt(Var_U * Var_V)

E_U, E_V, Var_U, Var_V, rho_UV

(0.0, 0.0, 0.27777777931630593, 0.2777565638412289, 0.5000186020617616)

モンテカルロ法による計算

# 2022 Q2(5)  2024.8.2

import numpy as np

# モンテカルロ法の設定
num_points = 10000

# UとVの一様乱数を生成
U = np.random.uniform(-1, 1, num_points)
V = np.random.uniform(-1, 1, num_points)

# 条件 U^2 - 2UV + V^2 <= 1 を満たす点をフィルタリング
condition = U**2 - 2*U*V + V**2 <= 1
U_filtered = U[condition]
V_filtered = V[condition]

# フィルタリングされた点数
num_filtered_points = len(U_filtered)

# 期待値 E[U] と E[V] の計算
E_U = np.mean(U_filtered)
E_V = np.mean(V_filtered)

# 分散 Var(U) と Var(V) の計算
Var_U = np.var(U_filtered, ddof=1)
Var_V = np.var(V_filtered, ddof=1)

# 共分散 Cov(U, V) の計算
Cov_UV = np.cov(U_filtered, V_filtered, ddof=1)[0, 1]

# 相関係数 ρ(U, V) の計算
rho_UV = Cov_UV / np.sqrt(Var_U * Var_V)

E_U, E_V, Var_U, Var_V, rho_UV, num_filtered_points

(-0.0030089627151707833,
 -0.0028658643404494717,
 0.2780688116636205,
 0.2779098545694732,
 0.4835460989080292,
 7488)

アルゴリズム

モンテカルロ法

  1. 初期設定
    • シミュレーションで使用する点の数を設定する。
  2. 乱数の生成
    • 範囲 [−1,1] でランダムに U と V の値を生成する。
  3. 条件の適用
    • 生成した点の中で、条件 [ U^2 - 2UV + V^2 \leq 1 ] を満たす点を選び出す。
  4. 期待値の計算
    • 条件を満たした点の U と V の平均値を計算する。
  5. 分散の計算
    • 条件を満たした点の U と V の分散を計算する。
  6. 相関係数の計算
    • 条件を満たした点の U と V の共分散を計算し、それを用いて相関係数を求める。
  7. 結果の表示
    • 計算した期待値、分散、相関係数、および条件を満たした点の数を表示する。

プロット

2022 Q2(4)

投稿日:2024年1月25日 　最終更新日:2024年8月1日

不等式を図示して面積から確率を求める問題をやりました。

コード

計算

#2022 Q2(4) 2024.8.1

import numpy as np
import scipy.integrate as integrate

# 積分の被積分関数を定義
def integrand(v, u):
    return 1/4 if (u**2 - 2*u*v + v**2 <= 1) else 0

# 積分の範囲を定義する関数
def bounds_u():
    return -1, 1

def bounds_v(u):
    return -1, 1

# 確率 P(U^2 - 2UV + V^2 <= 1) を計算
probability, _ = integrate.nquad(integrand, [bounds_v, bounds_u])

probability

0.7500047152195971

プロット

#2022 Q2(4) 2024.8.1

import matplotlib.pyplot as plt
import numpy as np

# UとVの範囲を設定
u = np.linspace(-1, 1, 400)
v = np.linspace(-1, 1, 400)
U, V = np.meshgrid(u, v)

# 条件 U^2 - 2UV + V^2 <= 1 を満たす領域を計算
condition = (U**2 - 2*U*V + V**2 <= 1)

# グラフの描画
plt.figure(figsize=(5, 5))
plt.contourf(U, V, condition, levels=[0.5, 1], colors=['lightblue'], alpha=0.5)
plt.title(r'条件 $U^2 - 2UV + V^2 \leq 1$')
plt.xlabel('U')
plt.ylabel('V')
plt.axhline(0, color='gray', lw=0.5)
plt.axvline(0, color='gray', lw=0.5)
plt.grid(True)

# 条件を満たす境界を描画
plt.contour(U, V, condition, levels=[0.5], colors=['blue'])

plt.show()

2022 Q2(2)

投稿日:2024年1月24日 　最終更新日:2024年7月31日

2変数の独立性の確認と、周辺累積分布関数から周辺確率密度関数を求めました。

コード (2)

#2022 Q2(2) 2024.7.30

from sympy import symbols, integrate

# 変数の定義
u, v = symbols('u v')

# 同時累積分布関数の定義
F_uv = (u*v + u + v + 1) / 4

# 周辺累積分布関数 F1(u) と F2(v) の計算
F1_u = integrate(F_uv, (v, -1, 1))
F2_v = integrate(F_uv, (u, -1, 1))

# 周辺確率密度関数 f1(u) と f2(v) の計算
f1_u = F1_u.diff(u)
f2_v = F2_v.diff(v)

# 同時確率密度関数 f(u,v) を求める
f_uv = F_uv.diff(u).diff(v)

# f(u,v) を f1(u) と f2(v) の積として表す
product_f1_f2 = f1_u * f2_v

# 結果の出力
results = {
    'F1_u': F1_u,
    'F2_v': F2_v,
    'f1_u': f1_u,
    'f2_v': f2_v,
    'f_uv': f_uv,
    'product_f1_f2': product_f1_f2
}

results
{'F1_u': u/2 + 1/2,
 'F2_v': v/2 + 1/2,
 'f1_u': 1/2,
 'f2_v': 1/2,
 'f_uv': 1/4,
 'product_f1_f2': 1/4}

コード(3)

計算

#2022 Q2(3) 2024.7.31

from sympy import symbols, pi
import numpy as np
from scipy.integrate import quad

# 符号変数の定義
u, v = symbols('u v')

# f(u,v) の定義
f_uv_value = 1/4

# 条件を満たす範囲でのf(u,v)（数値積分用）
def integrand(v, u):
    if u**2 + v**2 <= 1:
        return f_uv_value
    else:
        return 0

# uに関しての積分
def integrate_u(u):
    return quad(integrand, -1, 1, args=(u))[0]

# 数値積分による評価
probability_numeric = quad(integrate_u, -1, 1)[0]

# 手計算による値 (π/4) の計算
manual_calculation_value = pi / 4

# 結果の出力
print(f"数値積分による確率: {probability_numeric}\n")
print(f"手計算による値 (π/4): {manual_calculation_value.evalf()}\n")

数値積分による確率: 0.7855319594076084

手計算による値 (π/4): 0.785398163397448

3Dプロット

#2022 Q2(3) 2024.7.31

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# u, vの範囲と値を定義
u = np.linspace(-1, 1, 100)
v = np.linspace(-1, 1, 100)
U, V = np.meshgrid(u, v)

# f(u,v)の値を計算
f_uv_value = 1/4
F = np.where(U**2 + V**2 <= 1, f_uv_value, 0)

# 3Dプロットを作成
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(U, V, F, cmap='viridis')
ax.set_xlabel('U')
ax.set_ylabel('V')
ax.set_zlabel('f(u,v)')
ax.set_zlim(0, 1)  # z軸の範囲を設定
ax.set_title('領域 U^2 + V^2 <= 1 における f(u,v) の3Dプロット')
plt.show()

2022 Q2(1)

投稿日:2024年1月23日 　最終更新日:2024年7月29日

同時累積分布関数から周辺累積分布関数と周辺確率密度関数を求める問題をやりました。

コード（解答）

#2022 Q2(1)  2024.7.29

from sympy import symbols, integrate

# 変数の定義
u, v = symbols('u v')

# 同時累積分布関数の定義
F_uv = (u*v + u + v + 1) / 4

# 周辺累積分布関数 F1(u) と F2(v) の計算
F1_u = integrate(F_uv, (v, -1, 1))
F2_v = integrate(F_uv, (u, -1, 1))

F1_u, F2_v

(u/2 + 1/2, v/2 + 1/2)

コード（グラフ）

#2022 Q2(1)  2024.7.29

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 同時累積分布関数 F(u, v) を定義
def F(u, v):
    return (u * v + u + v + 1) / 4

# u, v の範囲を設定
u = np.linspace(-1, 1, 100)
v = np.linspace(-1, 1, 100)
u, v = np.meshgrid(u, v)
z = F(u, v)

# 3Dプロットの作成
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(u, v, z, cmap='viridis')

# グラフのラベル
ax.set_xlabel('U')
ax.set_ylabel('V')
ax.set_zlabel('F(U, V)')
ax.set_title('同時累積分布関数 F(U, V)')

plt.show()

2022 Q1(4)

投稿日:2024年1月22日 　最終更新日:2024年7月28日

更に条件が追加された場合のP(A∩B∩C)の取り得る範囲の計算。

コード

#2022 Q1(4)  2024.7.28

# 定数
total_points = 16
PA = 3 / 4
PB = 3 / 4
PC = 3 / 4
PAB = 9 / 16
PAC = 9 / 16
PBC = 9 / 16

# 結果を格納するリスト
valid_probabilities = []

# 最小値と最大値を取るときの組み合わせを格納する変数
min_combination = None
max_combination = None
min_value = float('inf')
max_value = float('-inf')

# 各エリアの点数の範囲を狭める
for n1 in range(total_points // 3 + 1):
    for n2 in range((total_points - 3 * n1) // 3 + 1):
        for n3 in range(total_points - 3 * n1 - 3 * n2 + 1):
            n4 = total_points - 3 * n1 - 3 * n2 - n3
            if n4 < 0:
                continue

            calc_PA = (n1 + 2 * n2 + n4) / total_points
            calc_PB = (n1 + 2 * n2 + n4) / total_points
            calc_PC = (n1 + 2 * n2 + n4) / total_points
            calc_PAB = (n2 + n4) / total_points
            calc_PAC = (n2 + n4) / total_points
            calc_PBC = (n2 + n4) / total_points

            # 条件を満たすか確認
            if (abs(calc_PA - PA) < 1e-9 and 
                abs(calc_PB - PB) < 1e-9 and 
                abs(calc_PC - PC) < 1e-9 and 
                abs(calc_PAB - PAB) < 1e-9 and 
                abs(calc_PAC - PAC) < 1e-9 and 
                abs(calc_PBC - PBC) < 1e-9):
                
                current_PABC = n4 / total_points
                valid_probabilities.append(current_PABC)
                
                if current_PABC < min_value:
                    min_value = current_PABC
                    min_combination = (n1, n2, n3, n4)
                
                if current_PABC > max_value:
                    max_value = current_PABC
                    max_combination = (n1, n2, n3, n4)

# 最小値と最大値を表示
min_PABC = min_value if valid_probabilities else None
max_PABC = max_value if valid_probabilities else None

# 結果の表示
print(f"最小値: {min_PABC}")
print(f"最大値: {max_PABC}")
print(f"最小値を取るときのN1, N2, N3, N4: {min_combination}")
print(f"最大値を取るときのN1, N2, N3, N4: {max_combination}")

最小値: 0.375
最大値: 0.4375
最小値を取るときのN1, N2, N3, N4: (0, 3, 1, 6)
最大値を取るときのN1, N2, N3, N4: (1, 2, 0, 7)

アルゴリズム

1.定数の設定

  • 総点数と目標とする確率を設定します。

2.初期化

  • 結果を保存するリストと、最小値・最大値、およびそのときの組み合わせを保存する変数を用意します。

3.探索

  • 3つのネストされたループを使って、各エリアに割り当てる点数のすべての可能な組み合わせを試します。
  • 各組み合わせについて、計算した確率が設定した目標確率と一致するかを確認します。

4.条件を満たす組み合わせの保存

  • 条件を満たす場合、その組み合わせをリストに追加します。
  • また、その組み合わせが現在の最小値または最大値である場合、それを更新します。

5.結果の表示

  • 最小値と最大値、およびそれぞれの組み合わせを表示します。

2022 Q1(3)

投稿日:2024年1月21日 　最終更新日:2024年7月28日

事象A,B,Cが独立とは限らない場合のP(A∩B∩C)の取り得る範囲を求めました。

コード

### 2022 Q1(3)   2024.7.27

import itertools

def generate_combinations(n, k):
    base_list = [1] * k + [0] * (n - k)
    all_combinations = set(itertools.permutations(base_list))
    return list(all_combinations)

def calculate_intersection(A, B, C):
    return np.mean(np.array(A) & np.array(B) & np.array(C))

def find_min_max_intersection(n, k):
    all_combinations = generate_combinations(n, k)
    min_intersection = 1.0
    max_intersection = 0.0

    for A in all_combinations:
        for B in all_combinations:
            for C in all_combinations:
                intersection_ratio = calculate_intersection(A, B, C)
                if intersection_ratio < min_intersection:
                    min_intersection = intersection_ratio
                if intersection_ratio > max_intersection:
                    max_intersection = intersection_ratio

    return min_intersection, max_intersection

# パラメータを設定
n = 8
k = 6

# 最小値と最大値を計算
min_intersection, max_intersection = find_min_max_intersection(n, k)

min_intersection, max_intersection

(0.25, 0.75)

アルゴリズム

  1. 組み合わせの生成
    • 配列の長さ n を 8、1の数 k を 6 と設定する。
    • 配列 11111100 を作成し、 itertools.permutations を使ってすべての並べ替えを生成する。
  2. 交差部分の確率計算
    • すべての組み合わせに対して、配列 A, B, C を選ぶ。
    • 例えば、 A = 11111100, B = 11011110, C = 11101110 の場合、交差部分 A & B & C は 11001100 となる。これにより、確率は \frac{4}{8} = 0.5 となる。
  3. 最小値と最大値の取得
    • すべての組み合わせに対して計算した確率を比較し、最小値と最大値を更新する。
    • 最終的に、最小値と最大値を返す。

2022 Q1(1)(2)

投稿日:2024年1月20日 　最終更新日:2024年7月28日

事象A,Bが独立なときと、そうとは限らない場合の確率の取り得る範囲の計算をしました。

(1)コード

### 2022 Q1(1)   2024.7.26

import numpy as np

# シミュレーションの回数
num_simulations = 1000000

# 与えられた確率
prob = 3/4

# 事象 A, B, C をシミュレーション
A = np.random.rand(num_simulations) < prob
B = np.random.rand(num_simulations) < prob
C = np.random.rand(num_simulations) < prob

# 交差確率の計算
P_A_and_B = np.mean(A & B)
P_A_and_B_and_C = np.mean(A & B & C)

P_A_and_B, P_A_and_B_and_C

(0.56311, 0.422594)
#検算
9/16,27/64
(0.5625, 0.421875)

(2)コード

### 2022 Q1(2)   2024.7.26

import numpy as np

def simulate_overlap_with_fixed_true(correlation):
    P_A = 3/4
    P_B = 3/4
    num_elements = 1000
    num_true_A = int(P_A * num_elements)
    num_true_B = int(P_B * num_elements)

    # A と B の配列を生成し、3/4がtrueになるように設定
    A = np.zeros(num_elements, dtype=bool)
    A[:num_true_A] = True
    np.random.shuffle(A)

    B = np.zeros(num_elements, dtype=bool)
    B[:num_true_B] = True
    np.random.shuffle(B)

    if correlation == 'positive':
        # 正の相関を持たせるために、スワップを行う
        A1_B0_indices = np.where(A & ~B)[0]
        A0_B1_indices = np.where(~A & B)[0]

        min_len = min(len(A1_B0_indices), len(A0_B1_indices))
        for i in range(min_len):
            B[A1_B0_indices[i]], B[A0_B1_indices[i]] = B[A0_B1_indices[i]], B[A1_B0_indices[i]]

    elif correlation == 'negative':
        # 負の相関を持たせるために、スワップを行う
        A1_B1_indices = np.where(A & B)[0]
        A0_B0_indices = np.where(~A & ~B)[0]

        min_len = min(len(A1_B1_indices), len(A0_B0_indices))
        for i in range(min_len):
            B[A1_B1_indices[i]], B[A0_B0_indices[i]] = B[A0_B0_indices[i]], B[A1_B1_indices[i]]

    else:
        raise ValueError("Invalid correlation type")

    # 重なりの割合を計算
    overlap_ratio = np.mean(A & B)

    return overlap_ratio

# 正の相関の場合
overlap_positive = simulate_overlap_with_fixed_true('positive')

# 負の相関の場合
overlap_negative = simulate_overlap_with_fixed_true('negative')

overlap_positive, overlap_negative
(0.75, 0.5)

アルゴリズム

  1. パラメータの設定
  • 事象 (A) の発生確率 (P_A) を設定する（ここでは 3/4）。
  • 事象 (B) の発生確率 (P_B) を設定する（ここでは 3/4）。
  • 配列の要素数 (num_elements) を設定する（ここでは 1000）。
  • シミュレーションの回数 (num_simulations) を設定する（ここでは 1000）。
  1. 配列の初期化
  • 配列 (A) と (B) をそれぞれ3/4がtrueになるように生成する。
    • 配列 (A) に対して、最初の (num_true_A) 個を true に設定し、残りを false に設定する。その後、配列をシャッフルする。
    • 配列 (B) に対しても同様の処理を行う。
  1. 相関の付与
  • 正の相関を持たせる場合：
    • 配列 (A) と (B) の要素をスワップして、 (A = 1, B = 0) と (A = 0, B = 1) の組を見つけたら、これらをスワップして値を揃える。
  • 負の相関を持たせる場合：
    • 配列 (A) と (B) の要素をスワップして、 (A = 1, B = 1) と (A = 0, B = 0) の組を見つけたら、これらをスワップして値が揃わないようにする。
  1. 重なりの割合の計算
  • 配列 (A) と (B) の共通部分（重なり部分）を計算し、その割合を求める。
  1. 結果の記録
  • 重なりの割合をリストに記録する。
  • シミュレーション回数分、上記のステップを繰り返す。
  1. 最小値と最大値の取得
  • 重なりの割合のリストから、最小値と最大値を取得する。

投稿ナビゲーション

← 過去の投稿 新しい投稿 →