114 lines
4.6 KiB
Python
114 lines
4.6 KiB
Python
import random
|
|
random.seed(100)
|
|
|
|
# Defining the local conditional probability distributions
|
|
P_F = {1: 0.01, 0: 0.99}
|
|
P_E = {1: 0.02, 0: 0.98}
|
|
|
|
# Problem 1a: Probabilities from Network (replace with your actual values)
|
|
P_A_given_FBE = {
|
|
(1, 1, 1): 0.99, # Alarm likely if fire, earthquake, and neighbor B calls
|
|
(1, 1, 0): 0.98, # Alarm likely if fire, earthquake, but no call from B
|
|
(1, 0, 1): 0.95, # Alarm likely if fire, no earthquake, but neighbor B calls
|
|
(1, 0, 0): 0.90, # Alarm likely if fire, no earthquake, and no call from B
|
|
(0, 1, 1): 0.70, # Alarm possible if earthquake and neighbor B calls, but no fire
|
|
(0, 1, 0): 0.60, # Alarm possible if earthquake, no fire, and no call from B
|
|
(0, 0, 1): 0.20, # Alarm less likely if no fire, no earthquake, but neighbor B calls
|
|
(0, 0, 0): 0.01 # Alarm very unlikely if no fire, no earthquake, and no call from B
|
|
}
|
|
P_B_given_A = {1: 0.90, 0: 0.05} # Probability of neighbor B calling given alarm state
|
|
P_C_given_A = {1: 0.70, 0.01: 0.01} # Probability of neighbor C calling given alarm state
|
|
|
|
|
|
def joint_probability(f, e, a, b, c):
|
|
return (P_F[f] * P_E[e] * P_A_given_FBE[(f, e, b)] *
|
|
P_B_given_A.get(a if b else 1-a, 0.5) * # Handle missing probabilities
|
|
P_C_given_A.get(a if c else 1-a, 0.5)) # Handle missing probabilities
|
|
|
|
|
|
# Problem 1a: Calculate P(F=1 | B=1) and P(F=1 | C=1)
|
|
def calculate_conditional_probs():
|
|
total_P_B_1 = sum([joint_probability(f, e, a, 1, c) for f in [True, False] for e in [True, False]
|
|
for a in [True, False] for c in [True, True]]) # Only consider B=1
|
|
total_P_C_1 = sum([joint_probability(f, e, a, b, 1) for f in [True, False] for e in [True, False]
|
|
for a in [True, False] for b in [True, True]]) # Only consider C=1
|
|
|
|
P_F_given_B = total_P_B_1 / sum([joint_probability(f, e, a, 1, c) for f in [True, False] for e in [True, False]
|
|
for a in [True, False] for c in [True, True]])
|
|
|
|
P_F_given_C = total_P_C_1 / sum([joint_probability(f, e, a, b, 1) for f in [True, False] for e in [True, False]
|
|
for a in [True, False] for b in [True, True]])
|
|
|
|
# Print the results for Problem 1a
|
|
print("Problem 1a:")
|
|
print(f"P(F=1 | B=1): {P_F_given_B:.6f}")
|
|
print(f"P(F=1 | C=1): {P_F_given_C:.6f}")
|
|
print("") # Empty line for separation
|
|
|
|
def rejection_sampling(N, evidence):
|
|
"""
|
|
Performs rejection sampling for a given evidence (B=1 or C=1) and number of samples (N).
|
|
|
|
Args:
|
|
N: Number of samples to draw.
|
|
evidence: Evidence variable (e.g., "B=1" or "C=1").
|
|
|
|
Returns:
|
|
A tuple containing the estimated probability and the number of rejected samples.
|
|
"""
|
|
count_accepted = 0
|
|
rejected_samples = 0
|
|
for _ in range(N):
|
|
f = random.choices([0, 1], weights=[P_F[v] for v in [0, 1]])[0] # Sample fire (f)
|
|
e = random.choices([0, 1], weights=[P_E[v] for v in [0, 1]])[0] # Sample earthquake (e)
|
|
|
|
# Handle missing probabilities in P_A_given_FBE (if needed)
|
|
a = P_A_given_FBE.get((f, e, 1), 0.5) # Use default 0.5 if (f, e) is missing
|
|
|
|
if evidence == "B=1":
|
|
b = 1 # Force neighbor B calling to be True based on evidence
|
|
else:
|
|
b = random.choices([0, 1], weights=[P_B_given_A.get(a, 0.5) for a in [0, 1]])[0] # Use default 0.5 if a is undefined
|
|
c = random.choices([0, 1], weights=[P_C_given_A.get(a, 0.5) for a in [0, 1]])[0] # Use default 0.5 if a is undefined
|
|
|
|
p = joint_probability(f, e, a, b, c)
|
|
u = random.random()
|
|
if u < p:
|
|
count_accepted += 1 if f == 1 else 0
|
|
rejected_samples = N - count_accepted
|
|
return count_accepted / N, rejected_samples
|
|
|
|
|
|
def main():
|
|
# True values from previous calculations (replace with your actual values)
|
|
true_prob_B = 0.5 # P(F=1 | B=1)
|
|
true_prob_C = 0.5 # P(F=1 | C=1)
|
|
|
|
# Sample sizes
|
|
sample_sizes = [10, 100, 1000, 10000]
|
|
|
|
# Estimated probabilities for B=1 and C=1
|
|
estimated_prob_B = []
|
|
estimated_prob_C = []
|
|
rejected_samples_B = []
|
|
rejected_samples_C = []
|
|
|
|
for N in sample_sizes:
|
|
est_prob_B, rejected_B = rejection_sampling(N, "B=1")
|
|
est_prob_C, rejected_C = rejection_sampling(N, "C=1")
|
|
|
|
estimated_prob_B.append(est_prob_B)
|
|
estimated_prob_C.append(est_prob_C)
|
|
rejected_samples_B.append(rejected_B)
|
|
rejected_samples_C.append(rejected_C)
|
|
|
|
# Print the results for Problem 1b
|
|
print("Problem 1b:")
|
|
print("Sample Size\tEstimated P(F=1 | B=1)\tEstimated P(F=1 | C=1)")
|
|
for i, N in enumerate(sample_sizes):
|
|
print(f"{N}\t\t{estimated_prob_B[i]:.6f}\t\t{estimated_prob_C[i]:.6f}")
|
|
print("") # Empty line for separation
|
|
|
|
if __name__ == "__main__":
|
|
main()
|
|
calculate_conditional_probs() |