44import uuid
55from collections .abc import Callable
66from datetime import datetime , timedelta
7+ from functools import partial
78from typing import TYPE_CHECKING
89
910import numpy as np
2324logger = logging .getLogger ('ABC' )
2425
2526
27+ def _simulate_one_from_prior (
28+ model_prior : RV ,
29+ parameter_priors : list [Distribution ],
30+ models : list [Model ],
31+ summary_statistics : Callable ,
32+ ):
33+ """Sample one particle from the prior."""
34+ from ..population import Particle
35+
36+ # sample model
37+ m = int (model_prior .rvs ())
38+ # sample parameter
39+ theta = parameter_priors [m ].rvs ()
40+ # simulate summary statistics
41+ model_result = models [m ].summary_statistics (0 , theta , summary_statistics )
42+ # sampled from prior, so all have uniform weight
43+ weight = 1.0
44+ # distance will be computed after initialization of the
45+ # distance function
46+ distance = np .inf
47+ # all are happy and accepted
48+ accepted = True
49+
50+ return Particle (
51+ m = m ,
52+ parameter = theta ,
53+ weight = weight ,
54+ sum_stat = model_result .sum_stat ,
55+ distance = distance ,
56+ accepted = accepted ,
57+ proposal_id = 0 ,
58+ preliminary = False ,
59+ )
60+
61+
62+ def _simulate_one (
63+ * ,
64+ t : int ,
65+ m : np .ndarray ,
66+ p : np .ndarray ,
67+ model_prior : RV ,
68+ parameter_priors : list [Distribution ],
69+ model_perturbation_kernel : ModelPerturbationKernel ,
70+ transitions : list [Transition ],
71+ models : list [Model ],
72+ summary_statistics : Callable ,
73+ x_0 : dict ,
74+ distance_function : Distance ,
75+ eps : Epsilon ,
76+ acceptor : Acceptor ,
77+ weight_function : Callable ,
78+ evaluate : bool ,
79+ proposal_id : int ,
80+ ):
81+ """Sample one parameter and evaluate/simulate one particle."""
82+ parameter = generate_valid_proposal (
83+ t = t ,
84+ m = m ,
85+ p = p ,
86+ model_prior = model_prior ,
87+ parameter_priors = parameter_priors ,
88+ model_perturbation_kernel = model_perturbation_kernel ,
89+ transitions = transitions ,
90+ )
91+ if evaluate :
92+ particle = evaluate_proposal (
93+ * parameter ,
94+ t = t ,
95+ models = models ,
96+ summary_statistics = summary_statistics ,
97+ distance_function = distance_function ,
98+ eps = eps ,
99+ acceptor = acceptor ,
100+ x_0 = x_0 ,
101+ weight_function = weight_function ,
102+ proposal_id = proposal_id ,
103+ )
104+ else :
105+ particle = only_simulate_data_for_proposal (
106+ * parameter ,
107+ t = t ,
108+ models = models ,
109+ summary_statistics = summary_statistics ,
110+ weight_function = weight_function ,
111+ proposal_id = proposal_id ,
112+ )
113+ return particle
114+
115+
116+ def _prior_pdf (
117+ m_ss : int ,
118+ theta_ss : Parameter ,
119+ model_prior : RV ,
120+ parameter_priors : list [Distribution ],
121+ ) -> float :
122+ """Evaluate the prior density for a proposed sample."""
123+ return model_prior .pmf (m_ss ) * parameter_priors [m_ss ].pdf (theta_ss )
124+
125+
126+ def _transition_pdf (
127+ m_ss : int ,
128+ theta_ss : Parameter ,
129+ transitions : list [Transition ],
130+ model_probabilities : pd .DataFrame ,
131+ model_perturbation_kernel : ModelPerturbationKernel ,
132+ ) -> float :
133+ """Evaluate the transition density for a proposed sample."""
134+ model_factor = sum (
135+ row .p * model_perturbation_kernel .pmf (m_ss , m )
136+ for m , row in model_probabilities .iterrows ()
137+ )
138+ particle_factor = transitions [m_ss ].pdf (theta_ss )
139+
140+ transition_pd = model_factor * particle_factor
141+ if transition_pd == 0 :
142+ logger .debug ('Transition density is zero!' )
143+ return transition_pd
144+
145+
146+ def _weight_function (
147+ m_ss : int ,
148+ theta_ss : Parameter ,
149+ acceptance_weight : float ,
150+ prior_pdf : Callable ,
151+ transition_pdf : Callable ,
152+ ) -> float :
153+ """Calculate total weight from sampling and acceptance weight."""
154+ prior_pd = prior_pdf (m_ss , theta_ss )
155+ transition_pd = transition_pdf (m_ss , theta_ss )
156+ return acceptance_weight * prior_pd / transition_pd
157+
158+
26159class AnalysisVars :
27160 """Contract object class for passing analysis variables.
28161
@@ -97,38 +230,13 @@ def create_simulate_from_prior_function(
97230 simulate_one:
98231 A function that returns a sampled particle.
99232 """
100- # simulation function, simplifying some parts compared to later
101- from ..population import Particle
102-
103- def simulate_one ():
104- # sample model
105- m = int (model_prior .rvs ())
106- # sample parameter
107- theta = parameter_priors [m ].rvs ()
108- # simulate summary statistics
109- model_result = models [m ].summary_statistics (
110- 0 , theta , summary_statistics
111- )
112- # sampled from prior, so all have uniform weight
113- weight = 1.0
114- # distance will be computed after initialization of the
115- # distance function
116- distance = np .inf
117- # all are happy and accepted
118- accepted = True
119-
120- return Particle (
121- m = m ,
122- parameter = theta ,
123- weight = weight ,
124- sum_stat = model_result .sum_stat ,
125- distance = distance ,
126- accepted = accepted ,
127- proposal_id = 0 ,
128- preliminary = False ,
129- )
130-
131- return simulate_one
233+ return partial (
234+ _simulate_one_from_prior ,
235+ model_prior = model_prior ,
236+ parameter_priors = parameter_priors ,
237+ models = models ,
238+ summary_statistics = summary_statistics ,
239+ )
132240
133241
134242def generate_valid_proposal (
@@ -273,11 +381,11 @@ def create_prior_pdf(
273381 prior_pdf: The prior density function.
274382 """
275383
276- def prior_pdf ( m_ss , theta_ss ):
277- prior_pd = model_prior . pmf ( m_ss ) * parameter_priors [ m_ss ]. pdf ( theta_ss )
278- return prior_pd
279-
280- return prior_pdf
384+ return partial (
385+ _prior_pdf ,
386+ model_prior = model_prior ,
387+ parameter_priors = parameter_priors ,
388+ )
281389
282390
283391def create_transition_pdf (
@@ -298,20 +406,12 @@ def create_transition_pdf(
298406 transition_pdf: The transition density function.
299407 """
300408
301- def transition_pdf (m_ss , theta_ss ):
302- model_factor = sum (
303- row .p * model_perturbation_kernel .pmf (m_ss , m )
304- for m , row in model_probabilities .iterrows ()
305- )
306- particle_factor = transitions [m_ss ].pdf (theta_ss )
307-
308- transition_pd = model_factor * particle_factor
309-
310- if transition_pd == 0 :
311- logger .debug ('Transition density is zero!' )
312- return transition_pd
313-
314- return transition_pdf
409+ return partial (
410+ _transition_pdf ,
411+ transitions = transitions ,
412+ model_probabilities = model_probabilities ,
413+ model_perturbation_kernel = model_perturbation_kernel ,
414+ )
315415
316416
317417def create_weight_function (
@@ -332,27 +432,11 @@ def create_weight_function(
332432 weight_function: The importance sample weight function.
333433 """
334434
335- def weight_function (m_ss , theta_ss , acceptance_weight : float ):
336- """Calculate total weight, from sampling and acceptance weight.
337-
338- Parameters
339- ----------
340- m_ss: The model sample.
341- theta_ss: The parameter sample.
342- acceptance_weight: The acceptance weight sample. In most cases 1.
343-
344- Returns
345- -------
346- weight: The total weight.
347- """
348- # prior and transition density (can be equal)
349- prior_pd = prior_pdf (m_ss , theta_ss )
350- transition_pd = transition_pdf (m_ss , theta_ss )
351- # calculate weight
352- weight = acceptance_weight * prior_pd / transition_pd
353- return weight
354-
355- return weight_function
435+ return partial (
436+ _weight_function ,
437+ prior_pdf = prior_pdf ,
438+ transition_pdf = transition_pdf ,
439+ )
356440
357441
358442def create_simulate_function (
@@ -431,42 +515,25 @@ def create_simulate_function(
431515 prior_pdf = prior_pdf , transition_pdf = transition_pdf
432516 )
433517
434- # simulation function
435- def simulate_one ():
436- parameter = generate_valid_proposal (
437- t = t ,
438- m = m ,
439- p = p ,
440- model_prior = model_prior ,
441- parameter_priors = parameter_priors ,
442- model_perturbation_kernel = model_perturbation_kernel ,
443- transitions = transitions ,
444- )
445- if evaluate :
446- particle = evaluate_proposal (
447- * parameter ,
448- t = t ,
449- models = models ,
450- summary_statistics = summary_statistics ,
451- distance_function = distance_function ,
452- eps = eps ,
453- acceptor = acceptor ,
454- x_0 = x_0 ,
455- weight_function = weight_function ,
456- proposal_id = proposal_id ,
457- )
458- else :
459- particle = only_simulate_data_for_proposal (
460- * parameter ,
461- t = t ,
462- models = models ,
463- summary_statistics = summary_statistics ,
464- weight_function = weight_function ,
465- proposal_id = proposal_id ,
466- )
467- return particle
468-
469- return simulate_one
518+ return partial (
519+ _simulate_one ,
520+ t = t ,
521+ m = m ,
522+ p = p ,
523+ model_prior = model_prior ,
524+ parameter_priors = parameter_priors ,
525+ model_perturbation_kernel = model_perturbation_kernel ,
526+ transitions = transitions ,
527+ models = models ,
528+ summary_statistics = summary_statistics ,
529+ x_0 = x_0 ,
530+ distance_function = distance_function ,
531+ eps = eps ,
532+ acceptor = acceptor ,
533+ weight_function = weight_function ,
534+ evaluate = evaluate ,
535+ proposal_id = proposal_id ,
536+ )
470537
471538
472539def only_simulate_data_for_proposal (
0 commit comments